1 #include "TStopwatch.h"
6 //----------------------------------------------------------
8 //flow analysis method can be: (set to kTRUE or kFALSE)
12 Bool_t LYZEP = kFALSE;
16 Bool_t MCEP = kFALSE; //does not work yet 24/12/08
17 //----------------------------------------------------------
19 //----------------------------------------------------------
21 //integrated selection
22 Double_t ptMaxInt = 10.;
23 Double_t ptMinInt = 0.;
24 Double_t etaMaxInt = 1.;
25 Double_t etaMinInt = -1.;
26 Double_t phiMaxInt = 7.5;
27 Double_t phiMinInt = 0.;
30 //differential selection
31 Double_t ptMaxDiff = 10.;
32 Double_t ptMinDiff = 0.;
33 Double_t etaMaxDiff = 1.;
34 Double_t etaMinDiff = -1.;
35 Double_t phiMaxDiff = 7.5;
36 Double_t phiMinDiff = 0.;
38 //----------------------------------------------------------
40 //----------------------------------------------------------
42 //to use or not to use the weights - that is a question!
43 Bool_t useWeightsPhi = kFALSE;//Phi
44 Bool_t useWeightsPt = kFALSE;//v'(pt)
45 Bool_t useWeightsEta = kFALSE;//v'(eta)
46 //----------------------------------------------------------
50 int runFlowAnalysis(Int_t aRuns = 44, const char*
51 dir="/data/alice1/kolk/KineOnly3/")
52 //dir="/Users/snelling/alice_data/KineOnly3/")
57 if (LYZ1 && LYZ2) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(); }
58 if (LYZ2 && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
59 if (LYZ1 && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
63 cout<<" ---- BEGIN ANALYSIS ---- "<<endl;
67 gSystem->AddIncludePath("-I$ROOTSYS/include");
70 gSystem->AddIncludePath("-I$ALICE_ROOT/include");
72 gSystem->Load("libANALYSIS.so");
73 gSystem->Load("libPWG2flow.so");
74 cerr<<"libPWG2flow.so loaded ..."<<endl;
77 // flow event in AliRoot
78 AliFlowEventSimpleMaker* fEventMaker = new AliFlowEventSimpleMaker();
81 //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
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)
89 fEventMaker->Init(file);
91 //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
99 gROOT->LoadMacro("code/AliFlowCommonConstants.cxx+");
100 gROOT->LoadMacro("code/AliFlowLYZConstants.cxx+");
101 gROOT->LoadMacro("code/AliFlowCumuConstants.cxx+");
104 gROOT->LoadMacro("code/AliFlowVector.cxx+");
105 gROOT->LoadMacro("code/AliFlowTrackSimple.cxx+");
106 gROOT->LoadMacro("code/AliFlowEventSimple.cxx+");
109 gROOT->LoadMacro("code/AliFlowTrackSimpleCuts.cxx+");
111 // output histosgrams
112 gROOT->LoadMacro("code/AliFlowCommonHist.cxx+");
113 gROOT->LoadMacro("code/AliFlowCommonHistResults.cxx+");
114 gROOT->LoadMacro("code/AliFlowLYZHist1.cxx+");fPhiWeights
115 gROOT->LoadMacro("code/AliFlowLYZHist2.cxx+");
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+");
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+");
132 // Class to fill the FlowEvent without aliroot dependence
133 // can be found in the directory other
134 gROOT->LoadMacro("code/FlowEventSimpleMaker.cxx+");
136 cout << "finished loading macros!" << endl;
139 FlowEventSimpleMaker* fEventMaker = new FlowEventSimpleMaker();
140 //------------------------------------------------------------------------
143 //load needed libraries
144 gSystem->Load("libTree.so");
146 //------------------------------------------------------------------------
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);
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);
167 //open the file with the weights:
169 TH1F *phiWeights = NULL;
170 TH1D *ptWeights = NULL;
171 TH1D *etaWeights = NULL;
173 if(useWeightsPhi||useWeightsPt||useWeightsEta)
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);
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;
189 AliFlowAnalysisWithMCEventPlane *mcep = NULL;
191 //MCEP = monte carlo event plane
193 AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
199 AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
201 if(useWeightsPhi) qc->SetPhiWeights(*phiWeights);
202 if(useWeightsPt) qc->SetPtWeights(*ptWeights);
203 if(useWeightsEta) qc->SetEtaWeights(*etaWeights);
206 //GFC = Generating Function Cumulants
208 AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
209 gfc->CreateOutputObjects();
212 //FQD = Fitting q-distribution
214 AliFittingQDistribution* fqd = new AliFittingQDistribution();
215 fqd->CreateOutputObjects();
218 //SP = Scalar Product
220 AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
224 //LYZ1 = Lee-Yang Zeroes first run
226 AliFlowAnalysisWithLeeYangZeros* lyz1 = new AliFlowAnalysisWithLeeYangZeros();
227 lyz1->SetFirstRun(kTRUE);
228 lyz1->SetUseSum(kTRUE);
232 //LYZ2 = Lee-Yang Zeroes second run
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 ;
242 TList* inputListLYZ2 = (TList*)inputFileLYZ2->Get("cobjLYZ1");
243 if (!inputListLYZ2) {cout<<"list is NULL pointer!"<<endl;}
245 cout<<"LYZ2 input file/list read..."<<endl;
246 lyz2->SetFirstRunList(inputListLYZ2);
247 lyz2->SetFirstRun(kFALSE);
248 lyz2->SetUseSum(kTRUE);
254 //LYZEP = Lee-Yang Zeroes event plane
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 ; }
264 TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2");
265 if (!inputListLYZEP) {cout<<"list is NULL pointer!"<<endl;}
267 cout<<"LYZEP input file/list read..."<<endl;
268 ep ->SetSecondRunList(inputListLYZEP);
269 lyzep->SetSecondRunList(inputListLYZEP);
275 //------------------------------------------------------------------------
280 TString execDir(gSystem->pwd());
281 TString targetDir(dir);
282 TSystemDirectory* baseDir = new TSystemDirectory(".", dir);
283 TList* dirList = baseDir->GetListOfFiles();
285 cout << endl << "No input files in: " << targetDir.Data() << endl;
288 Int_t nDirs = dirList->GetEntries();
290 cout<<"Int_t nDirs = "<<nDirs<<endl;
291 gSystem->cd(execDir);
293 for(Int_t iDir=0;iDir<nDirs;++iDir)
295 TSystemFile* presentDir = (TSystemFile*)dirList->At(iDir);
296 if(!presentDir || !presentDir->IsDirectory() ||
297 strcmp(presentDir->GetName(), ".") == 0 ||
298 strcmp(presentDir->GetName(), "..") == 0)
301 cout << "Directory (" << iDir << "): " << presentDir->GetName() <<
302 " - Skipping ... " << endl;
306 if(offset > 0) { --offset ; continue ; }
307 if((aRuns > 0) && (fCount >= aRuns)) { break ; }
309 TString presentDirName(dir); // aDataDir
310 presentDirName += presentDir->GetName();
311 presentDirName += "/";
312 //cerr<<" presentDirName = "<<presentDirName<<endl;
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))
319 cout << " File : " << fileName << " does NOT exist ! - Skipping ... "
323 // cout << endl ; cout << "Directory (" << iDir << "): " << presentDirName << " ... " << endl;
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();
330 //cout<<" Int_t nFiles = "<<nFiles<<endl;
331 gSystem->cd(execDir);
332 for(Int_t iFiles=0; iFiles<nFiles; ++iFiles)
334 TSystemFile* presentFile = (TSystemFile*) fileList->At(iFiles);
335 TString presentFileName(presentDirName);
336 presentFileName += presentFile->GetName();
338 if(!(presentFileName.Contains("Kinematics") &&
339 presentFileName.Contains("root"))) { continue ; }
341 //cout << " found: " << presentFileName.Data() << endl;
343 TFile* kineFile = new TFile(presentFileName.Data(), "READ");
345 Int_t nEvts = kineFile->GetNkeys() ;
346 //cout << " . found: " << nEvts << " KineTree(s) in " << presentFileName.Data() << endl;
347 TList* kineEventsList = (TList*)kineFile->GetListOfKeys();
349 TIter next(kineEventsList);
352 //loop over the events
353 while( key=(TKey *)next() )
355 TDirectory* tDir = (TDirectory*)key->ReadObj();
358 TString evtDir(tDir->GetName());
359 //cout << " . . found: " << tDir->GetName() << endl;
361 kTree = (TTree *)tDir->Get("TreeK");
364 Int_t nPart = kTree->GetEntries();
365 //cout << " . . . kTree " << fCount << " has " << nPart << " particles " << endl;
367 //-----------------------------------------------------------
368 //fill and save the flow event
369 AliFlowEventSimple *fEvent = fEventMaker->FillTracks(kTree, cutsInt, cutsDiff);
371 //pass the flow event to flow methods for analysis
372 Double_t fEP = 0.; // temporary number need true value
373 if(MCEP) mcep->Make(fEvent,fEP); //fix fEP
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);
381 //-----------------------------------------------------------
384 cout << "# " << fCount << " events processed" << endl;
393 //--------------------------------------------------------------
394 //calculating and storing the final results of flow analysis
395 if(MCEP) {mcep->Finish(); mcep->WriteHistograms("outputMCEPanalysis.root");}
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");}
403 //--------------------------------------------------------------
406 cout << " Finished ... " << endl;