more for weigths
[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
26c4cbb9 6//----------------------------------------------------------
334e3256 7//RUN SETTINGS
8//flow analysis method can be: (set to kTRUE or kFALSE)
26c4cbb9 9Bool_t SP = kFALSE;
334e3256 10Bool_t LYZ1 = kTRUE;
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
26c4cbb9 17//----------------------------------------------------------
334e3256 18
26c4cbb9 19//----------------------------------------------------------
334e3256 20//CUT SETTINGS
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;
26c4cbb9 38//----------------------------------------------------------
39
40//----------------------------------------------------------
41//WEIGHTS SETTINGS
42//to use or not to use the weights - that is a question!
43Bool_t useWeightsPhi = kFALSE;//Phi
44Bool_t useWeightsPt = kFALSE;//v'(pt)
45Bool_t useWeightsEta = kFALSE;//v'(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();
e085f1a9 79
80 /*
81 //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
82 // !!!! to be removed
26c4cbb9 83 // if weights are being used initialize the weight's histos in AliFlowEvenSimpleMaker:
84 fEventMaker->SetUseWeightsPhi(useWeightsPhi);
85 fEventMaker->SetUseWeightsPt(useWeightsPt);
86 fEventMaker->SetUseWeightsEta(useWeightsEta);
87 if(useWeightsPhi||useWeightsPt||useWeightsEta)
88 {
89 fEventMaker->Init(file);
90 }
e085f1a9 91 //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
92 */
93
26c4cbb9 94
2a413443 95 // In root
96
97 /*
334e3256 98 // constants
99 gROOT->LoadMacro("code/AliFlowCommonConstants.cxx+");
100 gROOT->LoadMacro("code/AliFlowLYZConstants.cxx+");
101 gROOT->LoadMacro("code/AliFlowCumuConstants.cxx+");
102
103 // flow event
104 gROOT->LoadMacro("code/AliFlowVector.cxx+");
105 gROOT->LoadMacro("code/AliFlowTrackSimple.cxx+");
106 gROOT->LoadMacro("code/AliFlowEventSimple.cxx+");
107
108 // cuts
109 gROOT->LoadMacro("code/AliFlowTrackSimpleCuts.cxx+");
110
111 // output histosgrams
112 gROOT->LoadMacro("code/AliFlowCommonHist.cxx+");
113 gROOT->LoadMacro("code/AliFlowCommonHistResults.cxx+");
e085f1a9 114 gROOT->LoadMacro("code/AliFlowLYZHist1.cxx+");fPhiWeights
334e3256 115 gROOT->LoadMacro("code/AliFlowLYZHist2.cxx+");
116
117 // functions needed for various methods
118 gROOT->LoadMacro("code/AliCumulantsFunctions.cxx+");
119 gROOT->LoadMacro("code/AliQCumulantsFunctions.cxx+");
120 gROOT->LoadMacro("code/AliFittingFunctionsForQDistribution.cxx+");
121 gROOT->LoadMacro("code/AliFlowLYZEventPlane.cxx+");
122
123 // Flow Analysis code for various methods
124 gROOT->LoadMacro("code/AliFlowAnalysisWithMCEventPlane.cxx+");
125 gROOT->LoadMacro("code/AliFlowAnalysisWithScalarProduct.cxx+");
126 gROOT->LoadMacro("code/AliFlowAnalysisWithLYZEventPlane.cxx+");
127 gROOT->LoadMacro("code/AliFlowAnalysisWithLeeYangZeros.cxx+");
128 gROOT->LoadMacro("code/AliFlowAnalysisWithCumulants.cxx+");
129 gROOT->LoadMacro("code/AliFlowAnalysisWithQCumulants.cxx+");
130 gROOT->LoadMacro("code/AliFittingQDistribution.cxx+");
131
2a413443 132 // Class to fill the FlowEvent without aliroot dependence
133 // can be found in the directory other
334e3256 134 gROOT->LoadMacro("code/FlowEventSimpleMaker.cxx+");
135
2a413443 136 cout << "finished loading macros!" << endl;
137
138 //flow event
139 FlowEventSimpleMaker* fEventMaker = new FlowEventSimpleMaker();
140 //------------------------------------------------------------------------
141 */
142
334e3256 143 //load needed libraries
144 gSystem->Load("libTree.so");
145
334e3256 146 //------------------------------------------------------------------------
e085f1a9 147 //cuts:
334e3256 148 AliFlowTrackSimpleCuts* cutsInt = new AliFlowTrackSimpleCuts();
149 cutsInt->SetPtMax(ptMaxInt);
150 cutsInt->SetPtMin(ptMinInt);
151 cutsInt->SetEtaMax(etaMaxInt);
152 cutsInt->SetEtaMin(etaMinInt);
153 cutsInt->SetPhiMax(phiMaxInt);
154 cutsInt->SetPhiMin(phiMinInt);
155 cutsInt->SetPID(PIDInt);
156
157 AliFlowTrackSimpleCuts* cutsDiff = new AliFlowTrackSimpleCuts();
158 cutsDiff->SetPtMax(ptMaxDiff);
159 cutsDiff->SetPtMin(ptMinDiff);
160 cutsDiff->SetEtaMax(etaMaxDiff);
161 cutsDiff->SetEtaMin(etaMinDiff);
162 cutsDiff->SetPhiMax(phiMaxDiff);
163 cutsDiff->SetPhiMin(phiMinDiff);
164 cutsDiff->SetPID(PIDDiff);
165
e085f1a9 166 //weights:
167 //open the file with the weights:
168 TFile *file = NULL;
169 TH1F *phiWeights = NULL;
170 TH1D *ptWeights = NULL;
171 TH1D *etaWeights = NULL;
334e3256 172
e085f1a9 173 if(useWeightsPhi||useWeightsPt||useWeightsEta)
174 {
175 file = TFile::Open("weightsForTheSecondRun.root", "READ");
176 if(useWeightsPhi) file->GetObject("phi_weights",phiWeights);
177 if(useWeightsPt) file->GetObject("pt_weights",ptWeights);
178 if(useWeightsEta) file->GetObject("eta_weights",etaWeights);
179 }
180
181 //flow methods:
334e3256 182 AliFlowAnalysisWithQCumulants *qc = NULL;
183 AliFlowAnalysisWithCumulants *gfc = NULL;
184 AliFittingQDistribution *fqd = NULL;
185 AliFlowAnalysisWithLeeYangZeros *lyz1 = NULL;
186 AliFlowAnalysisWithLeeYangZeros *lyz2 = NULL;
187 AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
188 AliFlowAnalysisWithScalarProduct *sp = NULL;
e085f1a9 189 AliFlowAnalysisWithMCEventPlane *mcep = NULL;
deebd72f 190
334e3256 191 //MCEP = monte carlo event plane
deebd72f 192 if (MCEP) {
193 AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
194 mcep->Init();
195 }
334e3256 196
197 //QC = Q-cumulants
deebd72f 198 if(QC) {
199 AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
e085f1a9 200 qc->Init();
201 if(useWeightsPhi) qc->SetPhiWeights(*phiWeights);
202 if(useWeightsPt) qc->SetPtWeights(*ptWeights);
203 if(useWeightsEta) qc->SetEtaWeights(*etaWeights);
deebd72f 204 }
334e3256 205
206 //GFC = Generating Function Cumulants
deebd72f 207 if(GFC) {
208 AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
209 gfc->CreateOutputObjects();
210 }
334e3256 211
212 //FQD = Fitting q-distribution
deebd72f 213 if(FQD) {
214 AliFittingQDistribution* fqd = new AliFittingQDistribution();
215 fqd->CreateOutputObjects();
216 }
60ec66fb 217
218 //SP = Scalar Product
219 if(SP) {
220 AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
221 sp->Init();
222 }
223
334e3256 224 //LYZ1 = Lee-Yang Zeroes first run
deebd72f 225 if(LYZ1) {
226 AliFlowAnalysisWithLeeYangZeros* lyz1 = new AliFlowAnalysisWithLeeYangZeros();
227 lyz1->SetFirstRun(kTRUE);
228 lyz1->SetUseSum(kTRUE);
229 lyz1->Init();
230 }
334e3256 231
232 //LYZ2 = Lee-Yang Zeroes second run
deebd72f 233 if(LYZ2) {
234 AliFlowAnalysisWithLeeYangZeros* lyz2 = new AliFlowAnalysisWithLeeYangZeros();
235 // read the input file from the first run
236 TString inputFileNameLYZ2 = "outputLYZ1analysis.root" ;
237 TFile* inputFileLYZ2 = new TFile(inputFileNameLYZ2.Data(),"READ");
238 if(!inputFileLYZ2 || inputFileLYZ2->IsZombie()) {
239 cerr << " ERROR: NO First Run file... " << endl ;
240 }
241 else {
242 TList* inputListLYZ2 = (TList*)inputFileLYZ2->Get("cobjLYZ1");
243 if (!inputListLYZ2) {cout<<"list is NULL pointer!"<<endl;}
244 else {
245 cout<<"LYZ2 input file/list read..."<<endl;
246 lyz2->SetFirstRunList(inputListLYZ2);
247 lyz2->SetFirstRun(kFALSE);
248 lyz2->SetUseSum(kTRUE);
249 lyz2->Init();
334e3256 250 }
251 }
deebd72f 252 }
253
334e3256 254 //LYZEP = Lee-Yang Zeroes event plane
deebd72f 255 if(LYZEP) {
256 AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
257 AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
258 // read the input file from the second lyz run
259 TString inputFileNameLYZEP = "outputLYZ2analysis.root" ;
260 TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
261 if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) {
262 cerr << " ERROR: NO Second Run file... " << endl ; }
263 else {
264 TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2");
265 if (!inputListLYZEP) {cout<<"list is NULL pointer!"<<endl;}
266 else {
267 cout<<"LYZEP input file/list read..."<<endl;
268 ep ->SetSecondRunList(inputListLYZEP);
269 lyzep->SetSecondRunList(inputListLYZEP);
270 ep ->Init();
271 lyzep->Init();
334e3256 272 }
273 }
deebd72f 274 }
334e3256 275 //------------------------------------------------------------------------
276
26c4cbb9 277
334e3256 278 //standard code
279 Int_t fCount = 0;
280 TString execDir(gSystem->pwd());
281 TString targetDir(dir);
282 TSystemDirectory* baseDir = new TSystemDirectory(".", dir);
283 TList* dirList = baseDir->GetListOfFiles();
284 if (!dirList) {
285 cout << endl << "No input files in: " << targetDir.Data() << endl;
286 break;
287 }
288 Int_t nDirs = dirList->GetEntries();
289 cout<<endl;
290 cout<<"Int_t nDirs = "<<nDirs<<endl;
291 gSystem->cd(execDir);
292
293 for(Int_t iDir=0;iDir<nDirs;++iDir)
294 {
295 TSystemFile* presentDir = (TSystemFile*)dirList->At(iDir);
296 if(!presentDir || !presentDir->IsDirectory() ||
297 strcmp(presentDir->GetName(), ".") == 0 ||
298 strcmp(presentDir->GetName(), "..") == 0)
299 {
300 cout << endl;
301 cout << "Directory (" << iDir << "): " << presentDir->GetName() <<
302 " - Skipping ... " << endl;
303 continue ;
304 }
305
306 if(offset > 0) { --offset ; continue ; }
307 if((aRuns > 0) && (fCount >= aRuns)) { break ; }
308
309 TString presentDirName(dir); // aDataDir
310 presentDirName += presentDir->GetName();
311 presentDirName += "/";
312 //cerr<<" presentDirName = "<<presentDirName<<endl;
313
314 TString fileName = presentDirName;
315 fileName += "galice.root";
316 Long_t *id, *size, *flags, *modtime;
317 if(gSystem->GetPathInfo(fileName.Data(),id,size,flags,modtime))
318 {
319 cout << " File : " << fileName << " does NOT exist ! - Skipping ... "
320 << endl;
321 continue;
322 }
deebd72f 323 // cout << endl ; cout << "Directory (" << iDir << "): " << presentDirName << " ... " << endl;
334e3256 324
325 //loop (simulations in the present dir)
326 TSystemDirectory* evtsDir = new TSystemDirectory(".",
327 presentDirName.Data());
328 TList* fileList = evtsDir->GetListOfFiles();
329 Int_t nFiles = fileList->GetEntries();
deebd72f 330 //cout<<" Int_t nFiles = "<<nFiles<<endl;
334e3256 331 gSystem->cd(execDir);
332 for(Int_t iFiles=0; iFiles<nFiles; ++iFiles)
333 {
334 TSystemFile* presentFile = (TSystemFile*) fileList->At(iFiles);
335 TString presentFileName(presentDirName);
336 presentFileName += presentFile->GetName();
337
338 if(!(presentFileName.Contains("Kinematics") &&
339 presentFileName.Contains("root"))) { continue ; }
340
deebd72f 341 //cout << " found: " << presentFileName.Data() << endl;
334e3256 342
343 TFile* kineFile = new TFile(presentFileName.Data(), "READ");
344 // kineFile->ls();
345 Int_t nEvts = kineFile->GetNkeys() ;
deebd72f 346 //cout << " . found: " << nEvts << " KineTree(s) in " << presentFileName.Data() << endl;
334e3256 347 TList* kineEventsList = (TList*)kineFile->GetListOfKeys();
348 TTree* kTree;
349 TIter next(kineEventsList);
350 TKey* key;
351
352 //loop over the events
353 while( key=(TKey *)next() )
354 {
355 TDirectory* tDir = (TDirectory*)key->ReadObj();
356 if(!tDir) break;
357
358 TString evtDir(tDir->GetName());
deebd72f 359 //cout << " . . found: " << tDir->GetName() << endl;
334e3256 360
361 kTree = (TTree *)tDir->Get("TreeK");
362 if(!kTree) break;
363
364 Int_t nPart = kTree->GetEntries();
deebd72f 365 //cout << " . . . kTree " << fCount << " has " << nPart << " particles " << endl;
334e3256 366
367 //-----------------------------------------------------------
26c4cbb9 368 //fill and save the flow event
369 AliFlowEventSimple *fEvent = fEventMaker->FillTracks(kTree, cutsInt, cutsDiff);
370
334e3256 371 //pass the flow event to flow methods for analysis
60ec66fb 372 Double_t fEP = 0.; // temporary number need true value
373 if(MCEP) mcep->Make(fEvent,fEP); //fix fEP
deebd72f 374 if(QC) qc->Make(fEvent);
375 if(GFC) gfc->Make(fEvent);
376 if(FQD) fqd->Make(fEvent);
377 if(LYZ1) lyz1->Make(fEvent);
378 if(LYZ2) lyz2->Make(fEvent);
379 if(LYZEP) lyzep->Make(fEvent,ep);
380 if(SP) sp->Make(fEvent);
334e3256 381 //-----------------------------------------------------------
382
383 fCount++;
deebd72f 384 cout << "# " << fCount << " events processed" << endl;
334e3256 385 delete kTree;
386 delete fEvent;
387 }
388 delete kineFile ;
389 }
390 delete evtsDir ;
391 }
392
393 //--------------------------------------------------------------
394 //calculating and storing the final results of flow analysis
60ec66fb 395 if(MCEP) {mcep->Finish(); mcep->WriteHistograms("outputMCEPanalysis.root");}
deebd72f 396 if(QC) {qc->Finish(); qc->WriteHistograms("outputQCanalysis.root");}
397 if(GFC) {gfc->Finish(); gfc->WriteHistograms("outputGFCanalysis.root");}
398 if(FQD) {fqd->Finish(); fqd->WriteHistograms("outputFQDanalysis.root");}
399 if(LYZ1) {lyz1->Finish(); lyz1->WriteHistograms("outputLYZ1analysis.root");}
400 if(LYZ2) {lyz2->Finish(); lyz2->WriteHistograms("outputLYZ2analysis.root");}
401 if(LYZEP) {lyzep->Finish(); lyzep->WriteHistograms("outputLYZEPanalysis.root");}
402 if(SP) {sp->Finish(); sp->WriteHistograms("outputSPanalysis.root");}
334e3256 403 //--------------------------------------------------------------
404
405 cout << endl;
406 cout << " Finished ... " << endl;
407 cout << endl;
408
409 timer.Stop();
410 cout << endl;
411 timer.Print();
412}