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