StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
combineHistograms3.C
1 void combineHistograms3(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 combineHistograms3.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 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 **ptsedp;
22  TH2D **ptsedpC;
23  TH2D **ptdedp;
24  TH2D **ptdedpC;
25  TH2D **sedp;
26  TH2D **sedpC;
27  TH2D **dedp;
28  TH2D **dedpC;
29  TH2D **ptetaeta;
30  TH2D **ptetaetaC;
31  TH2D **etaeta;
32  TH2D **etaetaC;
33  TH2D **ptetaetaSS;
34  TH2D **ptetaetaSSC;
35  TH2D **etaetaSS;
36  TH2D **etaetaSSC;
37  TH2D **ptetaetaAS;
38  TH2D **ptetaetaASC;
39  TH2D **etaetaAS;
40  TH2D **etaetaASC;
41  TH2D **ptphiphi;
42  TH2D **ptphiphiC;
43  TH2D **phiphi;
44  TH2D **phiphiC;
45 
46  const char* binName[]={"all","awayside","nearside","soft","softAS","softNS","neck","neckAS","neckNS","hard","hardAS","hardNS","softHard","softHardAS","softHardNS"};
47  const char* chargeName[] = {"_LS_", "_US_", "_CD_", "_CI_"};
48  const char* chargeType[] = {"_PP_", "_PM_", "_MP_", "_MM_"};
49  char outFileName[1024];
50  sprintf(outFileName,"%s/%s.root",dirName,outName);
51  TFile *out = new TFile(outFileName,"RECREATE");
52 
53  char inFileName[1024];
54 
55  // For axial space near/away side is not useful.
56  int axialBins[] = {0, 3, 6, 9, 12};
57  for (int ic=0;ic<nCent;ic++) {
58  printf("Centrality %2i: d2N/dEdP pHat_{A+} pHat_{A-} pHat_{B+} pHat_{B-}\n",ic);
59  for (int ibin=0;ibin<5;ibin++) {
60  sprintf(inFileName,"%s/%s%s.root",dirName,inNames[ic],binName[axialBins[ibin]]);
61  tf = new TFile(inFileName);
62  tf->cd();
63  ehelp = new StEStructSupport(tf,0);
64  ehelp->msilent = true;
65  ehelp->mapplyDEtaFix = false;
66  ehelp->mPairNormalization = true;
67  ehelp->mIdenticalPair = true;
68  ehelp->setBGMode(0);
69  int subtract = 1;
70 
71  ptsedpC = ehelp->buildPtCommon("SEtaDPhi",2,subtract);
72  ptdedpC = ehelp->buildPtCommon("DEtaDPhi",2,subtract);
73  ptetaetaC = ehelp->buildPtCommon("EtaEta",2,subtract);
74  ptetaetaSSC = ehelp->buildPtCommon("EtaEtaSS",2,subtract);
75  ptetaetaASC = ehelp->buildPtCommon("EtaEtaAS",2,subtract);
76  ptphiphiC = ehelp->buildPtCommon("PhiPhi",2,subtract);
77 
78  ptsedp = ehelp->buildPtChargeTypes("SEtaDPhi",2,subtract);
79  ptdedp = ehelp->buildPtChargeTypes("DEtaDPhi",2,subtract);
80  ptetaeta = ehelp->buildPtChargeTypes("EtaEta",2,subtract);
81  ptetaetaSS = ehelp->buildPtChargeTypes("EtaEtaSS",2,subtract);
82  ptetaetaAS = ehelp->buildPtChargeTypes("EtaEtaAS",2,subtract);
83  ptphiphi = ehelp->buildPtChargeTypes("PhiPhi",2,subtract);
84 
85  sedpC = ehelp->buildCommon("SEtaDPhi",2);
86  dedpC = ehelp->buildCommon("DEtaDPhi",2);
87  etaetaC = ehelp->buildCommon("EtaEta",2);
88  etaetaSSC = ehelp->buildCommon("EtaEtaSS",2);
89  etaetaASC = ehelp->buildCommon("EtaEtaAS",2);
90  phiphiC = ehelp->buildCommon("PhiPhi",2);
91 
92  sedp = ehelp->buildChargeTypes("SEtaDPhi",2);
93  dedp = ehelp->buildChargeTypes("DEtaDPhi",2);
94  etaeta = ehelp->buildChargeTypes("EtaEta",2);
95  etaetaSS = ehelp->buildChargeTypes("EtaEtaSS",2);
96  etaetaAS = ehelp->buildChargeTypes("EtaEtaAS",2);
97  phiphi = ehelp->buildChargeTypes("PhiPhi",2);
98 
99  out->cd();
100  for (int icharge=0;icharge<4;icharge++) {
101  TString name(binName[axialBins[ibin]]);
102  name += "_NDEtaDPhi"; name += chargeName[icharge]; name += ic;
103  if (dedp) {
104  dedp[icharge]->SetName(name.Data());
105  dedp[icharge]->SetTitle(name.Data());
106  dedp[icharge]->Write();
107  }
108  TString name(binName[axialBins[ibin]]);
109  name += "_NSEtaDPhi"; name += chargeName[icharge]; name += ic;
110  if (sedp) {
111  sedp[icharge]->SetName(name.Data());
112  sedp[icharge]->SetTitle(name.Data());
113  sedp[icharge]->Write();
114  }
115  TString name(binName[axialBins[ibin]]);
116  name += "_PtDEtaDPhi"; name += chargeName[icharge]; name += ic;
117  if (ptdedp) {
118  ptdedp[icharge]->SetName(name.Data());
119  ptdedp[icharge]->SetTitle(name.Data());
120  ptdedp[icharge]->Write();
121  }
122  TString name(binName[axialBins[ibin]]);
123  name += "_PtSEtaDPhi"; name += chargeName[icharge]; name += ic;
124  if (ptsedp) {
125  ptsedp[icharge]->SetName(name.Data());
126  ptsedp[icharge]->SetTitle(name.Data());
127  ptsedp[icharge]->Write();
128  }
129 
130  TString name(binName[axialBins[ibin]]);
131  name += "_NDEtaDPhi"; name += chargeType[icharge]; name += ic;
132  if (dedpC) {
133  dedpC[icharge]->SetName(name.Data());
134  dedpC[icharge]->SetTitle(name.Data());
135  dedpC[icharge]->Write();
136  }
137  TString name(binName[axialBins[ibin]]);
138  name += "_NSEtaDPhi"; name += chargeType[icharge]; name += ic;
139  if (sedpC) {
140  sedpC[icharge]->SetName(name.Data());
141  sedpC[icharge]->SetTitle(name.Data());
142  sedpC[icharge]->Write();
143  }
144  TString name(binName[axialBins[ibin]]);
145  name += "_PtDEtaDPhi"; name += chargeType[icharge]; name += ic;
146  if (ptdedpC) {
147  ptdedpC[icharge]->SetName(name.Data());
148  ptdedpC[icharge]->SetTitle(name.Data());
149  ptdedpC[icharge]->Write();
150  }
151  TString name(binName[axialBins[ibin]]);
152  name += "_PtSEtaDPhi"; name += chargeType[icharge]; name += ic;
153  if (ptsedpC) {
154  ptsedpC[icharge]->SetName(name.Data());
155  ptsedpC[icharge]->SetTitle(name.Data());
156  ptsedpC[icharge]->Write();
157  }
158 
159 
160  TString name(binName[axialBins[ibin]]);
161  name += "_NEtaEta"; name += chargeName[icharge]; name += ic;
162  if (etaeta) {
163  etaeta[icharge]->SetName(name.Data());
164  etaeta[icharge]->SetTitle(name.Data());
165  etaeta[icharge]->Write();
166  }
167  TString name(binName[axialBins[ibin]]);
168  name += "_PtEtaEta"; name += chargeName[icharge]; name += ic;
169  if (ptetaeta) {
170  ptetaeta[icharge]->SetName(name.Data());
171  ptetaeta[icharge]->SetTitle(name.Data());
172  ptetaeta[icharge]->Write();
173  }
174  TString name(binName[axialBins[ibin]]);
175  name += "_NEtaEtaSS"; name += chargeName[icharge]; name += ic;
176  if (etaetaSS) {
177  etaetaSS[icharge]->SetName(name.Data());
178  etaetaSS[icharge]->SetTitle(name.Data());
179  etaetaSS[icharge]->Write();
180  }
181  TString name(binName[axialBins[ibin]]);
182  name += "_PtEtaEtaSS"; name += chargeName[icharge]; name += ic;
183  if (ptetaetaSS) {
184  ptetaetaSS[icharge]->SetName(name.Data());
185  ptetaetaSS[icharge]->SetTitle(name.Data());
186  ptetaetaSS[icharge]->Write();
187  }
188  TString name(binName[axialBins[ibin]]);
189  name += "_NEtaEtaAS"; name += chargeName[icharge]; name += ic;
190  if (etaetaAS) {
191  etaetaAS[icharge]->SetName(name.Data());
192  etaetaAS[icharge]->SetTitle(name.Data());
193  etaetaAS[icharge]->Write();
194  }
195  TString name(binName[axialBins[ibin]]);
196  name += "_PtEtaEtaAS"; name += chargeName[icharge]; name += ic;
197  if (ptetaetaAS) {
198  ptetaetaAS[icharge]->SetName(name.Data());
199  ptetaetaAS[icharge]->SetTitle(name.Data());
200  ptetaetaAS[icharge]->Write();
201  }
202 
203  TString name(binName[axialBins[ibin]]);
204  name += "_NEtaEta"; name += chargeType[icharge]; name += ic;
205  if (etaetaC) {
206  etaetaC[icharge]->SetName(name.Data());
207  etaetaC[icharge]->SetTitle(name.Data());
208  etaetaC[icharge]->Write();
209  }
210  TString name(binName[axialBins[ibin]]);
211  name += "_PtEtaEta"; name += chargeType[icharge]; name += ic;
212  if (ptetaetaC) {
213  ptetaetaC[icharge]->SetName(name.Data());
214  ptetaetaC[icharge]->SetTitle(name.Data());
215  ptetaetaC[icharge]->Write();
216  }
217  TString name(binName[axialBins[ibin]]);
218  name += "_NEtaEtaSS"; name += chargeType[icharge]; name += ic;
219  if (etaetaSSC) {
220  etaetaSSC[icharge]->SetName(name.Data());
221  etaetaSSC[icharge]->SetTitle(name.Data());
222  etaetaSSC[icharge]->Write();
223  }
224  TString name(binName[axialBins[ibin]]);
225  name += "_PtEtaEtaSS"; name += chargeType[icharge]; name += ic;
226  if (ptetaetaSSC) {
227  ptetaetaSSC[icharge]->SetName(name.Data());
228  ptetaetaSSC[icharge]->SetTitle(name.Data());
229  ptetaetaSSC[icharge]->Write();
230  }
231  TString name(binName[axialBins[ibin]]);
232  name += "_NEtaEtaAS"; name += chargeType[icharge]; name += ic;
233  if (etaetaASC) {
234  etaetaASC[icharge]->SetName(name.Data());
235  etaetaASC[icharge]->SetTitle(name.Data());
236  etaetaASC[icharge]->Write();
237  }
238  TString name(binName[axialBins[ibin]]);
239  name += "_PtEtaEtaAS"; name += chargeType[icharge]; name += ic;
240  if (ptetaetaASC) {
241  ptetaetaASC[icharge]->SetName(name.Data());
242  ptetaetaASC[icharge]->SetTitle(name.Data());
243  ptetaetaASC[icharge]->Write();
244  }
245 
246 
247  TString name(binName[axialBins[ibin]]);
248  name += "_NPhiPhi"; name += chargeName[icharge]; name += ic;
249  if (phiphi) {
250  phiphi[icharge]->SetName(name.Data());
251  phiphi[icharge]->SetTitle(name.Data());
252  phiphi[icharge]->Write();
253  }
254  TString name(binName[axialBins[ibin]]);
255  name += "_PtPhiPhi"; name += chargeName[icharge]; name += ic;
256  if (ptphiphi) {
257  ptphiphi[icharge]->SetName(name.Data());
258  ptphiphi[icharge]->SetTitle(name.Data());
259  ptphiphi[icharge]->Write();
260  }
261 
262  TString name(binName[axialBins[ibin]]);
263  name += "_NPhiPhi"; name += chargeType[icharge]; name += ic;
264  if (phiphiC) {
265  phiphiC[icharge]->SetName(name.Data());
266  phiphiC[icharge]->SetTitle(name.Data());
267  phiphiC[icharge]->Write();
268  }
269  TString name(binName[axialBins[ibin]]);
270  name += "_PtPhiPhi"; name += chargeType[icharge]; name += ic;
271  if (ptphiphiC) {
272  ptphiphiC[icharge]->SetName(name.Data());
273  ptphiphiC[icharge]->SetTitle(name.Data());
274  ptphiphiC[icharge]->Write();
275  }
276  }
277  // Calculate and print scale factors
278  double scale = ehelp->getCIdNdEtadPhi();
279  double *ptHat = ehelp->getptHat();
280  printf(" bin %8s %7.3f %7.3f %7.3f %7.3f %7.3f\n",binName[axialBins[ibin]],scale,ptHat[0],ptHat[1],ptHat[2],ptHat[3]);
281  delete [] ptHat;
282 
283  delete tf;
284  delete ehelp;
285  }
286  }
287 
288  TH2D **ytyt = 0;
289  TH2D **ytytC = 0;
290  TH2D **sytdyt = 0;
291  TH2D **sytdytC = 0;
292  gROOT->LoadMacro("minimizeNegative.C");
293  double *sFactor[2][3], *eSFactor[2][3];
294  for (int it=0;it<2;it++) {
295  for (int ibin=0;ibin<3;ibin++) {
296  sFactor[it][ibin] = new double[nCent];
297  eSFactor[it][ibin] = new double[nCent];
298  }
299  }
300  float sf[2];
301  bool scaleYt = false;
302  double argList[10];
303  double start = 0.95;
304  double step = 0.001;
305  double bmin = 0.9;
306  double bmax = 1.1;
307  int errFlag = 0;
308  minData.mCorrType = 1;
309  minData.mLambda = 10;
310 
311  // For ytyt space soft/neck/hard is not useful.
312  int ytytBins[] = {0, 1, 2};
313  for (int ic=0;ic<nCent;ic++) {
314  for (int ibin=0;ibin<3;ibin++) {
315  sprintf(inFileName,"%s/%s%s.root",dirName,inNames[ic],binName[ytytBins[ibin]]);
316  tf = new TFile(inFileName);
317  tf->cd();
318  ehelp = new StEStructSupport(tf,0);
319  ehelp->msilent = true;
320  ehelp->mapplyDEtaFix = false;
321  ehelp->mPairNormalization = false;
322  ehelp->mIdenticalPair = true;
323  ehelp->mYtYtNormalization = true;
324 
325  // Make sure LS histogram for zBin 0 is in file.
326  if (0 == ehelp->histogramExists("YtYt", 0)) {
327  continue;
328  }
329 
330  // A lot of stuff so we can find a scaling factor for \rho_{ref}
331  // such that \Delta\rho is almost always positive.
332  minuit = new TMinuit(1);
333  if (!scaleYt) {
334  ehelp->setBGMode(0);
335  ehelp->mYtYtNormalization = true;
336  sf[0] = 1;
337  sf[1] = 1;
338  cout << ">>>>>Not fitting scale factors for centrality " << ic << " yt bin " << ibin << ", Using YtYtNormalization = " << ehelp->mYtYtNormalization << endl;
339  ytyt = ehelp->buildChargeTypes("YtYt",5,sf);
340  ytytC = ehelp->buildCommon("YtYt",5,sf);
341  sytdyt = ehelp->buildChargeTypes("SYtDYt",5,sf);
342  sytdytC = ehelp->buildCommon("SYtDYt",5,sf);
343  } else {
344  ehelp->setBGMode(1);
345  minData.mSupport = ehelp;
346  minData.mChargeType = 0;
347 
348  argList[0] = -1;
349  minuit->mnexcm("SET PRINT",argList,1,errFlag);
350  minuit->SetFCN(minimizeNegative);
351  minuit->mnparm( 0, "rho_ref scale factor", start, step, bmin, bmax, errFlag );
352  minuit->SetErrorDef(1);
353  minuit->SetPrintLevel(-1);
354  argList[0] = 1;
355  minuit->mnexcm("SET STR",argList,1,errFlag);
356  argList[0] = 500;
357  cout << ">>>>>Starting scale factor 0 fit for centrality " << ic << " yt bin " << ibin << endl;
358  minuit->mnexcm("MIGRAD",argList,1,errFlag);
359  double cVal, eVal;
360  minuit->GetParameter(0,cVal,eVal);
361  sFactor[0][ibin][ic] = cVal;
362  eSFactor[0][ibin][ic] = eVal;
363  if (4 == errFlag) {
364  minuit->mnexcm("MINOS",argList,1,errFlag);
365  minuit->GetParameter(0,cVal,eVal);
366  sFactor[0][ibin][ic] = cVal;
367  eSFactor[0][ibin][ic] = eVal;
368  if (4 == errFlag) {
369  cout << "++++Out of chesse error; Fit failed. Using scale factor = 1" << endl;
370  sFactor[0][ibin][ic] = 1;
371  eSFactor[0][ibin][ic] = 100;
372  }
373  }
374  delete minuit;
375 
376  // Seems like I should be able to reset the TMinuit object to do a
377  // different fit. Didn't work on my first attempts, so I just create
378  // a new one.
379  minData.mChargeType = 1;
380  minuit = new TMinuit(1);
381  argList[0] = -1;
382  minuit->mnexcm("SET PRINT",argList,1,errFlag);
383  minuit->SetFCN(minimizeNegative);
384  minuit->mnparm( 0, "rho_ref scale factor", start, step, bmin, bmax, errFlag );
385  minuit->SetErrorDef(1);
386  minuit->SetPrintLevel(-1);
387  argList[0] = 1;
388  minuit->mnexcm("SET STR",argList,1,errFlag);
389  argList[0] = 500;
390  cout << ">>>>>Starting scale factor 1 fit for centrality " << ic << " yt bin " << ibin << endl;
391  minuit->mnexcm("MIGRAD",argList,1,errFlag);
392  minuit->GetParameter(0,cVal,eVal);
393  sFactor[1][ibin][ic] = cVal;
394  eSFactor[1][ibin][ic] = eVal;
395  if (4 == errFlag) {
396  minuit->mnexcm("MINOS",argList,1,errFlag);
397  minuit->GetParameter(0,cVal,eVal);
398  sFactor[1][ibin][ic] = cVal;
399  eSFactor[1][ibin][ic] = eVal;
400  if (4 == errFlag) {
401  cout << "++++Out of chesse error; Fit failed. Using scale factor = 1" << endl;
402  sFactor[1][ibin][ic] = 1;
403  eSFactor[1][ibin][ic] = 100;
404  }
405  }
406 
407  sf[0] = sFactor[0][ibin][ic];
408  sf[1] = sFactor[1][ibin][ic];
409  ytyt = ehelp->buildChargeTypes("YtYt",5,sf);
410  ytytC = ehelp->buildCommon("YtYt",5,sf);
411  sytdyt = ehelp->buildChargeTypes("SYtDYt",5,sf);
412  sytdytC = ehelp->buildCommon("SYtDYt",5,sf);
413  }
414  delete minuit;
415 
416  out->cd();
417  for (int icharge=0;icharge<4;icharge++) {
418  if (ytyt) {
419  TString name(binName[ytytBins[ibin]]);
420  name += "_YtYt"; name += chargeName[icharge]; name += ic;
421  ytyt[icharge]->SetName(name.Data());
422  ytyt[icharge]->SetTitle(name.Data());
423  ytyt[icharge]->Write();
424  }
425  if (ytytC) {
426  TString name(binName[ytytBins[ibin]]);
427  name += "_YtYt"; name += chargeType[icharge]; name += ic;
428  ytytC[icharge]->SetName(name.Data());
429  ytytC[icharge]->SetTitle(name.Data());
430  ytytC[icharge]->Write();
431  }
432 
433  if (sytdyt) {
434  TString name(binName[ytytBins[ibin]]);
435  name += "_SYtDYt"; name += chargeName[icharge]; name += ic;
436  sytdyt[icharge]->SetName(name.Data());
437  sytdyt[icharge]->SetTitle(name.Data());
438  sytdyt[icharge]->Write();
439  }
440  if (sytdytC) {
441  TString name(binName[ytytBins[ibin]]);
442  name += "_SYtDYt"; name += chargeType[icharge]; name += ic;
443  sytdytC[icharge]->SetName(name.Data());
444  sytdytC[icharge]->SetTitle(name.Data());
445  sytdytC[icharge]->Write();
446  }
447  }
448  delete tf;
449  delete ehelp;
450  }
451  }
452  if (scaleYt) {
453  cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << endl;
454  cout << " Here are Yt scale factors for rho_{ref}" << endl;
455  printf(" Centrality all LS SS LS AS LS all US SS US AS US \n");
456  for (int ic=0;ic<nCent;ic++) {
457  printf("%6i %9.3f(%3.0f) %7.3f(%3.0f) %7.3f(%3.0f) %7.3f(%3.0f) %7.3f(%3.0f) %7.3f(%3.0f)\n",ic,
458  sFactor[0][0][ic],1000*eSFactor[0][0][ic],sFactor[0][1][ic],1000*eSFactor[0][1][ic],sFactor[0][2][ic],1000*eSFactor[0][2][ic],
459  sFactor[1][0][ic],1000*eSFactor[1][0][ic],sFactor[1][1][ic],1000*eSFactor[1][1][ic],sFactor[1][2][ic],1000*eSFactor[1][2][ic]);
460  }
461  for (int it=0;it<2;it++) {
462  for (int ibin=0;ibin<3;ibin++) {
463  delete sFactor[it][ibin];
464  delete eSFactor[it][ibin];
465  }
466  }
467  }
468  out->Close();
469 // delete out;
470 }
Generated by   doxygen 1.8.5