default: run all methods
[u/mrichter/AliRoot.git] / PWG2 / FLOW / macros / runFlowAnalysisOnTheFly.C
CommitLineData
93ff27bd 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)
06868ea2 9Bool_t MCEP = kTRUE;
c741f5d0 10Bool_t SP = kTRUE;
06868ea2 11Bool_t GFC = kTRUE;
12Bool_t QC = kTRUE;
13Bool_t FQD = kTRUE;
14Bool_t LYZ1SUM = kTRUE;
15Bool_t LYZ1PROD = kTRUE;
c741f5d0 16Bool_t LYZ2SUM = kFALSE;
17Bool_t LYZ2PROD = kFALSE;
5b40431d 18Bool_t LYZEP = kFALSE;
93ff27bd 19//--------------------------------------------------------------------------------------
20
7fb984f6 21// Weights
22// use weights for Q-vector:
23Bool_t usePhiWeights = kFALSE; // phi weights (correction for non-uniform azimuthal acceptance)
24Bool_t usePtWeights = kFALSE; // pt weights
25Bool_t useEtaWeights = kFALSE; // eta weights
26
ce4a88f5 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):
29Bool_t GFC_Additional_Analysis = kFALSE;
30Bool_t QC_Additional_Analysis = kFALSE;
31Bool_t FQD_Additional_Analysis = kFALSE;
32
7a01f4a7 33// Define the range for eta subevents
34Double_t minA = -0.9;
35Double_t maxA = -0.01;
36Double_t minB = 0.01;
37Double_t maxB = 0.9;
38
7fb984f6 39// Parameters for the simulation of events 'on the fly':
7a01f4a7 40Double_t dTemperatureOfRP = 0.44; // 'temperature' of RPs in GeV/c (increase this parameter to get more high pt RPs)
4d603348 41
7a01f4a7 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
47Bool_t bSameSeed = kTRUE;
4d603348 48
e23cbd39 49
7a01f4a7 50//===NONFLOW=============================================
51// number of times to use each track (to simulate nonflow)
52Int_t iLoops = 1;
678f8252 53
93ff27bd 54
7a01f4a7 55//===FLOW HARMONICS===============================================================
56// harmonics V1, V2, V4... are constant (kTRUE) or functions of pt and eta (kFALSE)
57Bool_t bConstantHarmonics = kTRUE;
58
4d603348 59// if you use (pt,eta) dependent harmonics (bConstantHarmonics = kFALSE):
60Double_t dPtCutOff = 2.0; // V2(pt) is linear up to pt = 2 GeV and for pt > 2 GeV it is constant: V2(pt) = dVRPMax
cf90787f 61Double_t dV2RPMax = 0.20; // maximum value of V2(pt) for pt >= 2GeV
4d603348 62
7a01f4a7 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
68Bool_t bV2DistrOfRPsIsGauss = kTRUE;
69
55381065 70// if you use constant harmonics (bConstantHarmonics = kTRUE) (i.e. no pt dependence):
71Double_t dV2RP = 0.05; // elliptic flow of RPs (if sampled from Gaussian)
72Double_t dV2SpreadRP = 0.0; // elliptic flow spread of RPs (if sampled from Gaussian)
7a01f4a7 73Double_t dMinV2RP = 0.02; // minimal elliptic flow of RPs (if sampled uniformly)
74Double_t dMaxV2RP = 0.08; // maximal elliptic flow of RPs (if sampled uniformly)
93ff27bd 75
e8c7e66c 76Double_t dV1RP = 0.0; // directed flow of RPs
77Double_t dV1SpreadRP = 0.0; // directed flow spread of RPs
93ff27bd 78
4d603348 79Double_t dV4RP = 0.0; // harmonic V4 of RPs (to be improved: name needed)
80Double_t dV4SpreadRP = 0.0; // harmonic V4's spread of RPs (to be improved: name needed)
4d603348 81
7a01f4a7 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
89Bool_t bMultDistrOfRPsIsGauss = kTRUE;
90
91Int_t iMultiplicityOfRP = 500; // mean multiplicity of RPs (if sampled from Gaussian)
92Double_t dMultiplicitySpreadOfRP = 0; // multiplicity spread of RPs (if sampled from Gaussian)
93Int_t iMinMultOfRP = 50; // minimal multiplicity of RPs (if sampled uniformly)
94Int_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.
107Bool_t uniformAcceptance = kTRUE;
108
e23cbd39 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:
113Double_t phimin1 = 60; // first non-uniform sector starts at this azimuth
114Double_t phimax1 = 120; // first non-uniform sector ends at this azimuth
7a01f4a7 115Double_t p1 = 0.6; // e.g. if p1 = 0 all particles emitted in phimin1 < phi < phimax1 are blocked
e23cbd39 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):
7a01f4a7 119Double_t phimin2 = 200.0; // second non-uniform sector starts at this azimuth (make sure phimin2 > phimax1 !!!!)
120Double_t phimax2 = 210.0; // second non-uniform sector ends at this azimuth
121Double_t p2 = 0.01;
122
123//=====================================================================================================
e23cbd39 124
4d603348 125enum anaModes {mLocal,mLocalSource,mLocalPAR};
93ff27bd 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
7a01f4a7 130int runFlowAnalysisOnTheFly( Int_t nEvts=44, Int_t mode=mLocal)
93ff27bd 131{
132 TStopwatch timer;
133 timer.Start();
134
c741f5d0 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(); }
e23cbd39 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 }
93ff27bd 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
93ff27bd 175 LoadLibraries(mode);
176
4d603348 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 }
5b8dccde 184
4d603348 185 if(!bSameSeed)
186 {
187 TTimeStamp dt;
188 sseed = dt.GetNanoSec()/1000;
189 }
190
e23cbd39 191 cout<<endl;
192 cout<<"Seed for the random generators is "<<sseed<<endl;
193 cout<<endl;
194
93ff27bd 195 //---------------------------------------------------------------------------------------
7fb984f6 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 //---------------------------------------------------------------------------------------
93ff27bd 212 // Initialize the flowevent maker
5b8dccde 213 AliFlowEventSimpleMakerOnTheFly* eventMakerOnTheFly = new AliFlowEventSimpleMakerOnTheFly(sseed);
4d603348 214 eventMakerOnTheFly->Init();
7a01f4a7 215 //set the range for eta subevents
216 eventMakerOnTheFly->SetSubeventEtaRange(minA,maxA,minB,maxB);
93ff27bd 217
93ff27bd 218 //---------------------------------------------------------------------------------------
ce4a88f5 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;
93ff27bd 230
231 // MCEP = monte carlo event plane
232 if (MCEP) {
233 AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
afa4af05 234 // mcep->SetHarmonic(2); // default is v2
93ff27bd 235 mcep->Init();
236 }
237
ce4a88f5 238 // QC = Q-cumulants
93ff27bd 239 if(QC) {
240 AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
7fb984f6 241 if(listWithWeights) qc->SetWeightsList(listWithWeights);
242 if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
243 if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
244 if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
afa4af05 245 // qc->SetHarmonic(2); // default is v2
246 // qc->SetApplyCorrectionForNUA(kTRUE); // default
a5b7efd0 247 // qc->SetCalculate2DFlow(kFALSE); // default
248 // qc->SetMultiplicityWeight("combinations"); // default
249 // qc->SetMultiplicityWeight("unit");
250 // qc->SetMultiplicityWeight("multiplicity");
ce4a88f5 251 qc->Init();
93ff27bd 252 }
253
254 // GFC = Generating Function Cumulants
255 if(GFC) {
256 AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
7fb984f6 257 if(listWithWeights) gfc->SetWeightsList(listWithWeights);
258 if(usePhiWeights) gfc->SetUsePhiWeights(usePhiWeights);
259 if(usePtWeights) gfc->SetUsePtWeights(usePtWeights);
260 if(useEtaWeights) gfc->SetUseEtaWeights(useEtaWeights);
ce4a88f5 261 gfc->Init();
93ff27bd 262 }
263
264 // FQD = Fitting q-distribution
265 if(FQD) {
ce4a88f5 266 AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
7fb984f6 267 if(listWithWeights) fqd->SetWeightsList(listWithWeights);
268 if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);
ce4a88f5 269 fqd->Init();
93ff27bd 270 }
271
272 // SP = Scalar Product
273 if(SP) {
274 AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
7a01f4a7 275 if(listWithWeights) sp->SetWeightsList(listWithWeights);
276 sp->SetUsePhiWeights(usePhiWeights);
93ff27bd 277 sp->Init();
278 }
279
280 // LYZ1 = Lee-Yang Zeroes first run
c741f5d0 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();
93ff27bd 292 }
93ff27bd 293 // LYZ2 = Lee-Yang Zeroes second run
c741f5d0 294 if(LYZ2SUM) {
295 AliFlowAnalysisWithLeeYangZeros* lyz2sum = new AliFlowAnalysisWithLeeYangZeros();
93ff27bd 296 // read the input file from the first run
c741f5d0 297 TString inputFileNameLYZ2SUM = "outputLYZ1SUManalysis.root" ;
298 TFile* inputFileLYZ2SUM = new TFile(inputFileNameLYZ2SUM.Data(),"READ");
299 if(!inputFileLYZ2SUM || inputFileLYZ2SUM->IsZombie()) {
5b40431d 300 cerr << " ERROR: To run LYZ2SUM you need the output file from LYZ1SUM. This file is not there! Please run LYZ1SUM first." << endl ;
93ff27bd 301 break;
302 }
303 else {
c741f5d0 304 TList* inputListLYZ2SUM = (TList*)inputFileLYZ2SUM->Get("cobjLYZ1SUM");
305 if (!inputListLYZ2SUM) {cout<<"Input list LYZ2SUM is NULL pointer!"<<endl; break;}
93ff27bd 306 else {
c741f5d0 307 cout<<"LYZ2SUM input file/list read..."<<endl;
308 lyz2sum->SetFirstRunList(inputListLYZ2SUM);
309 lyz2sum->SetFirstRun(kFALSE);
310 lyz2sum->SetUseSum(kTRUE);
311 lyz2sum->Init();
93ff27bd 312 }
313 }
314 }
c741f5d0 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()) {
5b40431d 321 cerr << " ERROR: To run LYZ2PROD you need the output file from LYZ1PROD. This file is not there! Please run LYZ1PROD first." << endl ;
c741f5d0 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
93ff27bd 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
c741f5d0 342 TString inputFileNameLYZEP = "outputLYZ2SUManalysis.root" ;
93ff27bd 343 TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
344 if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) {
5b40431d 345 cerr << " ERROR: To run LYZEP you need the output file from LYZ2SUM. This file is not there! Please run LYZ2SUM first." << endl ;
93ff27bd 346 break;
347 }
348 else {
c741f5d0 349 TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2SUM");
350 if (!inputListLYZEP) {cout<<"Input list LYZEP is NULL pointer!"<<endl; break;}
93ff27bd 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 //---------------------------------------------------------------------------------------
ce4a88f5 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 //---------------------------------------------------------------------------------------
93ff27bd 422
e8c7e66c 423 // set the global event parameters:
aaa62cf0 424 eventMakerOnTheFly->SetNoOfLoops(iLoops);
678f8252 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
cf90787f 438 eventMakerOnTheFly->SetTemperatureOfRP(dTemperatureOfRP);
4d603348 439
e8c7e66c 440 eventMakerOnTheFly->SetV1RP(dV1RP);
441 eventMakerOnTheFly->SetV1SpreadRP(dV1SpreadRP);
4d603348 442 eventMakerOnTheFly->SetV4RP(dV4RP);
443 eventMakerOnTheFly->SetV4SpreadRP(dV4SpreadRP);
444
445 // constant harmonic V2:
446 if(bConstantHarmonics)
447 {
448 eventMakerOnTheFly->SetUseConstantHarmonics(bConstantHarmonics);
55381065 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 }
4d603348 460 }
e23cbd39 461
4d603348 462 // (pt,eta) dependent harmonic V2:
463 if(!bConstantHarmonics)
464 {
465 eventMakerOnTheFly->SetUseConstantHarmonics(bConstantHarmonics);
466 eventMakerOnTheFly->SetV2RPMax(dV2RPMax);
467 eventMakerOnTheFly->SetPtCutOff(dPtCutOff);
468 }
e23cbd39 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 }
4d603348 480
93ff27bd 481 //---------------------------------------------------------------------------------------
482 // create and analyze events 'on the fly':
483
93ff27bd 484 for(Int_t i=0;i<nEvts;i++) {
485 // creating the event with above settings:
93ff27bd 486 AliFlowEventSimple *event = eventMakerOnTheFly->CreateEventOnTheFly();
93ff27bd 487
488 // analyzing the created event 'on the fly':
93ff27bd 489 // do flow analysis for various methods:
c741f5d0 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);
93ff27bd 500
ce4a88f5 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
93ff27bd 517 delete event;
518 } // end of for(Int_t i=0;i<nEvts;i++)
519 //---------------------------------------------------------------------------------------
520
521
522
523 //---------------------------------------------------------------------------------------
ad87ae62 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
ce4a88f5 541 // calculating and storing the final results of default flow analysis:
ad87ae62 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]);}
93ff27bd 552 //---------------------------------------------------------------------------------------
553
ad87ae62 554 outputFile->Close();
555 delete outputFile;
93ff27bd 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
567void 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
616void LoadLibraries(const anaModes mode) {
617
618 //--------------------------------------
619 // Load the needed libraries most of them already loaded by aliroot
620 //--------------------------------------
ad87ae62 621 //gSystem->Load("libTree");
5d040cf3 622 gSystem->Load("libGeom");
623 gSystem->Load("libVMC");
624 gSystem->Load("libXMLIO");
625 gSystem->Load("libPhysics");
93ff27bd 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");
5d040cf3 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;
93ff27bd 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+");
93ff27bd 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+");
ce4a88f5 705 gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithFittingQDistribution.cxx+");
93ff27bd 706
a47dc7e4 707 // Class to fill the FlowEvent on the fly (generate Monte Carlo events)
708 gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.cxx+");
93ff27bd 709
710 cout << "finished loading macros!" << endl;
711
712 }
713
714}
715
716