default: run all methods
[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 MCEP     = kTRUE;
10 Bool_t SP       = kTRUE;
11 Bool_t GFC      = kTRUE;
12 Bool_t QC       = kTRUE;
13 Bool_t FQD      = kTRUE;
14 Bool_t LYZ1SUM  = kTRUE;
15 Bool_t LYZ1PROD = kTRUE;
16 Bool_t LYZ2SUM  = kFALSE;
17 Bool_t LYZ2PROD = kFALSE;
18 Bool_t LYZEP    = kFALSE;
19 //--------------------------------------------------------------------------------------
20
21 // Weights 
22 // use weights for Q-vector:
23 Bool_t usePhiWeights = kFALSE; // phi weights (correction for non-uniform azimuthal acceptance)
24 Bool_t usePtWeights  = kFALSE; // pt weights 
25 Bool_t useEtaWeights = kFALSE; // eta weights
26
27 // Run same flow analysis method but with different settings/aims
28 // You will have to label each setting/aim with your own label (see examples bellow): 
29 Bool_t GFC_Additional_Analysis = kFALSE;
30 Bool_t QC_Additional_Analysis  = kFALSE;
31 Bool_t FQD_Additional_Analysis = kFALSE;
32
33 // Define the range for eta subevents
34 Double_t minA = -0.9;
35 Double_t maxA = -0.01;
36 Double_t minB = 0.01;
37 Double_t maxB = 0.9;
38
39 // Parameters for the simulation of events 'on the fly': 
40 Double_t dTemperatureOfRP = 0.44; // 'temperature' of RPs in GeV/c (increase this parameter to get more high pt RPs) 
41
42 //===SEED========================================================
43 // use always the same seed for random generators. 
44 // usage of same seed (kTRUE) is relevant in two cases:
45 // 1.) If you want to use LYZ method to calcualte differential flow;
46 // 2.) If you want to use phi weights for GFC, QC and FQD
47 Bool_t bSameSeed = kTRUE;  
48
49
50 //===NONFLOW=============================================
51 // number of times to use each track (to simulate nonflow)
52 Int_t iLoops = 1; 
53
54
55 //===FLOW HARMONICS===============================================================
56 // harmonics V1, V2, V4... are constant (kTRUE) or functions of pt and eta (kFALSE)                     
57 Bool_t bConstantHarmonics = kTRUE; 
58
59 // if you use (pt,eta) dependent harmonics (bConstantHarmonics = kFALSE):
60 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
61 Double_t dV2RPMax = 0.20; // maximum value of V2(pt) for pt >= 2GeV
62
63 // 1.) if kTRUE  = elliptic flow of RPs is sampled e-b-e from Gaussian distribution with
64 //                 mean = dV2RP and spread = dV2SpreadRP
65 // 2.) if kFALSE = elliptic flow of RPs is sampled e-b-e uniformly from 
66 //                 interval [dMinV2RP,dMaxV2RP]
67 // 3.) for a fixed elliptic flow use Gaussian with zero spread or use uniform with dMinV2RP=dMaxV2RP
68 Bool_t bV2DistrOfRPsIsGauss = kTRUE; 
69
70 // if you use constant harmonics (bConstantHarmonics = kTRUE) (i.e. no pt dependence):
71 Double_t dV2RP = 0.05;       // elliptic flow of RPs (if sampled from Gaussian)
72 Double_t dV2SpreadRP = 0.0;  // elliptic flow spread of RPs (if sampled from Gaussian)
73 Double_t dMinV2RP = 0.02;    // minimal elliptic flow of RPs (if sampled uniformly)
74 Double_t dMaxV2RP = 0.08;    // maximal elliptic flow of RPs (if sampled uniformly)
75
76 Double_t dV1RP = 0.0; // directed flow of RPs
77 Double_t dV1SpreadRP = 0.0; // directed flow spread of RPs
78
79 Double_t dV4RP = 0.0; // harmonic V4 of RPs (to be improved: name needed)
80 Double_t dV4SpreadRP = 0.0; // harmonic V4's spread of RPs (to be improved: name needed)
81
82
83 //===MULTIPLICITY===============================================================
84 // 1.) if kTRUE  = multiplicitiy of RPs is sampled e-b-e from Gaussian distribution with
85 //                 mean = iMultiplicityOfRP and spread = dMultiplicitySpreadOfRP
86 // 2.) if kFALSE = multiplicitiy of RPs is sampled e-b-e uniformly from 
87 //                 interval [iMinMultOfRP,iMaxMultOfRP]
88 // 3.) for a fixed multiplicity use Gaussian with zero spread or use uniform with iMinMult=iMaxMult
89 Bool_t bMultDistrOfRPsIsGauss = kTRUE; 
90                     
91 Int_t iMultiplicityOfRP = 500;        // mean multiplicity of RPs (if sampled from Gaussian)
92 Double_t dMultiplicitySpreadOfRP = 0; // multiplicity spread of RPs (if sampled from Gaussian)
93 Int_t iMinMultOfRP = 50;             // minimal multiplicity of RPs (if sampled uniformly)
94 Int_t iMaxMultOfRP = 500;             // maximal multiplicity of RPs (if sampled uniformly)
95
96
97                     
98 //===DETECTOR ACCEPTANCE===============================================================
99
100 // 1.) if kTRUE = detectors has uniform azimuthal acceptance
101 // 2.) if kFALSE = you will simulate detector with non-uniform acceptance in one or two sectors. 
102 //                 For each of two sectors you specify phi_min, phi_max and probability p. Then all particles 
103 //                 going in direction phi_min < phi < phi_max will be taken with probability p. If p = 0, that
104 //                 sector is blocked. Set bellow phimin1, phimax1, p1 for the first sector and phimin2, phimax2, p2 
105 //                 for the second sector. If you set phimin2 = phimax2 = p2 = 0, only first non-uniform sector is 
106 //                 simulated.
107 Bool_t uniformAcceptance = kTRUE;
108                                                                          
109 // settings for non-uniform acceptance:
110 // Remark: set the angles in degrees from interval [0,360] and probability from interval [0,1]
111
112 // 1st non-uniform sector:
113 Double_t phimin1 = 60;  // first non-uniform sector starts at this azimuth
114 Double_t phimax1 = 120; // first non-uniform sector ends at this azimuth
115 Double_t p1 = 0.6;     // e.g. if p1 = 0 all particles emitted in phimin1 < phi < phimax1 are blocked
116                         // e.g. if p1 = 0.5 half of the particles emitted in phimin1 < phi < phimax1 are blocked
117
118 // 2nd non-uniform sector (Remark: if you do NOT want to simulate this sector, set phimin2 = phimax2 = p2 = 0):                 
119 Double_t phimin2 = 200.0; // second non-uniform sector starts at this azimuth (make sure phimin2 > phimax1 !!!!)
120 Double_t phimax2 = 210.0; // second non-uniform sector ends at this azimuth
121 Double_t p2 = 0.01;
122
123 //=====================================================================================================
124
125 enum anaModes {mLocal,mLocalSource,mLocalPAR};
126 // mLocal: Analyze data on your computer using aliroot
127 // mLocalPAR: Analyze data on your computer using root + PAR files
128 // mLocalSource: Analyze data on your computer using root + source files
129                                           
130 int runFlowAnalysisOnTheFly( Int_t nEvts=44, Int_t mode=mLocal)
131 {
132  TStopwatch timer;
133  timer.Start();
134  
135  if (LYZ1SUM && LYZ2SUM)  {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1.  "<<endl; exit(); }
136  if (LYZ1PROD && LYZ2PROD)  {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1.  "<<endl; exit(); }
137  if (LYZ2SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
138  if (LYZ1SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
139
140  if(!uniformAcceptance && phimin1 > phimax1)
141  {
142   cout<<"WARNING: you must have phimin1 < phimax1 !!!!"<<endl;
143   break;
144  }
145
146  if (!uniformAcceptance && !((phimin2 == 0.) && (phimax2 == 0.) && (p2 == 0.)) && (phimin2 < phimax1 || phimin2 > phimax2))
147  {
148   cout<<"WARNING: you must have phimin2 > phimax1 and phimin2 < phimax2 !!!!"<<endl;
149   break;
150  }
151  
152  if((phimin1 < 0 || phimin1 > 360) || (phimax1 < 0 || phimax1 > 360) || 
153     (phimin2 < 0 || phimin2 > 360) || (phimax2 < 0 || phimax2 > 360) )
154  {
155   cout<<"WARNING: you must take azimuthal angles from interval [0,360] !!!!"<<endl;
156   break;
157  }
158  
159  if((p1 < 0 || p1 > 1) || (p2 < 0 || p2 > 1))
160  {
161   cout<<"WARNING: you must take p1 and p2 from interval [0,1] !!!!"<<endl;
162   break;
163  }
164  
165  cout<<endl;
166  cout<<endl;
167  cout<<"      ---- ARE YOU READY TO FLY ? ----      "<<endl;
168  cout<<endl;
169  
170  cout<<endl;
171  cout<<" ---- BEGIN FLOW ANALYSIS 'ON THE FLY' ---- "<<endl;
172  cout<<endl;
173  cout<<endl;
174  
175  LoadLibraries(mode);
176
177  // Initialize the seed for random generator
178  UInt_t sseed = 0;
179  
180  if(bSameSeed) 
181  {
182   sseed = 44; // the default constant value for seed for random generators
183  } 
184
185  if(!bSameSeed)
186  {
187   TTimeStamp dt;
188   sseed = dt.GetNanoSec()/1000;
189  }
190  
191  cout<<endl;
192  cout<<"Seed for the random generators is "<<sseed<<endl;
193  cout<<endl;
194  
195  //---------------------------------------------------------------------------------------
196  // If the weights are used: 
197  TFile *fileWithWeights = NULL;
198  TList *listWithWeights = NULL;
199   
200  if(usePhiWeights||usePtWeights||useEtaWeights) {
201    fileWithWeights = TFile::Open("weights.root","READ");
202    if(fileWithWeights) {
203      listWithWeights = (TList*)fileWithWeights->Get("weights");
204    }
205    else
206      {cout << " WARNING: the file <weights.root> with weights from the previous run was not found."<<endl;
207       break;
208      }    
209  }
210  
211  //---------------------------------------------------------------------------------------
212  // Initialize the flowevent maker
213  AliFlowEventSimpleMakerOnTheFly* eventMakerOnTheFly = new AliFlowEventSimpleMakerOnTheFly(sseed);
214  eventMakerOnTheFly->Init();
215  //set the range for eta subevents
216  eventMakerOnTheFly->SetSubeventEtaRange(minA,maxA,minB,maxB); 
217   
218  //---------------------------------------------------------------------------------------
219  // Initialize all the flow methods for default analysis:  
220  AliFlowAnalysisWithQCumulants *qc = NULL;
221  AliFlowAnalysisWithCumulants *gfc = NULL;
222  AliFlowAnalysisWithFittingQDistribution *fqd = NULL;
223  AliFlowAnalysisWithLeeYangZeros *lyz1sum  = NULL;
224  AliFlowAnalysisWithLeeYangZeros *lyz1prod = NULL;
225  AliFlowAnalysisWithLeeYangZeros *lyz2sum  = NULL;
226  AliFlowAnalysisWithLeeYangZeros *lyz2prod = NULL;
227  AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
228  AliFlowAnalysisWithScalarProduct *sp = NULL;
229  AliFlowAnalysisWithMCEventPlane *mcep = NULL;   
230
231  // MCEP = monte carlo event plane
232  if (MCEP) {
233    AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
234    // mcep->SetHarmonic(2); // default is v2
235    mcep->Init();
236  }
237
238  // QC = Q-cumulants  
239  if(QC) { 
240    AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
241    if(listWithWeights) qc->SetWeightsList(listWithWeights);
242    if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
243    if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
244    if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
245    // qc->SetHarmonic(2); // default is v2
246    // qc->SetApplyCorrectionForNUA(kTRUE); // default
247    // qc->SetCalculate2DFlow(kFALSE); // default
248    // qc->SetMultiplicityWeight("combinations"); // default
249    // qc->SetMultiplicityWeight("unit");
250    // qc->SetMultiplicityWeight("multiplicity");  
251    qc->Init();  
252  }
253   
254  // GFC = Generating Function Cumulants 
255  if(GFC) {
256    AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
257    if(listWithWeights) gfc->SetWeightsList(listWithWeights);
258    if(usePhiWeights) gfc->SetUsePhiWeights(usePhiWeights);
259    if(usePtWeights) gfc->SetUsePtWeights(usePtWeights);
260    if(useEtaWeights) gfc->SetUseEtaWeights(useEtaWeights);
261    gfc->Init();
262  }
263  
264  // FQD = Fitting q-distribution 
265  if(FQD) {
266    AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
267    if(listWithWeights) fqd->SetWeightsList(listWithWeights);
268    if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);  
269    fqd->Init();
270  }
271  
272  // SP = Scalar Product 
273  if(SP) {
274    AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
275    if(listWithWeights) sp->SetWeightsList(listWithWeights);
276    sp->SetUsePhiWeights(usePhiWeights);  
277    sp->Init();
278  }
279
280  // LYZ1 = Lee-Yang Zeroes first run
281  if(LYZ1SUM) {
282    AliFlowAnalysisWithLeeYangZeros* lyz1sum = new AliFlowAnalysisWithLeeYangZeros();
283    lyz1sum->SetFirstRun(kTRUE);
284    lyz1sum->SetUseSum(kTRUE);
285    lyz1sum->Init();
286  }
287  if(LYZ1PROD) {
288    AliFlowAnalysisWithLeeYangZeros* lyz1prod = new AliFlowAnalysisWithLeeYangZeros();
289    lyz1prod->SetFirstRun(kTRUE);
290    lyz1prod->SetUseSum(kFALSE);
291    lyz1prod->Init();
292  }
293  // LYZ2 = Lee-Yang Zeroes second run
294  if(LYZ2SUM) {
295    AliFlowAnalysisWithLeeYangZeros* lyz2sum = new AliFlowAnalysisWithLeeYangZeros();
296    // read the input file from the first run 
297    TString inputFileNameLYZ2SUM = "outputLYZ1SUManalysis.root" ;
298    TFile* inputFileLYZ2SUM = new TFile(inputFileNameLYZ2SUM.Data(),"READ");
299    if(!inputFileLYZ2SUM || inputFileLYZ2SUM->IsZombie()) { 
300      cerr << " ERROR: To run LYZ2SUM you need the output file from LYZ1SUM. This file is not there! Please run LYZ1SUM first." << endl ;
301      break; 
302    }
303    else { 
304      TList* inputListLYZ2SUM = (TList*)inputFileLYZ2SUM->Get("cobjLYZ1SUM");  
305      if (!inputListLYZ2SUM) {cout<<"Input list LYZ2SUM is NULL pointer!"<<endl; break;}
306      else {
307        cout<<"LYZ2SUM input file/list read..."<<endl;
308        lyz2sum->SetFirstRunList(inputListLYZ2SUM);
309        lyz2sum->SetFirstRun(kFALSE);
310        lyz2sum->SetUseSum(kTRUE);
311        lyz2sum->Init();
312      }
313    }
314  }
315  if(LYZ2PROD) {
316    AliFlowAnalysisWithLeeYangZeros* lyz2prod = new AliFlowAnalysisWithLeeYangZeros();
317    // read the input file from the first run 
318    TString inputFileNameLYZ2PROD = "outputLYZ1PRODanalysis.root" ;
319    TFile* inputFileLYZ2PROD = new TFile(inputFileNameLYZ2PROD.Data(),"READ");
320    if(!inputFileLYZ2PROD || inputFileLYZ2PROD->IsZombie()) { 
321      cerr << " ERROR: To run LYZ2PROD you need the output file from LYZ1PROD. This file is not there! Please run LYZ1PROD first." << endl ;
322      break; 
323    }
324    else { 
325      TList* inputListLYZ2PROD = (TList*)inputFileLYZ2PROD->Get("cobjLYZ1PROD");  
326      if (!inputListLYZ2PROD) {cout<<"Input list LYZ2PROD is NULL pointer!"<<endl; break;}
327      else {
328        cout<<"LYZ2PROD input file/list read..."<<endl;
329        lyz2prod->SetFirstRunList(inputListLYZ2PROD);
330        lyz2prod->SetFirstRun(kFALSE);
331        lyz2prod->SetUseSum(kFALSE);
332        lyz2prod->Init();
333      }
334    }
335  }
336
337  // LYZEP = Lee-Yang Zeroes event plane
338  if(LYZEP) {
339    AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
340    AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
341    // read the input file from the second lyz run 
342    TString inputFileNameLYZEP = "outputLYZ2SUManalysis.root" ;
343    TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
344    if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) { 
345      cerr << " ERROR: To run LYZEP you need the output file from LYZ2SUM. This file is not there! Please run LYZ2SUM first." << endl ; 
346      break;
347    }
348    else { 
349      TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2SUM");  
350      if (!inputListLYZEP) {cout<<"Input list LYZEP is NULL pointer!"<<endl; break;}
351      else {
352        cout<<"LYZEP input file/list read..."<<endl;
353        ep   ->SetSecondRunList(inputListLYZEP);
354        lyzep->SetSecondRunList(inputListLYZEP);
355        ep   ->Init();
356        lyzep->Init();
357      }
358    }
359  }
360  //---------------------------------------------------------------------------------------
361  
362  //---------------------------------------------------------------------------------------
363  // Initialize all the flow methods for additional analysis with different settings/aims
364  // Label each setting/aim with different label !!!!
365  
366  // GFC:    
367  TString gfcDefaultName = "outputGFCanalysis"; 
368  // 1.) GFC analysis for elliptic flow with r0 = 1.5: 
369  AliFlowAnalysisWithCumulants *gfc_1;
370  TString gfcAnalysisLabels_1 = "_r0_1.5"; // all histograms and output file name will have this label
371  TString gfcOutputFileName_1;
372  gfcOutputFileName_1 = gfcDefaultName.Data();
373  gfcOutputFileName_1 += gfcAnalysisLabels_1.Data();
374  gfcOutputFileName_1 += ".root";
375  if(GFC_Additional_Analysis)
376  {
377   gfc_1 = new AliFlowAnalysisWithCumulants();
378   //gfc_1->SetAnalysisLabel(gfcAnalysisLabels_1.Data());
379   if(listWithWeights) gfc_1->SetWeightsList(listWithWeights);
380   if(usePhiWeights) gfc_1->SetUsePhiWeights(usePhiWeights);
381   if(usePtWeights) gfc_1->SetUsePtWeights(usePtWeights);
382   if(useEtaWeights) gfc_1->SetUseEtaWeights(useEtaWeights);
383   gfc_1->Init();
384  }
385  
386  // QC:
387  TString qcDefaultName = "outputQCanalysis"; 
388  // 1.) QC analysis for directed flow: 
389  AliFlowAnalysisWithQCumulants *qc_1;
390  TString qcAnalysisLabels_1 = "_v1"; // all histograms and output file name will have this label
391  TString qcOutputFileName_1;
392  qcOutputFileName_1 = qcDefaultName.Data();
393  qcOutputFileName_1 += qcAnalysisLabels_1.Data();
394  qcOutputFileName_1 += ".root";
395  if(QC_Additional_Analysis)
396  {
397   qc_1 = new AliFlowAnalysisWithQCumulants();
398   //qc_1->SetAnalysisLabel(qcAnalysisLabels_1->Data());
399   if(listWithWeights) qc_1->SetWeightsList(listWithWeights);
400   if(usePhiWeights) qc_1->SetUsePhiWeights(usePhiWeights);
401   qc_1->Init();
402  }
403  
404  // FQD:
405  TString fqdDefaultName = "outputFQDanalysis";
406  // 1.) FQD fitting with fixed sigma:
407  AliFlowAnalysisWithFittingQDistribution *fqd_1;
408  TString fqdAnalysisLabels_1 = "_fixedSigma"; // all histograms and output file name will have this label
409  TString fqdOutputFileName_1;
410  fqdOutputFileName_1 = fqdDefaultName.Data();
411  fqdOutputFileName_1 += fqdAnalysisLabels_1.Data();
412  fqdOutputFileName_1 += ".root";
413  if(FQD_Additional_Analysis)
414  { 
415   fqd_1 = new AliFlowAnalysisWithFittingQDistribution();
416   //fqd_1->SetAnalysisLabel(fqdAnalysisLabels_1->Data());
417   if(listWithWeights) fqd_1->SetWeightsList(listWithWeights);
418   if(usePhiWeights) fqd_1->SetUsePhiWeights(usePhiWeights);
419   fqd_1->Init();
420  }
421  //---------------------------------------------------------------------------------------
422   
423  // set the global event parameters: 
424  eventMakerOnTheFly->SetNoOfLoops(iLoops);
425  
426  if(bMultDistrOfRPsIsGauss)
427  {
428   eventMakerOnTheFly->SetMultDistrOfRPsIsGauss(bMultDistrOfRPsIsGauss);
429   eventMakerOnTheFly->SetMultiplicityOfRP(iMultiplicityOfRP);
430   eventMakerOnTheFly->SetMultiplicitySpreadOfRP(dMultiplicitySpreadOfRP);
431  } else
432    {
433     eventMakerOnTheFly->SetMultDistrOfRPsIsGauss(bMultDistrOfRPsIsGauss);
434     eventMakerOnTheFly->SetMinMultOfRP(iMinMultOfRP);
435     eventMakerOnTheFly->SetMaxMultOfRP(iMaxMultOfRP); 
436    }
437   
438  eventMakerOnTheFly->SetTemperatureOfRP(dTemperatureOfRP);
439
440  eventMakerOnTheFly->SetV1RP(dV1RP);
441  eventMakerOnTheFly->SetV1SpreadRP(dV1SpreadRP);  
442  eventMakerOnTheFly->SetV4RP(dV4RP);
443  eventMakerOnTheFly->SetV4SpreadRP(dV4SpreadRP);  
444  
445  // constant harmonic V2:
446  if(bConstantHarmonics)
447  { 
448   eventMakerOnTheFly->SetUseConstantHarmonics(bConstantHarmonics);
449   if(bV2DistrOfRPsIsGauss)
450   {
451    eventMakerOnTheFly->SetV2DistrOfRPsIsGauss(bV2DistrOfRPsIsGauss);
452    eventMakerOnTheFly->SetV2RP(dV2RP);
453    eventMakerOnTheFly->SetV2SpreadRP(dV2SpreadRP);  
454   } else
455     {
456      eventMakerOnTheFly->SetV2DistrOfRPsIsGauss(bV2DistrOfRPsIsGauss);
457      eventMakerOnTheFly->SetMinV2RP(dMinV2RP);
458      eventMakerOnTheFly->SetMaxV2RP(dMaxV2RP);  
459     }
460  }
461  
462  // (pt,eta) dependent harmonic V2:
463  if(!bConstantHarmonics)
464  {
465   eventMakerOnTheFly->SetUseConstantHarmonics(bConstantHarmonics);
466   eventMakerOnTheFly->SetV2RPMax(dV2RPMax);
467   eventMakerOnTheFly->SetPtCutOff(dPtCutOff);  
468  }
469  
470  // non-uniform acceptance:
471  if(!uniformAcceptance)
472  {
473   eventMakerOnTheFly->SetFirstSectorPhiMin(phimin1);
474   eventMakerOnTheFly->SetFirstSectorPhiMax(phimax1);
475   eventMakerOnTheFly->SetFirstSectorProbability(p1);
476   eventMakerOnTheFly->SetSecondSectorPhiMin(phimin2);
477   eventMakerOnTheFly->SetSecondSectorPhiMax(phimax2);
478   eventMakerOnTheFly->SetSecondSectorProbability(p2);
479  }
480        
481  //---------------------------------------------------------------------------------------  
482  // create and analyze events 'on the fly':
483
484  for(Int_t i=0;i<nEvts;i++) {   
485    // creating the event with above settings:
486    AliFlowEventSimple *event = eventMakerOnTheFly->CreateEventOnTheFly(); 
487    
488    // analyzing the created event 'on the fly':
489    // do flow analysis for various methods:
490    if(MCEP)    mcep->Make(event);
491    if(QC)      qc->Make(event);
492    if(GFC)     gfc->Make(event);
493    if(FQD)     fqd->Make(event);
494    if(LYZ1SUM) lyz1sum->Make(event);
495    if(LYZ1PROD)lyz1prod->Make(event);
496    if(LYZ2SUM) lyz2sum->Make(event);
497    if(LYZ2PROD)lyz2prod->Make(event);
498    if(LYZEP)   lyzep->Make(event,ep);
499    if(SP)      sp->Make(event);
500    
501    if(GFC_Additional_Analysis)
502    {
503     // r0 = 1.5:
504     gfc_1->Make(event);
505    }
506    if(QC_Additional_Analysis)
507    {
508     // v1:
509     qc_1->Make(event);
510    }
511    if(FQD_Additional_Analysis)
512    {
513     // fixed sigma:
514     fqd_1->Make(event);
515    }
516    
517    delete event;
518  } // end of for(Int_t i=0;i<nEvts;i++)
519  //---------------------------------------------------------------------------------------  
520
521
522
523  //---------------------------------------------------------------------------------------  
524  // create a new file which will hold the final results of all methods:
525  TString outputFileName = "AnalysisResults.root";  
526  TFile *outputFile = new TFile(outputFileName.Data(),"RECREATE");
527  // create a new file for each method wich will hold list with final results:
528  const Int_t nMethods = 10;
529  TString method[nMethods] = {"MCEP","SP","GFC","QC","FQD","LYZ1SUM","LYZ1PROD","LYZ2SUM","LYZ2PROD","LYZEP"};
530  TDirectoryFile *dirFileFinal[nMethods] = {NULL};
531  TString fileName[nMethods]; 
532  for(Int_t i=0;i<nMethods;i++)
533  {
534   // form a file name for each method:
535   fileName[i]+="output";
536   fileName[i]+=method[i].Data();
537   fileName[i]+="analysis";
538   dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data());
539  } 
540  
541  // calculating and storing the final results of default flow analysis:
542  if(MCEP)    {mcep->Finish();    mcep->WriteHistograms(dirFileFinal[0]);}
543  if(SP)      {sp->Finish();      sp->WriteHistograms(dirFileFinal[1]);}
544  if(GFC)     {gfc->Finish();     gfc->WriteHistograms(dirFileFinal[2]);}
545  if(QC)      {qc->Finish();      qc->WriteHistograms(dirFileFinal[3]);}
546  if(FQD)     {fqd->Finish();     fqd->WriteHistograms(dirFileFinal[4]);}
547  if(LYZ1SUM) {lyz1sum->Finish(); lyz1sum->WriteHistograms(dirFileFinal[5]);}
548  if(LYZ1PROD){lyz1prod->Finish();lyz1prod->WriteHistograms(dirFileFinal[6]);}
549  if(LYZ2SUM) {lyz2sum->Finish(); lyz2sum->WriteHistograms(dirFileFinal[7]);}
550  if(LYZ2PROD){lyz2prod->Finish();lyz2prod->WriteHistograms(dirFileFinal[8]);}
551  if(LYZEP)   {lyzep->Finish();   lyzep->WriteHistograms(dirFileFinal[9]);}
552  //---------------------------------------------------------------------------------------  
553  
554  outputFile->Close();
555  delete outputFile;
556  
557  cout<<endl;
558  cout<<endl;
559  cout<<" ---- LANDED SUCCESSFULLY ---- "<<endl;
560  cout<<endl; 
561  
562  timer.Stop();
563  cout << endl;
564  timer.Print();
565 }
566
567 void SetupPar(char* pararchivename)
568 {
569   //Load par files, create analysis libraries
570   //For testing, if par file already decompressed and modified
571   //classes then do not decompress.
572  
573   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ; 
574   TString parpar(Form("%s.par", pararchivename)) ; 
575   if ( gSystem->AccessPathName(parpar.Data()) ) {
576     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
577     TString processline(Form(".! make %s", parpar.Data())) ; 
578     gROOT->ProcessLine(processline.Data()) ;
579     gSystem->ChangeDirectory(cdir) ; 
580     processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
581     gROOT->ProcessLine(processline.Data()) ;
582   } 
583   if ( gSystem->AccessPathName(pararchivename) ) {  
584     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
585     gROOT->ProcessLine(processline.Data());
586   }
587   
588   TString ocwd = gSystem->WorkingDirectory();
589   gSystem->ChangeDirectory(pararchivename);
590   
591   // check for BUILD.sh and execute
592   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
593     printf("*******************************\n");
594     printf("*** Building PAR archive    ***\n");
595     cout<<pararchivename<<endl;
596     printf("*******************************\n");
597     
598     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
599       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
600       return -1;
601     }
602   }
603   // check for SETUP.C and execute
604   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
605     printf("*******************************\n");
606     printf("*** Setup PAR archive       ***\n");
607     cout<<pararchivename<<endl;
608     printf("*******************************\n");
609     gROOT->Macro("PROOF-INF/SETUP.C");
610   }
611   
612   gSystem->ChangeDirectory(ocwd.Data());
613   printf("Current dir: %s\n", ocwd.Data());
614 }
615
616 void LoadLibraries(const anaModes mode) {
617   
618   //--------------------------------------
619   // Load the needed libraries most of them already loaded by aliroot
620   //--------------------------------------
621   //gSystem->Load("libTree");
622   gSystem->Load("libGeom");
623   gSystem->Load("libVMC");
624   gSystem->Load("libXMLIO");
625   gSystem->Load("libPhysics");
626   
627   //----------------------------------------------------------
628   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<< 
629   //----------------------------------------------------------
630   if (mode==mLocal) {
631     //--------------------------------------------------------
632     // If you want to use already compiled libraries 
633     // in the aliroot distribution
634     //--------------------------------------------------------
635     gSystem->Load("libSTEERBase");
636     gSystem->Load("libESD");
637     gSystem->Load("libAOD");
638     gSystem->Load("libANALYSIS");
639     gSystem->Load("libANALYSISalice");
640     gSystem->Load("libCORRFW");
641     cerr<<"libCORRFW loaded..."<<endl;
642     gSystem->Load("libPWG2flowCommon");
643     cerr<<"libPWG2flowCommon loaded..."<<endl;
644     gSystem->Load("libPWG2flowTasks");
645     cerr<<"libPWG2flowTasks loaded..."<<endl;
646   }
647   
648   else if (mode == mLocalPAR) {
649     //--------------------------------------------------------
650     //If you want to use root and par files from aliroot
651     //--------------------------------------------------------  
652      //If you want to use root and par files from aliroot
653     //--------------------------------------------------------  
654     SetupPar("STEERBase");
655     SetupPar("ESD");
656     SetupPar("AOD");
657     SetupPar("ANALYSIS");
658     SetupPar("ANALYSISalice");
659     SetupPar("PWG2AOD");
660     SetupPar("CORRFW");
661     SetupPar("PWG2flowCommon");
662     cerr<<"PWG2flowCommon.par loaded..."<<endl;
663     SetupPar("PWG2flowTasks");
664     cerr<<"PWG2flowTasks.par loaded..."<<endl;
665   }
666   
667   //---------------------------------------------------------
668   // <<<<<<<<<< Source mode >>>>>>>>>>>>
669   //---------------------------------------------------------
670   else if (mode==mLocalSource) {
671  
672     // In root inline compile
673
674    
675     // Constants  
676     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
677     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
678     gROOT->LoadMacro("AliFlowCommon/AliFlowCumuConstants.cxx+");
679     
680     // Flow event
681     gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+"); 
682     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");    
683     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
684     
685     // Cuts
686     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");    
687     
688     // Output histosgrams
689     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
690     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
691     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
692     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
693     
694     // Functions needed for various methods
695     gROOT->LoadMacro("AliFlowCommon/AliCumulantsFunctions.cxx+");
696     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZEventPlane.cxx+");
697     
698     // Flow Analysis code for various methods
699     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMCEventPlane.cxx+"); 
700     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx+");
701     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.cxx+");
702     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.cxx+");
703     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithCumulants.cxx+");
704     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithQCumulants.cxx+"); 
705     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithFittingQDistribution.cxx+");
706     
707     // Class to fill the FlowEvent on the fly (generate Monte Carlo events)
708     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.cxx+");   
709     
710     cout << "finished loading macros!" << endl;  
711     
712   }  
713   
714 }
715
716