more for weigths
[u/mrichter/AliRoot.git] / PWG2 / FLOW / macros / runFlowAnalysis.C
1 #include "TStopwatch.h"
2 #include "TObjArray"
3 #include "Riostream.h"
4 #include "TFile.h"
5
6 //----------------------------------------------------------
7 //RUN SETTINGS
8 //flow analysis method can be: (set to kTRUE or kFALSE)
9 Bool_t SP    = kFALSE;
10 Bool_t LYZ1  = kTRUE;
11 Bool_t LYZ2  = kFALSE;  
12 Bool_t LYZEP = kFALSE; 
13 Bool_t GFC   = kTRUE;
14 Bool_t QC    = kTRUE;
15 Bool_t FQD   = kFALSE;
16 Bool_t MCEP  = kFALSE; //does not work yet 24/12/08
17 //----------------------------------------------------------
18
19 //----------------------------------------------------------
20 //CUT SETTINGS
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.;
28 Int_t PIDInt       = 211;
29
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.;
37 Int_t PIDDiff       = 211;
38 //----------------------------------------------------------
39
40 //----------------------------------------------------------
41 //WEIGHTS SETTINGS 
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 //----------------------------------------------------------
47
48 Int_t offset = 0;
49                                           
50 int runFlowAnalysis(Int_t aRuns = 44, const char* 
51                                                   dir="/data/alice1/kolk/KineOnly3/")
52                           //dir="/Users/snelling/alice_data/KineOnly3/")
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(); }
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(); }
60
61
62   cout<<endl;
63   cout<<" ---- BEGIN ANALYSIS ---- "<<endl;
64   cout<<endl;
65   
66
67   gSystem->AddIncludePath("-I$ROOTSYS/include");
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; 
76   
77   // flow event in AliRoot
78   AliFlowEventSimpleMaker* fEventMaker = new AliFlowEventSimpleMaker();
79   
80   /*   
81   //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
82   //                       !!!! to be removed
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   }
91   //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
92   */  
93   
94   
95   // In root
96
97   /*    
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+");
114   gROOT->LoadMacro("code/AliFlowLYZHist1.cxx+");fPhiWeights
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
132   // Class to fill the FlowEvent without aliroot dependence
133   // can be found in the directory other
134   gROOT->LoadMacro("code/FlowEventSimpleMaker.cxx+");   
135
136   cout << "finished loading macros!" << endl;  
137
138   //flow event
139   FlowEventSimpleMaker* fEventMaker = new FlowEventSimpleMaker(); 
140   //------------------------------------------------------------------------
141   */
142
143   //load needed libraries
144   gSystem->Load("libTree.so");
145
146   //------------------------------------------------------------------------
147   //cuts:
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
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;
172   
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:  
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;   
190
191   //MCEP = monte carlo event plane
192   if (MCEP) {
193     AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
194     mcep->Init();
195   }
196
197   //QC = Q-cumulants  
198   if(QC) { 
199     AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
200     qc->Init();
201     if(useWeightsPhi) qc->SetPhiWeights(*phiWeights);
202     if(useWeightsPt) qc->SetPtWeights(*ptWeights);
203     if(useWeightsEta) qc->SetEtaWeights(*etaWeights);    
204   }
205   
206   //GFC = Generating Function Cumulants 
207   if(GFC) {
208     AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
209     gfc->CreateOutputObjects();
210   }
211   
212   //FQD = Fitting q-distribution 
213   if(FQD) {
214     AliFittingQDistribution* fqd = new AliFittingQDistribution();
215     fqd->CreateOutputObjects();
216   }
217
218   //SP = Scalar Product 
219   if(SP) {
220     AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
221     sp->Init();
222   }
223
224   //LYZ1 = Lee-Yang Zeroes first run
225   if(LYZ1) {
226     AliFlowAnalysisWithLeeYangZeros* lyz1 = new AliFlowAnalysisWithLeeYangZeros();
227     lyz1->SetFirstRun(kTRUE);
228     lyz1->SetUseSum(kTRUE);
229     lyz1->Init();
230   }
231
232   //LYZ2 = Lee-Yang Zeroes second run
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();
250       }
251     }
252   }
253   
254  //LYZEP = Lee-Yang Zeroes event plane
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();
272       }
273     }
274   }
275   //------------------------------------------------------------------------
276   
277   
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         }
323       //      cout << endl ; cout << "Directory (" << iDir << "): " << presentDirName << " ... " << endl;
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();
330       //cout<<" Int_t nFiles = "<<nFiles<<endl;
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           
341           //cout << " found: " << presentFileName.Data() << endl; 
342           
343           TFile* kineFile = new TFile(presentFileName.Data(), "READ"); 
344           // kineFile->ls();
345           Int_t nEvts = kineFile->GetNkeys() ; 
346           //cout << "  . found: " << nEvts << " KineTree(s) in " << presentFileName.Data() << endl;
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()); 
359               //cout << "  . . found: " << tDir->GetName() << endl;
360               
361               kTree = (TTree *)tDir->Get("TreeK");
362               if(!kTree) break;
363               
364               Int_t nPart = kTree->GetEntries();
365               //cout << "  . . . kTree " << fCount << " has " << nPart << " particles " << endl; 
366               
367               //-----------------------------------------------------------
368               //fill and save the flow event          
369               AliFlowEventSimple *fEvent = fEventMaker->FillTracks(kTree, cutsInt, cutsDiff); 
370                             
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               //-----------------------------------------------------------
382               
383               fCount++;
384               cout << "# " << fCount << " events processed" << endl;
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
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   //--------------------------------------------------------------
404   
405   cout << endl;
406   cout << " Finished ... " << endl;
407   cout << endl;
408   
409   timer.Stop();
410   cout << endl;
411   timer.Print();
412 }