StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
combineHistograms5.C
1 void combineHistograms5(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 combineHistograms5.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 **ptdedp;
22  TH2D **ptdedpC;
23  TH2D **dedp;
24  TH2D **dedpC;
25  TH2D **ptsedp;
26  TH2D **ptsedpC;
27  TH2D **sedp;
28  TH2D **sedpC;
29  TH2D **ptetaeta;
30  TH2D **ptetaetaC;
31  TH2D **etaeta;
32  TH2D **etaetaC;
33  TH2D **ptphiphi;
34  TH2D **ptphiphiC;
35  TH2D **phiphi;
36  TH2D **phiphiC;
37 
38  const char* pidName[] = {"all", "pipi", "piK", "pip", "KK", "Kp", "pp", "oo"};
39  const char* chargeName[] = {"_LS_", "_US_", "_CD_", "_CI_"};
40  const char* chargeType[] = {"_PP_", "_PM_", "_MP_", "_MM_"};
41  char outFileName[1024];
42  sprintf(outFileName,"%s/%s.root",dirName,outName);
43  TFile *out = new TFile(outFileName,"RECREATE");
44 
45  const char* oname[]={"all","pi_pi","pi_K","pi_p","K_K","K_p","p_p","o_o"};
46  char inFileName[1024];
47 
48  for (int ic=0;ic<nCent;ic++) {
49  printf("Centrality %2i: d2N/dEdP pHat_{A+} pHat_{A-} pHat_{B+} pHat_{B-}\n",ic);
50  for (int ipid=0;ipid<8;ipid++) {
51  sprintf(inFileName,"%s/%s%s.root",dirName,inNames[ic],oname[ipid]);
52  tf = new TFile(inFileName);
53  tf->cd();
54  ehelp = new StEStructSupport(tf,0);
55  ehelp->msilent = true;
56  ehelp->mapplyDEtaFix = false;
57  if ((0 ==ipid) || (1 == ipid) || (4 == ipid) || (6 == ipid) || (7 == ipid)) {
58  ehelp->mIdenticalPair = true;
59  } else {
60  ehelp->mIdenticalPair = false;
61  }
62  ehelp->setBGMode(0);
63  int subtract = 1;
64  ehelp->mPairNormalization = false;
65  ptdedpC = ehelp->buildPtCommon("DEtaDPhi",2,subtract);
66  ptsedpC = ehelp->buildPtCommon("SEtaDPhi",2,subtract);
67  ptetaetaC = ehelp->buildPtCommon("EtaEta",2,subtract);
68  ptphiphiC = ehelp->buildPtCommon("PhiPhi",2,subtract);
69 
70  ptdedp = ehelp->buildPtChargeTypes("DEtaDPhi",2,subtract);
71  ptsedp = ehelp->buildPtChargeTypes("SEtaDPhi",2,subtract);
72  ptetaeta = ehelp->buildPtChargeTypes("EtaEta",2,subtract);
73  ptphiphi = ehelp->buildPtChargeTypes("PhiPhi",2,subtract);
74 
75  ehelp->mPairNormalization = true;
76  dedpC = ehelp->buildCommon("DEtaDPhi",2);
77  sedpC = ehelp->buildCommon("SEtaDPhi",2);
78  etaetaC = ehelp->buildCommon("EtaEta",2);
79  phiphiC = ehelp->buildCommon("PhiPhi",2);
80 
81  dedp = ehelp->buildChargeTypes("DEtaDPhi",2);
82  sedp = ehelp->buildChargeTypes("SEtaDPhi",2);
83  etaeta = ehelp->buildChargeTypes("EtaEta",2);
84  phiphi = ehelp->buildChargeTypes("PhiPhi",2);
85 
86  out->cd();
87  for (int icharge=0;icharge<4;icharge++) {
88  TString name(pidName[ipid]);
89  name += "_NDEtaDPhi"; name += chargeName[icharge]; name += ic;
90  if (dedp) {
91  dedp[icharge]->SetName(name.Data());
92  dedp[icharge]->SetTitle(name.Data());
93  dedp[icharge]->Write();
94  }
95  TString name(pidName[ipid]);
96  name += "_PtDEtaDPhi"; name += chargeName[icharge]; name += ic;
97  if (ptdedp) {
98  ptdedp[icharge]->SetName(name.Data());
99  ptdedp[icharge]->SetTitle(name.Data());
100  ptdedp[icharge]->Write();
101  }
102  TString name(pidName[ipid]);
103  name += "_NSEtaDPhi"; name += chargeName[icharge]; name += ic;
104  if (sedp) {
105  sedp[icharge]->SetName(name.Data());
106  sedp[icharge]->SetTitle(name.Data());
107  sedp[icharge]->Write();
108  }
109  TString name(pidName[ipid]);
110  name += "_PtSEtaDPhi"; name += chargeName[icharge]; name += ic;
111  if (ptsedp) {
112  ptsedp[icharge]->SetName(name.Data());
113  ptsedp[icharge]->SetTitle(name.Data());
114  ptsedp[icharge]->Write();
115  }
116 
117  TString name(pidName[ipid]);
118  name += "_NDEtaDPhi"; name += chargeType[icharge]; name += ic;
119  if (dedpC) {
120  dedpC[icharge]->SetName(name.Data());
121  dedpC[icharge]->SetTitle(name.Data());
122  dedpC[icharge]->Write();
123  }
124  TString name(pidName[ipid]);
125  name += "_PtDEtaDPhi"; name += chargeType[icharge]; name += ic;
126  if (ptdedpC) {
127  ptdedpC[icharge]->SetName(name.Data());
128  ptdedpC[icharge]->SetTitle(name.Data());
129  ptdedpC[icharge]->Write();
130  }
131  TString name(pidName[ipid]);
132  name += "_NSEtaDPhi"; name += chargeType[icharge]; name += ic;
133  if (sedpC) {
134  sedpC[icharge]->SetName(name.Data());
135  sedpC[icharge]->SetTitle(name.Data());
136  sedpC[icharge]->Write();
137  }
138  TString name(pidName[ipid]);
139  name += "_PtSEtaDPhi"; name += chargeType[icharge]; name += ic;
140  if (ptsedpC) {
141  ptsedpC[icharge]->SetName(name.Data());
142  ptsedpC[icharge]->SetTitle(name.Data());
143  ptsedpC[icharge]->Write();
144  }
145 
146 
147  TString name(pidName[ipid]);
148  name += "_NEtaEta"; name += chargeName[icharge]; name += ic;
149  if (etaeta) {
150  etaeta[icharge]->SetName(name.Data());
151  etaeta[icharge]->SetTitle(name.Data());
152  etaeta[icharge]->Write();
153  }
154  TString name(pidName[ipid]);
155  name += "_PtEtaEta"; name += chargeName[icharge]; name += ic;
156  if (ptetaeta) {
157  ptetaeta[icharge]->SetName(name.Data());
158  ptetaeta[icharge]->SetTitle(name.Data());
159  ptetaeta[icharge]->Write();
160  }
161 
162  TString name(pidName[ipid]);
163  name += "_NEtaEta"; name += chargeType[icharge]; name += ic;
164  if (etaetaC) {
165  etaetaC[icharge]->SetName(name.Data());
166  etaetaC[icharge]->SetTitle(name.Data());
167  etaetaC[icharge]->Write();
168  }
169  TString name(pidName[ipid]);
170  name += "_PtEtaEta"; name += chargeType[icharge]; name += ic;
171  if (ptetaetaC) {
172  ptetaetaC[icharge]->SetName(name.Data());
173  ptetaetaC[icharge]->SetTitle(name.Data());
174  ptetaetaC[icharge]->Write();
175  }
176 
177 
178  TString name(pidName[ipid]);
179  name += "_NPhiPhi"; name += chargeName[icharge]; name += ic;
180  if (phiphi) {
181  phiphi[icharge]->SetName(name.Data());
182  phiphi[icharge]->SetTitle(name.Data());
183  phiphi[icharge]->Write();
184  }
185  TString name(pidName[ipid]);
186  name += "_PtPhiPhi"; name += chargeName[icharge]; name += ic;
187  if (ptphiphi) {
188  ptphiphi[icharge]->SetName(name.Data());
189  ptphiphi[icharge]->SetTitle(name.Data());
190  ptphiphi[icharge]->Write();
191  }
192 
193  TString name(pidName[ipid]);
194  name += "_NPhiPhi"; name += chargeType[icharge]; name += ic;
195  if (phiphiC) {
196  phiphiC[icharge]->SetName(name.Data());
197  phiphiC[icharge]->SetTitle(name.Data());
198  phiphiC[icharge]->Write();
199  }
200  TString name(pidName[ipid]);
201  name += "_PtPhiPhi"; name += chargeType[icharge]; name += ic;
202  if (ptphiphiC) {
203  ptphiphiC[icharge]->SetName(name.Data());
204  ptphiphiC[icharge]->SetTitle(name.Data());
205  ptphiphiC[icharge]->Write();
206  }
207  }
208  // Calculate and print scale factors
209  double scale = ehelp->getCIdNdEtadPhi();
210  double *ptHat = ehelp->getptHat();
211  printf(" bin %8s %7.3f %7.3f %7.3f %7.3f %7.3f\n",pidName[ipid],scale,ptHat[0],ptHat[1],ptHat[2],ptHat[3]);
212  delete [] ptHat;
213 
214  delete tf;
215  delete ehelp;
216  }
217  }
218 
219  TH2D **ytyt;
220  TH2D **ytytC;
221  TH2D **sytdyt;
222  TH2D **sytdytC;
223  gROOT->LoadMacro("minimizeNegative2.C");
224  double *sFactor[2][8];
225  for (int it=0;it<2;it++) {
226  for (int ipid=0;ipid<8;ipid++) {
227  sFactor[it][ipid] = new double[nCent];
228  }
229  }
230  float sf[2];
231  bool scaleYt = false;
232  double min = 0.85;
233  double max = 1.05;
234  double D = max - min;
235  int nBins = 51;
236  double d = D/(nBins-1.0);
237  int npar = 1;
238  double gin = 10;
239  double f;
240  double par;
241  int iflag = 1;
242  double aH, bH, aL, bL;
243  TH1F *hFunc = new TH1F("hFunc","hFunc",nBins,min-d/2,max+d/2);
244  TF1 *fHigh = new TF1("HighFit","[0]+[1]*x",1.01,1.05);
245  TF1 *fLow = new TF1("LowFit","[0]+[1]*x",0.85,0.95);
246 
247  minData.mCorrType = 1;
248  minData.mLambda = 10;
249 
250  // For ytyt space pid is not yet useful because of limited range due to dE/dx.
251  // Now we have ToF. Check if yt range is big enough that yty space is interesting.
252  for (int ic=0;ic<nCent;ic++) {
253  for (int ipid=0;ipid<8;ipid++) {
254  sprintf(inFileName,"%s/%s%s.root",dirName,inNames[ic],oname[ipid]);
255  tf = new TFile(inFileName);
256  tf->cd();
257  ehelp = new StEStructSupport(tf,0);
258  ehelp->msilent = true;
259  ehelp->mapplyDEtaFix = false;
260  ehelp->mPairNormalization = false;
261  ehelp->mYtYtNormalization = true;
262  if ((0 ==ipid) || (1 == ipid) || (4 == ipid) || (6 == ipid) || (7 == ipid)) {
263  ehelp->mIdenticalPair = true;
264  } else {
265  ehelp->mIdenticalPair = false;
266  }
267 
268  // Make sure LS histogram for zBin 0 is in file.
269  if (0 == ehelp->histogramExists("YtYt", 0)) {
270  continue;
271  }
272 
273  if (!scaleYt) {
274  ehelp->setBGMode(0);
275  ehelp->mYtYtVolumeNormalization = true;
276  ehelp->mYtYtNormalization = false;
277  ehelp->mPairNormalization = false;
278  sf[0] = 1;
279  sf[1] = 1;
280  cout << ">>>>>Not fitting scale factors for centrality " << ic << " pid bin " << ipid << ", Using YtYtNormalization = " << ehelp->mYtYtNormalization << endl;
281  ytyt = ehelp->buildChargeTypes("YtYt",5,sf);
282  ytytC = ehelp->buildCommon("YtYt",5,sf);
283  sytdyt = ehelp->buildChargeTypes("SYtDYt",5,sf);
284  sytdytC = ehelp->buildCommon("SYtDYt",5,sf);
285  } else {
286  // Here we try minimizing volume of \Delta\rho/\sqrt(\rho_{ref})
287  // constraining the quantity to be mostly positive.
288  ehelp->setBGMode(1);
289  minData.mSupport = ehelp;
290  minData.mChargeType = 0;
291  for (int ia=1;ia<=nBins;ia++) {
292  par = min + float(ia-1)*d;
293  minimizeNegative2(npar,&gin,f,&par,iflag);
294  hFunc->SetBinContent(ia,f);
295  }
296  hFunc->Fit("HighFit","NORQ");
297  aH = fHigh->GetParameter(0);
298  bH = fHigh->GetParameter(1);
299  hFunc->Fit("LowFit","NORQ");
300  aL = fLow->GetParameter(0);
301  bL = fLow->GetParameter(1);
302  if (bL > 0 || bH < 0) {
303  cout << "Suspect fit for ic = " << ic << ", ipid = " << ipid << ", LS " << endl;
304  }
305  sFactor[0][ipid][ic] = -(aH-aL)/(bH-bL);
306 
307  minData.mChargeType = 1;
308  for (int ia=1;ia<=nBins;ia++) {
309  par = min + float(ia-1)*d;
310  minimizeNegative2(npar,&gin,f,&par,iflag);
311  hFunc->SetBinContent(ia,f);
312  }
313  hFunc->Fit("HighFit","NORQ");
314  aH = fHigh->GetParameter(0);
315  bH = fHigh->GetParameter(1);
316  hFunc->Fit("LowFit","NORQ");
317  aL = fLow->GetParameter(0);
318  bL = fLow->GetParameter(1);
319  if (bL > 0 || bH < 0) {
320  cout << "Suspect fit for ic = " << ic << ", ipid = " << ipid << ", US " << endl;
321  }
322  sFactor[1][ipid][ic] = -(aH-aL)/(bH-bL);
323 
324 
325  sf[0] = sFactor[0][ipid][ic];
326  sf[1] = sFactor[1][ipid][ic];
327  ytyt = ehelp->buildChargeTypes("YtYt",5,sf);
328  ytytC = ehelp->buildCommon("YtYt",5,sf);
329  sytdyt = ehelp->buildChargeTypes("SYtDYt",5,sf);
330  sytdytC = ehelp->buildCommon("SYtDYt",5,sf);
331  }
332 
333  out->cd();
334  for (int icharge=0;icharge<4;icharge++) {
335  TString name(pidName[ipid]);
336  name += "_YtYt"; name += chargeName[icharge]; name += ic;
337  ytyt[icharge]->SetName(name.Data());
338  ytyt[icharge]->SetTitle(name.Data());
339  ytyt[icharge]->Write();
340  TString name(pidName[ipid]);
341  name += "_YtYt"; name += chargeType[icharge]; name += ic;
342  ytytC[icharge]->SetName(name.Data());
343  ytytC[icharge]->SetTitle(name.Data());
344  ytytC[icharge]->Write();
345 
346  TString name(pidName[ipid]);
347  name += "_SYtDYt"; name += chargeName[icharge]; name += ic;
348  sytdyt[icharge]->SetName(name.Data());
349  sytdyt[icharge]->SetTitle(name.Data());
350  sytdyt[icharge]->Write();
351  TString name(pidName[ipid]);
352  name += "_SYtDYt"; name += chargeType[icharge]; name += ic;
353  sytdytC[icharge]->SetName(name.Data());
354  sytdytC[icharge]->SetTitle(name.Data());
355  sytdytC[icharge]->Write();
356  }
357  delete tf;
358  delete ehelp;
359  }
360  }
361  if (scaleYt) {
362  cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << endl;
363  cout << " Here are Yt scale factors for rho_{ref}" << endl;
364  printf(" Centrality all pi_pi pi_K pi_p K_K K_p p_p o_o\n");
365  for (int ic=0;ic<nCent;ic++) {
366  printf("%4i US %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n",ic,
367  sFactor[0][0][ic],sFactor[0][1][ic],
368  sFactor[0][2][ic],sFactor[0][3][ic],
369  sFactor[0][4][ic],sFactor[0][5][ic],
370  sFactor[0][6][ic],sFactor[0][7][ic]);
371  printf("%4i LS %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n",ic,
372  sFactor[1][0][ic],sFactor[1][1][ic],
373  sFactor[1][2][ic],sFactor[1][3][ic],
374  sFactor[1][4][ic],sFactor[1][5][ic],
375  sFactor[1][6][ic],sFactor[1][7][ic]);
376  }
377  }
378  out->Close();
379 // delete out;
380  for (int it=0;it<2;it++) {
381  for (int ipid=0;ipid<8;ipid++) {
382  delete sFactor[it][ipid];
383  }
384  }
385 }
Generated by   doxygen 1.8.5