Implement weigths for Q cumulants in local run
[u/mrichter/AliRoot.git] / PWG2 / FLOW / macros / runFlowAnalysis.C
CommitLineData
334e3256 1#include "TStopwatch.h"
2#include "TObjArray"
3#include "Riostream.h"
26c4cbb9 4#include "TFile.h"
334e3256 5
03a02aca 6//--------------------------------------------------------------------------------------
7// RUN SETTINGS
334e3256 8//flow analysis method can be: (set to kTRUE or kFALSE)
26c4cbb9 9Bool_t SP = kFALSE;
03a02aca 10Bool_t LYZ1 = kFALSE;
334e3256 11Bool_t LYZ2 = kFALSE;
12Bool_t LYZEP = kFALSE;
13Bool_t GFC = kTRUE;
14Bool_t QC = kTRUE;
26c4cbb9 15Bool_t FQD = kFALSE;
334e3256 16Bool_t MCEP = kFALSE; //does not work yet 24/12/08
03a02aca 17//--------------------------------------------------------------------------------------
334e3256 18
03a02aca 19//--------------------------------------------------------------------------------------
20// CUT SETTINGS
334e3256 21//integrated selection
22Double_t ptMaxInt = 10.;
23Double_t ptMinInt = 0.;
24Double_t etaMaxInt = 1.;
25Double_t etaMinInt = -1.;
26Double_t phiMaxInt = 7.5;
27Double_t phiMinInt = 0.;
28Int_t PIDInt = 211;
29
30//differential selection
31Double_t ptMaxDiff = 10.;
32Double_t ptMinDiff = 0.;
33Double_t etaMaxDiff = 1.;
34Double_t etaMinDiff = -1.;
35Double_t phiMaxDiff = 7.5;
36Double_t phiMinDiff = 0.;
37Int_t PIDDiff = 211;
03a02aca 38//--------------------------------------------------------------------------------------
26c4cbb9 39
03a02aca 40//--------------------------------------------------------------------------------------
41// WEIGHTS SETTINGS
26c4cbb9 42//to use or not to use the weights - that is a question!
03a02aca 43Bool_t usePhiWeights = kFALSE; //Phi (correction for non-uniform azimuthal acceptance)
44Bool_t usePtWeights = kFALSE; //v'(pt) (differential flow in pt)
45Bool_t useEtaWeights = kFALSE; //v'(eta) (differential flow in eta)
46//--------------------------------------------------------------------------------------
334e3256 47
26c4cbb9 48Int_t offset = 0;
49
50int runFlowAnalysis(Int_t aRuns = 44, const char*
51 dir="/data/alice1/kolk/KineOnly3/")
52 //dir="/Users/snelling/alice_data/KineOnly3/")
334e3256 53{
54 TStopwatch timer;
55 timer.Start();
56
57 if (LYZ1 && LYZ2) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(); }
334e3256 58 if (LYZ2 && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
334e3256 59 if (LYZ1 && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
60
61
62 cout<<endl;
63 cout<<" ---- BEGIN ANALYSIS ---- "<<endl;
64 cout<<endl;
65
2a413443 66
334e3256 67 gSystem->AddIncludePath("-I$ROOTSYS/include");
2a413443 68
69 // In AliRoot
70 gSystem->AddIncludePath("-I$ALICE_ROOT/include");
71
72 gSystem->Load("libANALYSIS.so");
73 gSystem->Load("libPWG2flow.so");
74 cerr<<"libPWG2flow.so loaded ..."<<endl;
75 cout<<endl;
26c4cbb9 76
2a413443 77 // flow event in AliRoot
26c4cbb9 78 AliFlowEventSimpleMaker* fEventMaker = new AliFlowEventSimpleMaker();
03a02aca 79
2a413443 80 // In root
81
82 /*
334e3256 83 // constants
84 gROOT->LoadMacro("code/AliFlowCommonConstants.cxx+");
85 gROOT->LoadMacro("code/AliFlowLYZConstants.cxx+");
86 gROOT->LoadMacro("code/AliFlowCumuConstants.cxx+");
87
88 // flow event
89 gROOT->LoadMacro("code/AliFlowVector.cxx+");
90 gROOT->LoadMacro("code/AliFlowTrackSimple.cxx+");
91 gROOT->LoadMacro("code/AliFlowEventSimple.cxx+");
92
93 // cuts
94 gROOT->LoadMacro("code/AliFlowTrackSimpleCuts.cxx+");
95
96 // output histosgrams
97 gROOT->LoadMacro("code/AliFlowCommonHist.cxx+");
98 gROOT->LoadMacro("code/AliFlowCommonHistResults.cxx+");
03a02aca 99 gROOT->LoadMacro("code/AliFlowLYZHist1.cxx+");
334e3256 100 gROOT->LoadMacro("code/AliFlowLYZHist2.cxx+");
101
102 // functions needed for various methods
103 gROOT->LoadMacro("code/AliCumulantsFunctions.cxx+");
104 gROOT->LoadMacro("code/AliQCumulantsFunctions.cxx+");
105 gROOT->LoadMacro("code/AliFittingFunctionsForQDistribution.cxx+");
106 gROOT->LoadMacro("code/AliFlowLYZEventPlane.cxx+");
107
108 // Flow Analysis code for various methods
109 gROOT->LoadMacro("code/AliFlowAnalysisWithMCEventPlane.cxx+");
110 gROOT->LoadMacro("code/AliFlowAnalysisWithScalarProduct.cxx+");
111 gROOT->LoadMacro("code/AliFlowAnalysisWithLYZEventPlane.cxx+");
112 gROOT->LoadMacro("code/AliFlowAnalysisWithLeeYangZeros.cxx+");
113 gROOT->LoadMacro("code/AliFlowAnalysisWithCumulants.cxx+");
114 gROOT->LoadMacro("code/AliFlowAnalysisWithQCumulants.cxx+");
115 gROOT->LoadMacro("code/AliFittingQDistribution.cxx+");
116
2a413443 117 // Class to fill the FlowEvent without aliroot dependence
118 // can be found in the directory other
334e3256 119 gROOT->LoadMacro("code/FlowEventSimpleMaker.cxx+");
120
2a413443 121 cout << "finished loading macros!" << endl;
122
123 //flow event
124 FlowEventSimpleMaker* fEventMaker = new FlowEventSimpleMaker();
125 //------------------------------------------------------------------------
126 */
127
334e3256 128 //load needed libraries
129 gSystem->Load("libTree.so");
130
334e3256 131 //------------------------------------------------------------------------
e085f1a9 132 //cuts:
334e3256 133 AliFlowTrackSimpleCuts* cutsInt = new AliFlowTrackSimpleCuts();
134 cutsInt->SetPtMax(ptMaxInt);
135 cutsInt->SetPtMin(ptMinInt);
136 cutsInt->SetEtaMax(etaMaxInt);
137 cutsInt->SetEtaMin(etaMinInt);
138 cutsInt->SetPhiMax(phiMaxInt);
139 cutsInt->SetPhiMin(phiMinInt);
140 cutsInt->SetPID(PIDInt);
141
142 AliFlowTrackSimpleCuts* cutsDiff = new AliFlowTrackSimpleCuts();
143 cutsDiff->SetPtMax(ptMaxDiff);
144 cutsDiff->SetPtMin(ptMinDiff);
145 cutsDiff->SetEtaMax(etaMaxDiff);
146 cutsDiff->SetEtaMin(etaMinDiff);
147 cutsDiff->SetPhiMax(phiMaxDiff);
148 cutsDiff->SetPhiMin(phiMinDiff);
149 cutsDiff->SetPID(PIDDiff);
150
03a02aca 151 //if the weights are used:
152 TFile *fileWithWeights = NULL;
153 TList *listWithWeights = NULL;
334e3256 154
03a02aca 155 if(usePhiWeights||usePtWeights||useEtaWeights)
e085f1a9 156 {
03a02aca 157 fileWithWeights = TFile::Open("weights.root","READ");
158 if(fileWithWeights)
159 {
160 listWithWeights = (TList*)fileWithWeights->Get("weights");
161 }else{cout<<" WARNING: the file <weights.root> with weights from the previous run was not accessed."<<endl;}
e085f1a9 162 }
163
164 //flow methods:
334e3256 165 AliFlowAnalysisWithQCumulants *qc = NULL;
166 AliFlowAnalysisWithCumulants *gfc = NULL;
167 AliFittingQDistribution *fqd = NULL;
168 AliFlowAnalysisWithLeeYangZeros *lyz1 = NULL;
169 AliFlowAnalysisWithLeeYangZeros *lyz2 = NULL;
170 AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
171 AliFlowAnalysisWithScalarProduct *sp = NULL;
e085f1a9 172 AliFlowAnalysisWithMCEventPlane *mcep = NULL;
deebd72f 173
334e3256 174 //MCEP = monte carlo event plane
deebd72f 175 if (MCEP) {
176 AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
177 mcep->Init();
178 }
334e3256 179
180 //QC = Q-cumulants
deebd72f 181 if(QC) {
182 AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
e085f1a9 183 qc->Init();
03a02aca 184 if(listWithWeights) qc->SetWeightsList(listWithWeights);
185 if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
186 if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
187 if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
deebd72f 188 }
334e3256 189
190 //GFC = Generating Function Cumulants
deebd72f 191 if(GFC) {
192 AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
193 gfc->CreateOutputObjects();
194 }
334e3256 195
196 //FQD = Fitting q-distribution
deebd72f 197 if(FQD) {
198 AliFittingQDistribution* fqd = new AliFittingQDistribution();
199 fqd->CreateOutputObjects();
200 }
60ec66fb 201
202 //SP = Scalar Product
203 if(SP) {
204 AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
205 sp->Init();
206 }
207
334e3256 208 //LYZ1 = Lee-Yang Zeroes first run
deebd72f 209 if(LYZ1) {
210 AliFlowAnalysisWithLeeYangZeros* lyz1 = new AliFlowAnalysisWithLeeYangZeros();
211 lyz1->SetFirstRun(kTRUE);
212 lyz1->SetUseSum(kTRUE);
213 lyz1->Init();
214 }
334e3256 215
216 //LYZ2 = Lee-Yang Zeroes second run
deebd72f 217 if(LYZ2) {
218 AliFlowAnalysisWithLeeYangZeros* lyz2 = new AliFlowAnalysisWithLeeYangZeros();
219 // read the input file from the first run
220 TString inputFileNameLYZ2 = "outputLYZ1analysis.root" ;
221 TFile* inputFileLYZ2 = new TFile(inputFileNameLYZ2.Data(),"READ");
222 if(!inputFileLYZ2 || inputFileLYZ2->IsZombie()) {
223 cerr << " ERROR: NO First Run file... " << endl ;
224 }
225 else {
226 TList* inputListLYZ2 = (TList*)inputFileLYZ2->Get("cobjLYZ1");
227 if (!inputListLYZ2) {cout<<"list is NULL pointer!"<<endl;}
228 else {
229 cout<<"LYZ2 input file/list read..."<<endl;
230 lyz2->SetFirstRunList(inputListLYZ2);
231 lyz2->SetFirstRun(kFALSE);
232 lyz2->SetUseSum(kTRUE);
233 lyz2->Init();
334e3256 234 }
235 }
deebd72f 236 }
237
334e3256 238 //LYZEP = Lee-Yang Zeroes event plane
deebd72f 239 if(LYZEP) {
240 AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
241 AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
242 // read the input file from the second lyz run
243 TString inputFileNameLYZEP = "outputLYZ2analysis.root" ;
244 TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
245 if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) {
246 cerr << " ERROR: NO Second Run file... " << endl ; }
247 else {
248 TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2");
249 if (!inputListLYZEP) {cout<<"list is NULL pointer!"<<endl;}
250 else {
251 cout<<"LYZEP input file/list read..."<<endl;
252 ep ->SetSecondRunList(inputListLYZEP);
253 lyzep->SetSecondRunList(inputListLYZEP);
254 ep ->Init();
255 lyzep->Init();
334e3256 256 }
257 }
deebd72f 258 }
334e3256 259 //------------------------------------------------------------------------
260
26c4cbb9 261
334e3256 262 //standard code
263 Int_t fCount = 0;
264 TString execDir(gSystem->pwd());
265 TString targetDir(dir);
266 TSystemDirectory* baseDir = new TSystemDirectory(".", dir);
267 TList* dirList = baseDir->GetListOfFiles();
268 if (!dirList) {
269 cout << endl << "No input files in: " << targetDir.Data() << endl;
270 break;
271 }
272 Int_t nDirs = dirList->GetEntries();
273 cout<<endl;
274 cout<<"Int_t nDirs = "<<nDirs<<endl;
275 gSystem->cd(execDir);
276
277 for(Int_t iDir=0;iDir<nDirs;++iDir)
278 {
279 TSystemFile* presentDir = (TSystemFile*)dirList->At(iDir);
280 if(!presentDir || !presentDir->IsDirectory() ||
281 strcmp(presentDir->GetName(), ".") == 0 ||
282 strcmp(presentDir->GetName(), "..") == 0)
283 {
284 cout << endl;
285 cout << "Directory (" << iDir << "): " << presentDir->GetName() <<
286 " - Skipping ... " << endl;
287 continue ;
288 }
289
290 if(offset > 0) { --offset ; continue ; }
291 if((aRuns > 0) && (fCount >= aRuns)) { break ; }
292
293 TString presentDirName(dir); // aDataDir
294 presentDirName += presentDir->GetName();
295 presentDirName += "/";
296 //cerr<<" presentDirName = "<<presentDirName<<endl;
297
298 TString fileName = presentDirName;
299 fileName += "galice.root";
300 Long_t *id, *size, *flags, *modtime;
301 if(gSystem->GetPathInfo(fileName.Data(),id,size,flags,modtime))
302 {
303 cout << " File : " << fileName << " does NOT exist ! - Skipping ... "
304 << endl;
305 continue;
306 }
deebd72f 307 // cout << endl ; cout << "Directory (" << iDir << "): " << presentDirName << " ... " << endl;
334e3256 308
309 //loop (simulations in the present dir)
310 TSystemDirectory* evtsDir = new TSystemDirectory(".",
311 presentDirName.Data());
312 TList* fileList = evtsDir->GetListOfFiles();
313 Int_t nFiles = fileList->GetEntries();
deebd72f 314 //cout<<" Int_t nFiles = "<<nFiles<<endl;
334e3256 315 gSystem->cd(execDir);
316 for(Int_t iFiles=0; iFiles<nFiles; ++iFiles)
317 {
318 TSystemFile* presentFile = (TSystemFile*) fileList->At(iFiles);
319 TString presentFileName(presentDirName);
320 presentFileName += presentFile->GetName();
321
322 if(!(presentFileName.Contains("Kinematics") &&
323 presentFileName.Contains("root"))) { continue ; }
324
deebd72f 325 //cout << " found: " << presentFileName.Data() << endl;
334e3256 326
327 TFile* kineFile = new TFile(presentFileName.Data(), "READ");
328 // kineFile->ls();
329 Int_t nEvts = kineFile->GetNkeys() ;
deebd72f 330 //cout << " . found: " << nEvts << " KineTree(s) in " << presentFileName.Data() << endl;
334e3256 331 TList* kineEventsList = (TList*)kineFile->GetListOfKeys();
332 TTree* kTree;
333 TIter next(kineEventsList);
334 TKey* key;
335
336 //loop over the events
337 while( key=(TKey *)next() )
338 {
339 TDirectory* tDir = (TDirectory*)key->ReadObj();
340 if(!tDir) break;
341
342 TString evtDir(tDir->GetName());
deebd72f 343 //cout << " . . found: " << tDir->GetName() << endl;
334e3256 344
345 kTree = (TTree *)tDir->Get("TreeK");
346 if(!kTree) break;
347
348 Int_t nPart = kTree->GetEntries();
deebd72f 349 //cout << " . . . kTree " << fCount << " has " << nPart << " particles " << endl;
334e3256 350
351 //-----------------------------------------------------------
26c4cbb9 352 //fill and save the flow event
353 AliFlowEventSimple *fEvent = fEventMaker->FillTracks(kTree, cutsInt, cutsDiff);
354
334e3256 355 //pass the flow event to flow methods for analysis
60ec66fb 356 Double_t fEP = 0.; // temporary number need true value
357 if(MCEP) mcep->Make(fEvent,fEP); //fix fEP
deebd72f 358 if(QC) qc->Make(fEvent);
359 if(GFC) gfc->Make(fEvent);
360 if(FQD) fqd->Make(fEvent);
361 if(LYZ1) lyz1->Make(fEvent);
362 if(LYZ2) lyz2->Make(fEvent);
363 if(LYZEP) lyzep->Make(fEvent,ep);
364 if(SP) sp->Make(fEvent);
334e3256 365 //-----------------------------------------------------------
366
367 fCount++;
deebd72f 368 cout << "# " << fCount << " events processed" << endl;
334e3256 369 delete kTree;
370 delete fEvent;
371 }
372 delete kineFile ;
373 }
374 delete evtsDir ;
375 }
376
377 //--------------------------------------------------------------
378 //calculating and storing the final results of flow analysis
60ec66fb 379 if(MCEP) {mcep->Finish(); mcep->WriteHistograms("outputMCEPanalysis.root");}
deebd72f 380 if(QC) {qc->Finish(); qc->WriteHistograms("outputQCanalysis.root");}
381 if(GFC) {gfc->Finish(); gfc->WriteHistograms("outputGFCanalysis.root");}
382 if(FQD) {fqd->Finish(); fqd->WriteHistograms("outputFQDanalysis.root");}
383 if(LYZ1) {lyz1->Finish(); lyz1->WriteHistograms("outputLYZ1analysis.root");}
384 if(LYZ2) {lyz2->Finish(); lyz2->WriteHistograms("outputLYZ2analysis.root");}
385 if(LYZEP) {lyzep->Finish(); lyzep->WriteHistograms("outputLYZEPanalysis.root");}
386 if(SP) {sp->Finish(); sp->WriteHistograms("outputSPanalysis.root");}
334e3256 387 //--------------------------------------------------------------
388
389 cout << endl;
390 cout << " Finished ... " << endl;
391 cout << endl;
392
393 timer.Stop();
394 cout << endl;
395 timer.Print();
396}