PWG/FLOW/Base/AliFlowOnTheFlyEventGenerator.cxx

   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
  51 using std::endl;
  52 using std::cout;
  53
  54 static Int_t OnTheFlyEventGeneratorCounter = 0;
  55
  56 ClassImp(AliFlowOnTheFlyEventGenerator)
  57 //_____________________________________________________________________________
  58 AliFlowOnTheFlyEventGenerator::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 */}
  59 //_____________________________________________________________________________
  60 AliFlowOnTheFlyEventGenerator::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) {
  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)
  69     fAddV2Mothers = a;
  70     fAddV3Mothers = b;
  71     fAddV2Daughters = c;
  72     fAddV3Daughters = d;
  73     if(fAddV3Mothers||fAddV3Daughters) printf(" > WARNING - v3 is not impelmented yet ! \n <");
  74 }
  75 //_____________________________________________________________________________
  76 AliFlowOnTheFlyEventGenerator::~AliFlowOnTheFlyEventGenerator()
  77 {
  78     // destructor
  79     if(fGenerators)     delete fGenerators;
  80     if(fFlowEvent)      delete fFlowEvent;
  81 }
  82 //_____________________________________________________________________________
  83 AliFlowOnTheFlyEventGenerator::NaiveFlowAndSpectrumGenerator* AliFlowOnTheFlyEventGenerator::Find(Short_t pdg, Bool_t make)
  84 {
  85     // return instance of generator class for given fPDG value
  86     for(int i(0); i < fGenerators->GetEntriesFast(); i++) {
  87         if(((NaiveFlowAndSpectrumGenerator*)fGenerators->At(i))->GetPDGCode()==pdg) {
  88             return (NaiveFlowAndSpectrumGenerator*)(fGenerators->At(i));
  89         }
  90     }
  91     if(make) {
  92         AliFlowOnTheFlyEventGenerator::NaiveFlowAndSpectrumGenerator* _tmp = new AliFlowOnTheFlyEventGenerator::NaiveFlowAndSpectrumGenerator(pdg, fQA, fFF);
  93         fGenerators->Add(_tmp);
  94         return _tmp;
  95     }
  96     return 0x0;
  97 }
  98 //_____________________________________________________________________________
  99 void AliFlowOnTheFlyEventGenerator::SetPtSpectrum(const char* func, Short_t species)
 100 {
 101     // set pt spectrum for species, override default
 102     NaiveFlowAndSpectrumGenerator* gen = Find(species, kTRUE);
 103     TF1* spectrum = gen->GetPtSpectrum();
 104     TString name = "";
 105     if(spectrum) {
 106         name = spectrum->GetName();
 107         delete spectrum;
 108     }
 109     else name = Form("pt_%i", (int)species);
 110     gen->SetPtSpectrum(new TF1(name.Data(), func, 0., 20));
 111 }
 112 //_____________________________________________________________________________
 113 void AliFlowOnTheFlyEventGenerator::SetPtDependentV2(const char* func, Short_t species)
 114 {
 115     // set pt dependent v2 for species, override default
 116     NaiveFlowAndSpectrumGenerator* gen = Find(species, kTRUE);
 117     TF1* v2 = gen->GetDifferentialV2();
 118     TString name = "";
 119     if(v2) { // if v2 is already present, get its name and release its memory
 120         name = v2->GetName();
 121         delete v2;
 122     }
 123     else name = Form("v2_%i", (int)species);
 124     gen->SetDifferentialV2(new TF1(name.Data(),func, 0., 20));
 125 }
 126 //_____________________________________________________________________________
 127 void AliFlowOnTheFlyEventGenerator::SetPtDependentV3(const char* func, Short_t species)
 128 {
 129     // set pt dependent v2 for species, override default
 130     NaiveFlowAndSpectrumGenerator* gen = Find(species, kTRUE);
 131     TF1* v3 = gen->GetDifferentialV3();
 132     TString name = "";
 133     if(v3) { // if v3 is already present, get its name and release its memory
 134         name = v3->GetName();
 135         delete v3;
 136     }
 137     else name = Form("v3_%i", (int)species);
 138     gen->SetDifferentialV3(new TF1(name.Data(), func, 0., 20));
 139 }
 140 //_____________________________________________________________________________
 141 AliFlowEventSimple* AliFlowOnTheFlyEventGenerator::GenerateOnTheFlyEvent(TClonesArray* event, Int_t nSpecies, Int_t species[], Int_t mult[], Int_t bg, Bool_t fluc)
 142 {
 143     // generate a flow event according to settings provided in GENERATOR INTERFACE
 144     fPsi2 = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
 145     // generate a random number which will be used for flow fuctuations (if requested)
 146     Double_t fluc_poi(-1), fluc_rp(-1);
 147     if(fFF == 1 || fFF == 3 ) fluc_poi = gRandom->Uniform();
 148     if(fFF == 2) fluc_rp = gRandom->Uniform();
 149     if(fFF == 3) fluc_rp = fluc_poi;    // poi and rp fluctuations are fully correlated
 150     // calculate the total multiplicity for the event
 151     Int_t multPiKPr[] = {(int)(.7*bg/2.), (int)(.2*bg/2.), (int)(.1*bg/2.)};
 152     Int_t fPDGPiKPr[] = {211, 321, 2212};
 153     Int_t totalRPMultiplicity(0), totalPOIMultiplicity(0);
 154     for(int i(0); i < nSpecies; i++) totalPOIMultiplicity+=mult[i];
 155     for(int i(0); i < 3; i++)        totalRPMultiplicity+=multPiKPr[i];
 156     Int_t totalMultiplicity(totalRPMultiplicity+totalPOIMultiplicity);
 157     // generate the particles of interest. if requested, a vn boost is given to the primary particles
 158     for(int i(0); i < nSpecies; i++) {
 159         if(fluc) GenerateOnTheFlyTracks(gRandom->Uniform(mult[i]-TMath::Sqrt(mult[i]), mult[i]+TMath::Sqrt(mult[i])), species[i], event, fluc_poi);
 160         else GenerateOnTheFlyTracks(mult[i], species[i], event, fluc_poi);
 161     }
 162     // generate reference particles. if requested, a vn boost is given to the primary particles
 163     for(int i(0); i < 3; i++) {
 164         if(fluc) {
 165             GenerateOnTheFlyTracks(gRandom->Uniform(multPiKPr[i]-TMath::Sqrt(multPiKPr[i]), multPiKPr[i]+TMath::Sqrt(mult[i])), fPDGPiKPr[i], event, fluc_rp);
 166             GenerateOnTheFlyTracks(gRandom->Uniform(multPiKPr[i]-TMath::Sqrt(multPiKPr[i]), multPiKPr[i]+TMath::Sqrt(mult[i])), -1*fPDGPiKPr[i], event, fluc_rp);
 167         }
 168         else {
 169             GenerateOnTheFlyTracks(multPiKPr[i], fPDGPiKPr[i], event, fluc_rp);
 170             GenerateOnTheFlyTracks(multPiKPr[i], -1*fPDGPiKPr[i], event, fluc_rp);
 171         }
 172     }
 173     // if requested, decay the event
 174     if(fDecayer)        DecayOnTheFlyTracks(event);
 175     // if an embedded event is given to the generator at this point, embed it
 176     if(fEmbedMe) {
 177         event->AbsorbObjects(fEmbedMe); // event has ownership of embedded event
 178         fEmbedMe = 0x0;                 // reset to aviod embeddding more than once
 179     }
 180     // if requested, the secondaries are given a vn boost
 181     if(fAddV2Daughters) AddV2(event);
 182     // convert the event to an aliflowsimple event
 183     // the tagging (rp) is done in this step, all tracks are tagged as rp's.
 184     // daughters are made orphans (i.e. primaries with secondaries are rejected from the sample)
 185     // converting clears the event tclones array
 186     return ConvertTClonesToFlowEvent(event, totalMultiplicity);
 187 }
 188 //_____________________________________________________________________________
 189 AliFlowEventSimple* AliFlowOnTheFlyEventGenerator::ConvertTClonesToFlowEvent(TClonesArray* event, Int_t totalMultiplicity)
 190 {
 191     // convert the content of a tclones array to an aliflowevent simple for output and tag particles as poi's and rp's
 192     // first step, see if there's already a flow event available in memory. if not create one
 193     if(!fFlowEvent) fFlowEvent = new AliFlowEventSimple(totalMultiplicity, (AliFlowEventSimple::ConstructionMethod)0, 0x0, 0, TMath::TwoPi(), -.9, .9);
 194     // if there is a flow event, clear the members without releasing the allocated memory
 195     else fFlowEvent->ClearFast();
 196     fFlowEvent->SetMCReactionPlaneAngle(fPsi2);
 197     // prepare a track
 198     TParticle* pParticle = 0x0;
 199     Int_t iSelParticlesRP(0);
 200     for (Int_t i(0); i < event->GetEntries(); i++) {    // begin the loop over the input tracks
 201         pParticle = (TParticle*)event->At(i);           // get the particle
 202         if (!pParticle) continue;                       // skip if empty slot (no particle)
 203         if (pParticle->GetNDaughters()!=0) continue;    // see if the particle has daughters (if so, reject it)
 204         AliFlowTrackSimple* pTrack = new AliFlowTrackSimple(pParticle);         // allocate space for the new flow track
 205         pTrack->SetWeight(pParticle->Pz());                                     // ugly hack: store pz here ...
 206         pTrack->SetID(pParticle->GetPdgCode());                                 // set pid code as id
 207         pTrack->SetForRPSelection(kTRUE);                                       // tag ALL particles as RP's,
 208         // ugly hack 2: set charge to -1 for primaries and 1 for secondaries
 209         if(pParticle->GetFirstMother()==-1) pTrack->SetCharge(-1);
 210         else pTrack->SetCharge(1);
 211         iSelParticlesRP++;
 212         fFlowEvent->AddTrack(pTrack);
 213     }
 214     fFlowEvent->SetNumberOfRPs(iSelParticlesRP);
 215     // all trakcs have safely been copied so we can clear the event
 216     event->Clear();
 217     OnTheFlyEventGeneratorCounter = 0;
 218     return fFlowEvent;
 219 }
 220 //_____________________________________________________________________________
 221 void AliFlowOnTheFlyEventGenerator::AddV2(TParticle* part, Double_t v2, Double_t fluc)
 222 {
 223     // afterburner, adds v2, uses Newton-Raphson iteration
 224     Double_t phi = part->Phi();
 225     Double_t phi0=phi;
 226     Double_t f=0.;
 227     Double_t fp=0.;
 228     Double_t phiprev=0.;
 229     // introduce flow fluctuations (gaussian)
 230     if(fluc > -.5) {
 231         v2+=TMath::Sqrt(2*(v2*.25)*(v2*.25))*TMath::ErfInverse(2*fluc-1);
 232         if(fQA) Find(part->GetPdgCode(), kTRUE)->FillV2(part->Pt(), v2);
 233     }
 234     for (Int_t i=0; i!=fMaxNumberOfIterations; ++i) {
 235         phiprev=phi; //store last value for comparison
 236         f =  phi-phi0+v2*TMath::Sin(2.*(phi-fPsi2));
 237         fp = 1.0+2.0*v2*TMath::Cos(2.*(phi-fPsi2)); //first derivative
 238         phi -= f/fp;
 239         if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break;
 240     }
 241     part->SetMomentum( part->Pt()*TMath::Cos(phi), part->Pt()*TMath::Sin(phi), part->Pz(), part->Energy() );
 242 }
 243 //_____________________________________________________________________________
 244 void AliFlowOnTheFlyEventGenerator::AddV2(TClonesArray* event)
 245 {
 246     // afterburner, adds v2 for different species to all tracks in an event
 247     Double_t fluc(gRandom->Uniform());  // get a random number in case of fluctuations
 248     // FIXME at the moment no distincition between mothers and daughters
 249     TParticle *part;
 250     for(Int_t nTrack=0; nTrack!=event->GetEntriesFast(); ++nTrack) {
 251         part = (TParticle*) event->At(nTrack);
 252         // for each track, call overloaded addv2 function with the correct generator
 253         // create a generator in the case where the decayer has introduced a new particle species
 254         if(part)
 255           AddV2(part, Find(part->GetPdgCode(), kTRUE)->GetV2(part->Pt()), fluc);
 256     }
 257 }
 258 //_____________________________________________________________________________
 259 void AliFlowOnTheFlyEventGenerator::GenerateOnTheFlyTracks(Int_t mult, Int_t pid, TClonesArray* event, Double_t fluc)
 260 {
 261     // generate tracks for the event. the tracks that are generated here are the mothers
 262     Double_t ph, et, pt, mm, px, py, pz, ee;
 263     // access the descired generator through the static find method. if the particle species is unknown, we add it to the generator
 264     NaiveFlowAndSpectrumGenerator* generator = Find(pid, kTRUE);
 265     Int_t _tempCounter = OnTheFlyEventGeneratorCounter;
 266     for(int i=_tempCounter; i<=mult+_tempCounter; ++i) {
 267         OnTheFlyEventGeneratorCounter++;
 268         TParticle* part = (TParticle*)event->ConstructedAt(i);
 269         part->SetStatusCode(1);
 270         part->SetFirstMother(-1);
 271         part->SetLastMother(-1);
 272         part->SetFirstDaughter(-1);
 273         part->SetLastDaughter(-1);
 274         ph = gRandom->Uniform(0,TMath::TwoPi());
 275         et = gRandom->Uniform(-1.0,+1.0);
 276         pt = generator->GetPt();
 277         part->SetPdgCode( pid );
 278         mm =  part->GetMass();
 279         px = pt*TMath::Cos(ph);
 280         py = pt*TMath::Sin(ph);
 281         pz = pt*TMath::SinH(et);
 282         ee = TMath::Sqrt( mm*mm + pt*pt+pz*pz );
 283         part->SetMomentum( px, py, pz, ee );
 284         // if requested add v2 to the sample of mothers
 285         if(fAddV2Mothers) AddV2(part, generator->GetV2(pt), fluc);
 286     }
 287 }
 288 //_____________________________________________________________________________
 289 void AliFlowOnTheFlyEventGenerator::DecayOnTheFlyTracks(TClonesArray *event)
 290 {
 291     // decay the tracks using a decayer that is set in the GenerateEventsOnTheFly.C macro
 292     if(!fDecayer) return;       // shouldn't happen ...
 293     fDecayer->Init();
 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);
 336                 //else part->SetDaughter(nDaughter,nTracks+secondaries);
 337                 }
 338             }
 339         nStart = nTracks;
 340         nTracks = event->GetEntries();
 341     }
 342     delete arr;
 343     return;
 344 }
 345 //_____________________________________________________________________________
 346 void AliFlowOnTheFlyEventGenerator::ForceGammaDecay(TClonesArray* arr, TParticle* part)
 347 {
 348     // force gama decay
 349     const Float_t mass1(gRandom->Uniform(0.0000001, .1));   // FIXME randomly sample photon 'mass'
 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 //_____________________________________________________________________________
 396 void 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
 436    for(int i(0); i < 47; i++) fGenerators->Add(new NaiveFlowAndSpectrumGenerator(PDGcode[i], fQA, fFF));
 437 }
 438 //_____________________________________________________________________________
 439 void 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 //_____________________________________________________________________________
 452 void 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