small fix
[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
175   //SP = Scalar Product 
176   if(SP) {
177     AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
178     sp->Init();
179   }
180
181   //LYZ1 = Lee-Yang Zeroes first run
182   if(LYZ1) {
183     AliFlowAnalysisWithLeeYangZeros* lyz1 = new AliFlowAnalysisWithLeeYangZeros();
184     lyz1->SetFirstRun(kTRUE);
185     lyz1->SetUseSum(kTRUE);
186     lyz1->Init();
187   }
188
189   //LYZ2 = Lee-Yang Zeroes second run
190  
191   if(LYZ2) {
192     AliFlowAnalysisWithLeeYangZeros* lyz2 = new AliFlowAnalysisWithLeeYangZeros();
193     // read the input file from the first run 
194     TString inputFileNameLYZ2 = "outputLYZ1analysis.root" ;
195     TFile* inputFileLYZ2 = new TFile(inputFileNameLYZ2.Data(),"READ");
196     if(!inputFileLYZ2 || inputFileLYZ2->IsZombie()) { 
197       cerr << " ERROR: NO First Run file... " << endl ; 
198     }
199     else { 
200       TList* inputListLYZ2 = (TList*)inputFileLYZ2->Get("cobjLYZ1");  
201       if (!inputListLYZ2) {cout<<"list is NULL pointer!"<<endl;}
202       else {
203         cout<<"LYZ2 input file/list read..."<<endl;
204         lyz2->SetFirstRunList(inputListLYZ2);
205         lyz2->SetFirstRun(kFALSE);
206         lyz2->SetUseSum(kTRUE);
207         lyz2->Init();
208       }
209     }
210   }
211   
212  //LYZEP = Lee-Yang Zeroes event plane
213   if(LYZEP) {
214     AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
215     AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
216     // read the input file from the second lyz run 
217     TString inputFileNameLYZEP = "outputLYZ2analysis.root" ;
218     TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
219     if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) { 
220       cerr << " ERROR: NO Second Run file... " << endl ; }
221     else { 
222       TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2");  
223       if (!inputListLYZEP) {cout<<"list is NULL pointer!"<<endl;}
224       else {
225         cout<<"LYZEP input file/list read..."<<endl;
226         ep   ->SetSecondRunList(inputListLYZEP);
227         lyzep->SetSecondRunList(inputListLYZEP);
228         ep   ->Init();
229         lyzep->Init();
230       }
231     }
232   }
233   //------------------------------------------------------------------------
234   
235   //standard code
236   Int_t fCount = 0;
237   TString execDir(gSystem->pwd());
238   TString targetDir(dir);
239   TSystemDirectory* baseDir = new TSystemDirectory(".", dir);  
240   TList* dirList           = baseDir->GetListOfFiles();
241   if (!dirList) {
242     cout << endl << "No input files in: " << targetDir.Data() << endl;
243     break;
244   }
245   Int_t nDirs              = dirList->GetEntries();
246   cout<<endl;
247   cout<<"Int_t nDirs = "<<nDirs<<endl;
248   gSystem->cd(execDir);
249   
250   for(Int_t iDir=0;iDir<nDirs;++iDir)
251     {
252       TSystemFile* presentDir = (TSystemFile*)dirList->At(iDir);
253       if(!presentDir || !presentDir->IsDirectory() || 
254          strcmp(presentDir->GetName(), ".") == 0 || 
255          strcmp(presentDir->GetName(), "..") == 0) 
256         {
257           cout << endl; 
258           cout << "Directory (" << iDir << "): " << presentDir->GetName() << 
259             " - Skipping ... " << endl;
260           continue ;   
261         }
262       
263       if(offset > 0) { --offset ; continue ; }
264       if((aRuns > 0) && (fCount >= aRuns)) { break ; }
265       
266       TString presentDirName(dir); // aDataDir
267       presentDirName += presentDir->GetName();
268       presentDirName += "/";
269       //cerr<<" presentDirName = "<<presentDirName<<endl;
270       
271       TString fileName = presentDirName; 
272       fileName += "galice.root";
273       Long_t *id, *size, *flags, *modtime;
274       if(gSystem->GetPathInfo(fileName.Data(),id,size,flags,modtime)) 
275         { 
276           cout << " File : " << fileName << " does NOT exist ! - Skipping ... " 
277                << endl; 
278           continue; 
279         }
280       //      cout << endl ; cout << "Directory (" << iDir << "): " << presentDirName << " ... " << endl;
281       
282       //loop (simulations in the present dir) 
283       TSystemDirectory* evtsDir = new TSystemDirectory(".", 
284                                                        presentDirName.Data());
285       TList* fileList       = evtsDir->GetListOfFiles();
286       Int_t nFiles                  = fileList->GetEntries();
287       //cout<<" Int_t nFiles = "<<nFiles<<endl;
288       gSystem->cd(execDir);      
289       for(Int_t iFiles=0; iFiles<nFiles; ++iFiles)
290         {
291           TSystemFile* presentFile = (TSystemFile*) fileList->At(iFiles);
292           TString presentFileName(presentDirName);
293           presentFileName += presentFile->GetName();
294           
295           if(!(presentFileName.Contains("Kinematics") && 
296                presentFileName.Contains("root"))) { continue ; }
297           
298           //cout << " found: " << presentFileName.Data() << endl; 
299           
300           TFile* kineFile = new TFile(presentFileName.Data(), "READ"); 
301           // kineFile->ls();
302           Int_t nEvts = kineFile->GetNkeys() ; 
303           //cout << "  . found: " << nEvts << " KineTree(s) in " << presentFileName.Data() << endl;
304           TList* kineEventsList = (TList*)kineFile->GetListOfKeys(); 
305           TTree* kTree;
306           TIter next(kineEventsList); 
307           TKey* key;
308           
309           //loop over the events
310           while( key=(TKey *)next() ) 
311             {
312               TDirectory* tDir = (TDirectory*)key->ReadObj();
313               if(!tDir) break;
314               
315               TString evtDir(tDir->GetName()); 
316               //cout << "  . . found: " << tDir->GetName() << endl;
317               
318               kTree = (TTree *)tDir->Get("TreeK");
319               if(!kTree) break;
320               
321               Int_t nPart = kTree->GetEntries();
322               //cout << "  . . . kTree " << fCount << " has " << nPart << " particles " << endl; 
323               
324               //-----------------------------------------------------------
325               //fill and save the flow event
326               
327               AliFlowEventSimple* fEvent = fEventMaker->FillTracks(kTree, cutsInt, cutsDiff);
328               
329               //pass the flow event to flow methods for analysis 
330               Double_t fEP = 0.; // temporary number need true value
331               if(MCEP) mcep->Make(fEvent,fEP);  //fix fEP
332               if(QC) qc->Make(fEvent);
333               if(GFC) gfc->Make(fEvent);
334               if(FQD) fqd->Make(fEvent);
335               if(LYZ1) lyz1->Make(fEvent);
336               if(LYZ2) lyz2->Make(fEvent);
337               if(LYZEP) lyzep->Make(fEvent,ep);
338               if(SP) sp->Make(fEvent);
339               //-----------------------------------------------------------
340               
341               fCount++;
342               cout << "# " << fCount << " events processed" << endl;
343               delete kTree;
344               delete fEvent;
345             }
346           delete kineFile ;
347         }
348       delete evtsDir ;
349     }
350
351   //--------------------------------------------------------------
352   //calculating and storing the final results of flow analysis
353   if(MCEP) {mcep->Finish(); mcep->WriteHistograms("outputMCEPanalysis.root");}
354   if(QC) {qc->Finish(); qc->WriteHistograms("outputQCanalysis.root");}
355   if(GFC) {gfc->Finish(); gfc->WriteHistograms("outputGFCanalysis.root");}
356   if(FQD) {fqd->Finish(); fqd->WriteHistograms("outputFQDanalysis.root");}
357   if(LYZ1) {lyz1->Finish(); lyz1->WriteHistograms("outputLYZ1analysis.root");}
358   if(LYZ2) {lyz2->Finish(); lyz2->WriteHistograms("outputLYZ2analysis.root");}
359   if(LYZEP) {lyzep->Finish(); lyzep->WriteHistograms("outputLYZEPanalysis.root");}
360   if(SP) {sp->Finish(); sp->WriteHistograms("outputSPanalysis.root");}
361   //--------------------------------------------------------------
362   
363   cout << endl;
364   cout << " Finished ... " << endl;
365   cout << endl;
366   
367   timer.Stop();
368   cout << endl;
369   timer.Print();
370 }