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