hooks for PMD flow analysis
[u/mrichter/AliRoot.git] / PWG2 / FLOW / macros / plotCumulants.C
CommitLineData
a663ab5a 1// Macro plotCumulants.C is used to represent true correlations, a.k.a. cumulants,
2// as a function of reference multiplicity. Different order cumulants provide
3// independent estimates of the flow harmonics. In the ideal case when only flow
4// correlations are present we have:
5// c{2} = v^2, c{4} = -v^4, c{6} = 4v^6 and c{8} = -33v^8.
6// In practice cumulants can are calculated via Q-vectors (QC) or via formalism
7// of generating functions (GFC).
8
9dcf9913 9// Set how many output analysis files in total you want to access:
34f4fef7 10const Int_t nFiles = 2;
9dcf9913 11
12// Set how many of those output analysis files you want to represent with a mesh (usually used to represent results of simulations):
a9a2547e 13const Int_t nSim = 0;
dfb33a36 14
9dcf9913 15// Set paths of all output analysis files (first the ones to be represented with mesh (simulations), then the ones to be represented with markers (real data))
5c09ff70 16TString files[nFiles] = {"trackletsCorrectedNUA","default"}; // subdirectory names holding <commonOutputFileName>.root
17//TString files[nFiles] = {"outputCentrality0","outputCentrality1"}; // file names (use case: centrality train)
63cdf10a 18
5c09ff70 19// Set analysis types for all output analysis files (can be "ESD","AOD","MC","","MK", ....):
34f4fef7 20TString type[nFiles] = {"ESD","ESD"};
5c09ff70 21//TString type[nFiles] = {"MK","MK"};
9dcf9913 22
23// Set mesh color:
a9a2547e 24Int_t meshColor[nSim] = {};
9dcf9913 25
26// Set marker styles:
34f4fef7 27Int_t markerStyle[nFiles-nSim] = {kFullSquare,kOpenSquare};
9dcf9913 28
c109b8c7 29// Set marker colors:
34f4fef7 30Int_t markerColor[nFiles-nSim] = {kBlack,kRed};
c109b8c7 31
9dcf9913 32// Set legend entries:
34f4fef7 33TString legendEntry[nFiles] = {"",""};
63cdf10a 34
83cac266 35// Set if you want to rebin the histograms into wider multiplicity bins (set for each cumulant order separately):
34f4fef7 36Bool_t rebin = kFALSE;
37Int_t nMergedBins[4] = {10,10,10,10}; // set how many original multiplicity bins will be merged into 1 new one
9dcf9913 38
c109b8c7 39// Set if you whish to plot cumulants versus <reference multiplicity> (by default they are plotted versus # of RPs):
34f4fef7 40Bool_t plotCumulantsVsReferenceMultiplicity = kFALSE;
a9a2547e 41Bool_t showReferenceMultiplicityVsNoOfRPs = kFALSE;
c109b8c7 42
9dcf9913 43// Set flow values whose theoretical contribution to cumulants will be shown on the plots with the straight coloured lines:
a9a2547e 44Bool_t showTheoreticalLines = kFALSE;
63cdf10a 45const Int_t nFlowValues = 1;
a663ab5a 46Double_t v[nFlowValues] = {0.05};
63cdf10a 47Int_t lineColor[nFlowValues] = {kRed};
9dcf9913 48
49// If the statistical error of 6th and 8th order cumulant is huge you may prefer not to show them:
16eed38e 50Bool_t plotOnly2ndAnd4thOrderCumulant = kFALSE;
51
52// Set if you want independent canvas with results for reference flow vs multiplicity:
53Bool_t showRefFlowVsM = kTRUE;
34f4fef7 54Int_t refFlowVsMMarkerStyle[nFiles] = {kFullSquare,kOpenSquare}; // marker style is different for different file
55Int_t refFlowVsMMarkerColor[4] = {kBlack,kRed,kBlue,kGreen+2}; // marker color is different for different cumulant order
56Int_t refFlowVsMeshOrder = -1; // set results of which order will be plotted as mesh in all pads [0=2nd, 1=4th, 2=6th, 3=8th, -1=do not show]
57Int_t refFlowVsMMeshColor = kRed-10; // mesh color for above specified order
16eed38e 58Double_t refFlowVsMxRange[2] = {1.,11000.}; // x range on the plots for reference multiplicity vs M
59Double_t refFlowVsMyRange[2] = {0.0,0.194}; // x range on the plots for reference multiplicity vs M
9dcf9913 60
a9a2547e 61// Set if you want to show theoretical curves for the toy model:
62Bool_t showToyModel = kFALSE;
63const Int_t nToyModels = 2; // number of toy models with different parameters k, vn, v2n
64Double_t k[nToyModels] = {2.,4.};
65Double_t vn[nToyModels] = {0.25,0.25};
66Double_t v2n[nToyModels] = {0.0,0.0};
67Int_t lineColorToyModel[nToyModels] = {kBlack,kRed};
68
9dcf9913 69// For comparison sake show also GFC results with dotted line:
70Bool_t showAlsoGFCResults = kFALSE;
71Int_t gfcLineStyle = 3;
72
dfb33a36 73// For comparison sake show also Monte Carlo QC results with coloured mesh:
74Bool_t showAlsoMCEPResults = kFALSE;
75Bool_t showOnlyMCEPResults = kFALSE;
a9a2547e 76Int_t mcepMeshColor[nSim] = {};
dfb33a36 77
5c09ff70 78// Set the naming convention:
79Bool_t bLookInSubdirectories = kFALSE; // if kTRUE: Look in subdirectories <files[0]>, <files[1]>, ..., for output files <commonOutputFileName>
80 // if kFALSE: Look in <pwd> for files <files[0]>.root, <files[1]>.root, ....
81TString commonOutputFileName = "AnalysisResults.root";
82
dfb33a36 83// Set method names which calculate cumulants vs multiplicity (do not touch these settings unless you are looking for a trouble):
84const Int_t nMethods = 3;
85TString method[nMethods] = {"QC","GFC","MCEP"};
9dcf9913 86
87TFile *commonOutputFiles[nFiles] = {NULL}; // common output files "AnalysisResults.root"
88TList *lists[nFiles][nMethods] = {{NULL}}; // lists cobj<method> holding objects with results for each method
89TH1D *cumulantsVsM[nFiles][nMethods][4] = {{{NULL}}}; // histograms with results for cumulants vs multiplicity (4 stands for 4 cumulant orders)
16eed38e 90TH1D *refFlowVsM[nFiles][nMethods][4] = {{{NULL}}}; // histograms with results for reference flow vs multiplicity (4 stands for 4 cumulant orders)
63cdf10a 91TGraph *lines[nFlowValues][4] = {{NULL}}; // lines denoting theoretical flow contribution to cumulants
92TProfile *refMultVsNoOfRPs[nFiles] = {NULL}; // <reference multipicity> versus # of RPs
a9a2547e 93TF1 *toyModel[2][nToyModels] = {{NULL}}; // [cumulant order][number of toy models]
9dcf9913 94
95// Ranges for plots:
96Double_t xMin[4]={0.};
97Double_t xMax[4]={0.};
98Double_t yMin[4]={0.};
99Double_t yMax[4]={0.};
100
101enum libModes {mLocal,mLocalSource};
102
103void plotCumulants(Int_t analysisMode=mLocal)
104{
105 // analysisMode: if analysisMode = mLocal -> analyze data on your computer using aliroot
106 // if analysisMode = mLocalSource -> analyze data on your computer using root + source files
a663ab5a 107
108 // Cross-check user settings:
109 CrossCheckSettings();
63cdf10a 110
9dcf9913 111 // Load needed libraries:
112 LoadLibrariesPC(analysisMode);
113
dfb33a36 114 // Global settings which will affect all plots:
115 GlobalSettings();
5c09ff70 116
9dcf9913 117 // Access all common output files:
9dcf9913 118 AccessCommonOutputFiles(commonOutputFileName);
5c09ff70 119
9dcf9913 120 // Get from common output files the lists holding histograms for each method:
121 GetLists();
5c09ff70 122
9dcf9913 123 // Get histograms with results for cumulants vs multiplicity:
124 GetHistograms();
16eed38e 125
9dcf9913 126 // Determine ranges for plots:
127 DetermineMinMax();
dfb33a36 128
9dcf9913 129 // Make lines which indicate theoretical contributions of flow to cumulants:
130 Lines();
131
132 // Print number of events and average multiplicities for each common output file:
133 Print();
134
16eed38e 135 // Make plots for cumulants vs multiplicity: // to be improved
136 PlotCumulantsVsM();
137
138 // Make plots for reference flow vs multiplicity:
139 if(showRefFlowVsM){PlotRefFlowVsM();}
a663ab5a 140
34f4fef7 141} // end of void plotCumulants(Int_t analysisMode=mLocal)
142
9dcf9913 143// =====================================================================================
144
16eed38e 145void PlotRefFlowVsM()
9dcf9913 146{
34f4fef7 147 // Make plots for reference flow vs multiplicity.
148
149 // Calculate reference flow from cumulants vs M:
150 CalculateReferenceFlowFromCumulantsVsM();
151
16eed38e 152 TCanvas *cRefFlowVsM = new TCanvas("cRefFlowVsM","Reference Flow");
153 cRefFlowVsM->Divide(2,2);
154
155 TLegend *lRefFlowVsM = new TLegend(0.1,0.7,0.33,0.9);
156 lRefFlowVsM->SetFillStyle(0);
34f4fef7 157 lRefFlowVsM->SetHeader(" Therminator (20-30%)");
16eed38e 158
159 TString refFlowVsMFlag[4] = {"v_{2}{2,QC}","v_{2}{4,QC}","v_{2}{6,QC}","v_{2}{8,QC}"};
160
34f4fef7 161 for(Int_t f=0;f<nFiles;f++)
16eed38e 162 {
34f4fef7 163 for(Int_t m=0;m<1+(Int_t)(showAlsoGFCResults)+(Int_t)(showAlsoMCEPResults);m++)
16eed38e 164 {
34f4fef7 165 for(Int_t co=0;co<4;co++) // cumulant order
166 {
167 cRefFlowVsM->cd(co+1);
168 // Style histogram:
169 if(f==0) // superimpose histograms from other files on top of this one
170 {
171 TH1D *styleHist = (TH1D*)StyleHist(refFlowVsMFlag[co].Data(),co)->Clone();
172 if(styleHist)
173 {
174 styleHist->GetXaxis()->SetRangeUser(refFlowVsMxRange[0],refFlowVsMxRange[1]);
175 styleHist->GetYaxis()->SetRangeUser(refFlowVsMyRange[0],refFlowVsMyRange[1]);
176 styleHist->Draw();
177 }
178 } // end of if(f==0)
179 // Plot first the mesh for reference flow of specified order in all pads:
180 if(refFlowVsMeshOrder != -1)
181 {
182 TGraph *rfMeshOutOf = GetErrorMesh(refFlowVsM[f][m][refFlowVsMeshOrder]);
183 if(rfMeshOutOf)
184 {
185 rfMeshOutOf->SetFillColor(refFlowVsMMeshColor);
186 rfMeshOutOf->Draw("lfsame");
187 if(refFlowVsMeshOrder==co){lRefFlowVsM->AddEntry(rfMeshOutOf,Form("v_{2}{%i,QC} stat. error",2*(refFlowVsMeshOrder+1)),"f");}
188 }
189 } // end of if(refFlowVsMeshOrder != -1)
190 // Plot reference flow vs M:
191 refFlowVsM[f][m][co]->SetMarkerStyle(refFlowVsMMarkerStyle[f]);
192 refFlowVsM[f][m][co]->SetMarkerColor(refFlowVsMMarkerColor[co]);
193 refFlowVsM[f][m][co]->Draw("e1same");
194 lRefFlowVsM->AddEntry(refFlowVsM[f][m][co],Form("v_{2}{%i,QC}, %s",2*(co+1),files[f].Data()),"p");
195 } // end of for(Int_t co=0;co<4;co++) // cumulant order
196 } // end of for(Int_t m=0;m<nMethods;m++)
197 } // end of for(Int_t f=0;f<nFiles;f++)
16eed38e 198
34f4fef7 199 // Draw a common legend in the 1st pad:
16eed38e 200 cRefFlowVsM->cd(1);
201 lRefFlowVsM->Draw("same");
202
203} // end of void PlotRefFlowVsM()
204
205// =====================================================================================
206
207void PlotCumulantsVsM()
208{
209 // Make plots for cumulants vs multiplicity: // to be improved
9dcf9913 210
211 TCanvas *c = NULL;
212 Int_t coMax = 0;
213 if(!plotOnly2ndAnd4thOrderCumulant)
214 {
215 c = new TCanvas("c","cumulants");
216 c->Divide(2,2);
217 coMax = 4;
218 } else
219 {
220 c = new TCanvas("c","cumulants",1200,500);
221 c->Divide(2,1);
222 coMax = 2;
a9a2547e 223 }
9dcf9913 224
225 TLegend *legend = new TLegend(0.1,0.7,0.33,0.9);
226 legend->SetFillStyle(0);
63cdf10a 227 //legend->SetHeader(" minClustersTpcRP");
9dcf9913 228
229 TString qcFlag[4] = {"QC{2}","QC{4}","QC{6}","QC{8}"};
230
231 for(Int_t co=0;co<coMax;co++) // cumulant order
232 {
233 c->cd(co+1);
16eed38e 234 TH1D *styleHist = (TH1D*) StyleHist(qcFlag[co].Data(),co)->Clone();
235 if(styleHist){styleHist->Draw();}
a9a2547e 236 if(co==0)
237 {
238 if(showToyModel)
239 {
240 for(Int_t ntm=0;ntm<nToyModels;ntm++)
241 {
242 TF1 *tm = ToyModel(co,k[ntm],vn[ntm],v2n[ntm]);
243 if(tm)
244 {
245 tm->SetLineColor(lineColorToyModel[ntm]);
246 tm->Draw("same");
247 }
248 } // for(Int_t ntm=0;ntm<nToyModels;ntm++)
249 } // end of if(showToyModel)
250 } // end of if(co==0)
251 if(co==1)
252 {
253 if(showToyModel)
254 {
255 for(Int_t ntm=0;ntm<nToyModels;ntm++)
256 {
257 TF1 *tm = ToyModel(co,k[ntm],vn[ntm],v2n[ntm]);
258 if(tm)
259 {
260 tm->SetLineColor(lineColorToyModel[ntm]);
261 tm->Draw("same");
262 }
263 } // for(Int_t ntm=0;ntm<nToyModels;ntm++)
264 } // end of if(showToyModel)
265 } // end of if(co==1)
266
9dcf9913 267 // simulations:
268 for(Int_t s=0;s<nSim;s++)
269 {
dfb33a36 270 // Monte Carlo QC:
271 if(showAlsoMCEPResults)
272 {
273 TGraph *mcepQC = GetErrorMesh(cumulantsVsM[s][2][co]);
274 if(mcepQC)
275 {
276 mcepQC->SetFillColor(mcepMeshColor[s]);
277 mcepQC->Draw("lfsame");
278 }
279 if(co==0){legend->AddEntry(mcepQC,Form("%s (MC)",legendEntry[s].Data()),"f");}
280 } // end of if(showAlsoMCEPResults)
281 // QC:
9dcf9913 282 TGraph *errorMesh = GetErrorMesh(cumulantsVsM[s][0][co]);
283 if(errorMesh)
284 {
285 errorMesh->SetFillColor(meshColor[s]);
dfb33a36 286 if(!(showOnlyMCEPResults && showAlsoMCEPResults)){errorMesh->Draw("lfsame");}
9dcf9913 287 }
dfb33a36 288 if(co==0 && !(showOnlyMCEPResults && showAlsoMCEPResults)){legend->AddEntry(errorMesh,legendEntry[s].Data(),"f");}
9dcf9913 289 } // end of if(Int_t s=0;s<nSim;s++)
dfb33a36 290 // Theoretical lines:
291 if(showTheoreticalLines)
292 {
293 for(Int_t fv=0;fv<nFlowValues;fv++)
294 {
295 lines[fv][co]->Draw("lsame");
296 if(co==0){legend->AddEntry(lines[fv][co],Form("v_{2} = %g",v[fv]),"l");}
297 }
298 }
9dcf9913 299 // data:
300 for(Int_t f=nSim;f<nFiles;f++)
301 {
9dcf9913 302 // QC results:
303 if(cumulantsVsM[f][0][co])
304 {
305 cumulantsVsM[f][0][co]->Draw("e1same");
c109b8c7 306 cumulantsVsM[f][0][co]->SetMarkerStyle(markerStyle[f-nSim]);
307 cumulantsVsM[f][0][co]->SetMarkerColor(markerColor[f-nSim]);
34f4fef7 308 if(co==0)
9dcf9913 309 {
310 if(showAlsoGFCResults)
311 {
312 legend->AddEntry(cumulantsVsM[f][0][co],Form("%s (QC)",legendEntry[f].Data()),"p");
313 } else
314 {
315 legend->AddEntry(cumulantsVsM[f][0][co],legendEntry[f].Data(),"p");
316 }
317 }
318 }
319 // GFC results:
320 if(showAlsoGFCResults && cumulantsVsM[f][1][co])
321 {
322 cumulantsVsM[f][1][co]->Draw("lsame");
323 cumulantsVsM[f][1][co]->SetLineStyle(gfcLineStyle);
324 if(co==0){legend->AddEntry(cumulantsVsM[f][1][co],Form("%s (GFC)",legendEntry[f].Data()),"l");}
325 }
326 }
327 // Draw legend:
16eed38e 328 if(co==0){legend->Draw("same");}
9dcf9913 329 } // end of for(Int_t co=0;co<4;co++) // cumulant order
a663ab5a 330
c109b8c7 331 // Plot also <reference multiplicity> vs # of RPs:
a663ab5a 332 if(plotCumulantsVsReferenceMultiplicity && showReferenceMultiplicityVsNoOfRPs)
c109b8c7 333 {
63cdf10a 334 TCanvas *cRefMultVsNoOfRPs = new TCanvas("cRefMultVsNoOfRPs","#LTreference multiplicity#GT vs # of RPs",1200,600);
335 cRefMultVsNoOfRPs->Divide(nFiles,1);
336 for(Int_t f=0;f<nFiles;f++)
c109b8c7 337 {
63cdf10a 338 cRefMultVsNoOfRPs->cd(f+1);
339 if(refMultVsNoOfRPs[f])
340 {
341 refMultVsNoOfRPs[f]->SetTitle(legendEntry[f].Data());
342 refMultVsNoOfRPs[f]->GetXaxis()->SetRangeUser(0,refMultVsNoOfRPs[f]->FindLastBinAbove());
343 refMultVsNoOfRPs[f]->Draw();
344 } else
345 {
346 cout<<endl;
347 cout<<"WARNING: refMultVsNoOfRPs[f] is NULL in Plot(), f = "<<f<<" !!!!"<<endl;
348 cout<<endl;
349 }
350 }
c109b8c7 351 }
352
16eed38e 353} // end of void PlotCumulantsVsM()
9dcf9913 354
355// =====================================================================================
356
a663ab5a 357void CrossCheckSettings()
358{
359 // Cross-check user settings in this method.
360
361 if(showAlsoGFCResults && rebin)
362 {
363 cout<<endl;
364 cout<<" WARNING: Rebinning in M not supported for GFC yet !!!!"<<endl;
365 cout<<endl;
366 exit(0);
367 }
368
369 if(showAlsoGFCResults && plotCumulantsVsReferenceMultiplicity)
370 {
371 cout<<endl;
372 cout<<" WARNING: Showing GFC versus <reference multiplicity> not supported yet !!!!"<<endl;
373 cout<<endl;
374 exit(0);
375 }
376
377} // end of void CrossCheckSettings()
378
379// =====================================================================================
380
9dcf9913 381void Lines()
382{
383 // Make lines denoting theoretical contribution of flow to cumulants.
9dcf9913 384
385 for(Int_t fv=0;fv<nFlowValues;fv++)
386 {
63cdf10a 387 lines[fv][0] = new TGraph(2);
388 lines[fv][0]->SetPoint(0,xMin[0],pow(v[fv],2));
389 lines[fv][0]->SetPoint(1,xMax[0]+0.5,pow(v[fv],2));
9dcf9913 390 lines[fv][0]->SetLineColor(lineColor[fv]);
63cdf10a 391 lines[fv][1] = new TGraph(2);
392 lines[fv][1]->SetPoint(0,xMin[1],-pow(v[fv],4));
393 lines[fv][1]->SetPoint(1,xMax[1]+0.5,-pow(v[fv],4));
394 lines[fv][1]->SetLineColor(lineColor[fv]);
395 lines[fv][2] = new TGraph(2);
396 lines[fv][2]->SetPoint(0,xMin[2],4.*pow(v[fv],6));
397 lines[fv][2]->SetPoint(1,xMax[2]+0.5,4.*pow(v[fv],6));
398 lines[fv][2]->SetLineColor(lineColor[fv]);
399 lines[fv][3] = new TGraph(2);
400 lines[fv][3]->SetPoint(0,xMin[3],-33.*pow(v[fv],8));
401 lines[fv][3]->SetPoint(1,xMax[3]+0.5,-33.*pow(v[fv],8));
9dcf9913 402 lines[fv][3]->SetLineColor(lineColor[fv]);
403 }
404
405} // end of void Lines()
406
407// =====================================================================================
408
a9a2547e 409TF1* ToyModel(Int_t co, Double_t k, Double_t vn, Double_t v2n)
410{
411 // Make theoretical curves for the toy model.
412
413 TF1 *function = NULL;
414
415 if(co==0)
416 {
417 function = new TF1("function","([1]*[1]*(x-0.5-[0])+[0]-1)/((x-0.5)-1)",1.5,1000);
418 function->SetParameter(0,k);
419 function->SetParameter(1,vn);
420 }
421 else if(co==1)
422 {
423 function = new TF1("function","-1.*([1]*[1]*[1]*[1]*([0]-(x-0.5))*(12.*[0]-6.*[0]*[0]-12.*(x-0.5)-5.*[0]*(x-0.5)+6.*[0]*[0]*(x-0.5)+9.*(x-0.5)*(x-0.5)-3.*[0]*(x-0.5)*(x-0.5)-(x-0.5)*(x-0.5)*(x-0.5))+[1]*[1]*4.*([0]-1.)*(x-0.5-[0])*([0]*(x-0.5-1.)-2.*(x-0.5-2.))-[1]*[1]*[2]*2.*(x-0.5-1.)*([0]-1.)*(x-0.5-[0])*(x-0.5-2.*[0])-[2]*[2]*([0]-1.)*([0]-1.)*((x-0.5)-[0])*(x-0.5-1.)+([0]-1.)*(-6.+9*[0]-[0]*[0]+2.*(x-0.5)-5.*[0]*(x-0.5)+[0]*[0]*(x-0.5)))/((x-0.5-1.)*(x-0.5-1.)*(x-0.5-2.)*(x-0.5-3.))",3.5,100);
424 function->SetParameter(0,k);
425 function->SetParameter(1,vn);
426 function->SetParameter(2,v2n);
427 }
428
429 return function;
430
431} // end of TF1* ToyModel(Int_t co)
432
433// =====================================================================================
434
34f4fef7 435void CalculateReferenceFlowFromCumulantsVsM()
436{
437 // Calculate reference flow from cumulants vs M:
438
439 for(Int_t f=0;f<nFiles;f++)
440 {
441 for(Int_t m=0;m<1+(Int_t)(showAlsoGFCResults)+(Int_t)(showAlsoMCEPResults);m++)
442 {
443 for(Int_t co=0;co<4;co++) // cumulant order
444 {
445 if(cumulantsVsM[f][m][co])
446 {
447 refFlowVsM[f][m][co] = (TH1D*)cumulantsVsM[f][m][co]->Clone(Form("%i,%i,%i",f,m,co));
448 if(!refFlowVsM[f][m][co]){cout<<" WARNING: "<<Form("refFlowVsM[%i][%i][%i]",f,m,co)<<" is NULL in PlotRefFlowVsM() !!!!"<<exit(0)<<endl;}
449 refFlowVsM[f][m][co]->Reset(); // to have empty histogram with the same binning as the one with cumulants
450 Int_t nBins = refFlowVsM[f][m][co]->GetNbinsX();
451 for(Int_t b=1;b<=nBins;b++)
452 {
453 Double_t qcVsM = cumulantsVsM[f][m][co]->GetBinContent(b); // QC vs M
454 Double_t qcErrorVsM = cumulantsVsM[f][m][co]->GetBinError(b); // QC vs M stat. error
455 Double_t vVsM = 0.; // reference flow vs M
456 Double_t vErrorVsM = 0.; // reference flow vs M stat. error
457 if(co==0) // 2nd order
458 {
459 if(qcVsM>0.)
460 {
461 vVsM = pow(qcVsM,1./2.);
462 vErrorVsM = (1./2.)*pow(qcVsM,-1./2.)*qcErrorVsM;
463 }
464 } // end of if(co==0) 2nd order
465 else if(co==1) // 4th order
466 {
467 if(qcVsM<0.)
468 {
469 vVsM = pow(-1.*qcVsM,1./4.);
470 vErrorVsM = (1./4.)*pow(-qcVsM,-3./4.)*qcErrorVsM;
471 }
472 } // end of if(co==1) 4th order
473 else if(co==2) // 6th order
474 {
475 if(qcVsM>0.)
476 {
477 vVsM = pow((1./4.)*qcVsM,1./6.);
478 vErrorVsM = (1./6.)*pow(2.,-1./3.)*pow(qcVsM,-5./6.)*qcErrorVsM;
479 }
480 } // end of if(co==2) 6th order
481 else if(co==3) // 8th order
482 {
483 if(qcVsM<0.)
484 {
485 vVsM = pow((-1./33.)*qcVsM,1./8.);
486 vErrorVsM = (1./8.)*pow(33.,-1./8.)*pow(-qcVsM,-7./8.)*qcErrorVsM;
487 }
488 } // end of if(co==3) 8th order
489 // Store final results and statisticial errror for reference flow:
490 refFlowVsM[f][m][co]->SetBinContent(b,vVsM);
491 refFlowVsM[f][m][co]->SetBinError(b,vErrorVsM);
492 } // end of for(Int_t b=1;b<=nBins;b++)
493 } else
494 {
495 cout<<endl;
496 cout<<" WARNING: "<<Form("cumulantsVsM[%i][%i][%i]",f,m,co)<<" is NULL in CalculateReferenceFlowFromCumulantsVsM() !!!!"<<endl;
497 cout<<endl;
498 }
499 } // end of for(Int_t co=0;co<4;co++) // cumulant order
500 } // end of for(Int_t m=0;m<nMethods;m++)
501 } // end of for(Int_t f=0;f<nFiles;f++)
502
503} // end of void CalculateReferenceFlowFromCumulantsVsM()
504
505// =====================================================================================
506
9dcf9913 507void Print()
508{
509 // Print number of events and average multiplicities for each common output file.
510
511 cout<<endl;
512 cout<<"Accessed files:"<<endl;
513 cout<<endl;
514 for(Int_t f=0;f<nFiles;f++)
515 {
516 cout<<commonOutputFiles[f]->GetName()<<endl;
517 for(Int_t m=0;m<nMethods;m++)
518 {
63cdf10a 519 AliFlowCommonHist *commonHist = NULL;
520 if(lists[f][m])
521 {
522 commonHist = dynamic_cast<AliFlowCommonHist*> (lists[f][m]->FindObject(Form("AliFlowCommonHist%s",method[m].Data())));
523 }
524 Double_t nEvts = 0.;
525 Double_t AvM = 0.;
9dcf9913 526 if(commonHist && commonHist->GetHistMultRP())
527 {
528 nEvts = commonHist->GetHistMultRP()->GetEntries();
529 AvM = commonHist->GetHistMultRP()->GetMean();
530 }
531 if(!(strcmp(method[m].Data(),"QC")))
532 {
533 cout<<Form("%s:",method[m].Data())<<" <M> = "<<AvM<<", N = "<<nEvts<<endl;
534 }
535 if(!(strcmp(method[m].Data(),"GFC")) && showAlsoGFCResults)
536 {
537 cout<<Form("%s:",method[m].Data())<<" <M> = "<<AvM<<", N = "<<nEvts<<endl;
538 }
63cdf10a 539 } // end of for(Int_t m=0;m<nMethods;m++)
9dcf9913 540 cout<<endl;
541 } // end of for(Int_t f=0;f<nFiles;f++)
542
543} // end of void Print()
544
545// =====================================================================================
546
547void DetermineMinMax()
548{
549 // Determine ranges for plots.
550
551 for(Int_t co=0;co<4;co++)
552 {
553 xMin[co] = 0.; yMin[co] = 44.;
554 xMax[co] = -440000.; yMax[co] = -44.;
555 }
556
557 Double_t tfc[nFlowValues][4] = {{0.}}; // theoretical flow contribution
558 for(Int_t fv=0;fv<nFlowValues;fv++)
559 {
560 tfc[fv][0] = pow(v[fv],2);
561 tfc[fv][1] = -pow(v[fv],4);
562 tfc[fv][2] = 4.*pow(v[fv],6);
563 tfc[fv][3] = -33.*pow(v[fv],8);
564 }
565
566 for(Int_t f=0;f<nFiles;f++)
567 {
a663ab5a 568 for(Int_t m=0;m<1+(Int_t)(showAlsoGFCResults)+(Int_t)(showAlsoMCEPResults);m++)
9dcf9913 569 {
570 for(Int_t co=0;co<4;co++)
571 {
572 if(cumulantsVsM[f][m][co])
573 {
63cdf10a 574 Int_t nBins = cumulantsVsM[f][m][co]->GetXaxis()->GetNbins();
575 for(Int_t b=1;b<=nBins;b++)
9dcf9913 576 {
577 Double_t result = cumulantsVsM[f][m][co]->GetBinContent(b);
578 Double_t error = cumulantsVsM[f][m][co]->GetBinError(b);
579 if(TMath::Abs(result)>1.e-44 && TMath::Abs(error)>1.e-44)
580 {
581 // y-axis:
582 if(yMin[co] > result-error){yMin[co] = result-error;} // min value
a663ab5a 583 if(yMax[co] < result+error){yMax[co] = result+error;} // max value
9dcf9913 584 // x-axis:
a663ab5a 585 if(xMax[co] < cumulantsVsM[f][m][co]->GetBinLowEdge(b+1)){xMax[co] = cumulantsVsM[f][m][co]->GetBinLowEdge(b+1);}
9dcf9913 586 }
587 } // end of for(Int_t b=1;b<=cumulantsVsM[f][m][co]->GetXaxis()->GetNbins();b++)
588 // theoretical contributions:
63cdf10a 589 if(showTheoreticalLines)
9dcf9913 590 {
63cdf10a 591 for(Int_t fv=0;fv<nFlowValues;fv++)
592 {
593 if(yMin[co] > tfc[fv][co]) {yMin[co] = tfc[fv][co];} // min value
594 if(yMax[co] < tfc[fv][co]) {yMax[co] = tfc[fv][co];} // max value
595 } // end of for(Int_t fv=0;fv<nFlowValues;fv++)
596 } // end of if(showTheoreticalLines)
9dcf9913 597 } // end of if(cumulantsVsM[f][m][co])
598 } // end of for(Int_t co=0;co<4;co++)
599 } // end of for(Int_t m=0;m<nMethods;m++)
600 } // end of for(Int_t f=0;f<nFiles;f++)
601
602} // end of void DetermineMinMax()
603
604// =====================================================================================
605
606void GetHistograms()
607{
a663ab5a 608 // Get all histograms and profiles.
609
610 // a) Get profiles holding results for <refMult> vs number of Reference Particles (RPs);
16eed38e 611 // b) Get histograms holding results for cumulants and reference flow vs multiplicity.
c109b8c7 612
a663ab5a 613 // a) Get profiles holding results for <refMult> vs number of Reference Particles (RPs):
614 if(plotCumulantsVsReferenceMultiplicity){GetProfileRefMultVsNoOfRPs();}
9dcf9913 615
16eed38e 616 // b) Get histograms holding results for cumulants and reference flow vs multiplicity:
9dcf9913 617 TString qcFlag[4] = {"QC{2}","QC{4}","QC{6}","QC{8}"};
618 TString gfcFlag[4] = {"GFC{2}","GFC{4}","GFC{6}","GFC{8}"};
619 for(Int_t f=0;f<nFiles;f++)
620 {
621 for(Int_t m=0;m<nMethods;m++)
622 {
623 TList *temp = NULL;
63cdf10a 624 if(!(strcmp(method[m].Data(),"QC")) && lists[f][m])
9dcf9913 625 {
626 temp = dynamic_cast<TList*> (lists[f][m]->FindObject("Integrated Flow"));
627 if(temp) {temp = dynamic_cast<TList*> (temp->FindObject("Results"));}
628 if(temp)
629 {
630 for(Int_t co=0;co<4;co++)
631 {
16eed38e 632 // Cumulants vs multiplicity:
9dcf9913 633 cumulantsVsM[f][m][co] = dynamic_cast<TH1D*> (temp->FindObject(Form("fIntFlowQcumulantsVsM, %s",qcFlag[co].Data())));
dfb33a36 634 if(plotCumulantsVsReferenceMultiplicity && cumulantsVsM[f][m][co])
635 {
636 cumulantsVsM[f][m][co] = Map(cumulantsVsM[f][m][co],f);
637 }
6324dcfd 638 if(rebin && cumulantsVsM[f][m][co])
63cdf10a 639 {
83cac266 640 cumulantsVsM[f][m][co] = Rebin(cumulantsVsM[f][m][co],co);
63cdf10a 641 }
a663ab5a 642 } // end of for(Int_t co=0;co<4;co++)
643 } // end of if(temp)
9dcf9913 644 } // end of if(!(strcmp(method[m].Data(),"QC")))
63cdf10a 645 else if(!(strcmp(method[m].Data(),"GFC")) && lists[f][m])
9dcf9913 646 {
647 temp = dynamic_cast<TList*> (lists[f][m]->FindObject("Reference Flow"));
648 if(temp) {temp = dynamic_cast<TList*> (temp->FindObject("Results"));}
649 if(temp)
650 {
651 for(Int_t co=0;co<4;co++)
652 {
653 cumulantsVsM[f][m][co] = dynamic_cast<TH1D*> (temp->FindObject(Form("fReferenceFlowCumulantsVsM, %s",gfcFlag[co].Data())));
a663ab5a 654 if(showAlsoGFCResults && !cumulantsVsM[f][m][co])
c109b8c7 655 {
a663ab5a 656 cout<<endl;
657 cout<<Form(" WARNING: Couldn't access histogram fReferenceFlowCumulantsVsM, %s in the ",gfcFlag[co].Data())<<endl;
658 cout<<" file "<<commonOutputFiles[f]->GetName()<<" !!!!"<<endl;
659 cout<<" Did you enable calculation of GFC vs M in the analysis which produced this file?"<<endl;
660 cout<<endl;
c109b8c7 661 }
dfb33a36 662 if(plotCumulantsVsReferenceMultiplicity && cumulantsVsM[f][m][co])
663 {
664 cumulantsVsM[f][m][co] = Map(cumulantsVsM[f][m][co],f);
a663ab5a 665 }
666 if(rebin && cumulantsVsM[f][m][co])
667 {
668 cumulantsVsM[f][m][co] = Rebin(cumulantsVsM[f][m][co],co);
669 }
63cdf10a 670 } // end of for(Int_t co=0;co<4;co++)
671 } // end of if(temp)
dfb33a36 672 } // end of else if(!(strcmp(method[m].Data(),"GFC")))
673 else if(!(strcmp(method[m].Data(),"MCEP")) && lists[f][m])
674 {
675 TProfile *mcepVsM = dynamic_cast<TProfile*> (lists[f][m]->FindObject("FlowPro_VsM_MCEP"));
676 if(!mcepVsM)
677 {
678 cout<<endl;
679 cout<<" WARNING: Couldn't access TProfile FlowPro_VsM_MCEP in the file "<<commonOutputFiles[f]->GetName()<<" !!!!"<<endl;
680 cout<<endl;
681 return;
682 }
683 for(Int_t co=0;co<4;co++)
684 {
685 cumulantsVsM[f][m][co] = new TH1D("",Form("Monte Carlo QC{%i} #font[72]{vs} M",2*(co+1)),
686 mcepVsM->GetNbinsX(),mcepVsM->GetBinLowEdge(1),mcepVsM->GetBinLowEdge(1+mcepVsM->GetNbinsX()));
687 } // end of for(Int_t co=0;co<4;co++)
688 for(Int_t b=1;b<=mcepVsM->GetNbinsX();b++)
689 {
690 Double_t v = mcepVsM->GetBinContent(b);
691 Double_t vError = mcepVsM->GetBinError(b);
692 if(TMath::Abs(v)<1.e-44 || TMath::Abs(vError)<1.e-44){continue;}
693 // Monte Carlo QC{2}:
694 Double_t qc2 = pow(v,2.);
695 Double_t qc2Error = 2.*TMath::Abs(v)*vError; // error propagation for f(x) = x^2
696 cumulantsVsM[f][m][0]->SetBinContent(b,qc2);
697 cumulantsVsM[f][m][0]->SetBinError(b,qc2Error);
698 // Monte Carlo QC{4}:
699 Double_t qc4 = -pow(v,4.);
700 Double_t qc4Error = 4.*TMath::Abs(pow(v,3.))*vError; // error propagation for f(x) = -x^4
701 cumulantsVsM[f][m][1]->SetBinContent(b,qc4);
702 cumulantsVsM[f][m][1]->SetBinError(b,qc4Error);
703 // Monte Carlo QC{6}:
704 Double_t qc6 = 4.*pow(v,6.);
705 Double_t qc6Error = 24.*TMath::Abs(pow(v,5.))*vError; // error propagation for f(x) = 4x^6
706 cumulantsVsM[f][m][2]->SetBinContent(b,qc6);
707 cumulantsVsM[f][m][2]->SetBinError(b,qc6Error);
708 // Monte Carlo QC{8}:
709 Double_t qc8 = -33.*pow(v,8.);
710 Double_t qc8Error = 264.*TMath::Abs(pow(v,7.))*vError; // error propagation for f(x) = -33x^8
711 cumulantsVsM[f][m][3]->SetBinContent(b,qc8);
712 cumulantsVsM[f][m][3]->SetBinError(b,qc8Error);
713 } // end of for(Int_t b=1;b<=mcepVsM->GetNbinsX();b++)
714 for(Int_t co=0;co<4;co++)
715 {
716 if(plotCumulantsVsReferenceMultiplicity && cumulantsVsM[f][m][co])
717 {
718 cumulantsVsM[f][m][co] = Map(cumulantsVsM[f][m][co],f);
719 }
720 if(rebin && cumulantsVsM[f][m][co])
721 {
722 cumulantsVsM[f][m][co] = Rebin(cumulantsVsM[f][m][co],co);
723 }
724 } // end of for(Int_t co=0;co<4;co++)
725 } // end of else if(!(strcmp(method[m].Data(),"MCEP")))
9dcf9913 726 } // end of for(Int_t m=0;m<nMethods;m++)
727 } // end of for(Int_t f=0;f<nFiles;f++)
728
6324dcfd 729 Int_t counter = 0;
730 for(Int_t f=0;f<nFiles;f++)
731 {
732 for(Int_t m=0;m<nMethods;m++)
733 {
734 for(Int_t co=0;co<4;co++)
735 {
736 if(cumulantsVsM[f][m][co]){counter++;}
737 }
738 }
739 }
740
741 if(counter == 0)
742 {
743 cout<<endl;
744 cout<<" WARNING: Couldn't access a single histogram with results vs multiplicity !!!!"<<endl;
745 cout<<" Did you enable this calculation before running over data?"<<endl;
746 cout<<endl;
747 exit(0);
748 }
749
9dcf9913 750} // end of void GetHistograms()
751
752// =====================================================================================
753
a663ab5a 754void GetProfileRefMultVsNoOfRPs()
c109b8c7 755{
a663ab5a 756 // Get profiles holding results for <reference multiplicity> vs number of Reference Particles (RPs).
c109b8c7 757
a663ab5a 758 // Set here from which method's output file this profile will be accessed:
759 TString methodName = "QC"; // Alternatives are GFC and MCEP
760 Int_t i = -1;
761 for(Int_t m=0;m<nMethods;m++)
762 {
763 if(method[m] == methodName){i=m;}
764 }
765 if(i==-1)
766 {
767 cout<<endl;
768 cout<<" WARNING: Unknown method name in GetProfileRefMultVsNoOfRPs() !!!!"<<endl;
769 cout<<" Try something else for TString methodName in GetProfileRefMultVsNoOfRPs()."<<endl;
770 cout<<endl;exit(0);
771 }
772
63cdf10a 773 for(Int_t f=0;f<nFiles;f++)
c109b8c7 774 {
63cdf10a 775 AliFlowCommonHist *commonHist = NULL;
a663ab5a 776 if(lists[f][i])
63cdf10a 777 {
a663ab5a 778 commonHist = dynamic_cast<AliFlowCommonHist*> (lists[f][i]->FindObject(Form("AliFlowCommonHist%s",method[i].Data())));
63cdf10a 779 } else
780 {
781 cout<<endl;
a663ab5a 782 cout<<Form(" WARNING: lists[%i][%i] is NULL in GetProfileRefMultVsNoOfRPs() !!!!",f,i)<<endl;
783 cout<<Form(" Did you use method %s in the analysis which has produced output",method[i].Data())<<endl;
784 cout<<Form(" file %s?",commonOutputFiles[f]->GetName())<<endl;
63cdf10a 785 }
786 if(commonHist && commonHist->GetRefMultVsNoOfRPs())
787 {
788 refMultVsNoOfRPs[f] = commonHist->GetRefMultVsNoOfRPs();
789 } else
790 {
791 cout<<endl;
a663ab5a 792 cout<<" WARNING: commonHist && commonHist->GetRefMultVsNoOfRPs() is NULL in GetProfileRefMultVsNoOfRPs() !!!!"<<endl;
63cdf10a 793 cout<<endl;
794 }
a663ab5a 795 } // end of for(Int_t f=0;f<nFiles;f++)
dfb33a36 796
797 return;
c109b8c7 798
a663ab5a 799} // end of void GetProfileRefMultVsNoOfRPs()
c109b8c7 800
801// =====================================================================================
802
dfb33a36 803TH1D* Map(TH1D *hist, Int_t f)
c109b8c7 804{
805 // Map cumulant versus # of RPs into cumulants versus <reference multiplicity>.
806
63cdf10a 807 if(!refMultVsNoOfRPs[f])
808 {
809 cout<<endl;
a663ab5a 810 cout<<Form(" WARNING: refMultVsNoOfRPs[%i] is NULL in Map(...) !!!!",f)<<endl;
63cdf10a 811 cout<<endl;
812 }
813
a663ab5a 814 Int_t rpMinBin = refMultVsNoOfRPs[f]->FindFirstBinAbove(); // FindFirstBinAbove(Double_t threshold = 0, Int_t axis = 1)
dfb33a36 815 Int_t rpMaxBin = refMultVsNoOfRPs[f]->FindLastBinAbove(); // FindLastBinAbove(Double_t threshold = 0, Int_t axis = 1)
a663ab5a 816 Int_t rmMinBin = 440000;;
dfb33a36 817 Int_t rmMaxBin = (Int_t)TMath::Floor(refMultVsNoOfRPs[f]->GetMaximum());
a663ab5a 818 for(Int_t rpBin=rpMinBin;rpBin<=rpMaxBin;rpBin++) // non-empty # of RPs bins
819 {
820 if(refMultVsNoOfRPs[f]->GetBinContent(rpBin)>0. && refMultVsNoOfRPs[f]->GetBinContent(rpBin)<rmMinBin)
821 {
822 rmMinBin = (Int_t)TMath::Floor(refMultVsNoOfRPs[f]->GetBinContent(rpBin));
823 }
824 } // end of for(Int_t rpBin=rpMinBin;rpBin<=rpMaxBin;rpBin++) // non-empty # of RPs bins
dfb33a36 825
c109b8c7 826 if(hist)
827 {
828 temp = (TH1D*) hist->Clone();
dfb33a36 829 temp->Reset();
a663ab5a 830 for(Int_t rmBin=rmMinBin;rmBin<=rmMaxBin;rmBin++) // reference multiplicity bins
dfb33a36 831 {
832 Double_t value = 0.;
833 Double_t error = 0.;
834 Double_t dSum1 = 0.; // sum value_i/(error_i)^2
835 Double_t dSum2 = 0.; // sum 1/(error_i)^2
a663ab5a 836 for(Int_t rpBin=rpMinBin;rpBin<=rpMaxBin;rpBin++) // # of RPs bins
dfb33a36 837 {
838 if((Int_t)TMath::Floor(refMultVsNoOfRPs[f]->GetBinContent(rpBin)) >= temp->GetBinLowEdge(rmBin) &&
839 (Int_t)TMath::Floor(refMultVsNoOfRPs[f]->GetBinContent(rpBin)) < temp->GetBinLowEdge(rmBin+1))
840 {
841 value = hist->GetBinContent(rpBin);
842 error = hist->GetBinError(rpBin);
843 if(error>0.)
844 {
845 dSum1+=value/(error*error);
846 dSum2+=1./(error*error);
847 }
848 }
849 } // end of for(Int_t rpBin=1;rpBin<=nBins;rpBin++) // # of RPs bins
850 if(dSum2>0.)
851 {
852 temp->SetBinContent(rmBin,dSum1/dSum2);
853 temp->SetBinError(rmBin,pow(1./dSum2,0.5));
854 }
855 } // end of for(Int_t rmBin=1;rmBin<=nBins;rmBin++) // reference multiplicity bins
856 } // end of if(hist)
857
858 return temp;
63cdf10a 859
860} // end of Map()
861
862// =====================================================================================
863
83cac266 864TH1D* Rebin(TH1D *hist, Int_t co)
63cdf10a 865{
866 // Rebin original histograms.
867
83cac266 868 if(nMergedBins[co] == 0)
63cdf10a 869 {
870 cout<<endl;
83cac266 871 cout<<Form(" WARNING: nMergedBins[%i] == 0 !!!!",co)<<endl;
63cdf10a 872 cout<<endl;
873 exit(0);
874 }
875 if(!hist)
876 {
877 cout<<endl;
878 cout<<" WARNING: hist is NULL in Rebin() !!!!"<<endl;
879 cout<<endl;
880 exit(0);
881 }
882
883 Int_t nBinsOld = hist->GetXaxis()->GetNbins();
884 if(nBinsOld == 0){cout<<" WARNING: nBinsOld == 0 !!!!"<<endl;exit(0);}
885 Double_t xMinOld = hist->GetXaxis()->GetXmin();
886 Double_t xMaxOld = hist->GetXaxis()->GetXmax();
887 Double_t binWidthOld = (xMaxOld-xMinOld)/nBinsOld;
83cac266 888 Int_t nBinsNew = TMath::Floor(nBinsOld/nMergedBins[co]);
63cdf10a 889 Double_t xMinNew = xMinOld;
83cac266 890 Double_t xMaxNew = xMinOld + nBinsNew*nMergedBins[co]*binWidthOld;
c109b8c7 891
63cdf10a 892 TH1D *temp = new TH1D("","",nBinsNew,xMinNew,xMaxNew); // rebinned histogram
893 Int_t binNew = 1;
894 Double_t value = 0.;
895 Double_t error = 0.;
896 Double_t dSum1 = 0.; // sum value_i/(error_i)^2
897 Double_t dSum2 = 0.; // sum 1/(error_i)^2
898 for(Int_t b=1;b<=nBinsOld;b++)
899 {
900 value = hist->GetBinContent(b);
901 error = hist->GetBinError(b);
902 if(error>0.)
903 {
904 dSum1+=value/(error*error);
905 dSum2+=1./(error*error);
906 }
83cac266 907 if(b%nMergedBins[co] == 0)
63cdf10a 908 {
909 if(dSum2>0.)
910 {
911 temp->SetBinContent(binNew,dSum1/dSum2);
912 temp->SetBinError(binNew,pow(1./dSum2,0.5));
913 }
914 binNew++;
915 dSum1 = 0.;
916 dSum2 = 0.;
83cac266 917 } // end of if(b%nMergedBins[co] == 0)
63cdf10a 918 } // end of for(Int_t b=1;b<=nBinsOld;b++)
c109b8c7 919
63cdf10a 920 return temp;
921
83cac266 922} // end of TH1D* Rebin(TH1D *hist, Int_t co)
c109b8c7 923
924// =====================================================================================
925
9dcf9913 926TGraphErrors* GetGraphErrors(Int_t bin, Int_t nFiles, TH1D** qc)
927{
928 TGraphErrors *ge = new TGraphErrors(nFiles);
929 for(Int_t f=0;f<nFiles;f++)
930 {
931 ge->SetPoint(f,f+0.5,qc[f]->GetBinContent(bin+1));
932 ge->SetPointError(f,0,qc[f]->GetBinError(bin+1));
933 }
934
935 return ge;
936
937} // end of TGraphErrors* GetGraphErrors(Int_t bin, Int_t nFiles, TH1D** qc)
938
939// =====================================================================================
940
941TH1D* StyleHist(TString yAxisTitle, Int_t co)
942{
943 // Style histogram.
944
16eed38e 945 TH1D *styleHist = new TH1D(yAxisTitle.Data(),"",(Int_t)xMax[co],0,xMax[co]);
9dcf9913 946 // y-axis:
947 styleHist->GetYaxis()->SetRangeUser(yMin[co],yMax[co]);
63cdf10a 948
c109b8c7 949 if(plotCumulantsVsReferenceMultiplicity)
950 {
951 styleHist->GetXaxis()->SetTitle("#LTreference multiplicity#GT");
952 } else
953 {
954 styleHist->GetXaxis()->SetTitle("# of RPs");
955 }
956
9dcf9913 957 styleHist->GetYaxis()->SetTitle(yAxisTitle.Data());
958
959 return styleHist;
960
961} // end of TH1D* StyleHist(TString yAxisTitle, Int_t co)
962
963// ===========================================================================================
964
965TGraph* GetErrorMesh(TH1D *hist)
966{
967 // Error mesh.
968
9dcf9913 969 if(hist)
970 {
971 Int_t nBins = hist->GetNbinsX();
dfb33a36 972 Double_t value = 0.;
973 Double_t error = 0.;
9dcf9913 974 // Counting the non-empty bins:
dfb33a36 975 Int_t nNonEmptyBins = 0;
9dcf9913 976 for(Int_t i=1;i<nBins+1;i++)
977 {
dfb33a36 978 value = hist->GetBinContent(i);
979 error = hist->GetBinError(i);
980 if(TMath::Abs(value)>0.0 && error>0.0)
9dcf9913 981 {
982 nNonEmptyBins++;
983 }
dfb33a36 984 } // end of for(Int_t i=1;i<nBins+1;i++)
985 // Error mesh:
986 TGraph *errorMesh = new TGraph(2*nNonEmptyBins+1);
987 Int_t count = 0;
988 Double_t binCenter = 0.;
9dcf9913 989 for(Int_t i=1;i<nBins+1;i++)
990 {
9dcf9913 991 value = hist->GetBinContent(i);
dfb33a36 992 error = hist->GetBinError(i);
993 // Setting up the the mesh:
994 if(TMath::Abs(value)>0.0 && error>0.0)
9dcf9913 995 {
dfb33a36 996 binCenter = hist->GetBinCenter(i);
997 errorMesh->SetPoint(count,binCenter,value+error);
998 errorMesh->SetPoint(2*nNonEmptyBins-(count+1),binCenter,value-error);
999 count++;
1000 }
1001 } // end of for(Int_t i=1;i<nBins+1;i++)
9dcf9913 1002 // Closing the mesh area:
dfb33a36 1003 Double_t xStart = 0.;
1004 Double_t yStart = 0.;
1005 errorMesh->GetPoint(0,xStart,yStart);
1006 errorMesh->SetPoint(2*nNonEmptyBins,xStart,yStart);
9dcf9913 1007 } // end if(hist)
1008
9dcf9913 1009 return errorMesh;
1010
1011} // end of TGraph* GetErrorMesh(TH1D *hist)
1012
1013// ===========================================================================================
1014
1015void GlobalSettings()
1016{
1017 // Settings which will affect all plots.
1018
1019 gROOT->SetStyle("Plain"); // default color is white instead of gray
1020 gStyle->SetOptStat(0); // remove stat. box from all histos
5c09ff70 1021 TGaxis::SetMaxDigits(4); // prefer exp notation for 5 and more significant digits
9dcf9913 1022
1023} // end of void GlobalSettings()
1024
1025// ===========================================================================================
1026
1027void GetLists()
1028{
1029 // Get from common output files the lists holding histograms for each method.
1030
1031 TString fileName[nFiles][nMethods];
1032 TDirectoryFile *dirFile[nFiles][nMethods] = {{NULL}};
1033 TString listName[nFiles][nMethods];
1034 for(Int_t f=0;f<nFiles;f++)
1035 {
1036 for(Int_t i=0;i<nMethods;i++)
1037 {
1038 // Form a file name for each method:
1039 fileName[f][i]+="output";
1040 fileName[f][i]+=method[i].Data();
1041 fileName[f][i]+="analysis";
1042 fileName[f][i]+=type[f].Data();
1043 // Access this file:
1044 if(commonOutputFiles[f]){dirFile[f][i] = (TDirectoryFile*)commonOutputFiles[f]->FindObjectAny(fileName[f][i].Data());}
1045 // Form a list name for each method:
9dcf9913 1046 if(dirFile[f][i])
1047 {
5c09ff70 1048 TList* listTemp = dirFile[f][i]->GetListOfKeys();
1049 if(listTemp && listTemp->GetEntries() == 1)
63cdf10a 1050 {
5c09ff70 1051 listName[f][i] = listTemp->At(0)->GetName(); // to be improved - implemented better
1052 dirFile[f][i]->GetObject(listName[f][i].Data(),lists[f][i]);
1053 } else
1054 {
1055 cout<<endl;
1056 cout<<" WARNING: Couldn't find a list "<<listName[f][i].Data()<<" in "<<commonOutputFiles[f]->GetName()<<" !!!!"<<endl;
1057 cout<<" Did you use method "<<method[i].Data()<<" in the analysis which produced this file?"<<endl;
1058 }
9dcf9913 1059 } else
1060 {
dfb33a36 1061 cout<<" WARNING: Couldn't find a file "<<fileName[f][i].Data()<<".root in "<<commonOutputFiles[f]->GetName()<<" !!!!"<<endl;
5c09ff70 1062 }
9dcf9913 1063 } // end of for(Int_t i=0;i<nMethods;i++)
1064 } // end of for(Int_t f=0;f<nFiles;f++)
1065
1066} // end of void GetLists()
1067
1068// ===========================================================================================
1069
1070void AccessCommonOutputFiles(TString commonOutputFileName)
1071{
1072 // Access all output files.
a663ab5a 1073
9dcf9913 1074 for(Int_t f=0;f<nFiles;f++)
1075 {
5c09ff70 1076 if(bLookInSubdirectories)
9dcf9913 1077 {
5c09ff70 1078 if(!(gSystem->AccessPathName(Form("%s/%s/%s",gSystem->pwd(),files[f].Data(),commonOutputFileName.Data()),kFileExists)))
1079 {
1080 commonOutputFiles[f] = TFile::Open(Form("%s/%s/%s",gSystem->pwd(),files[f].Data(),commonOutputFileName.Data()),"READ");
1081 } else
1082 {
1083 cout<<endl;
1084 cout<<" WARNING: Couldn't find the file "<<Form("%s/%s/%s",gSystem->pwd(),files[f].Data(),commonOutputFileName.Data())<<" !!!!"<<endl;
1085 cout<<" Did you specify correctly all paths in TString files[nFiles]?"<<endl;
1086 cout<<" Or you are misusing 'Bool_t bLookInSubdirectories'? "<<endl;
1087 cout<<endl;
1088 exit(0);
1089 }
1090 } else // to if(bLookInSubdirectories)
1091 {
1092 if(!(gSystem->AccessPathName(Form("%s/%s%s",gSystem->pwd(),files[f].Data(),".root"),kFileExists)))
1093 {
1094 commonOutputFiles[f] = TFile::Open(Form("%s/%s%s",gSystem->pwd(),files[f].Data(),".root"),"READ");
1095 } else
1096 {
1097 cout<<endl;
1098 cout<<" WARNING: Couldn't find the file "<<Form("%s/%s%s",gSystem->pwd(),files[f].Data(),".root")<<" !!!!"<<endl;
1099 cout<<" Did you specify correctly all file names in TString files[nFiles]?"<<endl;
1100 cout<<endl;
1101 exit(0);
1102 }
9dcf9913 1103 }
1104 } // end of for(Int_t f=0;f<nFiles;f++)
1105
1106} // void AccessCommonOutputFiles(TString commonOutputFileName);
1107
1108// ===========================================================================================
1109
1110void LoadLibrariesPC(const libModes analysisMode) {
1111
1112 //--------------------------------------
1113 // Load the needed libraries most of them already loaded by aliroot
1114 //--------------------------------------
1115 //gSystem->Load("libTree");
1116 gSystem->Load("libGeom");
1117 gSystem->Load("libVMC");
1118 gSystem->Load("libXMLIO");
1119 gSystem->Load("libPhysics");
1120
1121 //----------------------------------------------------------
1122 // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
1123 //----------------------------------------------------------
1124 if (analysisMode==mLocal) {
1125 //--------------------------------------------------------
1126 // If you want to use already compiled libraries
1127 // in the aliroot distribution
1128 //--------------------------------------------------------
1129
1130 //==================================================================================
1131 //load needed libraries:
1132 gSystem->AddIncludePath("-I$ROOTSYS/include");
1133 //gSystem->Load("libTree");
1134
1135 // for AliRoot
1136 gSystem->AddIncludePath("-I$ALICE_ROOT/include");
1137 //gSystem->Load("libANALYSIS");
1138 gSystem->Load("libPWG2flowCommon");
1139 //cerr<<"libPWG2flowCommon loaded ..."<<endl;
1140
1141 }
1142
1143 else if (analysisMode==mLocalSource) {
1144
1145 // In root inline compile
1146
3b67e1f4 1147 // Constants
9dcf9913 1148 gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
1149 gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
9dcf9913 1150
1151 // Flow event
1152 gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+");
1153 gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");
9dcf9913 1154 gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");
3b67e1f4 1155 gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
1156
9dcf9913 1157 // Output histosgrams
1158 gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
1159 gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
1160 gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
1161 gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
1162
1163 cout << "finished loading macros!" << endl;
1164
1165 }
1166
1167} // end of void LoadLibraries(const libModes analysisMode)