1 //example script on what to do with the star events
3 // root runStarFlowAnalysis.C
4 // flow analysis method can be: (set to kTRUE or kFALSE)
10 Bool_t LYZ1SUM = kTRUE;
11 Bool_t LYZ1PROD = kTRUE;
12 Bool_t LYZ2SUM = kFALSE;
13 Bool_t LYZ2PROD = kFALSE;
14 Bool_t LYZEP = kFALSE;
15 Bool_t MH = kFALSE; // mixed harmonics
16 Bool_t NL = kFALSE; // nested loops
17 Bool_t MCEP_AH = kFALSE; // MCEP in another harmonic
18 //--------------------------------------------------------------------------------------
21 // use weights for Q-vector:
22 Bool_t usePhiWeights = kFALSE; // phi weights (correction for non-uniform azimuthal acceptance)
23 Bool_t usePtWeights = kFALSE; // pt weights
24 Bool_t useEtaWeights = kFALSE; // eta weights
26 // Define the range for eta subevents
28 Double_t maxA = -0.01;
33 void runStarFlowAnalysis(const char* inputDataFiles="/Users/snelling/alice_data/jthomas/testData/")
35 gSystem->Load("libTree.so");
36 gSystem->Load("libVMC.so");
37 gSystem->Load("libPhysics.so");
38 gSystem->Load("libPWGflowBase");
40 Int_t maxNumberOfEvents = 1000000;
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); }
52 //define reference particles
53 AliStarTrackCuts* rpCuts = AliStarTrackCuts::StandardCuts();
54 rpCuts->SetPtMin(0.05);
55 rpCuts->SetPtMax(10.);
57 //define particles of interest
58 AliStarTrackCuts* poiCuts = AliStarTrackCuts::StandardCuts();
59 poiCuts->SetPtMin(0.05);
60 poiCuts->SetPtMax(10.);
63 AliStarEventCuts* starEventCuts = AliStarEventCuts::StandardCuts();
64 starEventCuts-> SetCentralityIDMax(3);
65 starEventCuts-> SetCentralityIDMin(3);
67 //if the weights are used:
68 TFile *fileWithWeights = NULL;
69 TList *listWithWeights = NULL;
71 if(usePhiWeights||usePtWeights||useEtaWeights) {
72 fileWithWeights = TFile::Open("weights.root","READ");
74 listWithWeights = (TList*)fileWithWeights->Get("weights");
77 {cout << " WARNING: the file <weights.root> with weights from the previous run was not found."<<endl;
82 //---------------------------------------------------------------------------------------
83 // Initialize all the flow methods for default analysis:
84 AliFlowAnalysisWithQCumulants *qc = NULL;
85 AliFlowAnalysisWithCumulants *gfc = NULL;
86 AliFlowAnalysisWithFittingQDistribution *fqd = NULL;
87 AliFlowAnalysisWithLeeYangZeros *lyz1sum = NULL;
88 AliFlowAnalysisWithLeeYangZeros *lyz1prod = NULL;
89 AliFlowAnalysisWithLeeYangZeros *lyz2sum = NULL;
90 AliFlowAnalysisWithLeeYangZeros *lyz2prod = NULL;
91 AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
92 AliFlowAnalysisWithScalarProduct *sp = NULL;
93 AliFlowAnalysisWithMixedHarmonics *mh = NULL;
94 AliFlowAnalysisWithNestedLoops *nl = NULL;
95 AliFlowAnalysisWithMCEventPlane *mcep = NULL;
96 AliFlowAnalysisWithMCEventPlane *mcep_ah = NULL;
98 // MCEP = monte carlo event plane
100 AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
101 //mcep->SetHarmonic(2); // default is v2
102 //mcep->SetFlowOfResonances(kTRUE);
106 // MCEP = monte carlo event plane in another harmonic:
109 AliFlowAnalysisWithMCEventPlane *mcep_ah = new AliFlowAnalysisWithMCEventPlane();
110 mcep_ah->SetHarmonic(1);
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);
125 // NL = nested loops:
127 AliFlowAnalysisWithNestedLoops *nl = new AliFlowAnalysisWithNestedLoops();
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);
146 qc->SetMaxMult(10000);
150 // GFC = Generating Function Cumulants
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);
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);}
169 // FQD = Fitting q-distribution
171 AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
172 if(listWithWeights) fqd->SetWeightsList(listWithWeights);
173 if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);
177 // SP = Scalar Product
179 AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
180 if(listWithWeights) sp->SetWeightsList(listWithWeights);
181 sp->SetUsePhiWeights(usePhiWeights);
185 // LYZ1 = Lee-Yang Zeroes first run
187 AliFlowAnalysisWithLeeYangZeros* lyz1sum = new AliFlowAnalysisWithLeeYangZeros();
188 lyz1sum->SetFirstRun(kTRUE);
189 lyz1sum->SetUseSum(kTRUE);
193 AliFlowAnalysisWithLeeYangZeros* lyz1prod = new AliFlowAnalysisWithLeeYangZeros();
194 lyz1prod->SetFirstRun(kTRUE);
195 lyz1prod->SetUseSum(kFALSE);
198 // LYZ2 = Lee-Yang Zeroes second run
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 ;
209 TList* inputListLYZ2SUM = (TList*)inputFileLYZ2SUM->Get("cobjLYZ1SUM");
210 if (!inputListLYZ2SUM) {cout<<"Input list LYZ2SUM is NULL pointer!"<<endl; break;}
212 cout<<"LYZ2SUM input file/list read..."<<endl;
213 lyz2sum->SetFirstRunList(inputListLYZ2SUM);
214 lyz2sum->SetFirstRun(kFALSE);
215 lyz2sum->SetUseSum(kTRUE);
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 ;
230 TList* inputListLYZ2PROD = (TList*)inputFileLYZ2PROD->Get("cobjLYZ1PROD");
231 if (!inputListLYZ2PROD) {cout<<"Input list LYZ2PROD is NULL pointer!"<<endl; break;}
233 cout<<"LYZ2PROD input file/list read..."<<endl;
234 lyz2prod->SetFirstRunList(inputListLYZ2PROD);
235 lyz2prod->SetFirstRun(kFALSE);
236 lyz2prod->SetUseSum(kFALSE);
242 // LYZEP = Lee-Yang Zeroes event plane
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 ;
254 TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2SUM");
255 if (!inputListLYZEP) {cout<<"Input list LYZEP is NULL pointer!"<<endl; break;}
257 cout<<"LYZEP input file/list read..."<<endl;
258 ep ->SetSecondRunList(inputListLYZEP);
259 lyzep->SetSecondRunList(inputListLYZEP);
265 //---------------------------------------------------------------------------------------
269 AliStarEventReader starReader(inputDataFiles) ;
270 while ( starReader.GetNextEvent() ) // Get next event
272 AliStarEvent* starEvent = starReader.GetEvent();
273 if ( !starEventCuts->PassesCuts(starEvent) ) continue; // Test if the event is good
275 AliFlowEventSimple* flowEvent = new AliFlowEventStar(starEvent,rpCuts,poiCuts); // make a flow event from a star event (aka "the magic")
276 flowEvent->TagSubeventsInEta(minA, maxA, minB, maxB );
278 /////analysis here////////////////
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);
293 if(MCEP_AH) mcep_ah->Make(flowEvent);
295 //////////////////////////////////
302 if (i>maxNumberOfEvents) break;
305 //---------------------------------------------------------------------------------------
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:
310 const Int_t nMethods = 13;
311 TString method[nMethods] = {"MCEP","SP","GFC","QC","FQD","LYZ1SUM","LYZ1PROD","LYZ2SUM","LYZ2PROD","LYZEP","MH","NL","MCEP_AH"};
312 TDirectoryFile *dirFileFinal[nMethods] = {NULL};
313 TString fileName[nMethods];
314 for(Int_t i=0;i<nMethods;i++)
316 // form a file name for each method:
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());
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]);}
336 if(MCEP_AH) {mcep_ah->Finish(); mcep_ah->WriteHistograms(dirFileFinal[12]);}
337 //---------------------------------------------------------------------------------------
344 cout<<" ---- LANDED SUCCESSFULLY ---- "<<endl;
353 delete starEventCuts;