hooks for PMD flow analysis
[u/mrichter/AliRoot.git] / PWG2 / FLOW / macros / fqd.C
CommitLineData
f0e0fbcd 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:
9Double_t vStart = 0.05; // starting value for v fit
10Double_t vMin = 0.0; // lower bound for v fit
11Double_t vMax = 0.25; // upper bound for v fit
12Double_t sigma2Start = 0.75; // starting value for sigma^2 fit
13Double_t sigma2Min = 0.5; // lower bound for sigma^2 fit (according to theorists must be >= 0.5)
14Double_t sigma2Max = 2.5; // upper bound for sigma^2 fit
15Double_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
16Bool_t finalResultIsFromSigma2Fitted = kTRUE; // final saved result is obtained with sigma^2 fitted (kTRUE) or sigma^2 fixed (kFALSE)
17Bool_t printOnTheScreen = kTRUE; // print or not the final results on the screen
18Bool_t plotResults = kTRUE; // plot on the screen q-distribution and resulting fitting functions (for non-fixed sigma^2 and fixed sigma^2)
19
5c09ff70 20// Name of the common output file:
21TString outputFileName = "AnalysisResults.root";
22//TString outputFileName = "outputCentrality0.root";
23
f0e0fbcd 24enum libModes {mLocal,mLocalSource};
25
26TFile *commonOutputFile = NULL; // common output file "AnalysisResults.root"
27TDirectoryFile *dirFileFQD = NULL; // FQD's TDirectoryFile in "AnalysisResults.root"
28TList *outputListFQD = NULL; // output list holding all FQD objects
29
5c09ff70 30void fqd(TString analysisType="ESD", Int_t analysisMode=mLocal)
f0e0fbcd 31{
5c09ff70 32 // 1. analysisType: type of analysis can be ESD, AOD, MC, MK, ....
33 // (if type="" output files are from simulation 'on the fly')
f0e0fbcd 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
54void 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
79void GetOutputList(TString analysisType)
80{
81 // Get output list which holds all FQD objects.
82
83 // Access common output file:
f0e0fbcd 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
94void 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
5c09ff70 100 legend = new TLegend(0.6,0.55,0.85,0.7);
101 legend->SetFillStyle(0);
f0e0fbcd 102 // q-distribution:
103 TH1D *qDistribution = dynamic_cast<TH1D*>(outputListFQD->FindObject("fqDistribution"));
104 Cosmetics(qDistribution);
5c09ff70 105 legend->AddEntry(qDistribution,"q-distribution","f");
f0e0fbcd 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");
5c09ff70 116 legend->AddEntry(fittingFun[s2F]->GetName(),sigmaFlag[s2F].Data(),"l");
f0e0fbcd 117 }
5c09ff70 118 legend->Draw("SAME");
f0e0fbcd 119
120} // end of Plot()
121
122// ===========================================================================================
123
124void 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
136void 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
152TFile* 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
175void 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());
5c09ff70 189 // Form a list name for each method:
f0e0fbcd 190 if(dirFileFQD)
191 {
5c09ff70 192 TList* listTemp = dirFileFQD->GetListOfKeys();
193 if(listTemp && listTemp->GetEntries() == 1)
f0e0fbcd 194 {
5c09ff70 195 listName = listTemp->At(0)->GetName(); // to be improved - implemented better (use dynamic_cast)
f0e0fbcd 196 dirFileFQD->GetObject(listName.Data(),outputListFQD);
197 } else
198 {
5c09ff70 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 }
f0e0fbcd 204 } else
205 {
5c09ff70 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
f0e0fbcd 211} // end of void GetListWithHistograms(TFile *outputFile, TString analysisType, TString methodName)
212
213//===========================================================================================
214
215void 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
231void 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+");
f0e0fbcd 271
272 // Flow event
273 gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+");
274 gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");
f0e0fbcd 275 gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");
5062cb0f 276 gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
f0e0fbcd 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)