hooks for PMD flow analysis
[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 Bool_t MH       = kFALSE; // mixed harmonics 
20 Bool_t NL       = kFALSE; // nested loops
21 Bool_t MCEP_AH  = kFALSE; // MCEP in another harmonic 
22 //--------------------------------------------------------------------------------------
23
24 // Weights 
25 // use weights for Q-vector:
26 Bool_t usePhiWeights = kFALSE; // phi weights (correction for non-uniform azimuthal acceptance)
27 Bool_t usePtWeights  = kFALSE; // pt weights 
28 Bool_t useEtaWeights = kFALSE; // eta weights
29
30 // Define the range for eta subevents
31 Double_t minA = -0.9;
32 Double_t maxA = -0.01;
33 Double_t minB = 0.01;
34 Double_t maxB = 0.9;
35
36 // Define simple cuts for RP selection:
37 Double_t ptMinRP = 0.0; // in GeV
38 Double_t ptMaxRP = 10.0; // in GeV
39 Double_t etaMinRP  = -1.;
40 Double_t etaMaxRP  = 1.;
41 Double_t phiMinRP  = 0.0; // in degrees
42 Double_t phiMaxRP  = 360.0; // in degrees
43 Int_t pidRP = -1; // to be improved (supported eventually)
44
45 // Define simple cuts for POI selection:
46 Double_t ptMinPOI = 0.0; // in GeV
47 Double_t ptMaxPOI = 10.0; // in GeV
48 Double_t etaMinPOI  = -1.;
49 Double_t etaMaxPOI  = 1.;
50 Double_t phiMinPOI  = 0.0; // in degrees
51 Double_t phiMaxPOI  = 360.0; // in degrees
52 Int_t pidPOI = -1; // to be improved (supported eventually)
53
54 // Parameters for the simulation of events 'on the fly': 
55 //===SEED========================================================
56 // use always the same seed for random generators. 
57 // usage of same seed (kTRUE) is relevant in two cases:
58 // 1.) If you want to use LYZ method to calcualte differential flow;
59 // 2.) If you want to use phi weights for GFC, QC and FQD
60 Bool_t bSameSeed = kFALSE;  
61
62 //===NONFLOW=============================================
63 // number of times to use each track (to simulate nonflow)
64 Int_t iLoops = 1; 
65 Double_t phiRange = 0.0; // If the original track with azimuthal angle phi is splitted, 
66                          // the azimuthal angle of splitted track is sampled uniformly 
67                          // from (phi-phiRange, phi+phiRange). If phiRange = 0.0, the
68                          // azimuthal angle of splitted track is the same as azimuthal  
69                          // angle of original track. (Remark: phiRange's unit is degree.)
70 Double_t ptRange = 0.0; // If the original track with momentum pt is splitted, 
71                         // the momentum of splitted track is sampled uniformly 
72                         // from (pt-ptRange, pt+ptRange). If ptRange = 0.0, the
73                         // momentum of splitted track is the same as momentum 
74                         // of original track. (Remark: ptRange's unit is GeV.)
75 Double_t etaRange = 0.0; // If the original track with pseudorapidity eta is splitted, 
76                          // the pseudorapidity of splitted track is sampled uniformly 
77                          // from (eta-etaRange, eta+etaRange). If etaRange = 0.0, the
78                          // pseudorapidity of splitted track will be the same as the
79                          // pseudorapidity of original track.
80 // in addition one can simulate nonflow only in a certain detector's sector 
81 // ranging from nonflowSectorMin to nonflowSectorMax. Outside this sector the
82 // tracks will not be splitted. Use the following two settings only if iLoops>1:                         
83 Double_t nonflowSectorMin = 0.0; // detector's sector in which tracks will be splitted starts at this angle (in degrees)                         
84 Double_t nonflowSectorMax = 360.0; // detector's sector in which tracks will be splitted ends at this angle (in degrees)    
85 // to study nonflow fluctuations event-by-event one can create the jet-like structure in the following way:
86 Bool_t bCreateJets = kFALSE; // set kTRUE if you want to study nonflow fluctuations via jet creation
87 Double_t jetProbability = 0.2; // probability than event will have a jet
88 Double_t jetTracksFraction = 0.1; // percentage of total tracks in an event belonging to the jet
89 Double_t jetCone = 20; // angular size of jet cone (in degrees)                   
90
91 //===GLAUBER MODEL================================================================
92 Bool_t bUseGlauberModel = kFALSE; // 1.) if kTRUE = multiplicity and constant flow harmonics are 
93                                   //                determined event-by-event from Glauber model
94                                   //                (pt and/or eta dependence of flow harmonics not supported;
95                                   //                settings for multiplicity and flow harmonics bellow are irrelevant) 
96                                   // 2.) if kFALSE = multiplicity and flow harmonics are determined
97                                   //                 independently event-by-event with bellow settings
98                                   //                 (pt and/or eta dependence of flow harmonics supported)
99
100 //===FLOW HARMONICS===============================================================
101 // harmonics V1, V2, V4... are constants or functions of pt and eta:         
102 Bool_t bPtDependentHarmonicV1 = kFALSE; 
103 Bool_t bEtaDependentHarmonicV1 = kFALSE;            
104 Bool_t bPtDependentHarmonicV2 = kFALSE; 
105 Bool_t bEtaDependentHarmonicV2 = kFALSE; 
106 Bool_t bPtDependentHarmonicV4 = kFALSE; 
107 Bool_t bEtaDependentHarmonicV4 = kFALSE; 
108 // 1.) if you use constant harmonics (bPtDependentHarmonic* = kFALSE, bEtaDependentHarmonic* = kFALSE)
109 //     you can additionally select if V2 will be sampled e-b-e from Gaussian or from uniform distribution
110 //     (constant harmonics V1 and V4 will be always sampled from Gaussian):  
111 Bool_t bConstantV2IsSampledFromGauss = kTRUE;
112  //  1a.) if kTRUE = elliptic flow of RPs is sampled e-b-e from Gaussian distribution with
113  //                  mean = dV2RP and spread = dV2SpreadRP (analogously for V1 and V4).
114  //  1b.) if kFALSE = elliptic flow of RPs is sampled e-b-e uniformly from 
115  //                   interval [dMinV2RP,dMaxV2RP] (not supported for V1 and V4).
116  //  1c.) for a fixed elliptic flow e-b-e use Gaussian with zero spread or use uniform with dMinV2RP=dMaxV2RP 
117  // V2:
118  Double_t dV2RP = 0.05;       // elliptic flow of RPs (if sampled from Gaussian) 
119  Double_t dV2SpreadRP = 0.0;  // elliptic flow spread of RPs (if sampled from Gaussian)
120  Double_t dMinV2RP = 0.04;    // minimal elliptic flow of RPs (if sampled uniformly)
121  Double_t dMaxV2RP = 0.06;    // maximal elliptic flow of RPs (if sampled uniformly)
122  // V1:
123  Double_t dV1RP = 0.0; // directed flow of RPs 
124  Double_t dV1SpreadRP = 0.0; // directed flow spread of RPs
125  // V3:
126  Double_t dV3RP = 0.0; // triangular flow of RPs
127  Double_t dV3SpreadRP = 0.0; // triangular flow spread of RP
128  // V4:
129  Double_t dV4RP = 0.0; // harmonic V4 of RPs (to be improved: name needed) 
130  Double_t dV4SpreadRP = 0.0; // harmonic V4's spread of RPs (to be improved: name needed)
131 // 2.) if you use (pt,eta) dependent harmonic V1 (bPtDependentHarmonicV1 = kTRUE or bEtaDependentHarmonicV1 = kTRUE): 
132 //  2a.) V1(pt) is linear up to pt = dV1PtCutOff and for pt > dV1PtCutOff it is constant, V1(pt) = dV1vsPtEtaMax:
133 //  2b.) V1(eta) is determined from formula V1(eta) = -eta
134  Double_t dV1vsPtEtaMax = 0.10; // maximum value of V1(pt,eta) (for pt >= dV1PtCutOff and at eta = 0)
135  Double_t dV1PtCutOff = 2.0; // up to this pt value V1(pt) is growing linearly versus pt
136 // 3.) if you use (pt,eta) dependent harmonic V2 (bPtDependentHarmonicV2 = kTRUE or bEtaDependentHarmonicV2 = kTRUE): 
137 //  2a.) V2(pt) is linear up to pt = dV2PtCutOff and for pt > dV2PtCutOff it is constant, V2(pt) = dV2vsPtEtaMax:
138 //  2b.) V2(eta) is Gaussian centered at midrapidity (eta=0) with V2(0) = dV2vsPtEtaMax and spread = dV2vsEtaSpread:
139  Double_t dV2vsPtEtaMax = 0.20; // maximum value of V2(pt,eta) (for pt >= dV2PtCutOff and at eta = 0)
140  Double_t dV2PtCutOff = 2.0; // up to this pt value V2(pt) is growing linearly versus pt
141  Double_t dV2vsEtaSpread = 0.75; // V2(eta) is Gaussian centered at midrapidity (eta=0) with spread = dV2vsEtaSpread
142 // 4.) if you use (pt,eta) dependent harmonic V4 (bPtDependentHarmonicV4 = kTRUE or bEtaDependentHarmonicV4 = kTRUE):
143 //     V4(pt,eta) is determined as V4(pt,eta) = pow(V2(pt,eta),2.) 
144
145 //===MULTIPLICITY===============================================================
146 // 1.) if kTRUE  = multiplicitiy of RPs is sampled e-b-e from Gaussian distribution with
147 //                 mean = iMultiplicityOfRP and spread = dMultiplicitySpreadOfRP
148 // 2.) if kFALSE = multiplicitiy of RPs is sampled e-b-e uniformly from 
149 //                 interval [iMinMultOfRP,iMaxMultOfRP]
150 // 3.) for a fixed multiplicity use Gaussian with zero spread or use uniform with iMinMult=iMaxMult
151 Bool_t bMultDistrOfRPsIsGauss = kTRUE; 
152                     
153 Int_t iMultiplicityOfRP = 500;        // mean multiplicity of RPs (if sampled from Gaussian)
154 Double_t dMultiplicitySpreadOfRP = 0; // multiplicity spread of RPs (if sampled from Gaussian)
155 Int_t iMinMultOfRP = 400;              // minimal multiplicity of RPs(if sampled uniformly)
156 Int_t iMaxMultOfRP = 600;             // maximal multiplicity of RPs (if sampled uniformly)
157                     
158 //===DETECTOR ACCEPTANCE===============================================================
159
160 // 1.) if kTRUE = detectors has uniform azimuthal acceptance
161 // 2.) if kFALSE = you will simulate detector with non-uniform acceptance in one or two sectors. 
162 //                 For each of two sectors you specify phi_min, phi_max and probability p. Then all particles 
163 //                 going in direction phi_min < phi < phi_max will be taken with probability p. If p = 0, that
164 //                 sector is blocked. Set bellow phimin1, phimax1, p1 for the first sector and phimin2, phimax2, p2 
165 //                 for the second sector. If you set phimin2 = phimax2 = p2 = 0, only first non-uniform sector is 
166 //                 simulated.
167 Bool_t uniformAcceptance = kTRUE;
168                                                                          
169 // settings for non-uniform acceptance:
170 // Remark: set the angles in degrees from interval [0,360] and probability from interval [0,1]
171
172 // 1st non-uniform sector:
173 Double_t phimin1 = 60;  // first non-uniform sector starts at this azimuth
174 Double_t phimax1 = 120; // first non-uniform sector ends at this azimuth
175 Double_t p1 = 0.5;     // e.g. if p1 = 0 all particles emitted in phimin1 < phi < phimax1 are blocked
176                         // e.g. if p1 = 0.5 half of the particles emitted in phimin1 < phi < phimax1 are blocked
177
178 // 2nd non-uniform sector (Remark: if you do NOT want to simulate this sector, set phimin2 = phimax2 = p2 = 0):                 
179 Double_t phimin2 = 0.0; // second non-uniform sector starts at this azimuth (make sure phimin2 > phimax1 !!!!)
180 Double_t phimax2 = 0.0; // second non-uniform sector ends at this azimuth
181 Double_t p2 = 0.0;
182
183
184 //===MODIFYING Pt SPECTRA===============================================================
185 Double_t dTemperatureOfRP = 0.44; // 'temperature' of RPs in GeV/c (increase this parameter to get more high pt RPs) 
186
187 enum anaModes {mLocal,mLocalSource,mLocalPAR};
188 // mLocal: Analyze data on your computer using aliroot
189 // mLocalPAR: Analyze data on your computer using root + PAR files
190 // mLocalSource: Analyze data on your computer using root + source files
191                                           
192 int runFlowAnalysisOnTheFly(Int_t nEvts=1000, Int_t mode=mLocal)
193 {
194  CheckUserSettings();
195  TStopwatch timer;
196  timer.Start();
197  
198  cout<<endl;
199  cout<<endl;
200  cout<<"      ---- ARE YOU READY TO FLY ? ----      "<<endl;
201  cout<<endl;
202  
203  cout<<endl;
204  cout<<" ---- BEGIN FLOW ANALYSIS 'ON THE FLY' ---- "<<endl;
205  cout<<endl;
206  cout<<endl;
207  
208  LoadLibraries(mode);
209
210  // Initialize the seed for random generator
211  UInt_t sseed = 0;
212  
213  if(bSameSeed) 
214  {
215   sseed = 44; // the default constant value for seed for random generators
216  } 
217
218  if(!bSameSeed)
219  {
220   TTimeStamp dt;
221   sseed = dt.GetNanoSec()/1000;
222  }
223  
224  cout<<endl;
225  cout<<"Seed for the random generators is "<<sseed<<endl;
226  cout<<endl;
227  
228  //---------------------------------------------------------------------------------------
229  // If the weights are used: 
230  TFile *fileWithWeights = NULL;
231  TList *listWithWeights = NULL;
232   
233  if(usePhiWeights||usePtWeights||useEtaWeights) {
234    fileWithWeights = TFile::Open("weights.root","READ");
235    if(fileWithWeights) {
236      listWithWeights = (TList*)fileWithWeights->Get("weights");
237    }
238    else
239      {cout << " WARNING: the file <weights.root> with weights from the previous run was not found."<<endl;
240       break;
241      }    
242  }
243  
244  //---------------------------------------------------------------------------------------
245  // Initialize the flowevent maker
246  AliFlowEventSimpleMakerOnTheFly* eventMakerOnTheFly = new AliFlowEventSimpleMakerOnTheFly(sseed);
247  eventMakerOnTheFly->Init();
248  //set the range for eta subevents
249  eventMakerOnTheFly->SetSubeventEtaRange(minA,maxA,minB,maxB); 
250   
251  //---------------------------------------------------------------------------------------
252  // Initialize all the flow methods for default analysis:  
253  AliFlowAnalysisWithQCumulants *qc = NULL;
254  AliFlowAnalysisWithCumulants *gfc = NULL;
255  AliFlowAnalysisWithFittingQDistribution *fqd = NULL;
256  AliFlowAnalysisWithLeeYangZeros *lyz1sum  = NULL;
257  AliFlowAnalysisWithLeeYangZeros *lyz1prod = NULL;
258  AliFlowAnalysisWithLeeYangZeros *lyz2sum  = NULL;
259  AliFlowAnalysisWithLeeYangZeros *lyz2prod = NULL;
260  AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
261  AliFlowAnalysisWithScalarProduct *sp = NULL;
262  AliFlowAnalysisWithMixedHarmonics *mh = NULL;
263  AliFlowAnalysisWithNestedLoops *nl = NULL;
264  AliFlowAnalysisWithMCEventPlane *mcep = NULL;   
265  AliFlowAnalysisWithMCEventPlane *mcep_ah = NULL;   
266
267  // MCEP = monte carlo event plane
268  if (MCEP) {
269    AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
270    //mcep->SetHarmonic(2); // default is v2
271    //mcep->SetFlowOfResonances(kTRUE);
272    mcep->Init();
273  }
274
275  // MCEP = monte carlo event plane in another harmonic: 
276  if(MCEP_AH)
277  {
278   AliFlowAnalysisWithMCEventPlane *mcep_ah = new AliFlowAnalysisWithMCEventPlane();
279   mcep_ah->SetHarmonic(1);
280   mcep_ah->Init();
281  }
282  
283  // Mixed harmonics:
284  if(MH) 
285  {
286   AliFlowAnalysisWithMixedHarmonics *mh = new AliFlowAnalysisWithMixedHarmonics();
287   mh->SetHarmonic(1); // integer n in expression cos[n(2phi1-phi2-phi3)] = v2n*vn^2
288   mh->SetMinMultiplicity(100); 
289   mh->SetNoOfMultipicityBins(5);  
290   mh->SetMultipicityBinWidth(200);   
291   mh->Init(); 
292  }
293  
294  // NL = nested loops:
295  if(NL) {
296    AliFlowAnalysisWithNestedLoops *nl = new AliFlowAnalysisWithNestedLoops();
297    nl->Init();
298  }
299
300  // QC = Q-cumulants  
301  if(QC) { 
302    AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
303    if(listWithWeights) qc->SetWeightsList(listWithWeights);
304    if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
305    if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
306    if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
307    // qc->SetHarmonic(2); // default is v2
308    qc->SetApplyCorrectionForNUA(kFALSE);
309    qc->SetFillMultipleControlHistograms(kFALSE);     
310    // qc->SetCalculate2DFlow(kFALSE); // default
311    // qc->SetMultiplicityWeight("combinations"); // default
312    // qc->SetMultiplicityWeight("unit");
313    // qc->SetMultiplicityWeight("multiplicity");  
314    qc->SetCalculateCumulantsVsM(kFALSE);
315    qc->SetnBinsMult(10000);
316    qc->SetMinMult(0);
317    qc->SetMaxMult(10000);      
318    qc->Init();  
319  }
320   
321  // GFC = Generating Function Cumulants 
322  if(GFC) {
323    AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
324    if(listWithWeights) gfc->SetWeightsList(listWithWeights);
325    if(usePhiWeights) gfc->SetUsePhiWeights(usePhiWeights);
326    if(usePtWeights) gfc->SetUsePtWeights(usePtWeights);
327    if(useEtaWeights) gfc->SetUseEtaWeights(useEtaWeights);
328    // calculation vs multiplicity:
329    gfc->SetCalculateVsMultiplicity(kFALSE);   
330    gfc->SetnBinsMult(10000);
331    gfc->SetMinMult(0);
332    gfc->SetMaxMult(10000);   
333    // tuning of interpolating parameters:
334    gfc->SetTuneParameters(kFALSE);
335    Double_t r0[10] = {1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7}; // up to 10 values allowed
336    for(Int_t r=0;r<10;r++) {gfc->SetTuningR0(r0[r],r);}
337    gfc->Init();
338  }
339  
340  // FQD = Fitting q-distribution 
341  if(FQD) {
342    AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
343    if(listWithWeights) fqd->SetWeightsList(listWithWeights);
344    if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);  
345    fqd->Init();
346  }
347  
348  // SP = Scalar Product 
349  if(SP) {
350    AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
351    if(listWithWeights) sp->SetWeightsList(listWithWeights);
352    sp->SetUsePhiWeights(usePhiWeights);  
353    sp->Init();
354  }
355
356  // LYZ1 = Lee-Yang Zeroes first run
357  if(LYZ1SUM) {
358    AliFlowAnalysisWithLeeYangZeros* lyz1sum = new AliFlowAnalysisWithLeeYangZeros();
359    lyz1sum->SetFirstRun(kTRUE);
360    lyz1sum->SetUseSum(kTRUE);
361    lyz1sum->Init();
362  }
363  if(LYZ1PROD) {
364    AliFlowAnalysisWithLeeYangZeros* lyz1prod = new AliFlowAnalysisWithLeeYangZeros();
365    lyz1prod->SetFirstRun(kTRUE);
366    lyz1prod->SetUseSum(kFALSE);
367    lyz1prod->Init();
368  }
369  // LYZ2 = Lee-Yang Zeroes second run
370  if(LYZ2SUM) {
371    AliFlowAnalysisWithLeeYangZeros* lyz2sum = new AliFlowAnalysisWithLeeYangZeros();
372    // read the input file from the first run 
373    TString inputFileNameLYZ2SUM = "outputLYZ1SUManalysis.root" ;
374    TFile* inputFileLYZ2SUM = new TFile(inputFileNameLYZ2SUM.Data(),"READ");
375    if(!inputFileLYZ2SUM || inputFileLYZ2SUM->IsZombie()) { 
376      cerr << " ERROR: To run LYZ2SUM you need the output file from LYZ1SUM. This file is not there! Please run LYZ1SUM first." << endl ;
377      break; 
378    }
379    else { 
380      TList* inputListLYZ2SUM = (TList*)inputFileLYZ2SUM->Get("cobjLYZ1SUM");  
381      if (!inputListLYZ2SUM) {cout<<"Input list LYZ2SUM is NULL pointer!"<<endl; break;}
382      else {
383        cout<<"LYZ2SUM input file/list read..."<<endl;
384        lyz2sum->SetFirstRunList(inputListLYZ2SUM);
385        lyz2sum->SetFirstRun(kFALSE);
386        lyz2sum->SetUseSum(kTRUE);
387        lyz2sum->Init();
388      }
389    }
390  }
391  if(LYZ2PROD) {
392    AliFlowAnalysisWithLeeYangZeros* lyz2prod = new AliFlowAnalysisWithLeeYangZeros();
393    // read the input file from the first run 
394    TString inputFileNameLYZ2PROD = "outputLYZ1PRODanalysis.root" ;
395    TFile* inputFileLYZ2PROD = new TFile(inputFileNameLYZ2PROD.Data(),"READ");
396    if(!inputFileLYZ2PROD || inputFileLYZ2PROD->IsZombie()) { 
397      cerr << " ERROR: To run LYZ2PROD you need the output file from LYZ1PROD. This file is not there! Please run LYZ1PROD first." << endl ;
398      break; 
399    }
400    else { 
401      TList* inputListLYZ2PROD = (TList*)inputFileLYZ2PROD->Get("cobjLYZ1PROD");  
402      if (!inputListLYZ2PROD) {cout<<"Input list LYZ2PROD is NULL pointer!"<<endl; break;}
403      else {
404        cout<<"LYZ2PROD input file/list read..."<<endl;
405        lyz2prod->SetFirstRunList(inputListLYZ2PROD);
406        lyz2prod->SetFirstRun(kFALSE);
407        lyz2prod->SetUseSum(kFALSE);
408        lyz2prod->Init();
409      }
410    }
411  }
412
413  // LYZEP = Lee-Yang Zeroes event plane
414  if(LYZEP) {
415    AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
416    AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
417    // read the input file from the second lyz run 
418    TString inputFileNameLYZEP = "outputLYZ2SUManalysis.root" ;
419    TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
420    if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) { 
421      cerr << " ERROR: To run LYZEP you need the output file from LYZ2SUM. This file is not there! Please run LYZ2SUM first." << endl ; 
422      break;
423    }
424    else { 
425      TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2SUM");  
426      if (!inputListLYZEP) {cout<<"Input list LYZEP is NULL pointer!"<<endl; break;}
427      else {
428        cout<<"LYZEP input file/list read..."<<endl;
429        ep   ->SetSecondRunList(inputListLYZEP);
430        lyzep->SetSecondRunList(inputListLYZEP);
431        ep   ->Init();
432        lyzep->Init();
433      }
434    }
435  }
436  //---------------------------------------------------------------------------------------
437  
438  // set the global event parameters:
439  eventMakerOnTheFly->SetNoOfLoops(iLoops);
440  eventMakerOnTheFly->SetPhiRange(phiRange*TMath::Pi()/180.);
441  eventMakerOnTheFly->SetPtRange(ptRange);
442  eventMakerOnTheFly->SetEtaRange(etaRange);
443  eventMakerOnTheFly->SetNonflowSectorMin(nonflowSectorMin*TMath::Pi()/180.);
444  eventMakerOnTheFly->SetNonflowSectorMax(nonflowSectorMax*TMath::Pi()/180.);
445  eventMakerOnTheFly->SetUseGlauberModel(bUseGlauberModel);
446  eventMakerOnTheFly->SetCreateJets(bCreateJets);
447  eventMakerOnTheFly->SetJetProbability(jetProbability);
448  eventMakerOnTheFly->SetJetTracksFraction(jetTracksFraction);
449  eventMakerOnTheFly->SetJetCone(jetCone);
450  
451  if(!bUseGlauberModel)
452  {
453   if(bMultDistrOfRPsIsGauss)
454   {
455    eventMakerOnTheFly->SetMultDistrOfRPsIsGauss(bMultDistrOfRPsIsGauss);
456    eventMakerOnTheFly->SetMultiplicityOfRP(iMultiplicityOfRP);
457    eventMakerOnTheFly->SetMultiplicitySpreadOfRP(dMultiplicitySpreadOfRP);
458   } else
459     {
460      eventMakerOnTheFly->SetMultDistrOfRPsIsGauss(bMultDistrOfRPsIsGauss);
461      eventMakerOnTheFly->SetMinMultOfRP(iMinMultOfRP);
462      eventMakerOnTheFly->SetMaxMultOfRP(iMaxMultOfRP); 
463     }
464  } // end of if(!bUseGlauberModel)
465   
466  eventMakerOnTheFly->SetTemperatureOfRP(dTemperatureOfRP);
467  eventMakerOnTheFly->SetPtDependentHarmonicV1(bPtDependentHarmonicV1);
468  eventMakerOnTheFly->SetEtaDependentHarmonicV1(bEtaDependentHarmonicV1);
469  eventMakerOnTheFly->SetPtDependentHarmonicV2(bPtDependentHarmonicV2);
470  eventMakerOnTheFly->SetEtaDependentHarmonicV2(bEtaDependentHarmonicV2);
471  eventMakerOnTheFly->SetPtDependentHarmonicV4(bPtDependentHarmonicV4);
472  eventMakerOnTheFly->SetEtaDependentHarmonicV4(bEtaDependentHarmonicV4);
473  // V1:
474  if(!(bPtDependentHarmonicV1 || bEtaDependentHarmonicV1))
475  {
476   if(!bUseGlauberModel)
477   {
478    eventMakerOnTheFly->SetV1RP(dV1RP);
479    eventMakerOnTheFly->SetV1SpreadRP(dV1SpreadRP);  
480   } 
481  } else // (pt,eta) dependent V1
482    {
483     eventMakerOnTheFly->SetV1vsPtEtaMax(dV1vsPtEtaMax);
484     if(bPtDependentHarmonicV1)
485     {
486      eventMakerOnTheFly->SetV1PtCutOff(dV1PtCutOff);  
487     }
488    }
489  // V2:
490  if(!(bPtDependentHarmonicV2 || bEtaDependentHarmonicV2)) // constant V2
491  { 
492   if(!bUseGlauberModel)
493   {
494    eventMakerOnTheFly->SetConstantV2IsSampledFromGauss(bConstantV2IsSampledFromGauss);
495    if(bConstantV2IsSampledFromGauss) // Gauss
496    {
497     eventMakerOnTheFly->SetV2RP(dV2RP);
498     eventMakerOnTheFly->SetV2SpreadRP(dV2SpreadRP);  
499    } else // uniform
500      {
501       eventMakerOnTheFly->SetMinV2RP(dMinV2RP);
502       eventMakerOnTheFly->SetMaxV2RP(dMaxV2RP);  
503      }
504   } // end of if(!bUseGlauberModel)  
505  } else // (pt,eta) dependent V2
506    {
507     eventMakerOnTheFly->SetV2vsPtEtaMax(dV2vsPtEtaMax);
508     if(bPtDependentHarmonicV2)
509     {
510      eventMakerOnTheFly->SetV2PtCutOff(dV2PtCutOff);  
511     }
512     if(bEtaDependentHarmonicV2)
513     {
514      eventMakerOnTheFly->SetV2vsEtaSpread(dV2vsEtaSpread);      
515     }
516    }  
517  // V3:
518  if(!bUseGlauberModel)
519  {
520   eventMakerOnTheFly->SetV3RP(dV3RP);
521   eventMakerOnTheFly->SetV3SpreadRP(dV3SpreadRP);  
522  }    
523  // V4:
524  if(!(bPtDependentHarmonicV4 || bEtaDependentHarmonicV4))
525  {
526   if(!bUseGlauberModel)
527   {
528    eventMakerOnTheFly->SetV4RP(dV4RP);
529    eventMakerOnTheFly->SetV4SpreadRP(dV4SpreadRP);  
530   }
531  } 
532  
533  // non-uniform acceptance:
534  if(!uniformAcceptance)
535  {
536   eventMakerOnTheFly->SetFirstSectorPhiMin(phimin1);
537   eventMakerOnTheFly->SetFirstSectorPhiMax(phimax1);
538   eventMakerOnTheFly->SetFirstSectorProbability(p1);
539   eventMakerOnTheFly->SetSecondSectorPhiMin(phimin2);
540   eventMakerOnTheFly->SetSecondSectorPhiMax(phimax2);
541   eventMakerOnTheFly->SetSecondSectorProbability(p2);
542  }
543  
544  // simple cuts for RPs and POIs:
545  AliFlowTrackSimpleCuts *cutsRP = new AliFlowTrackSimpleCuts();
546  cutsRP->SetPtMax(ptMaxRP);
547  cutsRP->SetPtMin(ptMinRP);
548  cutsRP->SetEtaMax(etaMaxRP);
549  cutsRP->SetEtaMin(etaMinRP);
550  cutsRP->SetPhiMax(phiMaxRP*TMath::Pi()/180.);
551  cutsRP->SetPhiMin(phiMinRP*TMath::Pi()/180.);
552  cutsRP->SetPID(pidRP);
553   
554  AliFlowTrackSimpleCuts *cutsPOI = new AliFlowTrackSimpleCuts();
555  cutsPOI->SetPtMax(ptMaxPOI);
556  cutsPOI->SetPtMin(ptMinPOI);
557  cutsPOI->SetEtaMax(etaMaxPOI);
558  cutsPOI->SetEtaMin(etaMinPOI);
559  cutsPOI->SetPhiMax(phiMaxPOI*TMath::Pi()/180.);
560  cutsPOI->SetPhiMin(phiMinPOI*TMath::Pi()/180.);
561  cutsPOI->SetPID(pidPOI);
562                                        
563  //---------------------------------------------------------------------------------------  
564  // create and analyze events 'on the fly':
565
566  for(Int_t i=0;i<nEvts;i++) {   
567    // creating the event with above settings:
568    AliFlowEventSimple *event = eventMakerOnTheFly->CreateEventOnTheFly(cutsRP,cutsPOI); 
569    
570    // analyzing the created event 'on the fly':
571    // do flow analysis for various methods:
572    if(MCEP)    mcep->Make(event);
573    if(QC)      qc->Make(event);
574    if(GFC)     gfc->Make(event);
575    if(FQD)     fqd->Make(event);
576    if(LYZ1SUM) lyz1sum->Make(event);
577    if(LYZ1PROD)lyz1prod->Make(event);
578    if(LYZ2SUM) lyz2sum->Make(event);
579    if(LYZ2PROD)lyz2prod->Make(event);
580    if(LYZEP)   lyzep->Make(event,ep);
581    if(SP)      sp->Make(event);
582    if(MH)      mh->Make(event);
583    if(NL)      nl->Make(event);
584    if(MCEP_AH) mcep_ah->Make(event);
585    
586    delete event;
587  } // end of for(Int_t i=0;i<nEvts;i++)
588  //---------------------------------------------------------------------------------------  
589
590  //---------------------------------------------------------------------------------------  
591  // create a new file which will hold the final results of all methods:
592  TString outputFileName = "AnalysisResults.root";  
593  TFile *outputFile = new TFile(outputFileName.Data(),"RECREATE");
594  // create a new file for each method wich will hold list with final results:
595  const Int_t nMethods = 13;
596  TString method[nMethods] = {"MCEP","SP","GFC","QC","FQD","LYZ1SUM","LYZ1PROD","LYZ2SUM","LYZ2PROD","LYZEP","MH","NL","MCEP_AH"};
597  TDirectoryFile *dirFileFinal[nMethods] = {NULL};
598  TString fileName[nMethods]; 
599  for(Int_t i=0;i<nMethods;i++)
600  {
601   // form a file name for each method:
602   fileName[i]+="output";
603   fileName[i]+=method[i].Data();
604   fileName[i]+="analysis";
605   dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data());
606  } 
607  
608  // calculating and storing the final results of default flow analysis:
609  if(MCEP)    {mcep->Finish();    mcep->WriteHistograms(dirFileFinal[0]);}
610  if(SP)      {sp->Finish();      sp->WriteHistograms(dirFileFinal[1]);}
611  if(GFC)     {gfc->Finish();     gfc->WriteHistograms(dirFileFinal[2]);}
612  if(QC)      {qc->Finish();      qc->WriteHistograms(dirFileFinal[3]);}
613  if(FQD)     {fqd->Finish();     fqd->WriteHistograms(dirFileFinal[4]);}
614  if(LYZ1SUM) {lyz1sum->Finish(); lyz1sum->WriteHistograms(dirFileFinal[5]);}
615  if(LYZ1PROD){lyz1prod->Finish();lyz1prod->WriteHistograms(dirFileFinal[6]);}
616  if(LYZ2SUM) {lyz2sum->Finish(); lyz2sum->WriteHistograms(dirFileFinal[7]);}
617  if(LYZ2PROD){lyz2prod->Finish();lyz2prod->WriteHistograms(dirFileFinal[8]);}
618  if(LYZEP)   {lyzep->Finish();   lyzep->WriteHistograms(dirFileFinal[9]);}
619  if(MH)      {mh->Finish();      mh->WriteHistograms(dirFileFinal[10]);}
620  if(NL)      {nl->Finish();      nl->WriteHistograms(dirFileFinal[11]);}
621  if(MCEP_AH) {mcep_ah->Finish(); mcep_ah->WriteHistograms(dirFileFinal[12]);}
622  //---------------------------------------------------------------------------------------  
623  
624  outputFile->Close();
625  delete outputFile;
626  
627  cout<<endl;
628  cout<<endl;
629  cout<<" ---- LANDED SUCCESSFULLY ---- "<<endl;
630  cout<<endl; 
631  
632  timer.Stop();
633  cout << endl;
634  timer.Print();
635 }
636
637 void SetupPar(char* pararchivename)
638 {
639   //Load par files, create analysis libraries
640   //For testing, if par file already decompressed and modified
641   //classes then do not decompress.
642  
643   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ; 
644   TString parpar(Form("%s.par", pararchivename)) ; 
645   if ( gSystem->AccessPathName(parpar.Data()) ) {
646     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
647     TString processline(Form(".! make %s", parpar.Data())) ; 
648     gROOT->ProcessLine(processline.Data()) ;
649     gSystem->ChangeDirectory(cdir) ; 
650     processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
651     gROOT->ProcessLine(processline.Data()) ;
652   } 
653   if ( gSystem->AccessPathName(pararchivename) ) {  
654     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
655     gROOT->ProcessLine(processline.Data());
656   }
657   
658   TString ocwd = gSystem->WorkingDirectory();
659   gSystem->ChangeDirectory(pararchivename);
660   
661   // check for BUILD.sh and execute
662   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
663     printf("*******************************\n");
664     printf("*** Building PAR archive    ***\n");
665     cout<<pararchivename<<endl;
666     printf("*******************************\n");
667     
668     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
669       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
670       return -1;
671     }
672   }
673   // check for SETUP.C and execute
674   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
675     printf("*******************************\n");
676     printf("*** Setup PAR archive       ***\n");
677     cout<<pararchivename<<endl;
678     printf("*******************************\n");
679     gROOT->Macro("PROOF-INF/SETUP.C");
680   }
681   
682   gSystem->ChangeDirectory(ocwd.Data());
683   printf("Current dir: %s\n", ocwd.Data());
684 }
685
686 void CheckUserSettings()
687 {
688  // Check if user settings make sense before taking off:
689   
690  if (LYZ1SUM && LYZ2SUM)  {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1.  "<<endl; exit(); }
691  if (LYZ1PROD && LYZ2PROD)  {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1.  "<<endl; exit(); }
692  if (LYZ2SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
693  if (LYZ1SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
694
695  if(!uniformAcceptance && phimin1 > phimax1)
696  {
697   cout<<"WARNING: you must have phimin1 < phimax1 !!!!"<<endl;
698   break;
699  }
700
701  if (!uniformAcceptance && !((phimin2 == 0.) && (phimax2 == 0.) && (p2 == 0.)) && (phimin2 < phimax1 || phimin2 > phimax2))
702  {
703   cout<<"WARNING: you must have phimin2 > phimax1 and phimin2 < phimax2 !!!!"<<endl;
704   break;
705  }
706  
707  if((phimin1 < 0 || phimin1 > 360) || (phimax1 < 0 || phimax1 > 360) || 
708     (phimin2 < 0 || phimin2 > 360) || (phimax2 < 0 || phimax2 > 360) )
709  {
710   cout<<"WARNING: you must take azimuthal angles from interval [0,360] !!!!"<<endl;
711   break;
712  }
713  
714  if((p1 < 0 || p1 > 1) || (p2 < 0 || p2 > 1))
715  {
716   cout<<"WARNING: you must take p1 and p2 from interval [0,1] !!!!"<<endl;
717   break;
718  }
719  
720  if(bUseGlauberModel)
721  {
722   if(bPtDependentHarmonicV1||bEtaDependentHarmonicV1||
723      bPtDependentHarmonicV2||bEtaDependentHarmonicV2||
724      bPtDependentHarmonicV4||bEtaDependentHarmonicV4)
725   {   
726    cout<<"WARNING: When using Glauber model pt and/or eta dependence of flow harmonics is not supported !!!!"<<endl;
727    cout<<"         Set all booleans bPtDependentHarmonic* to kFALSE in the macro."<<endl;
728    exit(0); 
729   }
730  }
731
732  if(bPtDependentHarmonicV4 && !bPtDependentHarmonicV2))
733  {
734   cout<<"WARNING: V4(pt,eta) is determined as pow(V2(pt,eta),2.) !!!!"<<endl;
735   cout<<"         You must also set bPtDependentHarmonicV2 = kTRUE in the macro."<<endl;
736   exit(0);
737  }
738    
739  if(bEtaDependentHarmonicV4 && !bEtaDependentHarmonicV2))
740  {
741   cout<<"WARNING: V4(pt,eta) is determined as pow(V2(pt,eta),2.) !!!!"<<endl;
742   cout<<"         You must also set bEtaDependentHarmonicV2 = kTRUE in the macro."<<endl;
743   exit(0); 
744  }
745   
746  if(bPtDependentHarmonicV4 && (bEtaDependentHarmonicV4 != bEtaDependentHarmonicV2))
747  {
748   cout<<"WARNING: V4(pt,eta) is determined as pow(V2(pt,eta),2.) !!!!"<<endl;
749   cout<<"         You must also have bEtaDependentHarmonicV2 = bEtaDependentHarmonicV4 in the macro."<<endl;
750   exit(0); 
751  }
752  
753  if(bEtaDependentHarmonicV4 && (bPtDependentHarmonicV4 != bPtDependentHarmonicV2))
754  {
755   cout<<"WARNING: V4(pt,eta) is determined as pow(V2(pt,eta),2.) !!!!"<<endl;
756   cout<<"         You must also have bPtDependentHarmonicV2 = bPtDependentHarmonicV4 in the macro."<<endl;
757   exit(0); 
758  }
759  
760  if(!uniformAcceptance && bCreateJets)
761  {
762   cout<<"WARNING: Jet creation is supported only for uniform acceptance !!!!"<<endl;
763   exit(0); 
764  }
765   
766 } // end of void CheckUserSettings()
767
768 void LoadLibraries(const anaModes mode) {
769   
770   //--------------------------------------
771   // Load the needed libraries most of them already loaded by aliroot
772   //--------------------------------------
773   //gSystem->Load("libTree");
774   gSystem->Load("libGeom");
775   gSystem->Load("libVMC");
776   gSystem->Load("libXMLIO");
777   gSystem->Load("libPhysics");
778   
779   //----------------------------------------------------------
780   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<< 
781   //----------------------------------------------------------
782   if (mode==mLocal) {
783     //--------------------------------------------------------
784     // If you want to use already compiled libraries 
785     // in the aliroot distribution
786     //--------------------------------------------------------
787     gSystem->Load("libSTEERBase");
788     gSystem->Load("libESD");
789     gSystem->Load("libAOD");
790     gSystem->Load("libANALYSIS");
791     gSystem->Load("libANALYSISalice");
792     gSystem->Load("libCORRFW");
793     cerr<<"libCORRFW loaded..."<<endl;
794     gSystem->Load("libPWG2flowCommon");
795     cerr<<"libPWG2flowCommon loaded..."<<endl;
796     gSystem->Load("libPWG2flowTasks");
797     cerr<<"libPWG2flowTasks loaded..."<<endl;
798   }
799   
800   else if (mode == mLocalPAR) {
801     //--------------------------------------------------------
802     //If you want to use root and par files from aliroot
803     //--------------------------------------------------------  
804      //If you want to use root and par files from aliroot
805     //--------------------------------------------------------  
806     SetupPar("STEERBase");
807     SetupPar("ESD");
808     SetupPar("AOD");
809     SetupPar("ANALYSIS");
810     SetupPar("ANALYSISalice");
811     SetupPar("PWG2AOD");
812     SetupPar("CORRFW");
813     SetupPar("PWG2flowCommon");
814     cerr<<"PWG2flowCommon.par loaded..."<<endl;
815     SetupPar("PWG2flowTasks");
816     cerr<<"PWG2flowTasks.par loaded..."<<endl;
817   }
818   
819   //---------------------------------------------------------
820   // <<<<<<<<<< Source mode >>>>>>>>>>>>
821   //---------------------------------------------------------
822   else if (mode==mLocalSource) {
823  
824     // In root inline compile
825
826     // Constants  
827     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
828     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
829     
830     // Flow event
831     gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+"); 
832     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");    
833     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");    
834     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
835     
836     // Output histosgrams
837     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
838     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
839     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
840     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
841     
842     // Functions needed for various methods
843     gROOT->LoadMacro("AliFlowCommon/AliCumulantsFunctions.cxx+");
844     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZEventPlane.cxx+");
845     
846     // Flow Analysis code for various methods
847     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMCEventPlane.cxx+"); 
848     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx+");
849     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.cxx+");
850     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.cxx+");
851     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithCumulants.cxx+");
852     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithQCumulants.cxx+"); 
853     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithFittingQDistribution.cxx+");
854     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMixedHarmonics.cxx+"); 
855     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithNestedLoops.cxx+");
856     
857     // Class to fill the FlowEvent on the fly (generate Monte Carlo events)
858     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.cxx+");   
859     
860     cout << "finished loading macros!" << endl;  
861     
862   }  
863   
864 }
865
866