hooks for PMD flow analysis
[u/mrichter/AliRoot.git] / PWG2 / FLOW / macros / fqd.C
1 //==========================================================================//
2 // Macro to refit the q-distribution in FQD method. By using this macro the // 
3 // fitting parameters for q-distribution can be tuned and correspondingly   //
4 // the estimate for reference flow harmonic obtained from FQD method can be //
5 // significantly improved.                                                  // 
6 //==========================================================================//
7
8 // Fitting parameters for q-distribution:
9 Double_t vStart = 0.05; // starting value for v fit
10 Double_t vMin = 0.0; // lower bound for v fit 
11 Double_t vMax = 0.25; // upper bound for v fit  
12 Double_t sigma2Start = 0.75; // starting value for sigma^2 fit 
13 Double_t sigma2Min = 0.5; // lower bound for sigma^2 fit (according to theorists must be >= 0.5)   
14 Double_t sigma2Max = 2.5; // upper bound for sigma^2 fit
15 Double_t treshold = 5.; // first and last bin taken for fitting are determined as the first and last bin with more than 'treshold' number of entries
16 Bool_t finalResultIsFromSigma2Fitted = kTRUE; // final saved result is obtained with sigma^2 fitted (kTRUE) or sigma^2 fixed (kFALSE)
17 Bool_t printOnTheScreen = kTRUE; // print or not the final results on the screen
18 Bool_t plotResults = kTRUE; // plot on the screen q-distribution and resulting fitting functions (for non-fixed sigma^2 and fixed sigma^2)
19
20 // Name of the common output file:
21 TString outputFileName = "AnalysisResults.root";
22 //TString outputFileName = "outputCentrality0.root";
23
24 enum libModes {mLocal,mLocalSource};
25
26 TFile *commonOutputFile = NULL; // common output file "AnalysisResults.root"
27 TDirectoryFile *dirFileFQD = NULL; // FQD's TDirectoryFile in "AnalysisResults.root"
28 TList *outputListFQD = NULL; // output list holding all FQD objects
29
30 void fqd(TString analysisType="ESD", Int_t analysisMode=mLocal)
31 {
32  // 1. analysisType: type of analysis can be ESD, AOD, MC, MK, ....
33  //                  (if type="" output files are from simulation 'on the fly')
34  // 2. analysisMode: if analysisMode = mLocal -> analyze data on your computer using aliroot
35  //                  if analysisMode = mLocalSource -> analyze data on your computer using root + source files 
36   
37  // Cross-check if the user's settings make sense:
38  CrossCheckSettings();
39  // Load needed libraries:                       
40  LoadLibrariesFQD(analysisMode); 
41  // Get output list which holds all FQD objects:
42  GetOutputList(analysisType); 
43  // Redo fit of q-distribution:
44  RedoFit();
45  // Plot q-distribution and resulting fitting functions:
46  if(plotResults) Plot(); 
47  // Save the new results in the common output file: 
48  Save();
49   
50 } // end of void fqd(TString type="ESD", Int_t mode=mLocal)   
51
52 // =========================================================================================== 
53
54 void RedoFit()
55 {
56  // Redo the fit of q-distribution.
57
58  AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
59  fqd->GetOutputHistograms(outputListFQD);
60  // Set new fitting parameters:
61  fqd->SetvStart(vStart);  
62  fqd->SetvMin(vMin);
63  fqd->SetvMax(vMax);
64  fqd->SetSigma2Start(sigma2Start);  
65  fqd->SetSigma2Min(sigma2Min);
66  fqd->SetSigma2Max(sigma2Max);
67  fqd->SetTreshold(treshold);
68  fqd->SetFinalResultIsFromSigma2Fitted(finalResultIsFromSigma2Fitted);
69  fqd->SetPrintOnTheScreen(printOnTheScreen); 
70  // Save new fitting parameters:  
71  fqd->StoreFittingParameters();
72  // Redo fit:
73  fqd->Finish(kTRUE);
74  
75 } // end of void RedoFit(TList *outputList)
76  
77 // =========================================================================================== 
78  
79 void GetOutputList(TString analysisType)
80 {
81  // Get output list which holds all FQD objects.
82   
83  // Access common output file:
84  commonOutputFile = AccessOutputFile(outputFileName);
85  
86  // Access from common output file the TDirectoryFile for FQD method
87  // and from it the output list holding all objects:
88  GetListWithHistograms(commonOutputFile,analysisType,"FQD");
89    
90 } // end of TList* GetOutputList(TString analysisType)
91  
92 // ===========================================================================================
93
94 void Plot()
95 {
96  // Plot q-distribution and resulting fitting functions.
97  
98  gROOT->SetStyle("Plain"); // default color is white instead of gray
99  gStyle->SetOptStat(0); // remove stat. box from all histos
100  legend = new TLegend(0.6,0.55,0.85,0.7); 
101  legend->SetFillStyle(0);
102  // q-distribution:
103  TH1D *qDistribution = dynamic_cast<TH1D*>(outputListFQD->FindObject("fqDistribution"));
104  Cosmetics(qDistribution);
105  legend->AddEntry(qDistribution,"q-distribution","f");
106  qDistribution->Draw();
107  // resulting fitting functions:
108  TF1 *fittingFun[2] = {NULL};
109  TString sigmaFlag[2] = {"#sigma^{2} not fitted","#sigma^{2} fitted"};
110  Double_t lineColors[2] = {kBlue,kRed};
111  for(Int_t s2F=0;s2F<2;s2F++)
112  {
113   fittingFun[s2F] = dynamic_cast<TF1*>(outputListFQD->FindObject(Form("fFittingFunction, %s",sigmaFlag[s2F].Data())));
114   fittingFun[s2F]->SetLineColor(lineColors[s2F]);
115   fittingFun[s2F]->Draw("SAME");
116   legend->AddEntry(fittingFun[s2F]->GetName(),sigmaFlag[s2F].Data(),"l");
117  }
118  legend->Draw("SAME"); 
119   
120 } // end of Plot()
121
122 // ===========================================================================================
123
124 void Save()
125 {
126  // Save the new results of fitting for FQD method in the common output file.
127  
128  dirFileFQD->Add(outputListFQD,kTRUE);
129  dirFileFQD->Write(dirFileFQD->GetName(),TObject::kSingleKey+TObject::kOverwrite);
130  //delete commonOutputFile;
131
132 } // end of void Save()
133
134 // ===========================================================================================
135
136 void Cosmetics(TH1D *hist)
137 {
138  // Set cosmetics for the q-distribution.
139  
140  Int_t firstNonEmptyBin = hist->FindFirstBinAbove(0);
141  Double_t lowRange = hist->GetBinLowEdge(firstNonEmptyBin);
142  Int_t lastNonEmptyBin = hist->FindLastBinAbove(0);
143  Double_t upperRange = hist->GetBinLowEdge(lastNonEmptyBin+10);
144  hist->GetXaxis()->SetRangeUser(lowRange,upperRange);  
145  hist->SetFillColor(16);  
146  hist->SetTitle("Fitted q-distribution");
147  
148 } // end of void Cosmetics(TH1D *hist)
149
150 // ===========================================================================================
151
152 TFile* AccessOutputFile(TString outputFileName)
153 {
154  // Access the common output file.
155     
156  TFile *outputFile = NULL;
157  if(!(gSystem->AccessPathName(Form("%s%s%s",gSystem->pwd(),"/",outputFileName.Data()),kFileExists)))
158  {
159   outputFile = TFile::Open(outputFileName.Data(),"UPDATE");
160  } else
161    {
162     cout<<endl;
163     cout<<"WARNING: Couldn't find the file "<<outputFileName.Data()<<" in "<<endl;
164     cout<<"         directory "<<gSystem->pwd()<<" !!!!"<<endl;
165     cout<<endl;
166     exit(0);
167    }
168   
169  return outputFile;
170  
171 } // end of TFile* AccessOutputFile(TString outputFileName)
172   
173 // ===========================================================================================
174  
175 void GetListWithHistograms(TFile *outputFile, TString analysisType, TString methodName)
176 {
177  // Access from common output file the TDirectoryFile for FQD method and from it
178  // the output list holding all objects:
179
180  TString fileName = ""; 
181  TString listName = ""; 
182  // Form a file name: 
183  fileName+="output";
184  fileName+=methodName.Data();
185  fileName+="analysis";
186  fileName+=analysisType.Data();
187  // Access this file:
188  dirFileFQD = (TDirectoryFile*)outputFile->FindObjectAny(fileName.Data());
189  // Form a list name for each method:
190  if(dirFileFQD)
191  {
192   TList* listTemp = dirFileFQD->GetListOfKeys();
193   if(listTemp && listTemp->GetEntries() == 1)
194   {
195    listName = listTemp->At(0)->GetName(); // to be improved - implemented better (use dynamic_cast)
196    dirFileFQD->GetObject(listName.Data(),outputListFQD);
197   } else
198     {
199      cout<<" WARNING: Accessing TList from TDirectoryFile failed miserably for FQD method !!!!"<<endl;
200      cout<<"          Did you actually use FQD in the analysis?"<<endl;
201      cout<<endl;
202      exit(0);
203     }
204  } else 
205    {
206     cout<<" WARNING: Couldn't find a TDirectoryFile "<<fileName.Data()<<" in file "<<outputFileName.Data()<<" !!!!"<<endl;
207     cout<<"          Perhaps wrong analysis type? Can be \"ESD\",\"AOD\",\"\", etc."<<endl;
208     cout<<"          Or you didn't use FQD in the analysis which have produced "<<outputFileName.Data()<<"?"<<endl;
209    } 
210           
211 } // end of void GetListWithHistograms(TFile *outputFile, TString analysisType, TString methodName)
212  
213 //===========================================================================================
214
215 void CrossCheckSettings() 
216 {
217  // Check in this method if the settings make sense.
218  
219  if(sigma2Min<0.5)
220  {
221   cout<<endl; 
222   cout<<"WARNING: According to theorists sigma2Min >= 0.5 !!!!"<<endl;
223   cout<<endl; 
224   exit(0);
225  }
226  
227 } // end of void CrossCheckSettings()
228
229 //===========================================================================================
230
231 void LoadLibrariesFQD(const libModes mode) {
232   
233   //--------------------------------------
234   // Load the needed libraries most of them already loaded by aliroot
235   //--------------------------------------
236   gSystem->Load("libTree.so");
237   gSystem->Load("libGeom.so");
238   gSystem->Load("libVMC.so");
239   gSystem->Load("libXMLIO.so");
240   gSystem->Load("libPhysics.so");
241   
242   //----------------------------------------------------------
243   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<< 
244   //----------------------------------------------------------
245   if (mode==mLocal) {
246     //--------------------------------------------------------
247     // If you want to use already compiled libraries 
248     // in the aliroot distribution
249     //--------------------------------------------------------
250
251   //==================================================================================  
252   //load needed libraries:
253   gSystem->AddIncludePath("-I$ROOTSYS/include");
254   gSystem->Load("libTree.so");
255
256   // for AliRoot
257   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
258   gSystem->Load("libANALYSIS.so");
259   gSystem->Load("libPWG2flowCommon.so");
260   cerr<<"libPWG2flowCommon.so loaded ..."<<endl;
261   
262   }
263   
264   else if (mode==mLocalSource) {
265  
266     // In root inline compile
267   
268     // Constants  
269     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
270     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
271     
272     // Flow event
273     gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+"); 
274     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");    
275     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");    
276     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
277     
278     // Output histosgrams
279     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
280     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
281     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
282     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
283        
284     cout << "finished loading macros!" << endl;  
285     
286   } // end of else if (mode==mLocalSource) 
287   
288 } // end of void LoadLibrariesFQD(const libModes mode)