hooks for PMD flow analysis
[u/mrichter/AliRoot.git] / PWG2 / FLOW / macros / showSpread.C
1 // ... 
2
3 enum libModes {mLocal,mLocalSource};
4
5 const Int_t nFilesMax = 10; // number of files to be accessed to estimate spread
6
7 Bool_t showPlotForReferenceFlow = kTRUE;
8 Bool_t showPlotForIntegratedFlowRP = kTRUE;
9 Bool_t showPlotForIntegratedFlowPOI = kTRUE;
10 Bool_t showErrorOnMergedResult = kTRUE; 
11
12 const Int_t nMethods = 13;
13 TString method[nMethods] = {"MCEP","SP","2,GFC","2,QC","4,GFC","4,QC","6,GFC","6,QC","8,GFC","8,QC","FQD","LYZ1SUM","LYZ1PROD"};
14 Int_t methodMarkerStyle[nMethods] = {21,21,21,21,21,21,21,21,21,21,21,21,21};
15 Int_t methodMarkerColor[nMethods] = {kGray+1,kViolet-8,kBlue-9,kRed-7,kBlue-9,kRed-7,kBlue-9,kRed-7,kBlue-9,kRed-7,kOrange-8,kYellow-5,kYellow-2};
16 Int_t methodMeshColor[nMethods] = {kGray,kViolet-9,kBlue-10,kRed-10,kBlue-10,kRed-10,kBlue-10,kRed-10,kBlue-10,kRed-10,kOrange-9,kYellow-8,kYellow-6};
17
18 /*
19 const Int_t nMethods = 4;
20 TString method[nMethods] = {"2,QC","4,QC","6,QC","8,QC"};
21 Int_t methodMarkerStyle[nMethods] = {21,21,21,21};
22 Int_t methodMarkerColor[nMethods] = {kRed-7,kRed-7,kRed-7,kRed-7};
23 Int_t methodMeshColor[nMethods] = {kRed-10,kRed-10,kRed-10,kRed-10};
24 */
25
26 void showSpread(TString type="", Int_t mode=mLocal)
27 {
28  // type: type of analysis can be ESD, AOD, MC, ESDMCkineESD, ESDMCkineMC
29  //       (if type="" output files are from MC simulation (default))
30  // mode: if mode = mLocal: analyze data on your computer using aliroot
31  //       if mode = mLocalSource: analyze data on your computer using root + source files 
32
33  // Cross-check if the user's settings make sense:
34  CrossCheckUserSettings();
35  
36  // Load needed libraries:                       
37  LoadLibrariesSS(mode);  
38   
39  // Output file name:
40  TString outputFileName = "AnalysisResults.root"; 
41  
42  // Labels for reference flow, integrated flow of RPs and of POIs:
43  TString label[3] = {"","RP","POI"};
44   
45  // Standard magic:
46  TString *baseDirPath = new TString(gSystem->pwd());
47  TSystemDirectory *baseDir = new TSystemDirectory(".",baseDirPath->Data());          
48  TList *listOfFilesInBaseDir = baseDir->GetListOfFiles();
49  TStopwatch timer;
50  timer.Start();
51  // listOfFilesInBaseDir->Print();
52  Int_t nFiles = listOfFilesInBaseDir->GetEntries();
53  Int_t fileCounter = 0;
54  Double_t result[nMethods][3][nFilesMax] = {{{0.}}}; // [3 = "", "RP" or "POI"]
55  Double_t error[nMethods][3][nFilesMax] = {{{0.}}}; // [3 = "", "RP" or "POI"]
56  Double_t resultMinMax[nMethods][3][2] = {{{0.}}}; // [3 = "", "RP" or "POI"], [2 = min value, max value]
57  for(Int_t m=0;m<nMethods;m++)
58  {
59   for(Int_t l=0;l<3;l++)
60   {
61    resultMinMax[m][l][0] = 44.;
62    resultMinMax[m][l][1] = -44.;
63   }
64  } 
65  Double_t styleHistMinMax[3][2] = {{0.}}; // [3 = "", "RP" or "POI"], [2 = min value, max value]
66  for(Int_t l=0;l<3;l++)
67  {
68   styleHistMinMax[l][0] = 44.;
69   styleHistMinMax[l][1] = -44.;
70  }  
71  cout<<endl;
72  for(Int_t iFile=0;iFile<nFiles;iFile++)
73  {
74   TSystemFile *currentFile = (TSystemFile*)listOfFilesInBaseDir->At(iFile);
75   // Consider only subdirectories: 
76   if(!currentFile || 
77      !currentFile->IsDirectory() || 
78      strcmp(currentFile->GetName(), ".") == 0 || 
79     strcmp(currentFile->GetName(), "..") == 0) continue; 
80   // Accessing the output file "AnalysisResults.root" in current subdirectory: 
81   TString currentSubDirName = baseDirPath->Data();
82   (currentSubDirName+="/")+=currentFile->GetName();
83   currentSubDirName+="/";
84   TString fileName = currentSubDirName; 
85   fileName+=outputFileName.Data();
86   if(!(gSystem->AccessPathName(fileName.Data(),kFileExists)))
87   {
88    TFile *file = NULL; 
89    file = TFile::Open(fileName.Data(),"READ");
90    for(Int_t m=0;m<nMethods;m++)
91    {
92     // Access from common output file the output file for particular method:
93     TDirectoryFile *methodFile = NULL;
94     methodFile = GetMethodFile(file,method[m],type);  
95     for(Int_t l=0;l<3;l++)
96     {
97      TH1D *histResult = NULL;
98      histResult = GetHistogramWithResult(method[m],label[l],methodFile);
99      if(histResult)
100      {
101       // Access the results:
102       result[m][l][fileCounter] = histResult->GetBinContent(1); 
103       // Access the errors:
104       error[m][l][fileCounter] = histResult->GetBinError(1);
105       if(TMath::Abs(result[m][l][fileCounter])>pow(10.,-6.)) // take into account only if != 0 (to be improved - special care for < 0 is required)
106       {
107        // Establish min and max values for results:
108        if(resultMinMax[m][l][0] > result[m][l][fileCounter]) resultMinMax[m][l][0] = result[m][l][fileCounter]; // min value
109        if(resultMinMax[m][l][1] < result[m][l][fileCounter]) resultMinMax[m][l][1] = result[m][l][fileCounter]; // max value
110        // Establish min and max values for style histograms:
111        if(styleHistMinMax[l][0] > result[m][l][fileCounter]) styleHistMinMax[l][0] = result[m][l][fileCounter]; // min value
112        if(styleHistMinMax[l][1] < result[m][l][fileCounter]) styleHistMinMax[l][1] = result[m][l][fileCounter]; // max value    
113       }
114      }
115     } // end of for(Int_t l=0;l<3;l++)
116    } // end of for(Int_t m=0;m<nMethods;m++)
117    if(file) file->Close();
118    fileCounter++;
119    //if(fileCounter%10==0)
120    //{
121     cout<<Form("Accessed %d files \"AnalysisResults.root\" so far....",fileCounter)<<"\r"<<flush;
122    //}   
123   } 
124   if(fileCounter == nFilesMax) break;
125  } // end of for(Int_t iFile=0;iFile<nFiles;iFile++)   
126   
127  cout<<Form("Accessed %d files \"AnalysisResults.root\" in total to estimate spread. ",fileCounter)<<endl;
128  cout<<endl;
129  if(fileCounter==0){exit(0);}
130  const Int_t nFilesFinal = fileCounter;
131
132  // Make for each method graph holding results:
133  TGraph *methodGraph[nMethods][3] = {{NULL}}; // [3 = "", "RP" or "POI"] 
134  Double_t x[nMethods][nFilesFinal] = {{0.}};
135  for(Int_t m=0;m<nMethods;m++)
136  {
137   for(Int_t f=0;f<nFilesFinal;f++)
138   {
139    x[m][f]=m+0.5;
140   } 
141  }
142  for(Int_t m=0;m<nMethods;m++)
143  {
144   for(Int_t l=0;l<3;l++)
145   {
146    methodGraph[m][l] = new TGraph(nFilesFinal,x[m],result[m][l]);
147    methodGraph[m][l]->SetMarkerStyle(methodMarkerStyle[m]);
148    methodGraph[m][l]->SetMarkerColor(methodMarkerColor[m]);
149   } // end of for(Int_t l=0;l<3;l++)
150  } // for(Int_t m=0;m<nMethods;m++)
151  
152  // Make for each method coloured mesh out of min and max values:
153  Double_t meshWidth = 0.25;
154  TGraph *methodMesh[nMethods][3] = {{NULL}}; // [3 = "", "RP" or "POI"]
155  for(Int_t m=0;m<nMethods;m++)
156  {
157   for(Int_t l=0;l<3;l++)
158   {
159    if(resultMinMax[m][l][0]<44. && resultMinMax[m][l][1]>-44.)
160    {
161     methodMesh[m][l] = new TGraph(5);
162     methodMesh[m][l]->SetPoint(0,(m+1-0.5)-meshWidth,resultMinMax[m][l][0]);
163     methodMesh[m][l]->SetPoint(1,(m+1-0.5)+meshWidth,resultMinMax[m][l][0]);
164     methodMesh[m][l]->SetPoint(2,(m+1-0.5)+meshWidth,resultMinMax[m][l][1]);
165     methodMesh[m][l]->SetPoint(3,(m+1-0.5)-meshWidth,resultMinMax[m][l][1]);
166     methodMesh[m][l]->SetPoint(4,(m+1-0.5)-meshWidth,resultMinMax[m][l][0]);    
167     methodMesh[m][l]->SetFillStyle(1001);
168     methodMesh[m][l]->SetFillColor(methodMeshColor[m]);
169    } 
170   }
171  }
172   
173  // Access for each method the results from the merged, large statistics file:
174  Double_t resultMerged[nMethods][3] = {{0.}}; // [3 = "", "RP" or "POI"]
175  TString mergedFileName = Form("%s%s%s",gSystem->pwd(),"/",outputFileName.Data());
176  TFile *mergedFile = NULL; 
177  mergedFile = TFile::Open(mergedFileName.Data(),"READ"); 
178  for(Int_t m=0;m<nMethods;m++)
179  {
180   TDirectoryFile *methodFile = NULL;
181   if(!(gSystem->AccessPathName(fileName.Data(),kFileExists)))
182   {
183    if(mergedFile) methodFile = GetMethodFile(mergedFile,method[m],type);  
184   } else
185     {
186      cout<<"WARNING: Couldn't find the merged, large statistics file "<<endl;
187      cout<<"         "<<fileName.Data()<<endl;
188      cout<<"         in directory "<<gSystem->pwd()<<" !!!!"<<endl;     
189      cout<<"         Use macros mergeOuput.C and redoFinish.C to get it."<<endl;
190      cout<<endl;
191      break;
192     }
193   for(Int_t l=0;l<3;l++)
194   {
195    TH1D *histResult = NULL;
196    if(methodFile)
197    { 
198     histResult = GetHistogramWithResult(method[m],label[l],methodFile);
199    }
200    if(histResult)
201    {
202     // Access the results from the merged, large statistics file:
203     resultMerged[m][l] = histResult->GetBinContent(1);
204    } // end of for(Int_t l=0;l<3;l++)
205   }
206  } // end of for(Int_t m=0;m<nMethods;m++)  
207  if(mergedFile) mergedFile->Close();
208         
209  // Make for each method graph holding results from the merged, large statistics file
210  // and the errors from the randomly chosen small statistics file:
211  TGraphErrors *methodMergedGraph[nMethods][3] = {{NULL}}; // [3 = "", "RP" or "POI"] 
212  Double_t xMerged[nMethods] = {0.};
213  for(Int_t m=0;m<nMethods;m++)
214  {
215   xMerged[m]=m+0.5; 
216  }
217  // Select randomly small statistics file:
218  gRandom->SetSeed((UInt_t) (4400*timer.RealTime()/fileCounter));
219  Int_t randomFile = (Int_t)gRandom->Uniform(0,fileCounter); 
220  for(Int_t m=0;m<nMethods;m++)
221  {
222   for(Int_t l=0;l<3;l++)
223   {
224    methodMergedGraph[m][l] = new TGraphErrors(1);
225    methodMergedGraph[m][l]->SetPoint(0,xMerged[m],resultMerged[m][l]);
226    if(showErrorOnMergedResult) methodMergedGraph[m][l]->SetPointError(0,0.,error[m][l][randomFile]);
227    methodMergedGraph[m][l]->SetMarkerStyle(25);
228    methodMergedGraph[m][l]->SetMarkerColor(kBlack);
229   } // end of for(Int_t l=0;l<3;l++)
230  } // for(Int_t m=0;m<nMethods;m++)
231  
232  // Final drawing:
233  gROOT->SetStyle("Plain"); // removing default gray color and setting white instead
234  gStyle->SetOptStat(0); // removing statistics box from all histograms
235  Bool_t showPlot[3] = {showPlotForReferenceFlow,showPlotForIntegratedFlowRP,showPlotForIntegratedFlowPOI};
236  TString title[3] = {"Reference Flow","Integrated Flow (RP)","Integrated Flow (POI)"}
237  TCanvas *canvas[3] = {NULL}; // [3 = "", "RP" or "POI"]
238  for(Int_t l=0;l<3;l++)
239  {
240   if(showPlot[l])
241   {
242    canvas[l] = new TCanvas(Form("%s",title[l].Data()),Form("%s",title[l].Data()));
243    TH1D *styleHist = StyleHist(title[l]);
244    styleHist->SetMinimum(0.99*styleHistMinMax[l][0]);
245    styleHist->SetMaximum(1.01*styleHistMinMax[l][1]);
246    styleHist->Draw();
247    for(Int_t m=0;m<nMethods;m++)
248    {
249     if(methodMesh[m][l]) methodMesh[m][l]->Draw("lfsame");
250     if(methodGraph[m][l]) methodGraph[m][l]->Draw("psame");  
251     if(TMath::Abs(*(methodMergedGraph[m][l]->GetY()))>pow(10.,-6.)) // draw only if not == 0.
252     {
253      methodMergedGraph[m][l]->Draw("psame");
254     } 
255    } // end of for(Int_t m=0;m<nMethods;m++)
256   } // end of if(showPlot[l])
257  } // end of for(Int_t l=0;l<3;l++)
258  
259  timer.Stop();
260  timer.Print(); 
261  cout<<endl;
262  
263 } // end of void showSpread(TString type="", Int_t mode=mLocal)
264
265 // =============================================================================================
266
267 TH1D *StyleHist(TString histTitle)
268 {
269  // Make the style histogram.
270  Int_t n = 2; // harmonic (to be improved - access this from common control hist)
271  TH1D *styleHist = new TH1D(Form("%s",histTitle.Data()),Form("%s",histTitle.Data()),nMethods,0,nMethods);
272  styleHist->GetYaxis()->SetTickLength(0.01);
273  for(Int_t m=0;m<nMethods;m++)
274  {
275   styleHist->GetXaxis()->SetBinLabel(m+1,Form("v_{%d}{%s}",n,method[m].Data()));
276   if(method[m]=="LYZ1SUM" || method[m]=="LYZ2SUM")
277   {
278    styleHist->GetXaxis()->SetBinLabel(m+1,Form("v_{%d}{%s}",n,"LYZ,sum"));
279   } 
280   else if(method[m]=="LYZ1PROD" || method[m]=="LYZ2PROD")
281   {
282    styleHist->GetXaxis()->SetBinLabel(m+1,Form("v_{%d}{%s}",n,"LYZ,prod"));
283   } 
284  } // end of for(Int_t m=0;m<nMethods;m++)
285  
286  return styleHist;
287 }
288
289 // =============================================================================================
290
291 TDirectoryFile* GetMethodFile(TFile *commonFile, TString method, TString type)
292 {
293  // Form a file name for each method:
294  TString methodFileName = "output";
295  if(method.Contains("GFC"))
296  {
297   methodFileName+="GFC";
298  } else if(method.Contains("QC"))
299    {
300     methodFileName+="QC";
301    } else
302      {
303       methodFileName+=method.Data();
304      } 
305  methodFileName+="analysis";
306  methodFileName+=type.Data();
307  TDirectoryFile *methodFile = NULL;
308  if(commonFile)
309  {
310   methodFile = (TDirectoryFile*)commonFile->FindObjectAny(methodFileName.Data());
311  } 
312  
313  return methodFile;
314
315 } // end of TDirectoryFile* AccessMethodFile(TString commonFile, TString method, TString type)
316
317 // =============================================================================================
318
319 TH1D* GetHistogramWithResult(TString method, TString label, TDirectoryFile *methodFile)
320 {
321  // Access first the common list holding all output histograms:
322  TList *methodList = NULL;
323  if(method.Contains("GFC"))
324  {
325   methodFile->GetObject("cobjGFC",methodList);
326  } else if(method.Contains("QC"))
327    {
328     methodFile->GetObject("cobjQC",methodList);
329    } else
330      {
331       methodFile->GetObject(Form("cobj%s",method.Data()),methodList);
332      } 
333  // Access from the common list the needed histogram:
334  if(methodList)
335  {
336   AliFlowCommonHistResults *commonHistRes = NULL; 
337   if(!(method.Contains("GFC") || method.Contains("QC")))
338   {
339    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject(Form("AliFlowCommonHistResults%s",method.Data()));
340   } 
341   else if(method=="2,GFC")
342   {
343    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject("AliFlowCommonHistResults2ndOrderGFC");    
344   } 
345   else if(method=="4,GFC")
346   {
347    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject("AliFlowCommonHistResults4thOrderGFC");    
348   } 
349   else if(method=="6,GFC")
350   {
351    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject("AliFlowCommonHistResults6thOrderGFC");    
352   } 
353   else if(method=="8,GFC")
354   {
355    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject("AliFlowCommonHistResults8thOrderGFC");    
356   }    
357   else if(method=="2,QC")
358   {
359    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject("AliFlowCommonHistResults2ndOrderQC");    
360   } 
361   else if(method=="4,QC")
362   {
363    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject("AliFlowCommonHistResults4thOrderQC");    
364   } 
365   else if(method=="6,QC")
366   {
367    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject("AliFlowCommonHistResults6thOrderQC");    
368   } 
369   else if(method=="8,QC")
370   {
371    commonHistRes = dynamic_cast<AliFlowCommonHistResults*> methodList->FindObject("AliFlowCommonHistResults8thOrderQC");    
372   }  
373  } // end of if(methodList)
374  
375  // Access histogram with results for reference flow or integrated flow of RPs or POIs:
376  TH1D *hist = NULL;
377  if(label=="")
378  {
379   hist = commonHistRes->GetHistIntFlow();
380  } else if(label=="RP")
381    {
382     hist = commonHistRes->GetHistIntFlowRP(); 
383    } else if(label=="POI")
384      {
385       hist = commonHistRes->GetHistIntFlowPOI(); 
386      }
387
388  if(hist) return hist;
389
390 } // end of TH1D *GetHistogramWithResult(TString method, TString rf_rp_poi, TString file);
391
392 // =============================================================================================
393
394 void CrossCheckUserSettings() 
395 {
396  // Check in this method if the user settings make sense.
397  
398  if(nFilesMax<=0)
399  {
400   cout<<endl;
401   cout<<"WARNING: nFilesMax must be a positive integer (not too large, though) !!!!"<<endl;
402   cout<<endl;
403   exit(0);
404  }
405  if(nFilesMax>44)
406  { 
407   cout<<endl;
408   cout<<"WARNING: You may want to set nFilesMax to the smaller value."<<endl;
409   cout<<"         Otherwise you might wait forever to see the plots."<<endl;
410  } 
411  
412 } // end of void CrossCheckUserSettings()
413
414 // =============================================================================================
415
416 void LoadLibrariesSS(const libModes mode) {
417   
418   //--------------------------------------
419   // Load the needed libraries most of them already loaded by aliroot
420   //--------------------------------------
421   //gSystem->Load("libTree");
422   gSystem->Load("libGeom");
423   gSystem->Load("libVMC");
424   gSystem->Load("libXMLIO");
425   gSystem->Load("libPhysics");
426   
427   //----------------------------------------------------------
428   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<< 
429   //----------------------------------------------------------
430   if (mode==mLocal) {
431     //--------------------------------------------------------
432     // If you want to use already compiled libraries 
433     // in the aliroot distribution
434     //--------------------------------------------------------
435
436   //==================================================================================  
437   //load needed libraries:
438   gSystem->AddIncludePath("-I$ROOTSYS/include");
439   //gSystem->Load("libTree");
440
441   // for AliRoot
442   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
443   gSystem->Load("libANALYSIS");
444   gSystem->Load("libPWG2flowCommon");
445   //cerr<<"libPWG2flowCommon loaded ..."<<endl;
446   
447   }
448   
449   else if (mode==mLocalSource) {
450  
451     // In root inline compile
452   
453     // Constants  
454     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
455     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
456         
457     // Flow event
458     gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+"); 
459     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");    
460     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");    
461     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
462         
463     // Output histograms
464     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
465     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
466     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
467     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
468     
469     // Functions needed for various methods
470     gROOT->LoadMacro("AliFlowCommon/AliCumulantsFunctions.cxx+");
471     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZEventPlane.cxx+");
472     
473     // Flow Analysis code for various methods
474     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMCEventPlane.cxx+"); 
475     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx+");
476     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.cxx+");
477     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.cxx+");
478     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithCumulants.cxx+");
479     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithQCumulants.cxx+"); 
480     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithFittingQDistribution.cxx+");
481     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMixedHarmonics.cxx+");    
482     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithNestedLoops.cxx+");          
483     
484     cout << "finished loading macros!" << endl;  
485     
486   } // end of else if (mode==mLocalSource) 
487   
488 } // end of void LoadLibrariesSS(const libModes mode)