.so cleanup: removed from gSystem->Load()
[u/mrichter/AliRoot.git] / PWGCF / FLOW / macros / runStarFlowAnalysis.C
CommitLineData
6167f59c 1//example script on what to do with the star events
2//run e.g. like this:
3// root runStarFlowAnalysis.C
b80d67e1 4// flow analysis method can be: (set to kTRUE or kFALSE)
5Bool_t MCEP = kTRUE;
6167f59c 6Bool_t SP = kTRUE;
6167f59c 7Bool_t GFC = kTRUE;
8Bool_t QC = kTRUE;
9Bool_t FQD = kTRUE;
b80d67e1 10Bool_t LYZ1SUM = kTRUE;
11Bool_t LYZ1PROD = kTRUE;
12Bool_t LYZ2SUM = kFALSE;
13Bool_t LYZ2PROD = kFALSE;
14Bool_t LYZEP = kFALSE;
15Bool_t MH = kFALSE; // mixed harmonics
16Bool_t NL = kFALSE; // nested loops
17Bool_t MCEP_AH = kFALSE; // MCEP in another harmonic
6167f59c 18//--------------------------------------------------------------------------------------
19
20// Weights
b80d67e1 21// use weights for Q-vector:
22Bool_t usePhiWeights = kFALSE; // phi weights (correction for non-uniform azimuthal acceptance)
23Bool_t usePtWeights = kFALSE; // pt weights
24Bool_t useEtaWeights = kFALSE; // eta weights
25
26// Define the range for eta subevents
27Double_t minA = -0.9;
28Double_t maxA = -0.01;
29Double_t minB = 0.01;
30Double_t maxB = 0.9;
6167f59c 31
32
afc49a4b 33void runStarFlowAnalysis(const char* inputDataFiles="/Users/snelling/alice_data/jthomas/testData/")
6167f59c 34{
4070f709 35 gSystem->Load("libTree");
36 gSystem->Load("libVMC");
37 gSystem->Load("libPhysics");
2e311896 38 gSystem->Load("libPWGflowBase");
6167f59c 39
40 Int_t maxNumberOfEvents = 1000000;
41
42 TStopwatch timer;
43 timer.Start();
44
45 if (LYZ1SUM && LYZ2SUM) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(1); }
46 if (LYZ1PROD && LYZ2PROD) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(1); }
47 if (LYZ2SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(1); }
48 if (LYZ1SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(1); }
49
50
51
52 //define reference particles
53 AliStarTrackCuts* rpCuts = AliStarTrackCuts::StandardCuts();
cdfe7ebe 54 rpCuts->SetPtMin(0.05);
55 rpCuts->SetPtMax(10.);
6167f59c 56
57 //define particles of interest
58 AliStarTrackCuts* poiCuts = AliStarTrackCuts::StandardCuts();
59 poiCuts->SetPtMin(0.05);
cdfe7ebe 60 poiCuts->SetPtMax(10.);
6167f59c 61
62 //define event cuts
63 AliStarEventCuts* starEventCuts = AliStarEventCuts::StandardCuts();
cdfe7ebe 64 starEventCuts-> SetCentralityIDMax(3);
65 starEventCuts-> SetCentralityIDMin(3);
6167f59c 66
67 //if the weights are used:
68 TFile *fileWithWeights = NULL;
69 TList *listWithWeights = NULL;
70
71 if(usePhiWeights||usePtWeights||useEtaWeights) {
72 fileWithWeights = TFile::Open("weights.root","READ");
73 if(fileWithWeights) {
74 listWithWeights = (TList*)fileWithWeights->Get("weights");
75 }
76 else
77 {cout << " WARNING: the file <weights.root> with weights from the previous run was not found."<<endl;
78 break;
79 }
80 }
81
b80d67e1 82 //---------------------------------------------------------------------------------------
83 // Initialize all the flow methods for default analysis:
6167f59c 84 AliFlowAnalysisWithQCumulants *qc = NULL;
85 AliFlowAnalysisWithCumulants *gfc = NULL;
86 AliFlowAnalysisWithFittingQDistribution *fqd = NULL;
b80d67e1 87 AliFlowAnalysisWithLeeYangZeros *lyz1sum = NULL;
6167f59c 88 AliFlowAnalysisWithLeeYangZeros *lyz1prod = NULL;
b80d67e1 89 AliFlowAnalysisWithLeeYangZeros *lyz2sum = NULL;
6167f59c 90 AliFlowAnalysisWithLeeYangZeros *lyz2prod = NULL;
91 AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
92 AliFlowAnalysisWithScalarProduct *sp = NULL;
6167f59c 93 AliFlowAnalysisWithMixedHarmonics *mh = NULL;
94 AliFlowAnalysisWithNestedLoops *nl = NULL;
b80d67e1 95 AliFlowAnalysisWithMCEventPlane *mcep = NULL;
96 AliFlowAnalysisWithMCEventPlane *mcep_ah = NULL;
6167f59c 97
b80d67e1 98// MCEP = monte carlo event plane
99 if (MCEP) {
100 AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
101 //mcep->SetHarmonic(2); // default is v2
102 //mcep->SetFlowOfResonances(kTRUE);
103 mcep->Init();
104 }
6167f59c 105
b80d67e1 106 // MCEP = monte carlo event plane in another harmonic:
107 if(MCEP_AH)
108 {
109 AliFlowAnalysisWithMCEventPlane *mcep_ah = new AliFlowAnalysisWithMCEventPlane();
110 mcep_ah->SetHarmonic(1);
111 mcep_ah->Init();
112 }
113
114 // Mixed harmonics:
115 if(MH)
116 {
117 AliFlowAnalysisWithMixedHarmonics *mh = new AliFlowAnalysisWithMixedHarmonics();
118 mh->SetHarmonic(1); // integer n in expression cos[n(2phi1-phi2-phi3)] = v2n*vn^2
119 mh->SetMinMultiplicity(100);
120 mh->SetNoOfMultipicityBins(5);
121 mh->SetMultipicityBinWidth(200);
122 mh->Init();
123 }
124
125 // NL = nested loops:
126 if(NL) {
127 AliFlowAnalysisWithNestedLoops *nl = new AliFlowAnalysisWithNestedLoops();
128 nl->Init();
129 }
130
131 // QC = Q-cumulants
132 if(QC) {
133 AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
134 if(listWithWeights) qc->SetWeightsList(listWithWeights);
135 if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
136 if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
137 if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
138 // qc->SetHarmonic(2); // default is v2
139 // qc->SetApplyCorrectionForNUA(kTRUE); // default
140 // qc->SetCalculate2DFlow(kFALSE); // default
141 // qc->SetMultiplicityWeight("combinations"); // default
142 // qc->SetMultiplicityWeight("unit");
143 // qc->SetMultiplicityWeight("multiplicity");
144 qc->SetnBinsMult(10000);
145 qc->SetMinMult(0);
146 qc->SetMaxMult(10000);
147 qc->Init();
148 }
6167f59c 149
b80d67e1 150 // GFC = Generating Function Cumulants
151 if(GFC) {
152 AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
153 if(listWithWeights) gfc->SetWeightsList(listWithWeights);
154 if(usePhiWeights) gfc->SetUsePhiWeights(usePhiWeights);
155 if(usePtWeights) gfc->SetUsePtWeights(usePtWeights);
156 if(useEtaWeights) gfc->SetUseEtaWeights(useEtaWeights);
157 // calculation vs multiplicity:
158 gfc->SetCalculateVsMultiplicity(kFALSE);
159 gfc->SetnBinsMult(10000);
160 gfc->SetMinMult(0);
161 gfc->SetMaxMult(10000);
162 // tuning of interpolating parameters:
163 gfc->SetTuneParameters(kFALSE);
164 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
165 for(Int_t r=0;r<10;r++) {gfc->SetTuningR0(r0[r],r);}
166 gfc->Init();
167 }
168
169 // FQD = Fitting q-distribution
170 if(FQD) {
171 AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
172 if(listWithWeights) fqd->SetWeightsList(listWithWeights);
173 if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);
174 fqd->Init();
175 }
176
177 // SP = Scalar Product
178 if(SP) {
179 AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
180 if(listWithWeights) sp->SetWeightsList(listWithWeights);
181 sp->SetUsePhiWeights(usePhiWeights);
182 sp->Init();
183 }
6167f59c 184
b80d67e1 185 // LYZ1 = Lee-Yang Zeroes first run
186 if(LYZ1SUM) {
187 AliFlowAnalysisWithLeeYangZeros* lyz1sum = new AliFlowAnalysisWithLeeYangZeros();
188 lyz1sum->SetFirstRun(kTRUE);
189 lyz1sum->SetUseSum(kTRUE);
190 lyz1sum->Init();
191 }
192 if(LYZ1PROD) {
193 AliFlowAnalysisWithLeeYangZeros* lyz1prod = new AliFlowAnalysisWithLeeYangZeros();
194 lyz1prod->SetFirstRun(kTRUE);
195 lyz1prod->SetUseSum(kFALSE);
196 lyz1prod->Init();
197 }
198 // LYZ2 = Lee-Yang Zeroes second run
199 if(LYZ2SUM) {
200 AliFlowAnalysisWithLeeYangZeros* lyz2sum = new AliFlowAnalysisWithLeeYangZeros();
201 // read the input file from the first run
202 TString inputFileNameLYZ2SUM = "outputLYZ1SUManalysis.root" ;
203 TFile* inputFileLYZ2SUM = new TFile(inputFileNameLYZ2SUM.Data(),"READ");
204 if(!inputFileLYZ2SUM || inputFileLYZ2SUM->IsZombie()) {
205 cerr << " ERROR: To run LYZ2SUM you need the output file from LYZ1SUM. This file is not there! Please run LYZ1SUM first." << endl ;
206 break;
207 }
208 else {
209 TList* inputListLYZ2SUM = (TList*)inputFileLYZ2SUM->Get("cobjLYZ1SUM");
210 if (!inputListLYZ2SUM) {cout<<"Input list LYZ2SUM is NULL pointer!"<<endl; break;}
211 else {
212 cout<<"LYZ2SUM input file/list read..."<<endl;
213 lyz2sum->SetFirstRunList(inputListLYZ2SUM);
214 lyz2sum->SetFirstRun(kFALSE);
215 lyz2sum->SetUseSum(kTRUE);
216 lyz2sum->Init();
217 }
218 }
219 }
220 if(LYZ2PROD) {
221 AliFlowAnalysisWithLeeYangZeros* lyz2prod = new AliFlowAnalysisWithLeeYangZeros();
222 // read the input file from the first run
223 TString inputFileNameLYZ2PROD = "outputLYZ1PRODanalysis.root" ;
224 TFile* inputFileLYZ2PROD = new TFile(inputFileNameLYZ2PROD.Data(),"READ");
225 if(!inputFileLYZ2PROD || inputFileLYZ2PROD->IsZombie()) {
226 cerr << " ERROR: To run LYZ2PROD you need the output file from LYZ1PROD. This file is not there! Please run LYZ1PROD first." << endl ;
227 break;
228 }
229 else {
230 TList* inputListLYZ2PROD = (TList*)inputFileLYZ2PROD->Get("cobjLYZ1PROD");
231 if (!inputListLYZ2PROD) {cout<<"Input list LYZ2PROD is NULL pointer!"<<endl; break;}
232 else {
233 cout<<"LYZ2PROD input file/list read..."<<endl;
234 lyz2prod->SetFirstRunList(inputListLYZ2PROD);
235 lyz2prod->SetFirstRun(kFALSE);
236 lyz2prod->SetUseSum(kFALSE);
237 lyz2prod->Init();
238 }
239 }
240 }
6167f59c 241
b80d67e1 242 // LYZEP = Lee-Yang Zeroes event plane
243 if(LYZEP) {
244 AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
245 AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
246 // read the input file from the second lyz run
247 TString inputFileNameLYZEP = "outputLYZ2SUManalysis.root" ;
248 TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
249 if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) {
250 cerr << " ERROR: To run LYZEP you need the output file from LYZ2SUM. This file is not there! Please run LYZ2SUM first." << endl ;
251 break;
252 }
253 else {
254 TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2SUM");
255 if (!inputListLYZEP) {cout<<"Input list LYZEP is NULL pointer!"<<endl; break;}
256 else {
257 cout<<"LYZEP input file/list read..."<<endl;
258 ep ->SetSecondRunList(inputListLYZEP);
259 lyzep->SetSecondRunList(inputListLYZEP);
260 ep ->Init();
261 lyzep->Init();
262 }
263 }
264 }
265 //---------------------------------------------------------------------------------------
6167f59c 266
267
268 Int_t i=0;
afc49a4b 269 AliStarEventReader starReader(inputDataFiles) ;
6167f59c 270 while ( starReader.GetNextEvent() ) // Get next event
271 {
272 AliStarEvent* starEvent = starReader.GetEvent();
273 if ( !starEventCuts->PassesCuts(starEvent) ) continue; // Test if the event is good
274
275 AliFlowEventSimple* flowEvent = new AliFlowEventStar(starEvent,rpCuts,poiCuts); // make a flow event from a star event (aka "the magic")
b80d67e1 276 flowEvent->TagSubeventsInEta(minA, maxA, minB, maxB );
6167f59c 277
278 /////analysis here////////////////
279
280 // do flow analysis for various methods
281 if(MCEP) mcep->Make(flowEvent);
282 if(QC) qc->Make(flowEvent);
283 if(GFC) gfc->Make(flowEvent);
284 if(FQD) fqd->Make(flowEvent);
285 if(LYZ1SUM) lyz1sum->Make(flowEvent);
286 if(LYZ1PROD)lyz1prod->Make(flowEvent);
287 if(LYZ2SUM) lyz2sum->Make(flowEvent);
288 if(LYZ2PROD)lyz2prod->Make(flowEvent);
289 if(LYZEP) lyzep->Make(flowEvent,ep);
290 if(SP) sp->Make(flowEvent);
291 if(MH) mh->Make(flowEvent);
292 if(NL) nl->Make(flowEvent);
b80d67e1 293 if(MCEP_AH) mcep_ah->Make(flowEvent);
6167f59c 294
295 //////////////////////////////////
296
b80d67e1 297 flowEvent->Print();
6167f59c 298
299 delete flowEvent;
300
301 i++;
302 if (i>maxNumberOfEvents) break;
303 }
304
b80d67e1 305 //---------------------------------------------------------------------------------------
6167f59c 306 // create a new file which will hold the final results of all methods:
307 TString outputFileName = "AnalysisResults.root";
308 TFile *outputFile = new TFile(outputFileName.Data(),"RECREATE");
309 // create a new file for each method wich will hold list with final results:
b80d67e1 310 const Int_t nMethods = 13;
311 TString method[nMethods] = {"MCEP","SP","GFC","QC","FQD","LYZ1SUM","LYZ1PROD","LYZ2SUM","LYZ2PROD","LYZEP","MH","NL","MCEP_AH"};
6167f59c 312 TDirectoryFile *dirFileFinal[nMethods] = {NULL};
b80d67e1 313 TString fileName[nMethods];
6167f59c 314 for(Int_t i=0;i<nMethods;i++)
315 {
316 // form a file name for each method:
b80d67e1 317 fileName[i]+="output";
318 fileName[i]+=method[i].Data();
319 fileName[i]+="analysis";
320 dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data());
6167f59c 321 }
322
323 // calculating and storing the final results of default flow analysis:
324 if(MCEP) {mcep->Finish(); mcep->WriteHistograms(dirFileFinal[0]);}
325 if(SP) {sp->Finish(); sp->WriteHistograms(dirFileFinal[1]);}
326 if(GFC) {gfc->Finish(); gfc->WriteHistograms(dirFileFinal[2]);}
327 if(QC) {qc->Finish(); qc->WriteHistograms(dirFileFinal[3]);}
328 if(FQD) {fqd->Finish(); fqd->WriteHistograms(dirFileFinal[4]);}
329 if(LYZ1SUM) {lyz1sum->Finish(); lyz1sum->WriteHistograms(dirFileFinal[5]);}
330 if(LYZ1PROD){lyz1prod->Finish();lyz1prod->WriteHistograms(dirFileFinal[6]);}
331 if(LYZ2SUM) {lyz2sum->Finish(); lyz2sum->WriteHistograms(dirFileFinal[7]);}
332 if(LYZ2PROD){lyz2prod->Finish();lyz2prod->WriteHistograms(dirFileFinal[8]);}
333 if(LYZEP) {lyzep->Finish(); lyzep->WriteHistograms(dirFileFinal[9]);}
334 if(MH) {mh->Finish(); mh->WriteHistograms(dirFileFinal[10]);}
335 if(NL) {nl->Finish(); nl->WriteHistograms(dirFileFinal[11]);}
b80d67e1 336 if(MCEP_AH) {mcep_ah->Finish(); mcep_ah->WriteHistograms(dirFileFinal[12]);}
6167f59c 337 //---------------------------------------------------------------------------------------
338
339 outputFile->Close();
340 delete outputFile;
341
b80d67e1 342 cout<<endl;
343 cout<<endl;
344 cout<<" ---- LANDED SUCCESSFULLY ---- "<<endl;
345 cout<<endl;
6167f59c 346
347 timer.Stop();
348 cout << endl;
349 timer.Print();
350
351 delete rpCuts;
352 delete poiCuts;
353 delete starEventCuts;
354}