]>
Commit | Line | Data |
---|---|---|
09c93afd | 1 | #include "TStopwatch.h" |
2 | #include "Riostream.h" | |
3 | #include "TFile.h" | |
4 | #include "TTimer.h" | |
5 | #include "TPythia6.h" | |
6 | #include "TMath.h" | |
7 | ||
8 | /* | |
9 | * macro to generate flow events on the fly | |
10 | * if the macro is called with kFALSE as argument, | |
11 | * an example function of how to do analysis on the output is called | |
12 | * | |
90555064 | 13 | * at the end of the macro, an example function is given that can be used |
14 | * to connect the on the fly events to the AliAnalysistwoParticleResonanceFlowTask | |
15 | * | |
09c93afd | 16 | * Redmer Alexander Bertens (rbertens@cern.ch) |
17 | * Carlos Eugenio Perez Lara | |
18 | * Andrea Dubla | |
19 | */ | |
20 | ||
21 | class AliFlowOnTheFlyGenerator; | |
90555064 | 22 | class AliAnalysisTwoParticleResonanceFlowTask; |
23 | class AliFlowAnalysisWithQCumulants; | |
09c93afd | 24 | // main function |
25 | // | |
26 | // you can call this macro in different modes. | |
27 | // a) default mode will launch the event genetor and write the output to file | |
28 | // b) called with kFALSE the macro will read events from file and do an example analysis | |
29 | // this could be expanded to support more modes (e.g. not writing to file but immediately | |
30 | // doing the desired v2 analysis) | |
31 | void GenerateEventsOnTheFly(Bool_t generate = kTRUE) { | |
32 | if(!generate) { | |
90555064 | 33 | LoadLibraries(); |
09c93afd | 34 | AnalyseEventsFromFile(); |
35 | return; | |
36 | } | |
37 | TStopwatch aa; | |
38 | aa.Start(); | |
39 | // load libraries | |
40 | if(!LoadLibraries()) { | |
41 | printf(" > problem loading libraries, maybe you need to fix your includes or pythia path <"); | |
42 | exit(0); | |
43 | } | |
44 | ///========= GENERATOR INTERFACE ======/// | |
45 | // this is the generator interface. here you can specify what your event should look like in terms of | |
46 | // mother particles (define the available species and their multiplicities and define a background) | |
47 | // see the decayer interface for decayer specifics | |
48 | // in this example macro, pythia is used as a decayer. | |
49 | // however, any decayer that inherits from TVirualMCDecayer can be used | |
50 | // | |
51 | // by default, all tracks are made rp's, so 'offline' tagging is necessary. | |
52 | // decayed particles are not saved in the output events. | |
53 | // output is built up from files containing AliFlowEventSimple objects, which in turn consist of | |
54 | // AliFlowSimpleTracks, in which a few things should be noted: | |
55 | // 1) Charge() tells whether a particle is a primary (-1) or a secondary (1) | |
56 | // 2) Weight() stores the longitudinal momentum (Pz) | |
57 | // 3) GetID() gives the pdg value of the genreated track | |
58 | // ===================================/// | |
59 | // | |
60 | // a) select which species you want to have in your initial event (pdg code) | |
61 | Int_t nSpecies = 7; // number of poi species | |
62 | Int_t species[] = { 22, 333, 221, 311, 3122, 3334, 313}; // gamma, phi eta, k0, lambda, omega, kstar | |
63 | // b) specify their multiplicities per event | |
64 | Int_t multiplicities[] = { 10, 10, 10, 10, 10, 10, 10}; // so we get 10 particles of each species here | |
65 | // c) set the multiplicity for charged background (mixture of pi, k, pi) | |
66 | Int_t setChargedBackground = 300; | |
67 | // d) poissonian fluctuations for above multiplicities | |
68 | Bool_t fluctuate = kFALSE; | |
69 | // e) which tracks should get a flow boost ? | |
70 | Bool_t addMotherV2 = kTRUE; | |
71 | Bool_t addDaughterV2 = kFALSE; | |
72 | // f) for each species (poi or rp, or species that will be created by the decayer) you can specify 'custom' | |
73 | // functions for pt, v2 or v3. this will override defaults | |
74 | // exapmle: modify the v2 for species 310 | |
75 | TString myFunction = "0.05*TMath::Sin(x)+.07"; | |
76 | Int_t modifyV2 = 310; | |
77 | // and call - once an instance of the class is created (see below) | |
78 | // AliFlowOnTheFlyEventGenerator::SetPtDependentV2(myFunction.Data(), modifyV2) | |
79 | ||
80 | // == done with the flow stuff, proceed to the decayer, qa and writing to file == // | |
81 | ||
82 | // g) decayer modes. select which decays needs to be forced (in the example by pythia) | |
83 | // note that by default gamma's are 'decayed' to e+ e- pairs | |
84 | Int_t decayModes[] = { TPythia6Decayer::kHadronicD }; // some decay modes | |
85 | // h) number of events you want to generate. | |
1d89e44d | 86 | const int events = 50000; |
09c93afd | 87 | // i) write output to file |
88 | Bool_t writeOutput = kTRUE; | |
89 | // j) do qa. qa histograms will be created at the end of generation adn written to a rootfile in your working directory | |
90 | Bool_t qa = kTRUE; | |
91 | // k) write analysis to an output file | |
92 | TString trunkName = "OnTheFlyEvent"; // trunk of output file name | |
93 | // l) specify the max number of events per file (new files will be generated automatically) | |
1d89e44d | 94 | const int maxEventsPerFile = 1000; // specify the maximum number of events that is stored per file |
09c93afd | 95 | // note that events are stored temporarily in RAM memory before writing to file |
96 | // so setting this number to a very large value could lead to an unresponsive system | |
97 | ///====== END OF GENERATOR INTERFACE===/// | |
98 | // from here there's nothing that needs to be changed for the event generation | |
99 | ||
100 | int nFileCounter(0); | |
101 | TClonesArray *event = new TClonesArray("TParticle",5000); | |
102 | TObjArray *dataContainer = new TObjArray(maxEventsPerFile); | |
103 | dataContainer->SetOwner(kTRUE); | |
104 | TFile* _tmpfile = 0x0; | |
105 | // setup the decayer | |
106 | TPythia6Decayer* decayer = new TPythia6Decayer(); | |
107 | Int_t nDecayModes(sizeof(decayModes)/sizeof(Int_t)); | |
108 | // get an estimate of no of tracks per event | |
109 | Int_t nMult(setChargedBackground); | |
110 | for(Int_t i(0); i < nSpecies; i++) nMult += species[i]; | |
111 | for(Int_t i(0); i < nDecayModes; i++) decayer->SetForceDecay(decayModes[i]); | |
112 | // get an instance of the class that we'll use to generate events | |
113 | AliFlowOnTheFlyEventGenerator* eventGenerator = new AliFlowOnTheFlyEventGenerator( qa, // make some QA plots | |
1d89e44d | 114 | 3, // introduce flow fluctuations |
09c93afd | 115 | nMult, // set initial value for evnet multiplcity |
116 | decayer, // specify the decayer (NULL is no decay) | |
117 | addMotherV2, // add v2 for mothers | |
118 | NULL, // add v3 for mothers (not implemented yet) | |
119 | addDaughterV2, // add v2 for daughters | |
120 | NULL); // add v3 for daughters (not implemented yet) | |
121 | // example of how to pass a custom function for v2 for a certain species to the generator | |
122 | // eventGenerator->SetPtDependentV2(myFunction.Data(), modifyV2); | |
123 | for(int i(0); i < events; i++) { // event generator loop | |
124 | // prepare I/O necessities | |
125 | if(!_tmpfile) _tmpfile = new TFile(Form("%s_%i.root", trunkName.Data(), nFileCounter), "RECREATE"); // new file, overwrite if already exists | |
126 | // for each event, one can embed a TClonesArray of type TParticle*, e.g. do | |
127 | // eventGenerator->EmbedEvent(myEmbeddedEvent); | |
128 | // now that the generator is prepared, generate the actual event | |
129 | AliFlowEventSimple* _tmpevent = eventGenerator->GenerateOnTheFlyEvent(event, nSpecies, species, multiplicities, setChargedBackground, fluctuate); | |
130 | if(i==0) printf(" \n\n > Generating events on the fly, please be patient ... < \n"); | |
131 | _tmpevent->Write(); | |
132 | // check if we need to open a new file, or if we've reached the end of the generator cycle | |
133 | if((i-nFileCounter*maxEventsPerFile)==(maxEventsPerFile-1) || i == (events-1)) { | |
134 | printf(" - writing memory buffer %i to output file ... \n", nFileCounter); | |
135 | dataContainer->Write(); | |
136 | dataContainer->Clear(); | |
137 | _tmpfile->Close(); | |
138 | _tmpfile->Delete(); | |
139 | _tmpfile = 0x0; | |
140 | nFileCounter++; | |
141 | } | |
142 | } | |
143 | if(qa) eventGenerator->DoGeneratorQA(kTRUE, kFALSE); | |
144 | delete eventGenerator; | |
145 | delete decayer; | |
146 | delete event; | |
147 | aa.Stop(); | |
148 | printf("\n > events have been generated and written to file ! < \n"); | |
149 | aa.Print(); | |
150 | } | |
151 | //_____________________________________________________________________________ | |
152 | void AnalyseEventsFromFile() { | |
153 | // example function: read events from file and do some analysis. | |
09c93afd | 154 | Int_t nFiles(0); // file counter |
155 | // setup analysis methods | |
90555064 | 156 | AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants(); |
157 | PrepareQCumulants(qc); | |
09c93afd | 158 | |
159 | while(kTRUE) { | |
160 | TFile file(Form("OnTheFlyEvent_%i.root",nFiles)); // open the root file | |
161 | if(file.IsZombie()) break; | |
162 | TIter iter(file.GetListOfKeys()); // get a list of keys | |
163 | ||
164 | // loop over all the events in the file | |
165 | while(kTRUE) { // infinite loop ... | |
166 | TKey* key = iter(); // get the next key from the file | |
167 | if(!key) break; // ... and exit the loop if the key is empty | |
168 | AliFlowEventSimple* event = (AliFlowEventSimple*)key->ReadObj(); // cast to a flow event | |
169 | printf(" - read event with %i tracks\n", event->NumberOfTracks()); | |
170 | // track loop | |
171 | for(int i(0); i < event->NumberOfTracks(); i++) { | |
172 | // do some offline fun stuff with the events | |
173 | AliFlowTrackSimple* track = event->GetTrack(i); | |
174 | if(!track) continue; // break the loop when we don't find a track | |
175 | // do some next stuff with the tracks | |
176 | // example: retag the primary kaons as poi's and the rest as reference | |
177 | if(track->Charge()!=-1) continue; // reject secondary tracks UGLY !!! | |
178 | if(TMath::Abs(track->GetID()==321)) { // find the kaons | |
179 | if(track->InRPSelection()) { // if the kaon is an rp (it should be by default) | |
180 | track->SetForPOISelection(kTRUE); // set it as poi ... | |
181 | track->SetForRPSelection(kFALSE); // ... and 'detag' it as rp | |
182 | // make sure to correct the reference multiplicity | |
183 | event->SetReferenceMultiplicity(-1+event->GetReferenceMultiplicity()); | |
184 | } | |
185 | else track->SetForPOISelection(kTRUE); // tag as poi | |
186 | } | |
187 | else { // if track is not a kaon | |
188 | if(!track->InRPSelection()) { // check if it's an rp | |
189 | track->SetForRPSelection(kTRUE); // if not, tag as rp (it should be there, but just to be safe ... ) | |
190 | // and again correct the reference multiplicity | |
191 | event->SetReferenceMultiplicity(1+event->GetReferenceMultiplicity()); | |
192 | } | |
193 | } | |
194 | } | |
195 | // insert the poor guy into the flow anaysis method | |
196 | qc->Make(event); | |
09c93afd | 197 | delete event; |
198 | } | |
199 | nFiles++; | |
200 | } | |
201 | // prepare output for the flow package analyses | |
202 | TString outputFileName = "AnalysisResults.root"; | |
203 | TFile *outputFile = new TFile(outputFileName.Data(),"RECREATE"); | |
90555064 | 204 | const Int_t nMethods(1); |
205 | TString method[] = {"QC"}; | |
09c93afd | 206 | TDirectoryFile *dirFileFinal[nMethods] = {NULL}; |
207 | TString fileName[nMethods]; | |
208 | for(Int_t i(0); i < nMethods; i++) { | |
209 | fileName[i]+="output"; | |
210 | fileName[i]+=method[i].Data(); | |
211 | fileName[i]+="analysis"; | |
212 | dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data()); | |
213 | } | |
90555064 | 214 | qc->Finish();qc->WriteHistograms(dirFileFinal[0]); |
09c93afd | 215 | outputFile->Close(); |
216 | delete outputFile; | |
217 | } | |
90555064 | 218 | |
219 | /* ==================================================== | |
220 | * | |
221 | * the following two functions serve as an exmaple | |
222 | * of how to run the phi and kstar analysis on on the | |
223 | * fly events | |
224 | * | |
225 | * ==================================================== | |
226 | */ | |
227 | ||
228 | //_____________________________________________________________________________ | |
229 | void TwoParticleResonanceFlowOnTheFly() { | |
230 | ||
231 | // load additional libraries | |
232 | gSystem->Load("libGeom"); | |
233 | gSystem->Load("libVMC"); | |
234 | gSystem->Load("libXMLIO"); | |
235 | gSystem->Load("libPhysics"); | |
236 | gSystem->Load("libANALYSIS"); | |
237 | gSystem->Load("libANALYSISalice"); | |
238 | gSystem->Load("libPWGflowBase"); | |
239 | gSystem->Load("libPWGflowTasks"); | |
240 | ||
241 | TString dirName = "reconstruction"; | |
242 | Int_t nFiles(0); | |
243 | // define the poi cuts, see AddTwoParticelResonanceFlowTask.C for more info | |
244 | const Int_t mb(30); // no of massbands available for analysis | |
245 | Float_t minMass(.99), maxMass(1.092); // upper and lower bound | |
246 | AliFlowTrackSimpleCuts* POIfilterQC[mb]; // pointers to poi filters | |
247 | Double_t flowBands[2][mb]; // define the invm regins | |
248 | Double_t _inc = (maxMass-minMass)/(float)mb; | |
249 | for (Int_t _mb = 0; _mb < mb; _mb++) { // create the poi filters | |
250 | flowBands[0][_mb] = minMass + _mb * _inc; | |
251 | flowBands[1][_mb] = minMass + (_mb + 1) * _inc; | |
252 | POIfilterQC[_mb] = new AliFlowTrackSimpleCuts(Form("FilterPOIQC_MB%d", _mb)); | |
253 | POIfilterQC[_mb]->SetEtaMin(-0.8); // eta range | |
254 | POIfilterQC[_mb]->SetEtaMax(+0.8); | |
255 | POIfilterQC[_mb]->SetMassMin(flowBands[0][_mb]); // invm range lower bound | |
256 | POIfilterQC[_mb]->SetMassMax(flowBands[1][_mb]); // invm rnage upper bound | |
257 | } | |
258 | // do the flow analysis | |
259 | AliFlowAnalysisWithQCumulants* qc[mb]; | |
260 | for(int i(0); i < mb; i++) { // init the q-cumulant tasks, one for each invm bin | |
261 | qc[i] = new AliFlowAnalysisWithQCumulants(); | |
262 | PrepareQCumulants(qc[i]); | |
263 | printf(" > init qc task %i < \n", i); | |
264 | } | |
265 | AliAnalysisTwoParticleResonanceFlowTask* task = new AliAnalysisTwoParticleResonanceFlowTask("onthefly"); | |
266 | SetUpTask(task, minMass,maxMass);// setup the task | |
267 | // open files | |
268 | while(kTRUE) { // infinite loop which we will break when we've looked through all the files | |
269 | TFile file(Form("OnTheFlyEvent_%i.root", nFiles)); // open the root file | |
270 | if(nFiles==0&&file.IsZombie()) { // something went wrong ... | |
271 | printf(" > cannot get input files ... exiting < \n"); | |
272 | exit(0); | |
273 | } | |
274 | if(file.IsZombie()) break; // break the loop on the first empty file | |
275 | TIter iter(file.GetListOfKeys()); // get a list of keys | |
276 | while(kTRUE) { // infinite loop over events ... | |
277 | TKey* key = iter(); // get the next key from the file | |
278 | if(!key) break; // ... and exit the loop if the key is empty | |
279 | AliFlowEventSimple* event = (AliFlowEventSimple*)key->ReadObj(); // cast to a flow event | |
280 | printf(" - read event with %i tracks from file %i \n > task ", event->NumberOfTracks(), nFiles); | |
281 | task->DoAnalysisOnTheFly(event); // do the on the fly analysis | |
282 | AliFlowEventSimple* flowEvent = task->GetFlowEvent(); // retrieve the flow event | |
283 | for(int j(0); j < mb; j++) { // loop over all invm bands | |
284 | flowEvent->TagPOI(POIfilterQC[j]); // 'offline' tagging of poi's in certain mass range | |
285 | flowEvent->TagSubeventsInEta(-.8, 0, 0, .8); // setup subevents | |
286 | qc[j]->Make(flowEvent); // do qc analysis | |
287 | printf(" %i", j); | |
288 | } | |
289 | } | |
290 | nFiles++; | |
291 | } | |
292 | // prepare output for the flow package analyses | |
293 | TFile *outputFile = new TFile("AnalysisResults.root","RECREATE"); // common outputfile | |
294 | TDirectoryFile* dirFileFinal[mb]; // tdirectory files | |
295 | TString fileName[mb]; // dir names | |
296 | for(Int_t i(0); i < mb; i++) { // loop over all bands ... | |
297 | fileName[i]+=Form("QC_minv_%i", i); | |
298 | dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data()); | |
299 | qc[i]->Finish(); // finalize the method | |
300 | qc[i]->WriteHistograms(dirFileFinal[i]); // and write it to file | |
301 | } | |
302 | // write the output of the phi reconstruction to the same file | |
303 | TDirectoryFile* dir = new TDirectoryFile(dirName.Data(), dirName.Data()); | |
304 | task->DoAnalysisOnTheFly(dir); | |
305 | // end of analysis | |
306 | outputFile->Close(); | |
307 | delete outputFile; | |
308 | } | |
309 | //_____________________________________________________________________________ | |
310 | void SetUpTask(AliAnalysisTwoParticleResonanceFlowTask* task, Float_t minMass, Float_t maxMass) | |
311 | { | |
312 | // some magic which is necessary to 'trick' the analysis task into thinking | |
313 | // that we're actually doing an analysis on aod events. | |
314 | // some of these configurations are used (e.g. setupSpecies, common constants | |
315 | // and the binning in pt) | |
316 | // others have no meaning (PID, DCA, RP and POI cuts) | |
317 | // note that UserCreateOutputObjects is necessary since it initializes | |
318 | // the output of the analysis | |
319 | gSystem->Load("libPWGflowTasks"); | |
320 | Float_t PIDconfig[] = {0,0,0,0,0,0,0.3}; // not used | |
321 | task->SetPIDConfiguration(PIDconfig); // not used | |
322 | task->SetupSpeciesA(321, 1, 4.93676999e-01, 0.15, 5.); // pid and charge | |
323 | task->SetupSpeciesB(-321, -1, 4.93676999e-01, 0.15, 5.); // pid and charge | |
324 | task->SetRPCuts(new AliFlowTrackCuts("unused_rp_cuts")); // not used | |
325 | task->SetPOICuts(new AliFlowTrackCuts("unused_poi_cuts")); // not used | |
326 | Float_t POIDCA[] = {0,0,0,0,0}; // not used | |
327 | task->SetPOIDCAXYZ(POIDCA); // not used | |
328 | task->SetCommonConstants(30, minMass, maxMass); // necesssary | |
329 | Float_t _pt[] = {0., 0.6, 1.2, 1.8, 2.4, 3., 4., 5., 6., 7.}; // same | |
330 | task->SetPtBins(_pt, (Int_t)(sizeof(_pt)/sizeof(_pt[1]))-1); // same | |
331 | task->UserCreateOutputObjects(); // necessary | |
332 | } | |
333 | //_____________________________________________________________________________ | |
334 | void PrepareQCumulants(AliFlowAnalysisWithQCumulants* qc) | |
335 | { | |
336 | // prepare the cumulants task | |
337 | qc->SetHarmonic(2); | |
338 | qc->SetCalculateDiffFlow(kTRUE); | |
339 | qc->SetCalculate2DDiffFlow(kFALSE); // vs (pt,eta) | |
340 | qc->SetApplyCorrectionForNUA(kFALSE); | |
341 | qc->SetFillMultipleControlHistograms(kFALSE); | |
342 | qc->SetMultiplicityWeight("combinations"); // default (other supported options are "unit" and "multiplicity") | |
343 | qc->SetCalculateCumulantsVsM(kFALSE); | |
344 | qc->SetCalculateAllCorrelationsVsM(kFALSE); // calculate all correlations in mixed harmonics "vs M" | |
345 | qc->SetnBinsMult(10000); | |
346 | qc->SetMinMult(0); | |
347 | qc->SetMaxMult(10000); | |
348 | qc->SetBookOnlyBasicCCH(kFALSE); // book only basic common control histograms | |
349 | qc->SetCalculateDiffFlowVsEta(kTRUE); // if you set kFALSE only differential flow vs pt is calculated | |
350 | qc->Init(); | |
351 | } | |
09c93afd | 352 | //_____________________________________________________________________________ |
353 | Bool_t LoadLibraries() { | |
354 | // load libraries | |
355 | gSystem->SetIncludePath("-I. -I$ROOTSYS/include -I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TPC -I$ALICE_ROOT/CONTAINERS -I$ALICE_ROOT/STEER -I$ALICE_ROOT/TRD -I$ALICE_ROOT/macros -I$ALICE_ROOT/ANALYSIS -I$ALICE_ROOT/OADB -I$ALICE_ROOT/PWGHF -I$ALICE_ROOT/PWGHF/base -I$ALICE_ROOT/PWGHF/vertexingHF -I$ALICE_ROOT/PWG/FLOW/Base -I$ALICE_ROOT/PWG/FLOW/Tasks -g"); | |
356 | gSystem->Load("libGeom"); | |
357 | gSystem->Load("libVMC"); | |
358 | gSystem->Load("libXMLIO"); | |
359 | gSystem->Load("libPhysics"); | |
360 | gSystem->Load("libEG"); | |
361 | gSystem->Load("libPWGflowBase"); | |
362 | if(gSystem->Load("$ALICE_ROOT/lib/tgt_linuxx8664gcc/libpythia6")!=0) { | |
363 | printf(" \n\n\n *** fatal error, couldn't load pythia, exiting ! *** \n\n\n"); | |
364 | return kFALSE; | |
365 | } | |
366 | return kTRUE; | |
367 | } | |
368 | //_____________________________________________________________________________ |