]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FLOW/Tools/macros/plotCumulants.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
added a directory with user macros
[u/mrichter/AliRoot.git] / PWG2 / FLOW / Tools / macros / plotCumulants.C
CommitLineData
9dcf9913 1// add comment
2
3// Set how many output analysis files in total you want to access:
4const Int_t nFiles = 2;
5
6// Set how many of those output analysis files you want to represent with a mesh (usually used to represent results of simulations):
7const Int_t nSim = 1;
8
9// 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))
10TString files[nFiles] = {"sim/pythia/LHC10a18","data/mergedBC"};
11
12// Set analysis types for all output analysis files (can be "ESD","AOD","MC",""):
13TString type[nFiles] = {"ESD","ESD"};
14
15// Set mesh color:
16Int_t meshColor[nSim] = {kBlue-10};
17
18// Set marker styles:
19Int_t markerStyle[nFiles-nSim] = {kStar};
20
21// Set legend entries:
22TString legendEntry[nFiles] = {"Pythia (LHC10a18)","7 TeV data (LHC10b, LHC10c)"};
23
24// Set flow values whose theoretical contribution to cumulants will be shown on the plots with the straight coloured lines:
25Bool_t showTheoreticalLines = kFALSE;
26const Int_t nFlowValues = 2;
27Double_t v[nFlowValues] = {0.1,0.05};
28Int_t lineColor[nFlowValues] = {kRed,kBlue};
29
30// If the statistical error of 6th and 8th order cumulant is huge you may prefer not to show them:
31Bool_t plotOnly2ndAnd4thOrderCumulant = kFALSE;
32
33// For comparison sake show also GFC results with dotted line:
34Bool_t showAlsoGFCResults = kFALSE;
35Int_t gfcLineStyle = 3;
36
37// Set method names which calculate cumulants vs multiplicity:
38const Int_t nMethods = 2;
39TString method[nMethods] = {"QC","GFC"};
40
41TFile *commonOutputFiles[nFiles] = {NULL}; // common output files "AnalysisResults.root"
42TList *lists[nFiles][nMethods] = {{NULL}}; // lists cobj<method> holding objects with results for each method
43TH1D *cumulantsVsM[nFiles][nMethods][4] = {{{NULL}}}; // histograms with results for cumulants vs multiplicity (4 stands for 4 cumulant orders)
44TLine *lines[nFlowValues][4] = {{NULL}}; // lines denoting theoretical flow contribution to cumulants
45
46// Ranges for plots:
47Double_t xMin[4]={0.};
48Double_t xMax[4]={0.};
49Double_t yMin[4]={0.};
50Double_t yMax[4]={0.};
51
52enum libModes {mLocal,mLocalSource};
53
54void plotCumulants(Int_t analysisMode=mLocal)
55{
56 // analysisMode: if analysisMode = mLocal -> analyze data on your computer using aliroot
57 // if analysisMode = mLocalSource -> analyze data on your computer using root + source files
58
59 // Load needed libraries:
60 LoadLibrariesPC(analysisMode);
61
62 // Access all common output files:
63 TString commonOutputFileName = "AnalysisResults.root";
64 AccessCommonOutputFiles(commonOutputFileName);
65
66 // Get from common output files the lists holding histograms for each method:
67 GetLists();
68
69 // Get histograms with results for cumulants vs multiplicity:
70 GetHistograms();
71
72 // Determine ranges for plots:
73 DetermineMinMax();
74
75 // Make lines which indicate theoretical contributions of flow to cumulants:
76 Lines();
77
78 // Print number of events and average multiplicities for each common output file:
79 Print();
80
81 // Make plots:
82 Plot();
83
84 // Global settings which will affect all plots:
85 GlobalSettings();
86
87} // end of void plotCumulants(Int_t analysisMode=mLocal)
88
89// =====================================================================================
90
91void Plot()
92{
93 // Make all plots.
94
95 TCanvas *c = NULL;
96 Int_t coMax = 0;
97 if(!plotOnly2ndAnd4thOrderCumulant)
98 {
99 c = new TCanvas("c","cumulants");
100 c->Divide(2,2);
101 coMax = 4;
102 } else
103 {
104 c = new TCanvas("c","cumulants",1200,500);
105 c->Divide(2,1);
106 coMax = 2;
107 }
108
109 TLegend *legend = new TLegend(0.1,0.7,0.33,0.9);
110 legend->SetFillStyle(0);
111
112 TString qcFlag[4] = {"QC{2}","QC{4}","QC{6}","QC{8}"};
113
114 for(Int_t co=0;co<coMax;co++) // cumulant order
115 {
116 c->cd(co+1);
117 StyleHist(qcFlag[co].Data(),co)->Draw();
118 // simulations:
119 for(Int_t s=0;s<nSim;s++)
120 {
121 TGraph *errorMesh = GetErrorMesh(cumulantsVsM[s][0][co]);
122 if(errorMesh)
123 {
124 errorMesh->SetFillColor(meshColor[s]);
125 errorMesh->Draw("lfsame");
126 }
127 if(co==0){legend->AddEntry(errorMesh,legendEntry[s].Data(),"f");}
128 } // end of if(Int_t s=0;s<nSim;s++)
129 // data:
130 for(Int_t f=nSim;f<nFiles;f++)
131 {
132 // Theoretical lines:
133 if(showTheoreticalLines)
134 {
135 if(f==nSim) // plot them only once
136 {
137 for(Int_t fv=0;fv<nFlowValues;fv++)
138 {
139 lines[fv][co]->Draw("same");
140 if(co==0){legend->AddEntry(lines[fv][co],Form("v_{2} = %g",v[fv]),"l");}
141 }
142 }
143 }
144 // QC results:
145 if(cumulantsVsM[f][0][co])
146 {
147 cumulantsVsM[f][0][co]->Draw("e1same");
148 cumulantsVsM[f][0][co]->SetMarkerStyle(markerStyle[f-nSim]);
149 if(co==0)
150 {
151 if(showAlsoGFCResults)
152 {
153 legend->AddEntry(cumulantsVsM[f][0][co],Form("%s (QC)",legendEntry[f].Data()),"p");
154 } else
155 {
156 legend->AddEntry(cumulantsVsM[f][0][co],legendEntry[f].Data(),"p");
157 }
158 }
159 }
160 // GFC results:
161 if(showAlsoGFCResults && cumulantsVsM[f][1][co])
162 {
163 cumulantsVsM[f][1][co]->Draw("lsame");
164 cumulantsVsM[f][1][co]->SetLineStyle(gfcLineStyle);
165 if(co==0){legend->AddEntry(cumulantsVsM[f][1][co],Form("%s (GFC)",legendEntry[f].Data()),"l");}
166 }
167 }
168 // Draw legend:
169 if(co==0){legend->Draw("same");}
170 } // end of for(Int_t co=0;co<4;co++) // cumulant order
171
172} // end of void Plot()
173
174// =====================================================================================
175
176void Lines()
177{
178 // Make lines denoting theoretical contribution of flow to cumulants.
179
180 for(Int_t co=0;co<4;co++)
181 {
182 xMin[co] = 0.;
183 //xMax[co] = 59.5;
184 }
185
186 for(Int_t fv=0;fv<nFlowValues;fv++)
187 {
188 lines[fv][0] = new TLine(xMin[0],pow(v[fv],2),xMax[0]+0.5,pow(v[fv],2));
189 lines[fv][0]->SetLineColor(lineColor[fv]);
190 lines[fv][1] = new TLine(xMin[1],-pow(v[fv],4),xMax[1]+0.5,-pow(v[fv],4));
191 lines[fv][1]->SetLineColor(lineColor[fv]);
192 lines[fv][2] = new TLine(xMin[2],4.*pow(v[fv],6),xMax[2]+0.5,4.*pow(v[fv],6));
193 lines[fv][2]->SetLineColor(lineColor[fv]);
194 lines[fv][3] = new TLine(xMin[3],-33.*pow(v[fv],8),xMax[3]+0.5,-33.*pow(v[fv],8));
195 lines[fv][3]->SetLineColor(lineColor[fv]);
196 }
197
198} // end of void Lines()
199
200// =====================================================================================
201
202void Print()
203{
204 // Print number of events and average multiplicities for each common output file.
205
206 cout<<endl;
207 cout<<"Accessed files:"<<endl;
208 cout<<endl;
209 for(Int_t f=0;f<nFiles;f++)
210 {
211 cout<<commonOutputFiles[f]->GetName()<<endl;
212 for(Int_t m=0;m<nMethods;m++)
213 {
214 AliFlowCommonHist *commonHist = dynamic_cast<AliFlowCommonHist*> (lists[f][m]->FindObject(Form("AliFlowCommonHist%s",method[m].Data())));
215 Double_t nEvts = -1.;
216 Double_t AvM = -1.;
217 if(commonHist && commonHist->GetHistMultRP())
218 {
219 nEvts = commonHist->GetHistMultRP()->GetEntries();
220 AvM = commonHist->GetHistMultRP()->GetMean();
221 }
222 if(!(strcmp(method[m].Data(),"QC")))
223 {
224 cout<<Form("%s:",method[m].Data())<<" <M> = "<<AvM<<", N = "<<nEvts<<endl;
225 }
226 if(!(strcmp(method[m].Data(),"GFC")) && showAlsoGFCResults)
227 {
228 cout<<Form("%s:",method[m].Data())<<" <M> = "<<AvM<<", N = "<<nEvts<<endl;
229 }
230 }
231 cout<<endl;
232 } // end of for(Int_t f=0;f<nFiles;f++)
233
234} // end of void Print()
235
236// =====================================================================================
237
238void DetermineMinMax()
239{
240 // Determine ranges for plots.
241
242 for(Int_t co=0;co<4;co++)
243 {
244 xMin[co] = 0.; yMin[co] = 44.;
245 xMax[co] = -440000.; yMax[co] = -44.;
246 }
247
248 Double_t tfc[nFlowValues][4] = {{0.}}; // theoretical flow contribution
249 for(Int_t fv=0;fv<nFlowValues;fv++)
250 {
251 tfc[fv][0] = pow(v[fv],2);
252 tfc[fv][1] = -pow(v[fv],4);
253 tfc[fv][2] = 4.*pow(v[fv],6);
254 tfc[fv][3] = -33.*pow(v[fv],8);
255 }
256
257 for(Int_t f=0;f<nFiles;f++)
258 {
259 for(Int_t m=0;m<nMethods;m++)
260 {
261 for(Int_t co=0;co<4;co++)
262 {
263 if(cumulantsVsM[f][m][co])
264 {
265 for(Int_t b=1;b<=cumulantsVsM[f][m][co]->GetXaxis()->GetNbins();b++)
266 {
267 Double_t result = cumulantsVsM[f][m][co]->GetBinContent(b);
268 Double_t error = cumulantsVsM[f][m][co]->GetBinError(b);
269 if(TMath::Abs(result)>1.e-44 && TMath::Abs(error)>1.e-44)
270 {
271 // y-axis:
272 if(yMin[co] > result-error){yMin[co] = result-error;} // min value
273 if(yMax[co] < result+error) {yMax[co] = result+error;} // max value
274 // x-axis:
275 xMax[co] = b;
276 }
277 } // end of for(Int_t b=1;b<=cumulantsVsM[f][m][co]->GetXaxis()->GetNbins();b++)
278 // theoretical contributions:
279 for(Int_t fv=0;fv<nFlowValues;fv++)
280 {
281 //if(yMin[co] > tfc[fv][0]) {yMin[co] = tfc[fv][0];} // min value
282 //if(yMax[co] < tfc[fv][0]) {yMax[co] = tfc[fv][0];} // max value
283 } // end of for(Int_t fv=0;fv<nFlowValues;fv++)
284 } // end of if(cumulantsVsM[f][m][co])
285 } // end of for(Int_t co=0;co<4;co++)
286 } // end of for(Int_t m=0;m<nMethods;m++)
287 } // end of for(Int_t f=0;f<nFiles;f++)
288
289} // end of void DetermineMinMax()
290
291// =====================================================================================
292
293void GetHistograms()
294{
295 // Get histograms with results for cumulants vs multiplicity.
296
297 TString qcFlag[4] = {"QC{2}","QC{4}","QC{6}","QC{8}"};
298 TString gfcFlag[4] = {"GFC{2}","GFC{4}","GFC{6}","GFC{8}"};
299 for(Int_t f=0;f<nFiles;f++)
300 {
301 for(Int_t m=0;m<nMethods;m++)
302 {
303 TList *temp = NULL;
304 if(!(strcmp(method[m].Data(),"QC")))
305 {
306 temp = dynamic_cast<TList*> (lists[f][m]->FindObject("Integrated Flow"));
307 if(temp) {temp = dynamic_cast<TList*> (temp->FindObject("Results"));}
308 if(temp)
309 {
310 for(Int_t co=0;co<4;co++)
311 {
312 cumulantsVsM[f][m][co] = dynamic_cast<TH1D*> (temp->FindObject(Form("fIntFlowQcumulantsVsM, %s",qcFlag[co].Data())));
313 }
314 }
315 } // end of if(!(strcmp(method[m].Data(),"QC")))
316 else if(!(strcmp(method[m].Data(),"GFC")))
317 {
318 temp = dynamic_cast<TList*> (lists[f][m]->FindObject("Reference Flow"));
319 if(temp) {temp = dynamic_cast<TList*> (temp->FindObject("Results"));}
320 if(temp)
321 {
322 for(Int_t co=0;co<4;co++)
323 {
324 cumulantsVsM[f][m][co] = dynamic_cast<TH1D*> (temp->FindObject(Form("fReferenceFlowCumulantsVsM, %s",gfcFlag[co].Data())));
325 }
326 }
327 } // end of else if(!(strcmp(method[m].Data(),"QC")))
328 } // end of for(Int_t m=0;m<nMethods;m++)
329 } // end of for(Int_t f=0;f<nFiles;f++)
330
331} // end of void GetHistograms()
332
333// =====================================================================================
334
335TGraphErrors* GetGraphErrors(Int_t bin, Int_t nFiles, TH1D** qc)
336{
337 TGraphErrors *ge = new TGraphErrors(nFiles);
338 for(Int_t f=0;f<nFiles;f++)
339 {
340 ge->SetPoint(f,f+0.5,qc[f]->GetBinContent(bin+1));
341 ge->SetPointError(f,0,qc[f]->GetBinError(bin+1));
342 }
343
344 return ge;
345
346} // end of TGraphErrors* GetGraphErrors(Int_t bin, Int_t nFiles, TH1D** qc)
347
348// =====================================================================================
349
350TH1D* StyleHist(TString yAxisTitle, Int_t co)
351{
352 // Style histogram.
353
354 TH1D *styleHist = new TH1D("","",10000,0,10000); // to be improved (hardwired 10000)
355 // x-axis:
356 styleHist->GetXaxis()->SetRangeUser(xMin[co],xMax[co]);
357 // y-axis:
358 styleHist->GetYaxis()->SetRangeUser(yMin[co],yMax[co]);
359
360 styleHist->GetXaxis()->SetTitle("M");
361 styleHist->GetYaxis()->SetTitle(yAxisTitle.Data());
362
363 return styleHist;
364
365} // end of TH1D* StyleHist(TString yAxisTitle, Int_t co)
366
367// ===========================================================================================
368
369TGraph* GetErrorMesh(TH1D *hist)
370{
371 // Error mesh.
372
373 TGraph* errorMesh = NULL;
374 if(hist)
375 {
376 Int_t nBins = hist->GetNbinsX();
377 Double_t binWidth = hist->GetBinWidth(1); // assuming that all bins have the same width
378 // Counting the non-empty bins:
379 Int_t nNonEmptyBins=0;
380 for(Int_t i=1;i<nBins+1;i++)
381 {
382 if(!(hist)->GetBinError(i)==0.0))
383 {
384 nNonEmptyBins++;
385 }
386 }
387 errorMesh = new TGraph(2*nNonEmptyBins+1);
388 Double_t value=0.,error=0.;
389 Int_t count=1;
390 Double_t xFirst=0.,yUpFirst=0.; // needed to close up the mesh
391 for(Int_t i=1;i<nBins+1;i++)
392 {
393 // Setting up the upper limit of the mesh:
394 value = hist->GetBinContent(i);
395 error = hist->GetBinError(i);
396 if(!(error==0.0))
397 {
398 errorMesh->SetPoint(count++,(i-0.5)*binWidth,value+error);
399 if(xFirst==0.)
400 {
401 xFirst=(i-0.5)*binWidth;
402 yUpFirst=value+error;
403 }
404 }
405 }
406 for(Int_t i=nBins+1;i<2*nBins+1;i++)
407 {
408 // Setting up the lower limit of the mesh:
409 value = hist->GetBinContent(2*nBins+1-i);
410 error = hist->GetBinError(2*nBins+1-i);
411 if(!(error==0.0))
412 {
413 errorMesh->SetPoint(count++,(2*nBins-i+0.5)*binWidth,value-error);
414 }
415 }
416 // Closing the mesh area:
417 errorMesh->SetPoint(2*nNonEmptyBins+1,xFirst,yUpFirst);
418 } // end if(hist)
419
420 errorMesh->SetFillStyle(1001);
421
422 return errorMesh;
423
424} // end of TGraph* GetErrorMesh(TH1D *hist)
425
426// ===========================================================================================
427
428void GlobalSettings()
429{
430 // Settings which will affect all plots.
431
432 gROOT->SetStyle("Plain"); // default color is white instead of gray
433 gStyle->SetOptStat(0); // remove stat. box from all histos
434
435} // end of void GlobalSettings()
436
437// ===========================================================================================
438
439void GetLists()
440{
441 // Get from common output files the lists holding histograms for each method.
442
443 TString fileName[nFiles][nMethods];
444 TDirectoryFile *dirFile[nFiles][nMethods] = {{NULL}};
445 TString listName[nFiles][nMethods];
446 for(Int_t f=0;f<nFiles;f++)
447 {
448 for(Int_t i=0;i<nMethods;i++)
449 {
450 // Form a file name for each method:
451 fileName[f][i]+="output";
452 fileName[f][i]+=method[i].Data();
453 fileName[f][i]+="analysis";
454 fileName[f][i]+=type[f].Data();
455 // Access this file:
456 if(commonOutputFiles[f]){dirFile[f][i] = (TDirectoryFile*)commonOutputFiles[f]->FindObjectAny(fileName[f][i].Data());}
457 // Form a list name for each method:
458 listName[f][i]+="cobj";
459 listName[f][i]+=method[i].Data();
460 // Access this lists:
461 if(dirFile[f][i])
462 {
463 dirFile[f][i]->GetObject(listName[f][i].Data(),lists[f][i]);
464 } else
465 {
466 cout<<"WARNING: Couldn't find a file "<<fileName[f][i].Data()<<".root in "<<commonOutputFiles[f]->GetName()<<" !!!!"<<endl;exit(0);
467 }
468 } // end of for(Int_t i=0;i<nMethods;i++)
469 } // end of for(Int_t f=0;f<nFiles;f++)
470
471} // end of void GetLists()
472
473// ===========================================================================================
474
475void AccessCommonOutputFiles(TString commonOutputFileName)
476{
477 // Access all output files.
478 for(Int_t f=0;f<nFiles;f++)
479 {
480 if(!(gSystem->AccessPathName(Form("%s/%s/%s",gSystem->pwd(),files[f].Data(),commonOutputFileName.Data()),kFileExists)))
481 {
482 commonOutputFiles[f] = TFile::Open(Form("%s/%s/%s",gSystem->pwd(),files[f].Data(),commonOutputFileName.Data()),"READ");
483 } else
484 {
485 cout<<endl;
486 cout<<"WARNING: Couldn't find the file "<<Form("%s/%s/%s",gSystem->pwd(),files[f].Data(),commonOutputFileName.Data())<<" !!!!"<<endl;
487 cout<<endl;
488 exit(0);
489 }
490 } // end of for(Int_t f=0;f<nFiles;f++)
491
492} // void AccessCommonOutputFiles(TString commonOutputFileName);
493
494// ===========================================================================================
495
496void LoadLibrariesPC(const libModes analysisMode) {
497
498 //--------------------------------------
499 // Load the needed libraries most of them already loaded by aliroot
500 //--------------------------------------
501 //gSystem->Load("libTree");
502 gSystem->Load("libGeom");
503 gSystem->Load("libVMC");
504 gSystem->Load("libXMLIO");
505 gSystem->Load("libPhysics");
506
507 //----------------------------------------------------------
508 // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
509 //----------------------------------------------------------
510 if (analysisMode==mLocal) {
511 //--------------------------------------------------------
512 // If you want to use already compiled libraries
513 // in the aliroot distribution
514 //--------------------------------------------------------
515
516 //==================================================================================
517 //load needed libraries:
518 gSystem->AddIncludePath("-I$ROOTSYS/include");
519 //gSystem->Load("libTree");
520
521 // for AliRoot
522 gSystem->AddIncludePath("-I$ALICE_ROOT/include");
523 //gSystem->Load("libANALYSIS");
524 gSystem->Load("libPWG2flowCommon");
525 //cerr<<"libPWG2flowCommon loaded ..."<<endl;
526
527 }
528
529 else if (analysisMode==mLocalSource) {
530
531 // In root inline compile
532
533 // Constants
534 gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
535 gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
536 gROOT->LoadMacro("AliFlowCommon/AliFlowCumuConstants.cxx+");
537
538 // Flow event
539 gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+");
540 gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");
541 gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
542
543 // Cuts
544 gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");
545
546 // Output histosgrams
547 gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
548 gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
549 gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
550 gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
551
552 cout << "finished loading macros!" << endl;
553
554 }
555
556} // end of void LoadLibraries(const libModes analysisMode)