few extra options for momentum dependence v_n
[u/mrichter/AliRoot.git] / PWG2 / FLOW / macros / runFlowAnalysisOnTheFly.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    = kTRUE;
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   = kTRUE;
16 Bool_t MCEP  = kTRUE; 
17 //--------------------------------------------------------------------------------------
18
19 Bool_t bSameSeed = kFALSE; // use always the same seed for random generators
20 Bool_t bConstantHarmonics = kTRUE; // harmonics V1, V2, V4... are constant (kTRUE) or functions of pt and eta (kFALSE)
21
22 // Set the event parameters:
23 Int_t iLoops = 1; // number of times to use each track (to simulate nonflow)
24
25 Int_t iMultiplicityOfRP = 500; // multiplicity of RPs
26 Double_t dMultiplicitySpreadOfRP = 0; // multiplicity spread of RPs
27
28 //......................................................................................  
29 // if you use (pt,eta) dependent harmonics (bConstantHarmonics = kFALSE):
30 Double_t dPtCutOff = 2.0; // V2(pt) is linear up to pt = 2 GeV and for pt > 2 GeV it is constant: V2(pt) = dVRPMax
31 Double_t dV2RPMax = 0.2; // maximum value of V2(pt) for pt >= 2GeV
32 //...................................................................................... 
33
34 //......................................................................................  
35 // if you use constant harmonics (bConstantHarmonics = kTRUE):
36 Double_t dV2RP = 0.05; // elliptic flow of RPs
37 Double_t dV2SpreadRP = 0.; // elliptic flow spread of RPs
38
39 Double_t dV1RP = 0.0; // directed flow of RPs
40 Double_t dV1SpreadRP = 0.0; // directed flow spread of RPs
41
42 Double_t dV4RP = 0.0; // harmonic V4 of RPs (to be improved: name needed)
43 Double_t dV4SpreadRP = 0.0; // harmonic V4's spread of RPs (to be improved: name needed)
44 //......................................................................................  
45
46 enum anaModes {mLocal,mLocalSource,mLocalPAR};
47 // mLocal: Analyze data on your computer using aliroot
48 // mLocalPAR: Analyze data on your computer using root + PAR files
49 // mLocalSource: Analyze data on your computer using root + source files
50                                           
51 int runFlowAnalysisOnTheFly(Int_t mode=mLocal, Int_t nEvts=1000)
52 {
53  TStopwatch timer;
54  timer.Start();
55  
56  if (LYZ1 && LYZ2)  {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1.  "<<endl; exit(); }
57  if (LYZ2 && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
58  if (LYZ1 && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
59  
60  cout<<endl;
61  cout<<endl;
62  cout<<"      ---- ARE YOU READY TO FLY ? ----      "<<endl;
63  cout<<endl;
64  
65  cout<<endl;
66  cout<<" ---- BEGIN FLOW ANALYSIS 'ON THE FLY' ---- "<<endl;
67  cout<<endl;
68  cout<<endl;
69  
70  LoadLibraries(mode);
71
72  // Initialize the seed for random generator
73  UInt_t sseed = 0;
74  
75  if(bSameSeed) 
76  {
77   sseed = 44; // the default constant value for seed for random generators
78  } 
79
80  if(!bSameSeed)
81  {
82   TTimeStamp dt;
83   sseed = dt.GetNanoSec()/1000;
84  }
85  
86  //---------------------------------------------------------------------------------------
87  // Initialize the flowevent maker
88  AliFlowEventSimpleMakerOnTheFly* eventMakerOnTheFly = new AliFlowEventSimpleMakerOnTheFly(sseed);
89  eventMakerOnTheFly->Init();
90   
91  //---------------------------------------------------------------------------------------
92  // Initialize all the flow methods:  
93  AliFlowAnalysisWithQCumulants    *qc    = NULL;
94  AliFlowAnalysisWithCumulants     *gfc   = NULL;
95  AliFittingQDistribution          *fqd   = NULL;
96  AliFlowAnalysisWithLeeYangZeros  *lyz1  = NULL;
97  AliFlowAnalysisWithLeeYangZeros  *lyz2  = NULL;
98  AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
99  AliFlowAnalysisWithScalarProduct *sp    = NULL;
100  AliFlowAnalysisWithMCEventPlane  *mcep  = NULL;   
101
102  // MCEP = monte carlo event plane
103  if (MCEP) {
104    AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
105    mcep->Init();
106  }
107
108   // QC = Q-cumulants  
109  if(QC) { 
110    AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
111    qc->Init();
112  }
113   
114  // GFC = Generating Function Cumulants 
115  if(GFC) {
116    AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
117    gfc->Init();
118  }
119  
120  // FQD = Fitting q-distribution 
121  if(FQD) {
122    AliFittingQDistribution* fqd = new AliFittingQDistribution();
123    fqd->Init();
124  }
125  
126  // SP = Scalar Product 
127  if(SP) {
128    AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
129    sp->Init();
130  }
131
132  // LYZ1 = Lee-Yang Zeroes first run
133  if(LYZ1) {
134    AliFlowAnalysisWithLeeYangZeros* lyz1 = new AliFlowAnalysisWithLeeYangZeros();
135    lyz1->SetFirstRun(kTRUE);
136    lyz1->SetUseSum(kTRUE);
137    lyz1->Init();
138  }
139
140  // LYZ2 = Lee-Yang Zeroes second run
141  if(LYZ2) {
142    AliFlowAnalysisWithLeeYangZeros* lyz2 = new AliFlowAnalysisWithLeeYangZeros();
143    // read the input file from the first run 
144    TString inputFileNameLYZ2 = "outputLYZ1analysis.root" ;
145    TFile* inputFileLYZ2 = new TFile(inputFileNameLYZ2.Data(),"READ");
146    if(!inputFileLYZ2 || inputFileLYZ2->IsZombie()) { 
147      cerr << " ERROR: NO First Run file... " << endl ;
148      break; 
149    }
150    else { 
151      TList* inputListLYZ2 = (TList*)inputFileLYZ2->Get("cobjLYZ1");  
152      if (!inputListLYZ2) {cout<<"Input list is NULL pointer!"<<endl; break;}
153      else {
154        cout<<"LYZ2 input file/list read..."<<endl;
155        lyz2->SetFirstRunList(inputListLYZ2);
156        lyz2->SetFirstRun(kFALSE);
157        lyz2->SetUseSum(kTRUE);
158        lyz2->Init();
159      }
160    }
161  }
162  
163  // LYZEP = Lee-Yang Zeroes event plane
164  if(LYZEP) {
165    AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
166    AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
167    // read the input file from the second lyz run 
168    TString inputFileNameLYZEP = "outputLYZ2analysis.root" ;
169    TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
170    if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) { 
171      cerr << " ERROR: NO Second Run file... " << endl ; 
172      break;
173    }
174    else { 
175      TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2");  
176      if (!inputListLYZEP) {cout<<"Input list is NULL pointer!"<<endl; break;}
177      else {
178        cout<<"LYZEP input file/list read..."<<endl;
179        ep   ->SetSecondRunList(inputListLYZEP);
180        lyzep->SetSecondRunList(inputListLYZEP);
181        ep   ->Init();
182        lyzep->Init();
183      }
184    }
185  }
186  //---------------------------------------------------------------------------------------
187   
188  // set the global event parameters: 
189  eventMakerOnTheFly->SetNoOfLoops(iLoops);
190  eventMakerOnTheFly->SetMultiplicityOfRP(iMultiplicityOfRP);
191  eventMakerOnTheFly->SetMultiplicitySpreadOfRP(dMultiplicitySpreadOfRP);
192
193  eventMakerOnTheFly->SetV1RP(dV1RP);
194  eventMakerOnTheFly->SetV1SpreadRP(dV1SpreadRP);  
195  eventMakerOnTheFly->SetV4RP(dV4RP);
196  eventMakerOnTheFly->SetV4SpreadRP(dV4SpreadRP);  
197  
198  // constant harmonic V2:
199  if(bConstantHarmonics)
200  { 
201   eventMakerOnTheFly->SetUseConstantHarmonics(bConstantHarmonics);
202   eventMakerOnTheFly->SetV2RP(dV2RP);
203   eventMakerOnTheFly->SetV2SpreadRP(dV2SpreadRP);  
204  }
205  // (pt,eta) dependent harmonic V2:
206  if(!bConstantHarmonics)
207  {
208   eventMakerOnTheFly->SetUseConstantHarmonics(bConstantHarmonics);
209   eventMakerOnTheFly->SetV2RPMax(dV2RPMax);
210   eventMakerOnTheFly->SetPtCutOff(dPtCutOff);  
211  }
212        
213  //---------------------------------------------------------------------------------------  
214  // create and analyze events 'on the fly':
215
216  for(Int_t i=0;i<nEvts;i++) {   
217    // creating the event with above settings:
218    AliFlowEventSimple *event = eventMakerOnTheFly->CreateEventOnTheFly(); 
219    
220    // analyzing the created event 'on the fly':
221    // do flow analysis for various methods:
222    if(MCEP) mcep->Make(event);
223    if(QC) qc->Make(event);
224    if(GFC) gfc->Make(event);
225    if(FQD) fqd->Make(event);
226    if(LYZ1) lyz1->Make(event);
227    if(LYZ2) lyz2->Make(event);
228    if(LYZEP) lyzep->Make(event,ep);
229    if(SP) sp->Make(event);
230    
231    delete event;
232  } // end of for(Int_t i=0;i<nEvts;i++)
233  //---------------------------------------------------------------------------------------  
234
235
236
237  //---------------------------------------------------------------------------------------  
238  // calculating and storing the final results of flow analysis
239  if(MCEP) {mcep->Finish(); mcep->WriteHistograms("outputMCEPanalysis.root");}
240  if(SP) {sp->Finish(); sp->WriteHistograms("outputSPanalysis.root");}
241  if(QC) {qc->Finish(); qc->WriteHistograms("outputQCanalysis.root");}
242  if(GFC) {gfc->Finish(); gfc->WriteHistograms("outputGFCanalysis.root");}
243  if(FQD) {fqd->Finish(); fqd->WriteHistograms("outputFQDanalysis.root");}
244  if(LYZ1) {lyz1->Finish(); lyz1->WriteHistograms("outputLYZ1analysis.root");}
245  if(LYZ2) {lyz2->Finish(); lyz2->WriteHistograms("outputLYZ2analysis.root");}
246  if(LYZEP) {lyzep->Finish(); lyzep->WriteHistograms("outputLYZEPanalysis.root");}
247  //---------------------------------------------------------------------------------------  
248  
249  
250  
251  cout<<endl;
252  cout<<endl;
253  cout<<" ---- LANDED SUCCESSFULLY ---- "<<endl;
254  cout<<endl; 
255  
256  timer.Stop();
257  cout << endl;
258  timer.Print();
259 }
260
261 void SetupPar(char* pararchivename)
262 {
263   //Load par files, create analysis libraries
264   //For testing, if par file already decompressed and modified
265   //classes then do not decompress.
266  
267   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ; 
268   TString parpar(Form("%s.par", pararchivename)) ; 
269   if ( gSystem->AccessPathName(parpar.Data()) ) {
270     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
271     TString processline(Form(".! make %s", parpar.Data())) ; 
272     gROOT->ProcessLine(processline.Data()) ;
273     gSystem->ChangeDirectory(cdir) ; 
274     processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
275     gROOT->ProcessLine(processline.Data()) ;
276   } 
277   if ( gSystem->AccessPathName(pararchivename) ) {  
278     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
279     gROOT->ProcessLine(processline.Data());
280   }
281   
282   TString ocwd = gSystem->WorkingDirectory();
283   gSystem->ChangeDirectory(pararchivename);
284   
285   // check for BUILD.sh and execute
286   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
287     printf("*******************************\n");
288     printf("*** Building PAR archive    ***\n");
289     cout<<pararchivename<<endl;
290     printf("*******************************\n");
291     
292     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
293       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
294       return -1;
295     }
296   }
297   // check for SETUP.C and execute
298   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
299     printf("*******************************\n");
300     printf("*** Setup PAR archive       ***\n");
301     cout<<pararchivename<<endl;
302     printf("*******************************\n");
303     gROOT->Macro("PROOF-INF/SETUP.C");
304   }
305   
306   gSystem->ChangeDirectory(ocwd.Data());
307   printf("Current dir: %s\n", ocwd.Data());
308 }
309
310 void LoadLibraries(const anaModes mode) {
311   
312   //--------------------------------------
313   // Load the needed libraries most of them already loaded by aliroot
314   //--------------------------------------
315   gSystem->Load("libTree.so");
316   gSystem->Load("libGeom.so");
317   gSystem->Load("libVMC.so");
318   gSystem->Load("libXMLIO.so");
319   gSystem->Load("libPhysics.so");
320   
321   //----------------------------------------------------------
322   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<< 
323   //----------------------------------------------------------
324   if (mode==mLocal) {
325     //--------------------------------------------------------
326     // If you want to use already compiled libraries 
327     // in the aliroot distribution
328     //--------------------------------------------------------
329     gSystem->Load("libSTEERBase");
330     gSystem->Load("libESD");
331     gSystem->Load("libAOD");
332     gSystem->Load("libANALYSIS");
333     gSystem->Load("libANALYSISalice");
334     gSystem->Load("libCORRFW.so");
335     cerr<<"libCORRFW.so loaded..."<<endl;
336     gSystem->Load("libPWG2flowCommon.so");
337     cerr<<"libPWG2flowCommon.so loaded..."<<endl;
338     gSystem->Load("libPWG2flowTasks.so");
339     cerr<<"libPWG2flowTasks.so loaded..."<<endl;
340   }
341   
342   else if (mode == mLocalPAR) {
343     //--------------------------------------------------------
344     //If you want to use root and par files from aliroot
345     //--------------------------------------------------------  
346      //If you want to use root and par files from aliroot
347     //--------------------------------------------------------  
348     SetupPar("STEERBase");
349     SetupPar("ESD");
350     SetupPar("AOD");
351     SetupPar("ANALYSIS");
352     SetupPar("ANALYSISalice");
353     SetupPar("PWG2AOD");
354     SetupPar("CORRFW");
355     SetupPar("PWG2flowCommon");
356     cerr<<"PWG2flowCommon.par loaded..."<<endl;
357     SetupPar("PWG2flowTasks");
358     cerr<<"PWG2flowTasks.par loaded..."<<endl;
359   }
360   
361   //---------------------------------------------------------
362   // <<<<<<<<<< Source mode >>>>>>>>>>>>
363   //---------------------------------------------------------
364   else if (mode==mLocalSource) {
365  
366     // In root inline compile
367
368    
369     // Constants  
370     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
371     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
372     gROOT->LoadMacro("AliFlowCommon/AliFlowCumuConstants.cxx+");
373     
374     // Flow event
375     gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+"); 
376     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");    
377     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
378     
379     // Cuts
380     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");    
381     
382     // Output histosgrams
383     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
384     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
385     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
386     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
387     
388     // Functions needed for various methods
389     gROOT->LoadMacro("AliFlowCommon/AliCumulantsFunctions.cxx+");
390     gROOT->LoadMacro("AliFlowCommon/AliFittingFunctionsForQDistribution.cxx+");
391     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZEventPlane.cxx+");
392     
393     // Flow Analysis code for various methods
394     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMCEventPlane.cxx+"); 
395     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx+");
396     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.cxx+");
397     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.cxx+");
398     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithCumulants.cxx+");
399     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithQCumulants.cxx+"); 
400     gROOT->LoadMacro("AliFlowCommon/AliFittingQDistribution.cxx+");
401     
402     // Class to fill the FlowEvent on the fly (generate Monte Carlo events)
403     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.cxx+");   
404     
405     cout << "finished loading macros!" << endl;  
406     
407   }  
408   
409 }
410
411