StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
combineHistograms8.C
1 void combineHistograms8(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  TH2F **ptsedp;
22  TH2F **ptsedpC;
23  TH2F **ptdedp;
24  TH2F **ptdedpC;
25  TH2F **sedp;
26  TH2F **sedpC;
27  TH2F **dedp;
28  TH2F **dedpC;
29  TH2F **ptetaeta;
30  TH2F **ptetaetaC;
31  TH2F **etaeta;
32  TH2F **etaetaC;
33  TH2F **ptphiphi;
34  TH2F **ptphiphiC;
35  TH2F **phiphi;
36  TH2F **phiphiC;
37 
38  const char* binName[]={"all","soft","neck","hard","softNeck","softHard","neckHard"};
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  char inFileName[1024];
46 
47  for (int ic=0;ic<nCent;ic++) {
48  printf("Scale factors for centrality %2i, ++ +- -+ -- CI \n",ic);
49  for (int ibin=0;ibin<7;ibin++) {
50  sprintf(inFileName,"%s/%s%s.root",dirName,inNames[ic],binName[ibin]);
51  tf = new TFile(inFileName);
52  tf->cd();
53  ehelp = new StEStructSupport(tf,0);
54  ehelp->msilent = true;
55  ehelp->mapplyDEtaFix = false;
56  ehelp->mPairNormalization = true;
57  ehelp->mIdenticalPair = true;
58  ehelp->setBGMode(0);
59  int subtract = 1;
60  ptsedpC = (TH2F**) ehelp->buildPtCommon("SEtaDPhi",2,subtract);
61  ptdedpC = (TH2F**) ehelp->buildPtCommon("DEtaDPhi",2,subtract);
62  ptetaetaC = (TH2F**) ehelp->buildPtCommon("EtaEta",2,subtract);
63  ptphiphiC = (TH2F**) ehelp->buildPtCommon("PhiPhi",2,subtract);
64 
65  ptsedp = (TH2F**) ehelp->buildPtChargeTypes("SEtaDPhi",2,subtract);
66  ptdedp = (TH2F**) ehelp->buildPtChargeTypes("DEtaDPhi",2,subtract);
67  ptetaeta = (TH2F**) ehelp->buildPtChargeTypes("EtaEta",2,subtract);
68  ptphiphi = (TH2F**) ehelp->buildPtChargeTypes("PhiPhi",2,subtract);
69 
70  sedpC = (TH2F**) ehelp->buildCommon("SEtaDPhi",2);
71  dedpC = (TH2F**) ehelp->buildCommon("DEtaDPhi",2);
72  etaetaC = (TH2F**) ehelp->buildCommon("EtaEta",2);
73  phiphiC = (TH2F**) ehelp->buildCommon("PhiPhi",2);
74 
75  sedp = (TH2F**) ehelp->buildChargeTypes("SEtaDPhi",2);
76  dedp = (TH2F**) ehelp->buildChargeTypes("DEtaDPhi",2);
77  etaeta = (TH2F**) ehelp->buildChargeTypes("EtaEta",2);
78  phiphi = (TH2F**) ehelp->buildChargeTypes("PhiPhi",2);
79 
80  out->cd();
81  for (int icharge=0;icharge<4;icharge++) {
82  TString name(binName[ibin]);
83  name += "_NDEtaDPhi"; name += chargeName[icharge]; name += ic;
84  if (dedp) {
85  dedp[icharge]->SetName(name.Data());
86  dedp[icharge]->SetTitle(name.Data());
87  dedp[icharge]->Write();
88  }
89  TString name(binName[ibin]);
90  name += "_NSEtaDPhi"; name += chargeName[icharge]; name += ic;
91  if (sedp) {
92  sedp[icharge]->SetName(name.Data());
93  sedp[icharge]->SetTitle(name.Data());
94  sedp[icharge]->Write();
95  }
96  TString name(binName[ibin]);
97  name += "_PtDEtaDPhi"; name += chargeName[icharge]; name += ic;
98  if (ptdedp) {
99  ptdedp[icharge]->SetName(name.Data());
100  ptdedp[icharge]->SetTitle(name.Data());
101  ptdedp[icharge]->Write();
102  }
103  TString name(binName[ibin]);
104  name += "_PtSEtaDPhi"; name += chargeName[icharge]; name += ic;
105  if (ptsedp) {
106  ptsedp[icharge]->SetName(name.Data());
107  ptsedp[icharge]->SetTitle(name.Data());
108  ptsedp[icharge]->Write();
109  }
110 
111  TString name(binName[ibin]);
112  name += "_NDEtaDPhi"; name += chargeType[icharge]; name += ic;
113  if (dedpC) {
114  dedpC[icharge]->SetName(name.Data());
115  dedpC[icharge]->SetTitle(name.Data());
116  dedpC[icharge]->Write();
117  }
118  TString name(binName[ibin]);
119  name += "_NSEtaDPhi"; name += chargeType[icharge]; name += ic;
120  if (sedpC) {
121  sedpC[icharge]->SetName(name.Data());
122  sedpC[icharge]->SetTitle(name.Data());
123  sedpC[icharge]->Write();
124  }
125  TString name(binName[ibin]);
126  name += "_PtEtaDPhi"; name += chargeType[icharge]; name += ic;
127  if (ptdedpC) {
128  ptdedpC[icharge]->SetName(name.Data());
129  ptdedpC[icharge]->SetTitle(name.Data());
130  ptdedpC[icharge]->Write();
131  }
132  TString name(binName[ibin]);
133  if (name += ) {
134  name += "_PtEtaSPhi"; name += chargeType[icharge]; name += ic;
135  ptsedpC[icharge]->SetName(name.Data());
136  ptsedpC[icharge]->SetTitle(name.Data());
137  }
138  ptsedpC[icharge]->Write();
139 
140 
141  TString name(binName[ibin]);
142  name += "_NEtaEta"; name += chargeName[icharge]; name += ic;
143  if (etaeta) {
144  etaeta[icharge]->SetName(name.Data());
145  etaeta[icharge]->SetTitle(name.Data());
146  etaeta[icharge]->Write();
147  }
148  TString name(binName[ibin]);
149  name += "_PtEtaEta"; name += chargeName[icharge]; name += ic;
150  if (ptetaeta) {
151  ptetaeta[icharge]->SetName(name.Data());
152  ptetaeta[icharge]->SetTitle(name.Data());
153  ptetaeta[icharge]->Write();
154  }
155 
156  TString name(binName[ibin]);
157  name += "_NEtaEta"; name += chargeType[icharge]; name += ic;
158  if (etaetaC) {
159  etaetaC[icharge]->SetName(name.Data());
160  etaetaC[icharge]->SetTitle(name.Data());
161  etaetaC[icharge]->Write();
162  }
163  TString name(binName[ibin]);
164  name += "_PtEtaEta"; name += chargeType[icharge]; name += ic;
165  if (ptetaetaC) {
166  ptetaetaC[icharge]->SetName(name.Data());
167  ptetaetaC[icharge]->SetTitle(name.Data());
168  ptetaetaC[icharge]->Write();
169  }
170 
171 
172  TString name(binName[ibin]);
173  name += "_NPhiPhi"; name += chargeName[icharge]; name += ic;
174  if (phiphi) {
175  phiphi[icharge]->SetName(name.Data());
176  phiphi[icharge]->SetTitle(name.Data());
177  phiphi[icharge]->Write();
178  }
179  TString name(binName[ibin]);
180  name += "_PtPhiPhi"; name += chargeName[icharge]; name += ic;
181  if (ptphiphi) {
182  ptphiphi[icharge]->SetName(name.Data());
183  ptphiphi[icharge]->SetTitle(name.Data());
184  ptphiphi[icharge]->Write();
185  }
186 
187  TString name(binName[ibin]);
188  name += "_NPhiPhi"; name += chargeType[icharge]; name += ic;
189  if (phiphiC) {
190  phiphiC[icharge]->SetName(name.Data());
191  phiphiC[icharge]->SetTitle(name.Data());
192  phiphiC[icharge]->Write();
193  }
194  TString name(binName[ibin]);
195  name += "_PtPhiPhi"; name += chargeType[icharge]; name += ic;
196  if (ptphiphiC) {
197  ptphiphiC[icharge]->SetName(name.Data());
198  ptphiphiC[icharge]->SetTitle(name.Data());
199  ptphiphiC[icharge]->Write();
200  }
201  }
202  // Calculate and print scale factors
203  double *scale = ehelp->getScaleFactors();
204  printf(" bin %8s %7.2f %7.2f %7.2f %7.2f %7.2f\n",binName[ibin],scale[0],scale[1],scale[2],scale[3],scale[4]);
205  delete [] scale;
206  delete tf;
207  delete ehelp;
208  }
209  }
210 
211  TH2D **ytyt;
212  TH2D **ytytC;
213  TH2D **sytdyt;
214  TH2D **sytdytC;
215  gROOT->LoadMacro("minimizeNegative.C");
216  double *sFactor[2][3], *eSFactor[2][3];
217  for (int it=0;it<2;it++) {
218  for (int ibin=0;ibin<3;ibin++) {
219  sFactor[it][ibin] = new double[nCent];
220  eSFactor[it][ibin] = new double[nCent];
221  }
222  }
223  float sf[2];
224  double argList[10];
225  double start = 0.95;
226  double step = 0.001;
227  double bmin = 0.9;
228  double bmax = 1.1;
229  int errFlag = 0;
230  minData.mCorrType = 1;
231  minData.mLambda = 10;
232 
233  // For ytyt space soft/neck/hard is not useful.
234  int ytytBins[] = {0};
235  for (int ic=0;ic<nCent;ic++) {
236  for (int ibin=0;ibin<1;ibin++) {
237  sprintf(inFileName,"%s/%s%s.root",dirName,inNames[ic],binName[ytytBins[ibin]]);
238  tf = new TFile(inFileName);
239  tf->cd();
240  ehelp = new StEStructSupport(tf,0);
241  ehelp->msilent = true;
242  ehelp->mapplyDEtaFix = false;
243  ehelp->mPairNormalization = false;
244  ehelp->mIdenticalPair = true;
245  ehelp->setBGMode(1);
246 
247  // Make sure LS histogram for zBin 0 is in file.
248  if (0 == ehelp->histogramExists("YtYt", 0)) {
249  continue;
250  }
251 
252  // A lot of stuff so we can find a scaling factor for \rho_{ref}
253  // such that \Delta\rho is almost always positive.
254  minData.mSupport = ehelp;
255  minData.mChargeType = 0;
256 
257  minuit = new TMinuit(1);
258  minuit->SetFCN(minimizeNegative);
259  minuit->mnparm( 0, "rho_ref scale factor", start, step, bmin, bmax, errFlag );
260  minuit->SetErrorDef(1);
261  minuit->SetPrintLevel(-1);
262  argList[0] = 1;
263  minuit->mnexcm("SET STR",argList,1,errFlag);
264  argList[0] = 500;
265  cout << ">>>>>Starting scale factor 0 fit for centrality " << ic << " yt bin " << ibin << endl;
266  minuit->mnexcm("MIGRAD",argList,1,errFlag);
267  minuit->GetParameter(0,sFactor[0][ibin][ic],eSFactor[0][ibin][ic]);
268  if (0 != errFlag) {
269  cout << "++++Out of chesse error; Fit failed. Using scale factor = 1" << endl;
270  sFactor[0][ibin][ic] = 1;
271  eSFactor[0][ibin][ic] = 100;
272  }
273  delete minuit;
274 
275  // Seems like I should be able to reset the TMinuit object to do a
276  // different fit. Didn't work on my first attempts, so I just create
277  // a new one.
278  minData.mChargeType = 1;
279  minuit = new TMinuit(1);
280  minuit->SetFCN(minimizeNegative);
281  minuit->mnparm( 0, "rho_ref scale factor", start, step, bmin, bmax, errFlag );
282  minuit->SetErrorDef(1);
283  minuit->SetPrintLevel(-1);
284  argList[0] = 1;
285  minuit->mnexcm("SET STR",argList,1,errFlag);
286  argList[0] = 500;
287  cout << ">>>>>Starting scale factor 1 fit for centrality " << ic << " yt bin " << ibin << endl;
288  minuit->mnexcm("MIGRAD",argList,1,errFlag);
289  minuit->GetParameter(0,sFactor[1][ibin][ic],eSFactor[1][ibin][ic]);
290  if (0 != errFlag) {
291  cout << "++++Out of chesse error; Fit failed. Using scale factor = 1" << endl;
292  sFactor[1][ibin][ic] = 1;
293  eSFactor[1][ibin][ic] = 100;
294  }
295  delete minuit;
296 
297  sf[0] = sFactor[0][ibin][ic];
298  sf[1] = sFactor[1][ibin][ic];
299  ytyt = ehelp->buildChargeTypes("YtYt",5,sf);
300  ytytC = ehelp->buildCommon("YtYt",5,sf);
301  sytdyt = ehelp->buildChargeTypes("SYtDYt",5,sf);
302  sytdytC = ehelp->buildCommon("SYtDYt",5,sf);
303 
304  out->cd();
305  for (int icharge=0;icharge<4;icharge++) {
306  TString name(binName[ytytBins[ibin]]);
307  name += "_YtYt"; name += chargeName[icharge]; name += ic;
308  ytyt[icharge]->SetName(name.Data());
309  ytyt[icharge]->SetTitle(name.Data());
310  ytyt[icharge]->Write();
311  TString name(binName[ytytBins[ibin]]);
312  name += "_YtYt"; name += chargeType[icharge]; name += ic;
313  ytytC[icharge]->SetName(name.Data());
314  ytytC[icharge]->SetTitle(name.Data());
315  ytytC[icharge]->Write();
316 
317  TString name(binName[ytytBins[ibin]]);
318  name += "_SYtDYt"; name += chargeName[icharge]; name += ic;
319  sytdyt[icharge]->SetName(name.Data());
320  sytdyt[icharge]->SetTitle(name.Data());
321  sytdyt[icharge]->Write();
322  TString name(binName[ytytBins[ibin]]);
323  name += "_SYtDYt"; name += chargeType[icharge]; name += ic;
324  sytdytC[icharge]->SetName(name.Data());
325  sytdytC[icharge]->SetTitle(name.Data());
326  sytdytC[icharge]->Write();
327  }
328  delete tf;
329  delete ehelp;
330  }
331  }
332  cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << endl;
333  cout << " Here are Yt scale factors for rho_{ref}" << endl;
334  printf(" Centrality all LS SS LS AS LS all US SS US AS US \n");
335  for (int ic=0;ic<nCent;ic++) {
336  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,
337  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],
338  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]);
339  }
340  out->Close();
341 // delete out;
342  for (int it=0;it<2;it++) {
343  for (int ibin=0;ibin<3;ibin++) {
344  delete sFactor[it][ibin];
345  delete eSFactor[it][ibin];
346  }
347  }
348 }
Generated by   doxygen 1.8.5