StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
combineHistograms0.C
1 void combineHistograms0(const char *dirName, const char **inNames, const char *outName, int nCent) {
2 
3  // -- Calculate \Delta\rho/\sqrt{\rho_{ref}} for all quantities, i.e. Phi-Phi, Eta-Eta, DPhi-DEta etc.
4  //
5  // root.exe -q -b combineHistograms0.C("dir","inFileList","outFileName",numFiles)
6  //
7  // inFileList contains centrality information, so I can have a list
8  // like {"Data1", "Data2",...} or {"Sum0_1", "Sum2_3",...}
9  // Before getting here we use hadd to add histograms from different jobs
10  // creating files we usually call Data{n}.root (where n is a centrality bin)
11  // then we may combine centralities be merging the Data{n} files tagging
12  // histograms with z-vertex index and creating files Sum{n1}_{n2}.root where
13  // n1 and n2 are the limits of the centrality bin range.
14 
15  gROOT->LoadMacro("load2ptLibs.C");
16  load2ptLibs();
17  gSystem->Load("StEStructPoolSupport.so");
18 
19  TFile *tf;
20  StEStructSupport *ehelp;
21  TH2D **ptdedp;
22  TH2D **ptdedpC;
23  TH2D **dedp;
24  TH2D **dedpC;
25  TH2D **ptetaeta;
26  TH2D **ptetaetaC;
27  TH2D **etaeta;
28  TH2D **etaetaC;
29  TH2D **ptphiphi;
30  TH2D **ptphiphiC;
31  TH2D **phiphi;
32  TH2D **phiphiC;
33 
34  const char* binName[]={"Symm"};
35  const char* chargeName[] = {"_LS_", "_US_", "_CD_", "_CI_"};
36  const char* chargeType[] = {"_PP_", "_PM_", "_MP_", "_MM_"};
37  char outFileName[1024];
38  sprintf(outFileName,"%s/%s.root",dirName,outName);
39  TFile *out = new TFile(outFileName,"RECREATE");
40 
41  char inFileName[1024];
42 
43  for (int ic=0;ic<nCent;ic++) {
44  printf("Scale factors for centrality %2i, ++ +- -+ -- CI \n",ic);
45  for (int ibin=0;ibin<1;ibin++) {
46  sprintf(inFileName,"%s/%s%s.root",dirName,inNames[ic],binName[ibin]);
47  tf = new TFile(inFileName);
48  tf->cd();
49  ehelp = new StEStructSupport(tf,0);
50  ehelp->msilent = true;
51  ehelp->mapplyDEtaFix = false;
52  ehelp->mPairNormalization = true;
53  ehelp->setBGMode(0);
54  ehelp->mIdenticalPair = true;
55  int subtract = 1;
56 
57  ptdedpC = ehelp->buildPtCommon("DEtaDPhi",2,subtract);
58  ptetaetaC = ehelp->buildPtCommon("EtaEta",2,subtract);
59  ptphiphiC = ehelp->buildPtCommon("PhiPhi",2,subtract);
60 
61  ptdedp = ehelp->buildPtChargeTypes("DEtaDPhi",2,subtract);
62  ptetaeta = ehelp->buildPtChargeTypes("EtaEta",2,subtract);
63  ptphiphi = ehelp->buildPtChargeTypes("PhiPhi",2,subtract);
64 
65  dedpC = ehelp->buildCommon("DEtaDPhi",2);
66  etaetaC = ehelp->buildCommon("EtaEta",2);
67  phiphiC = ehelp->buildCommon("PhiPhi",2);
68 
69  dedp = ehelp->buildChargeTypes("DEtaDPhi",2);
70  etaeta = ehelp->buildChargeTypes("EtaEta",2);
71  phiphi = ehelp->buildChargeTypes("PhiPhi",2);
72 
73  out->cd();
74  for (int icharge=0;icharge<4;icharge++) {
75  TString name("NDEtaDPhi"); name += chargeName[icharge]; name += ic;
76  if (dedp) {
77  dedp[icharge]->SetName(name.Data());
78  dedp[icharge]->SetTitle(name.Data());
79  dedp[icharge]->Write();
80  }
81  TString name("PtDEtaDPhi"); name += chargeName[icharge]; name += ic;
82  if (ptdedp) {
83  ptdedp[icharge]->SetName(name.Data());
84  ptdedp[icharge]->SetTitle(name.Data());
85  ptdedp[icharge]->Write();
86  }
87 
88  TString name("NDEtaDPhi"); name += chargeType[icharge]; name += ic;
89  if (dedpC) {
90  dedpC[icharge]->SetName(name.Data());
91  dedpC[icharge]->SetTitle(name.Data());
92  dedpC[icharge]->Write();
93  }
94  TString name("PtDEtaDPhi"); name += chargeType[icharge]; name += ic;
95  if (ptdedpC) {
96  ptdedpC[icharge]->SetName(name.Data());
97  ptdedpC[icharge]->SetTitle(name.Data());
98  ptdedpC[icharge]->Write();
99  }
100 
101 
102  TString name("NEtaEta"); name += chargeName[icharge]; name += ic;
103  if (etaeta) {
104  etaeta[icharge]->SetName(name.Data());
105  etaeta[icharge]->SetTitle(name.Data());
106  etaeta[icharge]->Write();
107  }
108  TString name("PtEtaEta"); name += chargeName[icharge]; name += ic;
109  if (ptetaeta) {
110  ptetaeta[icharge]->SetName(name.Data());
111  ptetaeta[icharge]->SetTitle(name.Data());
112  ptetaeta[icharge]->Write();
113  }
114 
115  TString name("NEtaEta"); name += chargeType[icharge]; name += ic;
116  if (etaetaC) {
117  etaetaC[icharge]->SetName(name.Data());
118  etaetaC[icharge]->SetTitle(name.Data());
119  etaetaC[icharge]->Write();
120  }
121  TString name("PtEtaEta"); name += chargeType[icharge]; name += ic;
122  if (ptetaetaC) {
123  ptetaetaC[icharge]->SetName(name.Data());
124  ptetaetaC[icharge]->SetTitle(name.Data());
125  ptetaetaC[icharge]->Write();
126  }
127 
128 
129  TString name("NPhiPhi"); name += chargeName[icharge]; name += ic;
130  if (phiphi) {
131  phiphi[icharge]->SetName(name.Data());
132  phiphi[icharge]->SetTitle(name.Data());
133  phiphi[icharge]->Write();
134  }
135  TString name("PtPhiPhi"); name += chargeName[icharge]; name += ic;
136  if (ptphiphi) {
137  ptphiphi[icharge]->SetName(name.Data());
138  ptphiphi[icharge]->SetTitle(name.Data());
139  ptphiphi[icharge]->Write();
140  }
141 
142  TString name("NPhiPhi"); name += chargeType[icharge]; name += ic;
143  if (phiphiC) {
144  phiphiC[icharge]->SetName(name.Data());
145  phiphiC[icharge]->SetTitle(name.Data());
146  phiphiC[icharge]->Write();
147  }
148  TString name("PtPhiPhi"); name += chargeType[icharge]; name += ic;
149  if (ptphiphiC) {
150  ptphiphiC[icharge]->SetName(name.Data());
151  ptphiphiC[icharge]->SetTitle(name.Data());
152  ptphiphiC[icharge]->Write();
153  }
154  }
155  // Calculate and print scale factors
156  double *scale = ehelp->getScaleFactors();
157  printf(" %7.2f %7.2f %7.2f %7.2f %7.2f\n",scale[0],scale[1],scale[2],scale[3],scale[4]);
158  delete [] scale;
159  delete tf;
160  delete ehelp;
161  }
162  }
163 
164  TH2D **ytyt;
165  TH2D **ytytC;
166  TH2D **sytdyt;
167  TH2D **sytdytC;
168  gROOT->LoadMacro("minimizeNegative.C");
169  double *sFactor[2][1], *eSFactor[2][1];
170  for (int it=0;it<2;it++) {
171  for (int ibin=0;ibin<1;ibin++) {
172  sFactor[it][ibin] = new double[nCent];
173  eSFactor[it][ibin] = new double[nCent];
174  }
175  }
176  float sf[2];
177  double argList[10];
178  double start = 0.95;
179  double step = 0.001;
180  double bmin = 0.9;
181  double bmax = 1.1;
182  int errFlag = 0;
183  minData.mCorrType = 1;
184  minData.mLambda = 10;
185 
186  // For ytyt space soft/neck/hard is not useful.
187  int ytytBins[] = {0};
188  for (int ic=0;ic<nCent;ic++) {
189  for (int ibin=0;ibin<1;ibin++) {
190  sprintf(inFileName,"%s/%s%s.root",dirName,inNames[ic],binName[ytytBins[ibin]]);
191  tf = new TFile(inFileName);
192  tf->cd();
193  ehelp = new StEStructSupport(tf,0);
194  ehelp->msilent = true;
195  ehelp->mapplyDEtaFix = false;
196  ehelp->mPairNormalization = true;
197  ehelp->mIdenticalPair = true;
198  ehelp->setBGMode(1);
199 
200  // Make sure LS histogram for zBin 0 is in file.
201  if (0 == ehelp->histogramExists("YtYt", 0)) {
202  continue;
203  }
204 
205  // A lot of stuff so we can find a scaling factor for \rho_{ref}
206  // such that \Delta\rho is almost always positive.
207  minData.mSupport = ehelp;
208  minData.mChargeType = 0;
209 
210  minuit = new TMinuit(1);
211  minuit->SetFCN(minimizeNegative);
212  minuit->mnparm( 0, "rho_ref scale factor", start, step, bmin, bmax, errFlag );
213  minuit->SetErrorDef(1);
214  minuit->SetPrintLevel(-1);
215  argList[0] = 1;
216  minuit->mnexcm("SET STR",argList,1,errFlag);
217  argList[0] = 500;
218  cout << ">>>>>Starting scale factor 0 fit for centrality " << ic << " yt bin " << ibin << endl;
219  minuit->mnexcm("MIGRAD",argList,1,errFlag);
220  minuit->GetParameter(0,sFactor[0][ibin][ic],eSFactor[0][ibin][ic]);
221  if (0 != errFlag) {
222  cout << "++++Out of chesse error; Fit failed. Using scale factor = 1" << endl;
223  sFactor[0][ibin][ic] = 1;
224  eSFactor[0][ibin][ic] = 100;
225  }
226  delete minuit;
227 
228  // Seems like I should be able to reset the TMinuit object to do a
229  // different fit. Didn't work on my first attempts, so I just create
230  // a new one.
231  minData.mChargeType = 1;
232  minuit = new TMinuit(1);
233  minuit->SetFCN(minimizeNegative);
234  int errFlag = 0;
235  minuit->mnparm( 0, "rho_ref scale factor", start, step, bmin, bmax, errFlag );
236  minuit->SetErrorDef(1);
237  minuit->SetPrintLevel(-1);
238  argList[0] = 1;
239  minuit->mnexcm("SET STR",argList,1,errFlag);
240  argList[0] = 500;
241  cout << ">>>>>Starting scale factor 1 fit for centrality " << ic << " yt bin " << ibin << endl;
242  minuit->mnexcm("MIGRAD",argList,1,errFlag);
243  minuit->GetParameter(0,sFactor[1][ibin][ic],eSFactor[1][ibin][ic]);
244  if (0 != errFlag) {
245  cout << "++++Out of chesse error; Fit failed. Using scale factor = 1" << endl;
246  sFactor[1][ibin][ic] = 1;
247  eSFactor[1][ibin][ic] = 100;
248  }
249  delete minuit;
250 
251  sf[0] = sFactor[0][ibin][ic];
252  sf[1] = sFactor[1][ibin][ic];
253  ytyt = ehelp->buildChargeTypes("YtYt",5,sf);
254  ytytC = ehelp->buildCommon("YtYt",5,sf);
255  sytdyt = ehelp->buildChargeTypes("SYtDYt",5,sf);
256  sytdytC = ehelp->buildCommon("SYtDYt",5,sf);
257 
258  out->cd();
259  for (int icharge=0;icharge<4;icharge++) {
260  TString name("YtYt"); name += chargeName[icharge]; name += ic;
261  ytyt[icharge]->SetName(name.Data());
262  ytyt[icharge]->SetTitle(name.Data());
263  ytyt[icharge]->Write();
264  TString name("YtYt"); name += chargeType[icharge]; name += ic;
265  ytytC[icharge]->SetName(name.Data());
266  ytytC[icharge]->SetTitle(name.Data());
267  ytytC[icharge]->Write();
268 
269  TString name("SYtDYt"); name += chargeName[icharge]; name += ic;
270  sytdyt[icharge]->SetName(name.Data());
271  sytdyt[icharge]->SetTitle(name.Data());
272  sytdyt[icharge]->Write();
273  TString name("SYtDYt"); name += chargeType[icharge]; name += ic;
274  sytdytC[icharge]->SetName(name.Data());
275  sytdytC[icharge]->SetTitle(name.Data());
276  sytdytC[icharge]->Write();
277  }
278  delete tf;
279  delete ehelp;
280  }
281  }
282  cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << endl;
283  cout << " Here are Yt scale factors for rho_{ref}" << endl;
284  printf(" Centrality all LS all US \n");
285  for (int ic=0;ic<nCent;ic++) {
286  printf("%6i %9.3f(%3.0f) %7.3f(%3.0f)\n",ic,
287  sFactor[0][0][ic],1000*eSFactor[0][0][ic],sFactor[1][0][ic],1000*eSFactor[1][0][ic]);
288  }
289  out->Close();
290 // delete out;
291  for (int it=0;it<2;it++) {
292  for (int ibin=0;ibin<1;ibin++) {
293  delete sFactor[it][ibin];
294  delete eSFactor[it][ibin];
295  }
296  }
297 }
Generated by   doxygen 1.8.5