EVGEN/AliGenParam.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 /* $Id$ */
  17
  18 // Class to generate particles from using paramtrized pT and y distributions.
  19 // Distributions are obtained from pointer to object of type AliGenLib.
  20 // (For example AliGenMUONlib)
  21 // Decays are performed using Pythia.
  22 // andreas.morsch@cern.ch
  23
  24 #include <TCanvas.h>
  25 #include <TClonesArray.h>
  26 #include <TDatabasePDG.h>
  27 #include <TF1.h>
  28 #include <TH1F.h>
  29 #include <TLorentzVector.h>
  30 #include <TMath.h>
  31 #include <TParticle.h>
  32 #include <TParticlePDG.h>
  33 #include <TROOT.h>
  34 #include <TVirtualMC.h>
  35
  36 #include "AliDecayer.h"
  37 #include "AliGenMUONlib.h"
  38 #include "AliGenParam.h"
  39 #include "AliMC.h"
  40 #include "AliRun.h"
  41 #include "AliGenEventHeader.h"
  42
  43 ClassImp(AliGenParam)
  44
  45 //------------------------------------------------------------
  46
  47 //Begin_Html
  48 /*
  49     <img src="picts/AliGenParam.gif">
  50 */
  51 //End_Html
  52
  53 //____________________________________________________________
  54 AliGenParam::AliGenParam()
  55 : fPtParaFunc(0),
  56         fYParaFunc(0),
  57         fIpParaFunc(0),
  58   fV2ParaFunc(0),
  59         fPtPara(0),
  60         fYPara(0),
  61   fV2Para(0),
  62   fdNdPhi(0),
  63         fParam(0),
  64         fdNdy0(0.),
  65         fYWgt(0.),
  66         fPtWgt(0.),
  67         fBias(0.),
  68         fTrials(0),
  69         fDeltaPt(0.01),
  70         fSelectAll(kFALSE),
  71   fDecayer(0),
  72   fForceConv(kFALSE)
  73 {
  74   // Default constructor
  75 }
  76 //____________________________________________________________
  77 AliGenParam::AliGenParam(Int_t npart, const AliGenLib * Library,  Int_t param, const char* tname)
  78     :AliGenMC(npart),
  79      fPtParaFunc(Library->GetPt(param, tname)),
  80      fYParaFunc (Library->GetY (param, tname)),
  81      fIpParaFunc(Library->GetIp(param, tname)),
  82    fV2ParaFunc(Library->GetV2(param, tname)),
  83      fPtPara(0),
  84      fYPara(0),
  85    fV2Para(0),
  86    fdNdPhi(0),
  87      fParam(param),
  88      fdNdy0(0.),
  89      fYWgt(0.),
  90      fPtWgt(0.),
  91      fBias(0.),
  92      fTrials(0),
  93      fDeltaPt(0.01),
  94      fSelectAll(kFALSE),
  95    fDecayer(0),
  96    fForceConv(kFALSE)
  97 {
  98   // Constructor using number of particles parameterisation id and library
  99     fName = "Param";
 100     fTitle= "Particle Generator using pT and y parameterisation";
 101     fAnalog = kAnalog;
 102     SetForceDecay();
 103 }
 104 //____________________________________________________________
 105 AliGenParam::AliGenParam(Int_t npart, Int_t param, const char* tname, const char* name):
 106     AliGenMC(npart),
 107     fPtParaFunc(0),
 108     fYParaFunc (0),
 109     fIpParaFunc(0),
 110   fV2ParaFunc(0),
 111     fPtPara(0),
 112     fYPara(0),
 113   fV2Para(0),
 114   fdNdPhi(0),
 115     fParam(param),
 116     fdNdy0(0.),
 117     fYWgt(0.),
 118     fPtWgt(0.),
 119     fBias(0.),
 120     fTrials(0),
 121     fDeltaPt(0.01),
 122     fSelectAll(kFALSE),
 123   fDecayer(0),
 124   fForceConv(kFALSE)
 125 {
 126   // Constructor using parameterisation id and number of particles
 127   //
 128     fName = name;
 129     fTitle= "Particle Generator using pT and y parameterisation";
 130
 131     AliGenLib* pLibrary = new AliGenMUONlib();
 132     fPtParaFunc = pLibrary->GetPt(param, tname);
 133     fYParaFunc  = pLibrary->GetY (param, tname);
 134     fIpParaFunc = pLibrary->GetIp(param, tname);
 135   fV2ParaFunc = pLibrary->GetV2(param, tname);
 136
 137     fAnalog = kAnalog;
 138     fChildSelect.Set(5);
 139     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
 140     SetForceDecay();
 141     SetCutOnChild();
 142     SetChildMomentumRange();
 143     SetChildPtRange();
 144     SetChildPhiRange();
 145     SetChildThetaRange();
 146 }
 147 //____________________________________________________________
 148
 149 AliGenParam::AliGenParam(Int_t npart, Int_t param,
 150                          Double_t (*PtPara) (const Double_t*, const Double_t*),
 151                          Double_t (*YPara ) (const Double_t* ,const Double_t*),
 152                          Double_t (*V2Para) (const Double_t* ,const Double_t*),
 153                          Int_t    (*IpPara) (TRandom *))
 154     :AliGenMC(npart),
 155
 156      fPtParaFunc(PtPara),
 157      fYParaFunc(YPara),
 158      fIpParaFunc(IpPara),
 159    fV2ParaFunc(V2Para),
 160      fPtPara(0),
 161      fYPara(0),
 162    fV2Para(0),
 163    fdNdPhi(0),
 164      fParam(param),
 165      fdNdy0(0.),
 166      fYWgt(0.),
 167      fPtWgt(0.),
 168      fBias(0.),
 169      fTrials(0),
 170      fDeltaPt(0.01),
 171      fSelectAll(kFALSE),
 172   fDecayer(0),
 173   fForceConv(kFALSE)
 174 {
 175   // Constructor
 176   // Gines Martinez 1/10/99
 177     fName = "Param";
 178     fTitle= "Particle Generator using pT and y parameterisation";
 179
 180     fAnalog = kAnalog;
 181     fChildSelect.Set(5);
 182     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
 183     SetForceDecay();
 184     SetCutOnChild();
 185     SetChildMomentumRange();
 186     SetChildPtRange();
 187     SetChildPhiRange();
 188     SetChildThetaRange();
 189 }
 190
 191 //____________________________________________________________
 192 AliGenParam::~AliGenParam()
 193 {
 194   // Destructor
 195     delete  fPtPara;
 196     delete  fYPara;
 197   delete  fV2Para;
 198   delete  fdNdPhi;
 199 }
 200
 201 //-------------------------------------------------------------------
 202 TVector3 AliGenParam::OrthogonalVector(TVector3 &inVec){
 203   double abc[]={inVec.x(), inVec.y(), inVec.z()};
 204   double xyz[]={1,1,1};
 205   int solvDim=0;
 206   double tmp=abc[0];
 207   for(int i=0; i<3; i++)
 208     if(abs(abc[i])>tmp){
 209       solvDim=i;
 210       tmp=abs(abc[i]);
 211     }
 212   xyz[solvDim]=(-abc[(1+solvDim)%3]-abc[(2+solvDim)%3])/abc[(0+solvDim)%3];
 213
 214   TVector3 res(xyz[0],xyz[1],xyz[2]);
 215   return res;
 216 }
 217
 218 double AliGenParam::ScreenFunc1(double d){
 219   if(d>1)
 220     return 21.12-4.184*log(d+0.952);
 221   else
 222     return 20.867-3.242*d+0.652*d*d;
 223 }
 224
 225 double AliGenParam::ScreenFunc2(double d){
 226   if(d>1)
 227     return 21.12-4.184*log(d+0.952);
 228   else
 229     return 20.209-1.93*d-0.086*d*d;
 230 }
 231
 232 double AliGenParam::EnergyFraction(double Z, double E){
 233   double e0=0.000511/E;
 234   double aZ=Z/137.036;
 235   double dmin=ScreenVar(Z,e0,0.5);
 236   double fcZ=(aZ*aZ)*(1/(1+aZ*aZ)+0.20206-0.0368*aZ*aZ+0.0083*aZ*aZ*aZ);
 237   double Fz=8*log(Z)/3;
 238   if(E>0.05)
 239     Fz+=8*fcZ;
 240   double dmax=exp((42.24-Fz)/8.368)-0.952;
 241   if(!(dmax>dmin))
 242     return fRandom->Uniform(e0,0.5);
 243
 244   double e1=0.5-0.5*sqrt(1-dmin/dmax);
 245   double emin=TMath::Max(e0,e1);
 246
 247   double normval=1/(0.5*(ScreenFunc1(dmin)-0.5*Fz)+0.1666667*(ScreenFunc2(dmin)-0.5*Fz));
 248   while(true){
 249     double y=fRandom->Uniform(emin,1);
 250     double eps=y*(0.38453+y*(0.10234+y*(0.026072+y*(0.014367-0.027313*y)))); //inverse to the enveloping cumulative probability distribution
 251     double val=fRandom->Uniform(0,2.12*(eps*eps-eps)+1.53); //enveloping probability density
 252     double d=ScreenVar(Z,e0,eps);
 253     double bh=((eps*eps)+(1-eps)*(1-eps))*(ScreenFunc1(d)-0.5*Fz)+0.6666667*eps*(1-eps)*(ScreenFunc2(d)-0.5*Fz);
 254     bh*=normval;
 255     if(val<bh)
 256       return eps;
 257   }
 258 }
 259
 260 double AliGenParam::PolarAngle(double E){
 261   float rand[3];
 262   AliRndm rndm;
 263   rndm.Rndm(rand,3);
 264   double u=-8*log(rand[1]*rand[2])/5;
 265   if(!(rand[0]<9.0/36))
 266     u/=3;
 267   return u*0.000511/E;
 268 }
 269
 270 Int_t AliGenParam::ForceGammaConversion(TClonesArray *particles, Int_t nPart)
 271 {
 272   //based on: http://geant4.cern.ch/G4UsersDocuments/UsersGuides/PhysicsReferenceManual/html/node27.html
 273   //     and: http://geant4.cern.ch/G4UsersDocuments/UsersGuides/PhysicsReferenceManual/html/node58.html
 274   Int_t nPartNew=nPart;
 275   for(int iPart=0; iPart<nPart; iPart++){
 276     TParticle *gamma = (TParticle *) particles->At(iPart);
 277     if(gamma->GetPdgCode()!=22) continue;
 278
 279     TVector3 gammaV3(gamma->Px(),gamma->Py(),gamma->Pz());
 280     Float_t az=fRandom->Uniform(TMath::Pi()*2);
 281     double frac=EnergyFraction(1,gamma->Energy());
 282     double Ee1=frac*gamma->Energy();
 283     double Ee2=(1-frac)*gamma->Energy();
 284     double Pe1=Ee1;//sqrt(Ee1*Ee1-0.000511*0.000511);
 285     double Pe2=Ee2;//sqrt(Ee2*Ee2-0.000511*0.000511);
 286
 287     TVector3 rotAxis(OrthogonalVector(gammaV3));
 288     rotAxis.Rotate(az,gammaV3);
 289     TVector3 e1V3(gammaV3);
 290     e1V3.Rotate(PolarAngle(Ee1),rotAxis);
 291     e1V3=e1V3.Unit();
 292     e1V3*=Pe1;
 293     TVector3 e2V3(gammaV3);
 294     e2V3.Rotate(-PolarAngle(Ee2),rotAxis);
 295     e2V3=e2V3.Unit();
 296     e2V3*=Pe2;
 297     // gamma = new TParticle(*gamma);
 298     // particles->RemoveAt(iPart);
 299     gamma->SetFirstDaughter(nPartNew+1);
 300     gamma->SetLastDaughter(nPartNew+2);
 301     // new((*particles)[iPart]) TParticle(*gamma);
 302     // delete gamma;
 303
 304     TLorentzVector vtx;
 305     gamma->ProductionVertex(vtx);
 306     new((*particles)[nPartNew]) TParticle(11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, TLorentzVector(e1V3,Ee1), vtx);
 307     nPartNew++;
 308     new((*particles)[nPartNew]) TParticle(-11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, TLorentzVector(e2V3,Ee2), vtx);
 309     nPartNew++;
 310   }
 311   // particles->Compress();
 312   return particles->GetEntriesFast();
 313 }
 314
 315 //____________________________________________________________
 316 void AliGenParam::Init()
 317 {
 318   // Initialisation
 319
 320     if (gMC) fDecayer = gMC->GetDecayer();
 321   //Begin_Html
 322   /*
 323     <img src="picts/AliGenParam.gif">
 324   */
 325   //End_Html
 326     char name[256];
 327     snprintf(name, 256, "pt-parameterisation for %s", GetName());
 328
 329     if (fPtPara) fPtPara->Delete();
 330     fPtPara = new TF1(name, fPtParaFunc, fPtMin, fPtMax,0);
 331     gROOT->GetListOfFunctions()->Remove(fPtPara);
 332   //  Set representation precision to 10 MeV
 333     Int_t npx= Int_t((fPtMax - fPtMin) / fDeltaPt);
 334
 335     fPtPara->SetNpx(npx);
 336
 337     snprintf(name, 256, "y-parameterisation  for %s", GetName());
 338     if (fYPara) fYPara->Delete();
 339     fYPara  = new TF1(name, fYParaFunc, fYMin, fYMax, 0);
 340     gROOT->GetListOfFunctions()->Remove(fYPara);
 341
 342   snprintf(name, 256, "v2-parameterisation for %s", GetName());
 343   if (fV2Para) fV2Para->Delete();
 344   fV2Para  = new TF1(name, fV2ParaFunc, fPtMin, fPtMax, 0);
 345   // fV2Para  = new TF1(name, "2*[0]/(1+TMath::Exp([1]*([2]-x)))-[0]", fPtMin, fPtMax);
 346   // fV2Para->SetParameter(0, 0.236910);
 347   // fV2Para->SetParameter(1, 1.71122);
 348   // fV2Para->SetParameter(2, 0.0827617);
 349   //gROOT->GetListOfFunctions()->Remove(fV2Para);  //TR: necessary?
 350
 351   snprintf(name, 256, "dNdPhi for %s", GetName());
 352   if (fdNdPhi) fdNdPhi->Delete();
 353   fdNdPhi  = new TF1(name, "1+2*[0]*TMath::Cos(2*(x-[1]))", fPhiMin, fPhiMax);
 354   //gROOT->GetListOfFunctions()->Remove(fdNdPhi);  //TR: necessary?
 355
 356     snprintf(name, 256, "pt-for-%s", GetName());
 357     TF1 ptPara(name ,fPtParaFunc, 0, 15, 0);
 358     snprintf(name, 256, "y-for-%s", GetName());
 359     TF1 yPara(name, fYParaFunc, -6, 6, 0);
 360
 361   //
 362   // dN/dy| y=0
 363     Double_t y1=0;
 364     Double_t y2=0;
 365
 366     fdNdy0=fYParaFunc(&y1,&y2);
 367   //
 368   // Integral over generation region
 369 #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
 370     Float_t intYS  = yPara.Integral(fYMin, fYMax,(Double_t*) 0x0,1.e-6);
 371     Float_t intPt0 = ptPara.Integral(0,15,(Double_t *) 0x0,1.e-6);
 372     Float_t intPtS = ptPara.Integral(fPtMin,fPtMax,(Double_t*) 0x0,1.e-6);
 373 #else
 374     Float_t intYS  = yPara.Integral(fYMin, fYMax,1.e-6);
 375     Float_t intPt0 = ptPara.Integral(0,15,1.e-6);
 376     Float_t intPtS = ptPara.Integral(fPtMin,fPtMax,1.e-6);
 377 #endif
 378   Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();    //TR: should probably be done differently in case of anisotropic phi...
 379     if (fAnalog == kAnalog) {
 380         fYWgt  = intYS/fdNdy0;
 381         fPtWgt = intPtS/intPt0;
 382         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
 383     } else {
 384         fYWgt = intYS/fdNdy0;
 385         fPtWgt = (fPtMax-fPtMin)/intPt0;
 386         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
 387     }
 388   //
 389   // particle decay related initialization
 390     fDecayer->SetForceDecay(fForceDecay);
 391     fDecayer->Init();
 392
 393   //
 394     AliGenMC::Init();
 395 }
 396
 397 //____________________________________________________________
 398 void AliGenParam::Generate()
 399 {
 400   //
 401   // Generate 1 event (see Generate(Int_t ntimes) for details
 402   //
 403   GenerateN(1);
 404 }
 405 //____________________________________________________________
 406 void AliGenParam::GenerateN(Int_t ntimes)
 407 {
 408   //
 409   // Generate ntimes*'npart' light and heavy mesons (J/Psi, upsilon or phi, Pion,
 410   // Kaons, Etas, Omegas) and Baryons (proton, antiprotons, neutrons and
 411   // antineutrons in the the desired theta, phi and momentum windows;
 412   // Gaussian smearing on the vertex is done if selected.
 413   // The decay of heavy mesons is done using lujet,
 414   //    and the childern particle are tracked by GEANT
 415   // However, light mesons are directly tracked by GEANT
 416   // setting fForceDecay = nodecay (SetForceDecay(nodecay))
 417   //
 418   //
 419   //  Reinitialize decayer
 420   fDecayer->SetForceDecay(fForceDecay);
 421   fDecayer->Init();
 422
 423   //
 424   Float_t polar[3]= {0,0,0};  // Polarisation of the parent particle (for GEANT tracking)
 425   Float_t origin0[3];         // Origin of the generated parent particle (for GEANT tracking)
 426   Float_t time0;              // Time0 of the generated parent particle
 427   Float_t pt, pl, ptot;       // Transverse, logitudinal and total momenta of the parent particle
 428   Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
 429   Float_t p[3], pc[3],
 430           och[3];             // Momentum, polarisation and origin of the children particles from lujet
 431   Double_t ty, xmt;
 432   Int_t nt, i, j;
 433   Float_t  wgtp, wgtch;
 434   Double_t dummy;
 435   static TClonesArray *particles;
 436   //
 437   if(!particles) particles = new TClonesArray("TParticle",1000);
 438
 439   TDatabasePDG *pDataBase = TDatabasePDG::Instance();
 440   //
 441   Float_t random[6];
 442
 443   // Calculating vertex position per event
 444   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
 445   time0 = fTimeOrigin;
 446   if(fVertexSmear==kPerEvent) {
 447       Vertex();
 448       for (j=0;j<3;j++) origin0[j]=fVertex[j];
 449       time0 = fTime;
 450   }
 451
 452   Int_t ipa=0;
 453
 454   // Generating fNpart particles
 455   fNprimaries = 0;
 456
 457   Int_t nGen = fNpart*ntimes;
 458   while (ipa<nGen) {
 459       while(1) {
 460       //
 461       // particle type
 462           Int_t iPart = fIpParaFunc(fRandom);
 463           fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;
 464           TParticlePDG *particle = pDataBase->GetParticle(iPart);
 465           Float_t am = particle->Mass();
 466
 467           Rndm(random,2);
 468       //
 469       // y
 470           ty = TMath::TanH(fYPara->GetRandom());
 471       //
 472       // pT
 473           if (fAnalog == kAnalog) {
 474               pt=fPtPara->GetRandom();
 475               wgtp=fParentWeight;
 476               wgtch=fChildWeight;
 477           } else {
 478               pt=fPtMin+random[1]*(fPtMax-fPtMin);
 479               Double_t ptd=pt;
 480               wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
 481               wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
 482           }
 483           xmt=sqrt(pt*pt+am*am);
 484           if (TMath::Abs(ty)==1.) {
 485               ty=0.;
 486               Fatal("AliGenParam",
 487                     "Division by 0: Please check you rapidity range !");
 488           }
 489       //
 490       // phi
 491           //      if(!ipa)
 492           //phi=fEvPlane; //align first particle of each event with event plane
 493           //else{
 494         double v2 = fV2Para->Eval(pt);
 495         fdNdPhi->SetParameter(0,v2);
 496         fdNdPhi->SetParameter(1,fEvPlane);
 497         phi=fdNdPhi->GetRandom();
 498         //     }
 499
 500           pl=xmt*ty/sqrt((1.-ty)*(1.+ty));
 501           theta=TMath::ATan2(pt,pl);
 502       // Cut on theta
 503           if(theta<fThetaMin || theta>fThetaMax) continue;
 504           ptot=TMath::Sqrt(pt*pt+pl*pl);
 505       // Cut on momentum
 506           if(ptot<fPMin || ptot>fPMax) continue;
 507       //
 508           p[0]=pt*TMath::Cos(phi);
 509           p[1]=pt*TMath::Sin(phi);
 510           p[2]=pl;
 511           if(fVertexSmear==kPerTrack) {
 512               Rndm(random,6);
 513               for (j=0;j<3;j++) {
 514                   origin0[j]=
 515                       fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 516                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 517               }
 518               Rndm(random,2);
 519               time0 = fTimeOrigin + fOsigma[2]/TMath::Ccgs()*
 520                 TMath::Cos(2*random[0]*TMath::Pi())*
 521                 TMath::Sqrt(-2*TMath::Log(random[1]));
 522           }
 523
 524       // Looking at fForceDecay :
 525       // if fForceDecay != none Primary particle decays using
 526       // AliPythia and children are tracked by GEANT
 527       //
 528       // if fForceDecay == none Primary particle is tracked by GEANT
 529       // (In the latest, make sure that GEANT actually does all the decays you want)
 530       //
 531           Bool_t decayed = kFALSE;
 532
 533
 534           if (fForceDecay != kNoDecay) {
 535         // Using lujet to decay particle
 536               Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
 537               TLorentzVector pmom(p[0], p[1], p[2], energy);
 538               fDecayer->Decay(iPart,&pmom);
 539         //
 540         // select decay particles
 541               Int_t np=fDecayer->ImportParticles(particles);
 542
 543         if(fForceConv) np=ForceGammaConversion(particles,np);
 544
 545         // for(int iPart=0; iPart<np; iPart++){
 546         //   TParticle *gamma = (TParticle *) particles->At(iPart);
 547         //   printf("%i %i:", iPart, gamma->GetPdgCode());
 548         //   printf("%i %i %i|",gamma->GetFirstMother(),gamma->GetFirstDaughter(),gamma->GetLastDaughter());
 549         // }
 550
 551               //  Selecting  GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut;
 552               if (fGeometryAcceptance)
 553                 if (!CheckAcceptanceGeometry(np,particles)) continue;
 554               Int_t ncsel=0;
 555               Int_t* pFlag      = new Int_t[np];
 556               Int_t* pParent    = new Int_t[np];
 557               Int_t* pSelected  = new Int_t[np];
 558               Int_t* trackIt    = new Int_t[np];
 559
 560               for (i=0; i<np; i++) {
 561                   pFlag[i]     =  0;
 562                   pSelected[i] =  0;
 563                   pParent[i]   = -1;
 564               }
 565
 566               if (np >1) {
 567                   decayed = kTRUE;
 568                   TParticle* iparticle =  0;
 569                   Int_t ipF, ipL;
 570                   for (i = 1; i<np ; i++) {
 571                       trackIt[i] = 1;
 572                       iparticle = (TParticle *) particles->At(i);
 573                       Int_t kf = iparticle->GetPdgCode();
 574                       Int_t ks = iparticle->GetStatusCode();
 575             // flagged particle
 576
 577                       if (pFlag[i] == 1) {
 578                           ipF = iparticle->GetFirstDaughter();
 579                           ipL = iparticle->GetLastDaughter();
 580                           if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
 581                           continue;
 582                       }
 583
 584             // flag decay products of particles with long life-time (c tau > .3 mum)
 585
 586                       if (ks != 1) {
 587               //                          TParticlePDG *particle = pDataBase->GetParticle(kf);
 588
 589                           Double_t lifeTime = fDecayer->GetLifetime(kf);
 590               //                          Double_t mass     = particle->Mass();
 591               //                          Double_t width    = particle->Width();
 592                           if (lifeTime > (Double_t) fMaxLifeTime) {
 593                               ipF = iparticle->GetFirstDaughter();
 594                               ipL = iparticle->GetLastDaughter();
 595                               if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
 596                           } else{
 597                               trackIt[i]     = 0;
 598                               pSelected[i]   = 1;
 599                           }
 600                       } // ks==1 ?
 601             //
 602             // children
 603
 604                       if ((ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll || fSelectAll) && trackIt[i])
 605                       {
 606                           if (fCutOnChild) {
 607                               pc[0]=iparticle->Px();
 608                               pc[1]=iparticle->Py();
 609                               pc[2]=iparticle->Pz();
 610                               Bool_t  childok = KinematicSelection(iparticle, 1);
 611                               if(childok) {
 612                                   pSelected[i]  = 1;
 613                                   ncsel++;
 614                               } else {
 615                                   ncsel=-1;
 616                                   break;
 617                               } // child kine cuts
 618                           } else {
 619                               pSelected[i]  = 1;
 620                               ncsel++;
 621                           } // if child selection
 622                       } // select muon
 623                   } // decay particle loop
 624               } // if decay products
 625
 626               Int_t iparent;
 627               if ((fCutOnChild && ncsel >0) || !fCutOnChild){
 628                   ipa++;
 629           //
 630           // Parent
 631
 632
 633                   PushTrack(0, -1, iPart, p, origin0, polar, time0, kPPrimary, nt, wgtp, ((decayed)? 11 : 1));
 634                   pParent[0] = nt;
 635                   KeepTrack(nt);
 636                   fNprimaries++;
 637
 638           //
 639           // Decay Products
 640           //
 641                   for (i = 1; i < np; i++) {
 642                       if (pSelected[i]) {
 643                           TParticle* iparticle = (TParticle *) particles->At(i);
 644                           Int_t kf   = iparticle->GetPdgCode();
 645                           Int_t ksc  = iparticle->GetStatusCode();
 646                           Int_t jpa  = iparticle->GetFirstMother()-1;
 647
 648                           och[0] = origin0[0]+iparticle->Vx();
 649                           och[1] = origin0[1]+iparticle->Vy();
 650                           och[2] = origin0[2]+iparticle->Vz();
 651                           pc[0]  = iparticle->Px();
 652                           pc[1]  = iparticle->Py();
 653                           pc[2]  = iparticle->Pz();
 654
 655                           if (jpa > -1) {
 656                               iparent = pParent[jpa];
 657                           } else {
 658                               iparent = -1;
 659                           }
 660
 661                           PushTrack(fTrackIt * trackIt[i], iparent, kf,
 662                                            pc, och, polar,
 663                                            time0 + iparticle->T(), kPDecay, nt, wgtch, ksc);
 664                           pParent[i] = nt;
 665                           KeepTrack(nt);
 666                           fNprimaries++;
 667                       } // Selected
 668                   } // Particle loop
 669               }  // Decays by Lujet
 670               particles->Clear();
 671               if (pFlag)      delete[] pFlag;
 672               if (pParent)    delete[] pParent;
 673               if (pSelected)  delete[] pSelected;
 674               if (trackIt)    delete[] trackIt;
 675           } // kinematic selection
 676           else  // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
 677           {
 678             gAlice->GetMCApp()->
 679                 PushTrack(fTrackIt,-1,iPart,p,origin0,polar,time0,kPPrimary,nt,wgtp, 1);
 680             ipa++;
 681             fNprimaries++;
 682           }
 683           break;
 684     } // while
 685   } // event loop
 686
 687   SetHighWaterMark(nt);
 688
 689   AliGenEventHeader* header = new AliGenEventHeader("PARAM");
 690   header->SetPrimaryVertex(fVertex);
 691   header->SetInteractionTime(fTime);
 692   header->SetNProduced(fNprimaries);
 693   AddHeader(header);
 694 }
 695 //____________________________________________________________________________________
 696 Float_t AliGenParam::GetRelativeArea(Float_t ptMin, Float_t ptMax, Float_t yMin, Float_t yMax, Float_t phiMin, Float_t phiMax)
 697 {
 698   //
 699   // Normalisation for selected kinematic region
 700   //
 701 #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
 702   Float_t ratio =
 703     fPtPara->Integral(ptMin,ptMax,(Double_t *)0,1.e-6) / fPtPara->Integral( fPtPara->GetXmin(), fPtPara->GetXmax(),(Double_t *)0,1.e-6) *
 704     fYPara->Integral(yMin,yMax,(Double_t *)0,1.e-6)/fYPara->Integral(fYPara->GetXmin(),fYPara->GetXmax(),(Double_t *)0,1.e-6)   *
 705     (phiMax-phiMin)/360.;
 706 #else
 707   Float_t ratio =
 708     fPtPara->Integral(ptMin,ptMax,1.e-6) / fPtPara->Integral( fPtPara->GetXmin(), fPtPara->GetXmax(),1.e-6) *
 709     fYPara->Integral(yMin,yMax,1.e-6)/fYPara->Integral(fYPara->GetXmin(),fYPara->GetXmax(),1.e-6)   *
 710     (phiMax-phiMin)/360.;
 711 #endif
 712   return TMath::Abs(ratio);
 713 }
 714
 715 //____________________________________________________________________________________
 716
 717 void AliGenParam::Draw( const char * /*opt*/)
 718 {
 719     //
 720     // Draw the pT and y Distributions
 721     //
 722      TCanvas *c0 = new TCanvas("c0","Canvas 0",400,10,600,700);
 723      c0->Divide(2,1);
 724      c0->cd(1);
 725      fPtPara->Draw();
 726      fPtPara->GetHistogram()->SetXTitle("p_{T} (GeV)");
 727      c0->cd(2);
 728      fYPara->Draw();
 729      fYPara->GetHistogram()->SetXTitle("y");
 730 }
 731
 732
 733
 734