]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/FLOW/Base/AliFlowOnTheFlyEventGenerator.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWG / FLOW / Base / AliFlowOnTheFlyEventGenerator.cxx
CommitLineData
09c93afd 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/*****************************************************************
17 AliFlowOnTheFlyEventGenerator
18
19 Create events with realstic spectra and v_n values
20 store as AliFlowEventSimple
21
22 usage: see $ALICE_ROOT/PWGCF/FLOW/macros/GenerateOnTheFlyEvents.C
23
24 origin: Redmer Alexander Bertens, rbertens@cern.ch, rbertens@nikhef.nl, rbertens@uu.nl
25 Carlos Eugenio Perez Lara
26 Andrea Dubla
27******************************************************************/
28
29#include "Riostream.h"
30#include "TObjArray.h"
31#include "TClonesArray.h"
32#include "TFile.h"
33#include "TList.h"
34#include "TTree.h"
35#include "TParticle.h"
36#include "TMath.h"
37#include "TH1.h"
38#include "TH1F.h"
39#include "TH1D.h"
40#include "TF1.h"
41#include "TString.h"
42#include "TProfile.h"
43#include "TParameter.h"
44#include "TBrowser.h"
45#include "AliFlowTrackSimple.h"
46#include "AliFlowEventSimple.h"
47#include "TRandom.h"
48#include "TVirtualMCDecayer.h"
49#include "AliFlowOnTheFlyEventGenerator.h"
50
51using std::endl;
52using std::cout;
53
54static Int_t OnTheFlyEventGeneratorCounter = 0;
55
56ClassImp(AliFlowOnTheFlyEventGenerator)
57//_____________________________________________________________________________
1164b0ff 58AliFlowOnTheFlyEventGenerator::AliFlowOnTheFlyEventGenerator() : fGenerators(0), fEmbedMe(0), fFlowEvent(0), fDecayer(0), fAddV2Mothers(0), fAddV3Mothers(0), fAddV2Daughters(0), fAddV3Daughters(0), fPsi2(0), fPsi3(0), fPrecisionPhi(.001), fMaxNumberOfIterations(100), fQA(0), fFF(0) {/* constructor used for IO by root */}
09c93afd 59//_____________________________________________________________________________
1d89e44d 60AliFlowOnTheFlyEventGenerator::AliFlowOnTheFlyEventGenerator(Bool_t qa, Int_t ff, Int_t mult, TVirtualMCDecayer* decayer, Bool_t a, Bool_t b, Bool_t c, Bool_t d) : fGenerators(0), fEmbedMe(0), fFlowEvent(0), fDecayer(0), fAddV2Mothers(0), fAddV3Mothers(0), fAddV2Daughters(0), fAddV3Daughters(0), fPsi2(0), fPsi3(0), fPrecisionPhi(.001), fMaxNumberOfIterations(100), fQA(qa), fFF(ff) {
09c93afd 61 // contructor
62 if(!fFlowEvent) fFlowEvent = new AliFlowEventSimple(mult, (AliFlowEventSimple::ConstructionMethod)0, 0x0, 0, TMath::TwoPi(), -.9, .9);
63 if(!fGenerators) {
64 fGenerators = new TObjArray();
65 fGenerators->SetOwner(kTRUE);
66 InitGenerators();
67 }
68 fDecayer = decayer; // decayer: user has ownership of decayer (see dtor)
ca33e334 69 if(fDecayer) fDecayer->Init();
09c93afd 70 fAddV2Mothers = a;
71 fAddV3Mothers = b;
72 fAddV2Daughters = c;
73 fAddV3Daughters = d;
74 if(fAddV3Mothers||fAddV3Daughters) printf(" > WARNING - v3 is not impelmented yet ! \n <");
75}
76//_____________________________________________________________________________
77AliFlowOnTheFlyEventGenerator::~AliFlowOnTheFlyEventGenerator()
78{
79 // destructor
80 if(fGenerators) delete fGenerators;
81 if(fFlowEvent) delete fFlowEvent;
82}
83//_____________________________________________________________________________
84AliFlowOnTheFlyEventGenerator::NaiveFlowAndSpectrumGenerator* AliFlowOnTheFlyEventGenerator::Find(Short_t pdg, Bool_t make)
85{
86 // return instance of generator class for given fPDG value
87 for(int i(0); i < fGenerators->GetEntriesFast(); i++) {
88 if(((NaiveFlowAndSpectrumGenerator*)fGenerators->At(i))->GetPDGCode()==pdg) {
89 return (NaiveFlowAndSpectrumGenerator*)(fGenerators->At(i));
90 }
91 }
92 if(make) {
1164b0ff 93 AliFlowOnTheFlyEventGenerator::NaiveFlowAndSpectrumGenerator* _tmp = new AliFlowOnTheFlyEventGenerator::NaiveFlowAndSpectrumGenerator(pdg, fQA, fFF);
09c93afd 94 fGenerators->Add(_tmp);
95 return _tmp;
96 }
97 return 0x0;
98}
99//_____________________________________________________________________________
100void AliFlowOnTheFlyEventGenerator::SetPtSpectrum(const char* func, Short_t species)
101{
102 // set pt spectrum for species, override default
103 NaiveFlowAndSpectrumGenerator* gen = Find(species, kTRUE);
104 TF1* spectrum = gen->GetPtSpectrum();
105 TString name = "";
106 if(spectrum) {
107 name = spectrum->GetName();
108 delete spectrum;
109 }
110 else name = Form("pt_%i", (int)species);
111 gen->SetPtSpectrum(new TF1(name.Data(), func, 0., 20));
112}
113//_____________________________________________________________________________
114void AliFlowOnTheFlyEventGenerator::SetPtDependentV2(const char* func, Short_t species)
115{
116 // set pt dependent v2 for species, override default
117 NaiveFlowAndSpectrumGenerator* gen = Find(species, kTRUE);
118 TF1* v2 = gen->GetDifferentialV2();
119 TString name = "";
120 if(v2) { // if v2 is already present, get its name and release its memory
121 name = v2->GetName();
122 delete v2;
123 }
124 else name = Form("v2_%i", (int)species);
125 gen->SetDifferentialV2(new TF1(name.Data(),func, 0., 20));
126}
127//_____________________________________________________________________________
128void AliFlowOnTheFlyEventGenerator::SetPtDependentV3(const char* func, Short_t species)
129{
130 // set pt dependent v2 for species, override default
131 NaiveFlowAndSpectrumGenerator* gen = Find(species, kTRUE);
132 TF1* v3 = gen->GetDifferentialV3();
133 TString name = "";
134 if(v3) { // if v3 is already present, get its name and release its memory
135 name = v3->GetName();
136 delete v3;
137 }
138 else name = Form("v3_%i", (int)species);
139 gen->SetDifferentialV3(new TF1(name.Data(), func, 0., 20));
140}
141//_____________________________________________________________________________
142AliFlowEventSimple* AliFlowOnTheFlyEventGenerator::GenerateOnTheFlyEvent(TClonesArray* event, Int_t nSpecies, Int_t species[], Int_t mult[], Int_t bg, Bool_t fluc)
143{
144 // generate a flow event according to settings provided in GENERATOR INTERFACE
09c93afd 145 fPsi2 = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
1d89e44d 146 // generate a random number which will be used for flow fuctuations (if requested)
147 Double_t fluc_poi(-1), fluc_rp(-1);
148 if(fFF == 1 || fFF == 3 ) fluc_poi = gRandom->Uniform();
149 if(fFF == 2) fluc_rp = gRandom->Uniform();
150 if(fFF == 3) fluc_rp = fluc_poi; // poi and rp fluctuations are fully correlated
09c93afd 151 // calculate the total multiplicity for the event
152 Int_t multPiKPr[] = {(int)(.7*bg/2.), (int)(.2*bg/2.), (int)(.1*bg/2.)};
153 Int_t fPDGPiKPr[] = {211, 321, 2212};
154 Int_t totalRPMultiplicity(0), totalPOIMultiplicity(0);
155 for(int i(0); i < nSpecies; i++) totalPOIMultiplicity+=mult[i];
156 for(int i(0); i < 3; i++) totalRPMultiplicity+=multPiKPr[i];
157 Int_t totalMultiplicity(totalRPMultiplicity+totalPOIMultiplicity);
158 // generate the particles of interest. if requested, a vn boost is given to the primary particles
159 for(int i(0); i < nSpecies; i++) {
8d1c89ad 160 if(fluc) GenerateOnTheFlyTracks((Int_t)(gRandom->Uniform(mult[i]-TMath::Sqrt(mult[i]), mult[i]+TMath::Sqrt(mult[i]))), species[i], event, fluc_poi);
1d89e44d 161 else GenerateOnTheFlyTracks(mult[i], species[i], event, fluc_poi);
09c93afd 162 }
163 // generate reference particles. if requested, a vn boost is given to the primary particles
164 for(int i(0); i < 3; i++) {
165 if(fluc) {
8d1c89ad 166 GenerateOnTheFlyTracks((Int_t)(gRandom->Uniform(multPiKPr[i]-TMath::Sqrt(multPiKPr[i]), multPiKPr[i]+TMath::Sqrt(mult[i]))), fPDGPiKPr[i], event, fluc_rp);
167 GenerateOnTheFlyTracks((Int_t)(gRandom->Uniform(multPiKPr[i]-TMath::Sqrt(multPiKPr[i]), multPiKPr[i]+TMath::Sqrt(mult[i]))), -1*fPDGPiKPr[i], event, fluc_rp);
09c93afd 168 }
169 else {
1d89e44d 170 GenerateOnTheFlyTracks(multPiKPr[i], fPDGPiKPr[i], event, fluc_rp);
171 GenerateOnTheFlyTracks(multPiKPr[i], -1*fPDGPiKPr[i], event, fluc_rp);
09c93afd 172 }
173 }
174 // if requested, decay the event
175 if(fDecayer) DecayOnTheFlyTracks(event);
176 // if an embedded event is given to the generator at this point, embed it
177 if(fEmbedMe) {
178 event->AbsorbObjects(fEmbedMe); // event has ownership of embedded event
179 fEmbedMe = 0x0; // reset to aviod embeddding more than once
180 }
181 // if requested, the secondaries are given a vn boost
182 if(fAddV2Daughters) AddV2(event);
183 // convert the event to an aliflowsimple event
184 // the tagging (rp) is done in this step, all tracks are tagged as rp's.
185 // daughters are made orphans (i.e. primaries with secondaries are rejected from the sample)
186 // converting clears the event tclones array
187 return ConvertTClonesToFlowEvent(event, totalMultiplicity);
188}
189//_____________________________________________________________________________
190AliFlowEventSimple* AliFlowOnTheFlyEventGenerator::ConvertTClonesToFlowEvent(TClonesArray* event, Int_t totalMultiplicity)
191{
192 // convert the content of a tclones array to an aliflowevent simple for output and tag particles as poi's and rp's
193 // first step, see if there's already a flow event available in memory. if not create one
194 if(!fFlowEvent) fFlowEvent = new AliFlowEventSimple(totalMultiplicity, (AliFlowEventSimple::ConstructionMethod)0, 0x0, 0, TMath::TwoPi(), -.9, .9);
195 // if there is a flow event, clear the members without releasing the allocated memory
196 else fFlowEvent->ClearFast();
26f043a1 197 fFlowEvent->SetMCReactionPlaneAngle(fPsi2);
09c93afd 198 // prepare a track
199 TParticle* pParticle = 0x0;
200 Int_t iSelParticlesRP(0);
201 for (Int_t i(0); i < event->GetEntries(); i++) { // begin the loop over the input tracks
202 pParticle = (TParticle*)event->At(i); // get the particle
203 if (!pParticle) continue; // skip if empty slot (no particle)
204 if (pParticle->GetNDaughters()!=0) continue; // see if the particle has daughters (if so, reject it)
205 AliFlowTrackSimple* pTrack = new AliFlowTrackSimple(pParticle); // allocate space for the new flow track
206 pTrack->SetWeight(pParticle->Pz()); // ugly hack: store pz here ...
207 pTrack->SetID(pParticle->GetPdgCode()); // set pid code as id
208 pTrack->SetForRPSelection(kTRUE); // tag ALL particles as RP's,
209 // ugly hack 2: set charge to -1 for primaries and 1 for secondaries
210 if(pParticle->GetFirstMother()==-1) pTrack->SetCharge(-1);
211 else pTrack->SetCharge(1);
212 iSelParticlesRP++;
213 fFlowEvent->AddTrack(pTrack);
214 }
215 fFlowEvent->SetNumberOfRPs(iSelParticlesRP);
216 // all trakcs have safely been copied so we can clear the event
217 event->Clear();
218 OnTheFlyEventGeneratorCounter = 0;
219 return fFlowEvent;
220}
221//_____________________________________________________________________________
1d89e44d 222void AliFlowOnTheFlyEventGenerator::AddV2(TParticle* part, Double_t v2, Double_t fluc)
09c93afd 223{
224 // afterburner, adds v2, uses Newton-Raphson iteration
225 Double_t phi = part->Phi();
226 Double_t phi0=phi;
227 Double_t f=0.;
228 Double_t fp=0.;
229 Double_t phiprev=0.;
1d89e44d 230 // introduce flow fluctuations (gaussian)
231 if(fluc > -.5) {
232 v2+=TMath::Sqrt(2*(v2*.25)*(v2*.25))*TMath::ErfInverse(2*fluc-1);
233 if(fQA) Find(part->GetPdgCode(), kTRUE)->FillV2(part->Pt(), v2);
234 }
09c93afd 235 for (Int_t i=0; i!=fMaxNumberOfIterations; ++i) {
236 phiprev=phi; //store last value for comparison
237 f = phi-phi0+v2*TMath::Sin(2.*(phi-fPsi2));
238 fp = 1.0+2.0*v2*TMath::Cos(2.*(phi-fPsi2)); //first derivative
239 phi -= f/fp;
240 if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break;
241 }
242 part->SetMomentum( part->Pt()*TMath::Cos(phi), part->Pt()*TMath::Sin(phi), part->Pz(), part->Energy() );
243}
244//_____________________________________________________________________________
245void AliFlowOnTheFlyEventGenerator::AddV2(TClonesArray* event)
246{
247 // afterburner, adds v2 for different species to all tracks in an event
1d89e44d 248 Double_t fluc(gRandom->Uniform()); // get a random number in case of fluctuations
09c93afd 249 // FIXME at the moment no distincition between mothers and daughters
5ac1a8e1 250 TParticle *part;
09c93afd 251 for(Int_t nTrack=0; nTrack!=event->GetEntriesFast(); ++nTrack) {
252 part = (TParticle*) event->At(nTrack);
253 // for each track, call overloaded addv2 function with the correct generator
254 // create a generator in the case where the decayer has introduced a new particle species
5ac1a8e1 255 if(part)
256 AddV2(part, Find(part->GetPdgCode(), kTRUE)->GetV2(part->Pt()), fluc);
09c93afd 257 }
258}
259//_____________________________________________________________________________
1d89e44d 260void AliFlowOnTheFlyEventGenerator::GenerateOnTheFlyTracks(Int_t mult, Int_t pid, TClonesArray* event, Double_t fluc)
09c93afd 261{
262 // generate tracks for the event. the tracks that are generated here are the mothers
263 Double_t ph, et, pt, mm, px, py, pz, ee;
264 // access the descired generator through the static find method. if the particle species is unknown, we add it to the generator
265 NaiveFlowAndSpectrumGenerator* generator = Find(pid, kTRUE);
266 Int_t _tempCounter = OnTheFlyEventGeneratorCounter;
267 for(int i=_tempCounter; i<=mult+_tempCounter; ++i) {
268 OnTheFlyEventGeneratorCounter++;
269 TParticle* part = (TParticle*)event->ConstructedAt(i);
270 part->SetStatusCode(1);
271 part->SetFirstMother(-1);
272 part->SetLastMother(-1);
273 part->SetFirstDaughter(-1);
274 part->SetLastDaughter(-1);
275 ph = gRandom->Uniform(0,TMath::TwoPi());
276 et = gRandom->Uniform(-1.0,+1.0);
277 pt = generator->GetPt();
278 part->SetPdgCode( pid );
279 mm = part->GetMass();
280 px = pt*TMath::Cos(ph);
281 py = pt*TMath::Sin(ph);
282 pz = pt*TMath::SinH(et);
283 ee = TMath::Sqrt( mm*mm + pt*pt+pz*pz );
284 part->SetMomentum( px, py, pz, ee );
285 // if requested add v2 to the sample of mothers
1d89e44d 286 if(fAddV2Mothers) AddV2(part, generator->GetV2(pt), fluc);
09c93afd 287 }
288}
289//_____________________________________________________________________________
290void AliFlowOnTheFlyEventGenerator::DecayOnTheFlyTracks(TClonesArray *event)
291{
292 // decay the tracks using a decayer that is set in the GenerateEventsOnTheFly.C macro
293 if(!fDecayer) return; // shouldn't happen ...
09c93afd 294 Int_t kf;
295 TClonesArray* arr = new TClonesArray("TParticle",10);
296 Int_t secondaries=0;
297 Int_t nStart=0;
298 Int_t nTracks = event->GetEntriesFast();
299 TParticle *part, *part1;
300 TClonesArray &in = *event;
301 for(Int_t loop=0; loop!=3; ++loop) {
302 secondaries=0;
303 for(Int_t nTrack=nStart; nTrack!=nTracks; ++nTrack) {
304 arr->Clear();
305 part = (TParticle*) event->At( nTrack );
306 if(!part) continue;
307 kf = TMath::Abs(part->GetPdgCode());
308 if(kf==22 ) ForceGammaDecay(arr, part); // if gamma, we decay it by hand
309 else { // else we let pythia do it
310 TLorentzVector pmom( part->Px(), part->Py(), part->Pz(), part->Energy() );
311 fDecayer->Decay(kf,&pmom);
312 fDecayer->ImportParticles(arr);
313 }
314 Int_t nDaughters = arr->GetEntries();
315 if(nDaughters<2) continue;
316 for(Int_t nDaughter=1; nDaughter!=nDaughters; ++nDaughter) {
317 part1 = (TParticle*) (arr->At(nDaughter));
318 if(!part1) continue; // was part1, mistake?
319 new(in[nTracks+secondaries]) TParticle( part1->GetPdgCode(),
320 part1->GetStatusCode(),
321 part1->GetFirstMother(),
322 part1->GetSecondMother(),
323 part1->GetFirstDaughter(),
324 part1->GetLastDaughter(),
325 part1->Px(),
326 part1->Py(),
327 part1->Pz(),
328 part1->Energy(),
329 part1->Vx(),
330 part1->Vy(),
331 part1->Vz(),
332 part1->T());
333 secondaries++;
334 if(nDaughter==1) part->SetFirstDaughter(nTracks+secondaries);
335 else if ((nDaughters-1)==nDaughter) part->SetLastDaughter(nTracks+secondaries);
5ac1a8e1 336 //else part->SetDaughter(nDaughter,nTracks+secondaries);
26f043a1 337 }
09c93afd 338 }
09c93afd 339 nStart = nTracks;
340 nTracks = event->GetEntries();
341 }
342 delete arr;
343 return;
344}
345//_____________________________________________________________________________
346void AliFlowOnTheFlyEventGenerator::ForceGammaDecay(TClonesArray* arr, TParticle* part)
347{
348 // force gama decay
1d89e44d 349 const Float_t mass1(gRandom->Uniform(0.0000001, .1)); // FIXME randomly sample photon 'mass'
09c93afd 350 const Float_t mass_d1(.000511); // daughter mass (electron)
351 const Float_t mass_d2(.000511);
352 Double_t p_dec1=sqrt((TMath::Abs(mass1*mass1-(mass_d1+mass_d2)*(mass_d1+mass_d2))*(mass1*mass1-(mass_d1-mass_d2)*(mass_d1-mass_d2))))/2/mass1; // energy
353 TLorentzVector *p_d1=new TLorentzVector(); // lorentz vectors for daughters
354 TLorentzVector *p_d2=new TLorentzVector();
355 Double_t yy_m=part->Y(); // pseudo rapidity of mother
356 Double_t pt_m=part->Pt(); // pt of mother
357 Double_t mt_m=sqrt(mass1*mass1+pt_m*pt_m); // transverse mass of mother
358 Double_t pz_m=mt_m*TMath::SinH(yy_m); // MAGIC
359 Double_t phi_m=gRandom->Rndm()*2*TMath::Pi();
360 Double_t px_m=pt_m*TMath::Cos(phi_m);
361 Double_t py_m=pt_m*TMath::Sin(phi_m);
362 Double_t costh_d=2*gRandom->Rndm()-1;
363 Double_t phi_d=gRandom->Rndm()*2*TMath::Pi();
364 Double_t pz_d=p_dec1*costh_d;
365 Double_t pt_d=p_dec1*sqrt(1-costh_d*costh_d);
366 Double_t px_d=pt_d*TMath::Cos(phi_d);
367 Double_t py_d=pt_d*TMath::Sin(phi_d);
368 p_d1->SetPxPyPzE(px_d,py_d,pz_d,sqrt(mass_d1*mass_d1+p_dec1*p_dec1));
369 p_d2->SetPxPyPzE(-px_d,-py_d,-pz_d,sqrt(mass_d2*mass_d2+p_dec1*p_dec1));
370 Double_t gamma_b=sqrt(mass1*mass1+pz_m*pz_m+pt_m*pt_m)/mass1;
371 Double_t bx=px_m/gamma_b/mass1;
372 Double_t by=py_m/gamma_b/mass1;
373 Double_t bz=pz_m/gamma_b/mass1;
374 p_d1->Boost(bx,by,bz);
375 p_d2->Boost(bx,by,bz);
376 TParticle* daughter_1 = (TParticle*)arr->ConstructedAt(0);
377 TParticle* daughter_2 = (TParticle*)arr->ConstructedAt(1);
378 daughter_1->SetStatusCode(1);
379 daughter_1->SetFirstMother(-1);
380 daughter_1->SetLastMother(-1);
381 daughter_1->SetFirstDaughter(-1);
382 daughter_1->SetLastDaughter(-1);
383 daughter_1->SetPdgCode(11);
384 TLorentzVector& ref_p_d1 = *p_d1;
385 daughter_1->SetMomentum(ref_p_d1);
386 daughter_2->SetStatusCode(1);
387 daughter_2->SetFirstMother(-1);
388 daughter_2->SetLastMother(-1);
389 daughter_2->SetFirstDaughter(-1);
390 daughter_2->SetLastDaughter(-1);
391 daughter_2->SetPdgCode(-11);
392 TLorentzVector& pp = *p_d2;
393 daughter_2->SetMomentum(pp);
394}
395//_____________________________________________________________________________
396void AliFlowOnTheFlyEventGenerator::InitGenerators()
397{
398 // initializes generators for a list of common particles
399 // new generators will be created 'on the fly' if they are encountered
400 // if a generator is requested for a non-existent pdg value, things will get unstable ...
401 Short_t PDGcode[] = {11, -11, // e-
402 12, // v_e
403 13, -13, // mu-
404 14, // v_mu
405 15, -15, // t-
406 16, // v_t
407 21, // g
408 22, // gamma
409 111, // pi0
410 211, -211, // pi+
411 113, // rho0
412 213, -213, // rho+
413 221, // eta
414 331, // eta prime
415 130, // k0l
416 310, // k0s
417 311, // k0
418 313, // k*
419 321, -321, // k+
420 323, -323, // k*+
421 223, // omega
422 411, -411, // d+
423 413, -413, // d*+
424 421, // d0
425 423, // d*0
426 333, // phi
427 443, // jpsi
428 2112, // neutron
429 2212, -2212, // proton
430 3122, // lambda0
431 3312, -3312, // xi- (cascade)
432 3314, -3314, // xi*-
433 3322, // xi0
434 3324, // xi*0
435 3334}; // Omeg
1164b0ff 436 for(int i(0); i < 47; i++) fGenerators->Add(new NaiveFlowAndSpectrumGenerator(PDGcode[i], fQA, fFF));
09c93afd 437}
438//_____________________________________________________________________________
439void AliFlowOnTheFlyEventGenerator::PrintGenerators()
440{
441 // make sure all went ok
442 if(!fGenerators) {
443 printf(" > PANIC <\n");
444 return;
445 }
446 printf(" > %i species available in generator\n", fGenerators->GetEntriesFast());
447 for(int i(0); i < fGenerators->GetEntriesFast(); i++)
448 printf(" - PDG: %i \n", ((NaiveFlowAndSpectrumGenerator*)fGenerators->At(i))->GetPDGCode());
449 printf(" more can be created 'on the fly', but beware of non-existent pdg values !");
450}
451//_____________________________________________________________________________
452void AliFlowOnTheFlyEventGenerator::DoGeneratorQA(Bool_t v2, Bool_t v3)
453{
454 // this function loops over all the species that are available in the fGenerators, which are
455 // a) created by the user (in the event generator) as mother particles or
456 // b) introduced by pythia (the decayer) as daughter particles
457 // c) made by the InitGenerators() function
458 // 'empty species' (i.e. species for which a generator exists but which were not actually sampled) are omitted
459 // the output histograms are written to a file named "GeneratorfQA.root"
460 if(!fQA) {
461 printf(" > Request has been made for QA plots but QA histograms have been disabled !\n");
462 return;
463 }
464 TFile *QAfile = new TFile("GeneratorfQA.root","RECREATE");
465 for(int i(0); i < fGenerators->GetEntriesFast(); i++) {
466 NaiveFlowAndSpectrumGenerator* gen = (NaiveFlowAndSpectrumGenerator*)fGenerators->At(i);
467 if(!gen) continue; // skip if we failed to get a generator
468 TH1F* QApt = (TH1F*)gen->GetQAType(0);
469 TH2F* QAv2 = (TH2F*)gen->GetQAType(1);
470 TH2F* QAv3 = (TH2F*)gen->GetQAType(2);
471 if((!QApt)||(!QAv2)||(!QAv3)) {
472 printf(" > couldn't read qa histogrmas for species %i <\n", gen->GetPDGCode());
473 continue;
474 }
475 if(QApt->GetEntries()==0&&QAv2->GetEntries()==0&&QAv3->GetEntries()==0) continue; // skip if no tracks of this species have been sampled
476 printf(" > saving QA histograms for sampled species %i <\n", gen->GetPDGCode());
477 QAfile->mkdir(Form("fPDG_%i", gen->GetPDGCode())); // create a directory in the output file
478 QAfile->cd(Form("fPDG_%i", gen->GetPDGCode())); // cd into this directory
479 if(!QApt->GetEntries()==0) { // only plot the pt fSpectrum if this guy was generated
480 if(!gen->GetPtSpectrum()->Integral(0,20)==0) {
481 QApt->Scale(gen->GetPtSpectrum()->Integral(0,20)/QApt->Integral(),"width");
482 QApt->Write();
483 }
484 }
485 gen->GetPtSpectrum()->SetNpx(10000); // otherwise tf1 plot is very ugly
486 gen->GetPtSpectrum()->Draw();
487 gen->GetPtSpectrum()->Write();
488 if(v2) {
489 QAv2->Draw();
490 gen->GetDifferentialV2()->SetNpx(10000);
491 gen->GetDifferentialV2()->Draw();
492 gen->GetDifferentialV2()->Write();
493 QAv2->Write();
494 }
495 if(v3) {
496 QAv3->Draw();
497 gen->GetDifferentialV3()->Draw();
498 gen->GetDifferentialV3()->SetNpx(10000);
499 gen->GetDifferentialV3()->Write();
500 QAv3->Write();
501 }
502 }
503 QAfile->Close();
504}
505//_____________________________________________________________________________
506