hooks for PMD flow analysis
[u/mrichter/AliRoot.git] / PWG2 / FLOW / macros / runFlowAnalysis.C
1 #include <Riostream.h>
2 #include <TStopwatch.h>
3 #include <TObjArray.h>
4 #include <TFile.h>
5 #include <TRandom3.h>
6 #include <TTimeStamp.h>
7
8 #include <AliFlowCommon/AliFlowCommonConstants.h>
9 #include <AliFlowCommon/AliFlowLYZConstants.h>
10 #include <AliFlowCommon/AliFlowCumuConstants.h>
11 #include <AliFlowCommon/AliFlowVector.h>
12 #include <AliFlowCommon/AliFlowTrackSimple.h>
13 #include <AliFlowCommon/AliFlowEvent.h>
14 #include <AliFlowCommon/AliFlowEventSimple.h>
15 #include <AliFlowCommon/AliFlowTrackSimpleCuts.h>
16 #include <AliFlowCommon/AliFlowCommonHist.h>
17 #include <AliFlowCommon/AliFlowCommonHistResults.h>
18 #include <AliFlowCommon/AliFlowLYZHist1.h>
19 #include <AliFlowCommon/AliFlowLYZHist2.h>
20 #include <AliFlowCommon/AliCumulantsFunctions.h>
21 #include <AliFlowCommon/AliFlowLYZEventPlane.h>
22 #include <AliFlowCommon/AliFlowAnalysisWithMCEventPlane.h>
23 #include <AliFlowCommon/AliFlowAnalysisWithScalarProduct.h>
24 #include <AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.h>
25 #include <AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.h>
26 #include <AliFlowCommon/AliFlowAnalysisWithCumulants.h>
27 #include <AliFlowCommon/AliFlowAnalysisWithQCumulants.h>
28 #include <AliFlowCommon/AliFlowAnalysisWithFittingQDistribution.h>
29 #include <AliFlowCommon/AliFlowAnalysisWithMixedHarmonics.h>
30 #include <AliFlowCommon/AliFlowAnalysisWithNestedLoops.h>
31
32 //--------------------------------------------------------------------------------------
33 // Run flow analysis on local data with custom FlowEvent maker
34 // RUN SETTINGS
35 //flow analysis method can be: (set to kTRUE or kFALSE)
36 Bool_t SP       = kTRUE;
37 Bool_t LYZ1SUM  = kTRUE;
38 Bool_t LYZ1PROD = kTRUE;
39 Bool_t LYZ2SUM  = kFALSE; 
40 Bool_t LYZ2PROD = kFALSE;
41 Bool_t LYZEP    = kFALSE; 
42 Bool_t GFC      = kTRUE;
43 Bool_t QC       = kTRUE;
44 Bool_t FQD      = kTRUE;
45 Bool_t MH       = kTRUE; 
46 Bool_t NL       = kFALSE; 
47 Bool_t MCEP     = kFALSE; //does not work yet 24/12/08
48 //--------------------------------------------------------------------------------------
49
50 // Weights 
51 // Use weights for Q vector
52 Bool_t usePhiWeights = kFALSE; //Phi (correction for non-uniform azimuthal acceptance)
53 Bool_t usePtWeights  = kFALSE; //v'(pt) (differential flow in pt)
54 Bool_t useEtaWeights = kFALSE; //v'(eta) (differential flow in eta)
55
56 //--------------------------------------------------------------------------------------
57 // CUT SETTINGS
58 //integrated selection
59 Double_t ptMaxInt  = 10.;
60 Double_t ptMinInt  = 0.;
61 Double_t etaMaxInt = 1.;
62 Double_t etaMinInt = -1.;
63 Double_t phiMaxInt = 7.5;
64 Double_t phiMinInt = 0.;
65 Int_t PIDInt       = 211;
66
67 //differential selection
68 Double_t ptMaxDiff  = 10.;
69 Double_t ptMinDiff  = 0.;
70 Double_t etaMaxDiff = 1.;
71 Double_t etaMinDiff = -1.;
72 Double_t phiMaxDiff = 7.5;
73 Double_t phiMinDiff = 0.;
74 Int_t PIDDiff       = 211;
75 /*
76 //--------------------------------------------------------------------------------------
77 // FLOW SETTINGS (R.Rietkerk)
78 Int_t nLoops=1;                 // Number of times to use the same particle (nonflow).
79 Double_t xEllipticFlowValue=0.1;// Add Elliptic Flow. Must be in range [0,1].
80 Int_t nMultiplicityOfEvent=500; // Set Average Multiplicity.
81 Double_t xSigmaFlow=0.00;       // Add Elliptic Flow. Must be in range [0,1].
82 Int_t nSigmaMult=50;            // Set Average Multiplicity.
83 //--------------------------------------------------------------------------------------
84 */
85
86 enum anaModes {mLocal,mLocalSource,mLocalPAR,};
87 //mLocal: Analyze data on your computer using aliroot
88 //mLocalPAR: Analyze data on your computer using root + PAR files
89 //mLocalSource: Analyze data on your computer using root + source files
90
91 void LoadLibraries(const anaModes mode);
92
93 Int_t offset = 0;
94                                           
95 int runFlowAnalysis(const anaModes mode=mLocal, Int_t aRuns = 100, const char* 
96                     dir="/data/alice1/kolk/KineOnly3/")
97                     //              dir="/Users/snelling/alice_data/KineOnly3/")
98                     //              dir="/Users/snelling/alice_data/stoomboot/5b/")
99 {
100   TStopwatch timer;
101   timer.Start();
102   
103   if (LYZ1SUM && LYZ2SUM) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(1); }
104   if (LYZ1PROD && LYZ2PROD) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(1); }
105   if (LYZ2SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(1); }
106   if (LYZ1SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(1); }
107
108
109   cout<<endl;
110   cout<<" ---- BEGIN ANALYSIS ---- "<<endl;
111   cout<<endl;
112   
113   LoadLibraries(mode);
114
115   TRandom3 random3Temp; //init for manual settings (R.Rietkerk)
116   TTimeStamp dt;
117   Int_t sseed = dt.GetNanoSec()/1000;
118   random3Temp.SetSeed(sseed);
119
120   if (mode == mLocal || mode == mLocalPAR) {
121   }
122   else if (mode == mLocalSource) {
123   }
124   else{
125     cout << "No supported running mode selected!" << endl;
126     break;
127   }
128
129
130   //------------------------------------------------------------------------
131   //cuts:
132   AliFlowTrackSimpleCuts* cutsInt = new AliFlowTrackSimpleCuts();
133   cutsInt->SetPtMax(ptMaxInt);
134   cutsInt->SetPtMin(ptMinInt);
135   cutsInt->SetEtaMax(etaMaxInt);
136   cutsInt->SetEtaMin(etaMinInt);
137   cutsInt->SetPhiMax(phiMaxInt);
138   cutsInt->SetPhiMin(phiMinInt);
139   cutsInt->SetPID(PIDInt);
140
141   AliFlowTrackSimpleCuts* cutsDiff = new AliFlowTrackSimpleCuts();
142   cutsDiff->SetPtMax(ptMaxDiff);
143   cutsDiff->SetPtMin(ptMinDiff);
144   cutsDiff->SetEtaMax(etaMaxDiff);
145   cutsDiff->SetEtaMin(etaMinDiff);
146   cutsDiff->SetPhiMax(phiMaxDiff);
147   cutsDiff->SetPhiMin(phiMinDiff);
148   cutsDiff->SetPID(PIDDiff);
149
150   //if the weights are used: 
151   TFile *fileWithWeights = NULL;
152   TList *listWithWeights = NULL;
153   
154   if(usePhiWeights||usePtWeights||useEtaWeights) {
155     fileWithWeights = TFile::Open("weights.root","READ");
156     if(fileWithWeights) {
157       listWithWeights = (TList*)fileWithWeights->Get("weights");
158     }
159     else
160       {cout << " WARNING: the file <weights.root> with weights from the previous run was not found."<<endl;
161         break;
162       }    
163   }
164
165   //flow methods:  
166   AliFlowAnalysisWithQCumulants *qc = NULL;
167   AliFlowAnalysisWithCumulants *gfc = NULL;
168   AliFlowAnalysisWithFittingQDistribution *fqd = NULL;
169   AliFlowAnalysisWithLeeYangZeros *lyz1sum = NULL;
170   AliFlowAnalysisWithLeeYangZeros *lyz1prod = NULL;
171   AliFlowAnalysisWithLeeYangZeros *lyz2sum = NULL;
172   AliFlowAnalysisWithLeeYangZeros *lyz2prod = NULL;
173   AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
174   AliFlowAnalysisWithScalarProduct *sp = NULL;
175   AliFlowAnalysisWithMCEventPlane *mcep = NULL;     
176   AliFlowAnalysisWithMixedHarmonics *mh = NULL;
177   AliFlowAnalysisWithNestedLoops *nl = NULL;
178
179   //MCEP = monte carlo event plane
180   if (MCEP) {
181     AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
182     mcep->Init();
183   }
184
185   //QC = Q-cumulants  
186   if(QC) { 
187     AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
188     if(listWithWeights) qc->SetWeightsList(listWithWeights);
189     if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
190     if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
191     if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
192     qc->Init();
193   }
194   
195   //GFC = Generating Function Cumulants 
196   if(GFC) {
197     AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
198     if(listWithWeights) gfc->SetWeightsList(listWithWeights);
199     if(usePhiWeights) gfc->SetUsePhiWeights(usePhiWeights);
200     if(usePtWeights) gfc->SetUsePtWeights(usePtWeights);
201     if(useEtaWeights) gfc->SetUseEtaWeights(useEtaWeights);
202     gfc->Init();
203   }
204   
205   //FQD = Fitting q-distribution 
206   if(FQD) {
207     AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
208     if(listWithWeights) fqd->SetWeightsList(listWithWeights);
209     if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);
210     fqd->Init();
211   }
212
213   //SP = Scalar Product 
214   if(SP) {
215     AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
216     if(usePhiWeights) sp->SetUsePhiWeights(usePhiWeights);
217     sp->Init();
218   }
219
220   //LYZ1 = Lee-Yang Zeroes first run
221   if(LYZ1SUM) {
222     AliFlowAnalysisWithLeeYangZeros* lyz1sum = new AliFlowAnalysisWithLeeYangZeros();
223     lyz1sum->SetFirstRun(kTRUE);
224     lyz1sum->SetUseSum(kTRUE);
225     lyz1sum->Init();
226   }
227   if(LYZ1PROD) {
228     AliFlowAnalysisWithLeeYangZeros* lyz1prod = new AliFlowAnalysisWithLeeYangZeros();
229     lyz1prod->SetFirstRun(kTRUE);
230     lyz1prod->SetUseSum(kFALSE);
231     lyz1prod->Init();
232   }
233   //LYZ2 = Lee-Yang Zeroes second run
234   if(LYZ2SUM) {
235     AliFlowAnalysisWithLeeYangZeros* lyz2sum = new AliFlowAnalysisWithLeeYangZeros();
236     // read the input file from the first run 
237     TString inputFileNameLYZ2SUM = "outputLYZ1SUManalysis.root" ;
238     TFile* inputFileLYZ2SUM = new TFile(inputFileNameLYZ2SUM.Data(),"READ");
239     if(!inputFileLYZ2SUM || inputFileLYZ2SUM->IsZombie()) { 
240       cerr << " ERROR: To run LYZ2SUM you need the output file from LYZ1SUM. This file is not there! Please run LYZ1SUM first." << endl ;
241       break; 
242     }
243     else { 
244       TList* inputListLYZ2SUM = (TList*)inputFileLYZ2SUM->Get("cobjLYZ1SUM");  
245       if (!inputListLYZ2SUM) {cout<<"SUM Input list is NULL pointer!"<<endl; break;}
246       else {
247         cout<<"LYZ2SUM input file/list read..."<<endl;
248         lyz2sum->SetFirstRunList(inputListLYZ2SUM);
249         lyz2sum->SetFirstRun(kFALSE);
250         lyz2sum->SetUseSum(kTRUE);
251         lyz2sum->Init();
252       }
253     }
254   }
255   if(LYZ2PROD) {
256     AliFlowAnalysisWithLeeYangZeros* lyz2prod = new AliFlowAnalysisWithLeeYangZeros();
257     // read the input file from the first run 
258     TString inputFileNameLYZ2PROD = "outputLYZ1PRODanalysis.root" ;
259     TFile* inputFileLYZ2PROD = new TFile(inputFileNameLYZ2PROD.Data(),"READ");
260     if(!inputFileLYZ2PROD || inputFileLYZ2PROD->IsZombie()) { 
261       cerr << " ERROR: To run LYZ2PROD you need the output file from LYZ1PROD. This file is not there! Please run LYZ1PROD first." << endl ;
262       break; 
263     }
264     else { 
265       TList* inputListLYZ2PROD = (TList*)inputFileLYZ2PROD->Get("cobjLYZ1PROD");  
266       if (!inputListLYZ2PROD) {cout<<"PROD Input list is NULL pointer!"<<endl; break;}
267       else {
268         cout<<"LYZ2PROD input file/list read..."<<endl;
269         lyz2prod->SetFirstRunList(inputListLYZ2PROD);
270         lyz2prod->SetFirstRun(kFALSE);
271         lyz2prod->SetUseSum(kTRUE);
272         lyz2prod->Init();
273       }
274     }
275   }
276  //LYZEP = Lee-Yang Zeroes event plane
277   if(LYZEP) {
278     AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
279     AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
280     // read the input file from the second lyz run 
281     TString inputFileNameLYZEP = "outputLYZ2SUManalysis.root" ;
282     TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
283     if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) { 
284       cerr << " ERROR: To run LYZEP you need the output file from LYZ2SUM. This file is not there! Please run LYZ2SUM first." << endl ; 
285       break;
286     }
287     else { 
288       TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2SUM");  
289       if (!inputListLYZEP) {cout<<"Input list is NULL pointer!"<<endl; break;}
290       else {
291         cout<<"LYZEP input file/list read..."<<endl;
292         ep   ->SetSecondRunList(inputListLYZEP);
293         lyzep->SetSecondRunList(inputListLYZEP);
294         ep   ->Init();
295         lyzep->Init();
296       }
297     }
298   }
299   // MH = Mixed Harmonics:  
300   if(MH) { 
301     AliFlowAnalysisWithMixedHarmonics* mh = new AliFlowAnalysisWithMixedHarmonics();
302     if(listWithWeights) mh->SetWeightsList(listWithWeights);
303     //if(usePhiWeights) mh->SetUsePhiWeights(usePhiWeights); // to be improved (enabled)
304     //if(usePtWeights) mh->SetUsePtWeights(usePtWeights); // to be improved (enabled)
305     //if(useEtaWeights) mh->SetUseEtaWeights(useEtaWeights); // to be improved (enabled)
306     mh->Init();
307   }
308   // NL = Nested Loops:  
309   if(NL) { 
310     AliFlowAnalysisWithNestedLoops* nl = new AliFlowAnalysisWithNestedLoops();
311     if(listWithWeights) nl->SetWeightsList(listWithWeights);
312     //if(usePhiWeights) nl->SetUsePhiWeights(usePhiWeights); // to be improved (enabled)
313     //if(usePtWeights) nl->SetUsePtWeights(usePtWeights); // to be improved (enabled)
314     //if(useEtaWeights) nl->SetUseEtaWeights(useEtaWeights); // to be improved (enabled)
315     nl->Init();
316   }
317
318   //------------------------------------------------------------------------
319   
320   
321   //standard code to read files in directory
322   Int_t fCount = 0;
323   TString execDir(gSystem->pwd());
324   TString targetDir(dir);
325   TSystemDirectory* baseDir = new TSystemDirectory(".", dir);  
326   TList* dirList           = baseDir->GetListOfFiles();
327   if (!dirList) {
328     cout << endl << "No input files in: " << targetDir.Data() << endl;
329     break;
330   }
331   Int_t nDirs              = dirList->GetEntries();
332   cout<<endl;
333   cout<<"Int_t nDirs = "<<nDirs<<endl;
334   gSystem->cd(execDir);
335   
336   for(Int_t iDir=0;iDir<nDirs;++iDir)
337     {
338       TSystemFile* presentDir = (TSystemFile*)dirList->At(iDir);
339       if(!presentDir || !presentDir->IsDirectory() || 
340          strcmp(presentDir->GetName(), ".") == 0 || 
341          strcmp(presentDir->GetName(), "..") == 0) 
342         {
343           cout << endl; 
344           cout << "Directory (" << iDir << "): " << presentDir->GetName() << 
345             " - Skipping ... " << endl;
346           continue ;   
347         }
348       
349       if(offset > 0) { --offset ; continue ; }
350       if((aRuns > 0) && (fCount >= aRuns)) { break ; }
351       
352       TString presentDirName(dir); // aDataDir
353       presentDirName += presentDir->GetName();
354       presentDirName += "/";
355       //cerr<<" presentDirName = "<<presentDirName<<endl;
356       
357       TString fileName = presentDirName; 
358       fileName += "galice.root";
359       Long_t *id, *size, *flags, *modtime;
360       if(gSystem->GetPathInfo(fileName.Data(),id,size,flags,modtime)) 
361         { 
362           cout << " File : " << fileName << " does NOT exist ! - Skipping ... " 
363                << endl; 
364           continue; 
365         }
366       //      cout << endl ; cout << "Directory (" << iDir << "): " << presentDirName << " ... " << endl;
367       
368       //loop (simulations in the present dir) 
369       TSystemDirectory* evtsDir = new TSystemDirectory(".", 
370                                                        presentDirName.Data());
371       TList* fileList       = evtsDir->GetListOfFiles();
372       Int_t nFiles                  = fileList->GetEntries();
373       //cout<<" Int_t nFiles = "<<nFiles<<endl;
374       gSystem->cd(execDir);      
375       for(Int_t iFiles=0; iFiles<nFiles; ++iFiles)
376         {
377           TSystemFile* presentFile = (TSystemFile*) fileList->At(iFiles);
378           TString presentFileName(presentDirName);
379           presentFileName += presentFile->GetName();
380           
381           if(!(presentFileName.Contains("Kinematics") && 
382                presentFileName.Contains("root"))) { continue ; }
383           
384           //cout << " found: " << presentFileName.Data() << endl; 
385           
386           TFile* kineFile = new TFile(presentFileName.Data(), "READ"); 
387           // kineFile->ls();
388           Int_t nEvts = kineFile->GetNkeys() ; 
389           //cout << "  . found: " << nEvts << " KineTree(s) in " << presentFileName.Data() << endl;
390           TList* kineEventsList = (TList*)kineFile->GetListOfKeys(); 
391           TTree* kTree;
392           TIter next(kineEventsList); 
393           TKey* key;
394
395           Double_t xRPAngle;      
396           //loop over the events
397           while( key=(TKey *)next() ) 
398             {
399               TDirectory* tDir = (TDirectory*)key->ReadObj();
400               if(!tDir) break;
401               
402               TString evtDir(tDir->GetName()); 
403               //cout << "  . . found: " << tDir->GetName() << endl;
404               
405               kTree = (TTree *)tDir->Get("TreeK");
406               if(!kTree) break;
407               
408               Int_t nPart = kTree->GetEntries();
409               //cout << "  . . . kTree " << fCount << " has " << nPart << " particles " << endl; 
410               
411               //-----------------------------------------------------------
412               //fill and save the flow event          
413
414               /*
415               Int_t nNewMultOfEvent = random3Temp.Gaus(nMultiplicityOfEvent,nSigmaMult);
416               cout << "new multiplicity: " << nNewMultOfEvent << endl;
417               Double_t xNewFlowValue = random3Temp.Gaus(xEllipticFlowValue,xSigmaFlow);
418               if ( (fCount % 100) == 0) {
419                 cout << "new multiplicity: " << nNewMultOfEvent << endl;
420                 cout << "new flow value: " << xNewFlowValue << endl;
421               }
422               */
423
424               AliFlowEventSimple *fEvent = new AliFlowEventSimple(kTree, cutsInt, cutsDiff); 
425               //xRPAngle=random3Temp.Uniform(0.0,TMath::TwoPi());
426         //fEvent->SetMCReactionPlaneAngle(xRPAngle);
427         //fEvent->SetV2(xNewFlowValue);
428         //fEvent->CloneTracks(nLoops);
429                             
430               // do flow analysis for various methods
431               if(MCEP)    mcep->Make(fEvent);
432               if(QC)      qc->Make(fEvent);
433               if(GFC)     gfc->Make(fEvent);
434               if(FQD)     fqd->Make(fEvent);
435               if(LYZ1SUM) lyz1sum->Make(fEvent);
436               if(LYZ1PROD)lyz1prod->Make(fEvent);
437               if(LYZ2SUM) lyz2sum->Make(fEvent);
438               if(LYZ2PROD)lyz2prod->Make(fEvent);
439               if(LYZEP)   lyzep->Make(fEvent,ep);
440               if(SP)      sp->Make(fEvent);           
441               if(MH)      mh->Make(fEvent);
442               if(NL)      nl->Make(fEvent);
443               //-----------------------------------------------------------
444               fCount++;
445               //cout << "# " << fCount << " events processed" << endl;
446               delete kTree;
447               delete fEvent;
448             }
449           delete kineFile ;
450         }
451       delete evtsDir ;
452     }
453
454  //---------------------------------------------------------------------------------------  
455  // create a new file which will hold the final results of all methods:
456  TString outputFileName = "AnalysisResults.root";  
457  TFile *outputFile = new TFile(outputFileName.Data(),"RECREATE");
458  // create a new file for each method wich will hold list with final results:
459  const Int_t nMethods = 12;
460  TString method[nMethods] = {"MCEP","SP","GFC","QC","FQD","LYZ1SUM","LYZ1PROD","LYZ2SUM","LYZ2PROD","LYZEP","MH","NL"};
461  TDirectoryFile *dirFileFinal[nMethods] = {NULL};
462  TString fileNameMethod[nMethods]; 
463  for(Int_t i=0;i<nMethods;i++)
464  {
465   // form a file name for each method:
466   fileNameMethod[i]+="output";
467   fileNameMethod[i]+=method[i].Data();
468   fileNameMethod[i]+="analysis";
469   dirFileFinal[i] = new TDirectoryFile(fileNameMethod[i].Data(),fileNameMethod[i].Data());
470  } 
471  
472  // calculating and storing the final results of default flow analysis:
473  if(MCEP)    {mcep->Finish();    mcep->WriteHistograms(dirFileFinal[0]);}
474  if(SP)      {sp->Finish();      sp->WriteHistograms(dirFileFinal[1]);}
475  if(GFC)     {gfc->Finish();     gfc->WriteHistograms(dirFileFinal[2]);}
476  if(QC)      {qc->Finish();      qc->WriteHistograms(dirFileFinal[3]);}
477  if(FQD)     {fqd->Finish();     fqd->WriteHistograms(dirFileFinal[4]);}
478  if(LYZ1SUM) {lyz1sum->Finish(); lyz1sum->WriteHistograms(dirFileFinal[5]);}
479  if(LYZ1PROD){lyz1prod->Finish();lyz1prod->WriteHistograms(dirFileFinal[6]);}
480  if(LYZ2SUM) {lyz2sum->Finish(); lyz2sum->WriteHistograms(dirFileFinal[7]);}
481  if(LYZ2PROD){lyz2prod->Finish();lyz2prod->WriteHistograms(dirFileFinal[8]);}
482  if(LYZEP)   {lyzep->Finish();   lyzep->WriteHistograms(dirFileFinal[9]);}
483  if(MH)      {mh->Finish();      mh->WriteHistograms(dirFileFinal[10]);}
484  if(NL)      {nl->Finish();      nl->WriteHistograms(dirFileFinal[11]);}
485  //---------------------------------------------------------------------------------------  
486  
487  outputFile->Close();
488  delete outputFile;
489   
490   cout << endl;
491   cout << " Fini ... " << endl;
492   cout << endl;
493   
494   timer.Stop();
495   cout << endl;
496   timer.Print();
497 }
498
499 void SetupPar(char* pararchivename)
500 {
501   //Load par files, create analysis libraries
502   //For testing, if par file already decompressed and modified
503   //classes then do not decompress.
504  
505   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ; 
506   TString parpar(Form("%s.par", pararchivename)) ; 
507   if ( gSystem->AccessPathName(parpar.Data()) ) {
508     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
509     TString processline(Form(".! make %s", parpar.Data())) ; 
510     gROOT->ProcessLine(processline.Data()) ;
511     gSystem->ChangeDirectory(cdir) ; 
512     processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
513     gROOT->ProcessLine(processline.Data()) ;
514   } 
515   if ( gSystem->AccessPathName(pararchivename) ) {  
516     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
517     gROOT->ProcessLine(processline.Data());
518   }
519   
520   TString ocwd = gSystem->WorkingDirectory();
521   gSystem->ChangeDirectory(pararchivename);
522   
523   // check for BUILD.sh and execute
524   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
525     printf("*******************************\n");
526     printf("*** Building PAR archive    ***\n");
527     cout<<pararchivename<<endl;
528     printf("*******************************\n");
529     
530     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
531       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
532       return -1;
533     }
534   }
535   // check for SETUP.C and execute
536   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
537     printf("*******************************\n");
538     printf("*** Setup PAR archive       ***\n");
539     cout<<pararchivename<<endl;
540     printf("*******************************\n");
541     gROOT->Macro("PROOF-INF/SETUP.C");
542   }
543   
544   gSystem->ChangeDirectory(ocwd.Data());
545   printf("Current dir: %s\n", ocwd.Data());
546 }
547
548 void LoadLibraries(const anaModes mode) {
549   
550   //--------------------------------------
551   // Load the needed libraries most of them already loaded by aliroot
552   //--------------------------------------
553   //gSystem->Load("libTree");
554   gSystem->Load("libGeom");
555   gSystem->Load("libVMC");
556   gSystem->Load("libXMLIO");
557   gSystem->Load("libPhysics");
558   
559   //----------------------------------------------------------
560   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<< 
561   //----------------------------------------------------------
562   if (mode==mLocal) {
563     //--------------------------------------------------------
564     // If you want to use already compiled libraries 
565     // in the aliroot distribution
566     //--------------------------------------------------------
567     gSystem->Load("libSTEERBase");
568     gSystem->Load("libESD");
569     gSystem->Load("libAOD");
570     gSystem->Load("libANALYSIS");
571     gSystem->Load("libANALYSISalice");
572     gSystem->Load("libCORRFW");
573     cerr<<"libCORRFW loaded..."<<endl;
574     gSystem->Load("libPWG2flowCommon");
575     cerr<<"libPWG2flowCommon loaded..."<<endl;
576     gSystem->Load("libPWG2flowTasks");
577     cerr<<"libPWG2flowTasks loaded..."<<endl;
578   }
579   
580   else if (mode == mLocalPAR) {
581     //--------------------------------------------------------
582     //If you want to use root and par files from aliroot
583     //--------------------------------------------------------  
584      //If you want to use root and par files from aliroot
585     //--------------------------------------------------------  
586     SetupPar("STEERBase");
587     SetupPar("ESD");
588     SetupPar("AOD");
589     SetupPar("ANALYSIS");
590     SetupPar("ANALYSISalice");
591     SetupPar("PWG2AOD");
592     SetupPar("CORRFW");
593     SetupPar("PWG2flowCommon");
594     cerr<<"PWG2flowCommon.par loaded..."<<endl;
595     SetupPar("PWG2flowTasks");
596     cerr<<"PWG2flowTasks.par loaded..."<<endl;
597   }
598   
599   //---------------------------------------------------------
600   // <<<<<<<<<< Source mode >>>>>>>>>>>>
601   //---------------------------------------------------------
602   else if (mode==mLocalSource) {
603  
604     // In root inline compile
605
606    
607     // Constants  
608     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
609     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
610     
611     // Flow event
612     gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+"); 
613     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");    
614     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");    
615     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
616       
617     // Output histosgrams
618     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
619     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
620     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
621     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
622     
623     // Functions needed for various methods
624     gROOT->LoadMacro("AliFlowCommon/AliCumulantsFunctions.cxx+");
625     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZEventPlane.cxx+");
626     
627     // Flow Analysis code for various methods
628     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMCEventPlane.cxx+"); 
629     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx+");
630     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.cxx+");
631     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.cxx+");
632     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithCumulants.cxx+");
633     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithQCumulants.cxx+"); 
634     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithFittingQDistribution.cxx+");    
635     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMixedHarmonics.cxx+");
636     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithNestedLoops.cxx+");
637     
638     cout << "finished loading macros!" << endl;  
639     
640   }  
641   
642 }
643
644