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