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 void AliGenParam::RotateVector(Double_t *pin, Double_t *pout, Double_t costheta, Double_t sintheta,
 219     Double_t cosphi, Double_t sinphi)
 220 {
 221   // Perform rotation
 222   pout[0] = pin[0]*costheta*cosphi-pin[1]*sinphi+pin[2]*sintheta*cosphi;
 223   pout[1] = pin[0]*costheta*sinphi+pin[1]*cosphi+pin[2]*sintheta*sinphi;
 224   pout[2] = -1.0  * pin[0] * sintheta + pin[2] * costheta;
 225   return;
 226 }
 227
 228 double AliGenParam::ScreenFunction1(double screenVariable){
 229   if(screenVariable>1)
 230     return 42.24 - 8.368 * log(screenVariable + 0.952);
 231   else
 232     return 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
 233 }
 234
 235 double AliGenParam::ScreenFunction2(double screenVariable){
 236   if(screenVariable>1)
 237     return 42.24 - 8.368 * log(screenVariable + 0.952);
 238   else
 239     return 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
 240 }
 241
 242 double AliGenParam::RandomEnergyFraction(double Z, double photonEnergy){
 243   double aZ=Z/137.036;
 244   double epsilon ;
 245   double epsilon0Local = 0.000511 / photonEnergy ;
 246
 247   // Do it fast if photon energy < 2. MeV
 248   if (photonEnergy < 0.002 )
 249     {
 250       epsilon = epsilon0Local + (0.5 - epsilon0Local) * fRandom->Rndm();
 251     }
 252   else
 253     {
 254       double fZ = 8*log(Z)/3;
 255   double fcZ=(aZ*aZ)*(1/(1+aZ*aZ)+0.20206-0.0368*aZ*aZ+0.0083*aZ*aZ*aZ);
 256       if (photonEnergy > 0.050) fZ += 8*fcZ;
 257
 258       // Limits of the screening variable
 259       double screenFactor = 136. * epsilon0Local / pow (Z,1/3);
 260       double screenMax = exp ((42.24 - fZ)/8.368) - 0.952 ;
 261       double screenMin = std::min(4.*screenFactor,screenMax) ;
 262
 263       // Limits of the energy sampling
 264       double epsilon1 = 0.5 - 0.5 * sqrt(1. - screenMin / screenMax) ;
 265       double epsilonMin = std::max(epsilon0Local,epsilon1);
 266       double epsilonRange = 0.5 - epsilonMin ;
 267
 268       // Sample the energy rate of the created electron (or positron)
 269       double screen;
 270       double gReject ;
 271
 272       double f10 = ScreenFunction1(screenMin) - fZ;
 273       double f20 = ScreenFunction2(screenMin) - fZ;
 274       double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
 275       double normF2 = std::max(1.5 * f20,0.);
 276
 277       do
 278         {
 279           if (normF1 / (normF1 + normF2) > fRandom->Rndm() )
 280             {
 281               epsilon = 0.5 - epsilonRange * pow(fRandom->Rndm(), 0.333333) ;
 282               screen = screenFactor / (epsilon * (1. - epsilon));
 283               gReject = (ScreenFunction1(screen) - fZ) / f10 ;
 284             }
 285           else
 286             {
 287               epsilon = epsilonMin + epsilonRange * fRandom->Rndm();
 288               screen = screenFactor / (epsilon * (1 - epsilon));
 289               gReject = (ScreenFunction2(screen) - fZ) / f20 ;
 290             }
 291         } while ( gReject < fRandom->Rndm() );
 292
 293     }   //  End of epsilon sampling
 294   return epsilon;
 295 }
 296
 297 double AliGenParam::RandomPolarAngle(){
 298   double u;
 299   const double a1 = 0.625;
 300   double a2 = 3. * a1;
 301   //  double d = 27. ;
 302
 303   //  if (9. / (9. + d) > fRandom->Rndm())
 304   if (0.25 > fRandom->Rndm())
 305     {
 306       u = - log(fRandom->Rndm() * fRandom->Rndm()) / a1 ;
 307     }
 308   else
 309     {
 310       u = - log(fRandom->Rndm() * fRandom->Rndm()) / a2 ;
 311     }
 312   return u*0.000511;
 313 }
 314
 315 Double_t AliGenParam::RandomMass(Double_t mh){
 316   while(true){
 317     double y=fRandom->Rndm();
 318     double mee=2*0.000511*TMath::Power(2*0.000511/mh,-y); //inverse of the enveloping cumulative distribution
 319     double apxkw=2.0/3.0/137.036/TMath::Pi()/mee; //enveloping probability density
 320     double val=fRandom->Uniform(0,apxkw);
 321     double kw=apxkw*sqrt(1-4*0.000511*0.000511/mee/mee)*(1+2*0.000511*0.000511/mee/mee)*1*1*TMath::Power(1-mee*mee/mh/mh,3);
 322     if(val<kw)
 323       return mee;
 324   }
 325 }
 326
 327 Int_t AliGenParam::VirtualGammaPairProduction(TClonesArray *particles, Int_t nPart)
 328 {
 329   Int_t nPartNew=nPart;
 330   for(int iPart=0; iPart<nPart; iPart++){
 331     TParticle *gamma = (TParticle *) particles->At(iPart);
 332     if(gamma->GetPdgCode()!=220000) continue;
 333     if(gamma->Pt()<0.002941) continue;  //approximation of kw in AliGenEMlib is 0 below 0.002941
 334     double mass=RandomMass(gamma->Pt());
 335
 336     // lepton pair kinematics in virtual photon rest frame
 337     double Ee=mass/2;
 338     double Pe=TMath::Sqrt((Ee+0.000511)*(Ee-0.000511));
 339
 340     double costheta = (2.0 * gRandom->Rndm()) - 1.;
 341     double sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta));
 342     double phi      = 2.0 * TMath::ACos(-1.) * gRandom->Rndm();
 343     double sinphi   = TMath::Sin(phi);
 344     double cosphi   = TMath::Cos(phi);
 345
 346     // momentum vectors of leptons in virtual photon rest frame
 347     Double_t pProd1[3] = {Pe * sintheta * cosphi,
 348                           Pe * sintheta * sinphi,
 349                           Pe * costheta};
 350
 351     Double_t pProd2[3] = {-1.0 * Pe * sintheta * cosphi,
 352                           -1.0 * Pe * sintheta * sinphi,
 353                           -1.0 * Pe * costheta};
 354
 355     // lepton 4-vectors in properly rotated virtual photon rest frame
 356     Double_t pRot1[3] = {0.};
 357     RotateVector(pProd1, pRot1, costheta, -sintheta, -cosphi, -sinphi);
 358     Double_t pRot2[3] = {0.};
 359     RotateVector(pProd2, pRot2, costheta, -sintheta, -cosphi, -sinphi);
 360
 361     TLorentzVector e1V4(pRot1[0],pRot1[1],pRot1[2],Ee);
 362     TLorentzVector e2V4(pRot2[0],pRot2[1],pRot2[2],Ee);
 363
 364     TVector3 boost(gamma->Px(),gamma->Py(),gamma->Pz());
 365     boost*=1/sqrt(gamma->P()*gamma->P()+mass*mass);
 366     e1V4.Boost(boost);
 367     e2V4.Boost(boost);
 368
 369     TLorentzVector vtx;
 370     gamma->ProductionVertex(vtx);
 371     new((*particles)[nPartNew]) TParticle(11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, e1V4, vtx);
 372     nPartNew++;
 373     new((*particles)[nPartNew]) TParticle(-11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, e2V4, vtx);
 374     nPartNew++;
 375   }
 376   return nPartNew;
 377 }
 378
 379 Int_t AliGenParam::ForceGammaConversion(TClonesArray *particles, Int_t nPart)
 380 {
 381   //based on: http://geant4.cern.ch/G4UsersDocuments/UsersGuides/PhysicsReferenceManual/html/node27.html
 382   //     and: http://geant4.cern.ch/G4UsersDocuments/UsersGuides/PhysicsReferenceManual/html/node58.html
 383   //     and: G4LivermoreGammaConversionModel.cc
 384   Int_t nPartNew=nPart;
 385   for(int iPart=0; iPart<nPart; iPart++){
 386     TParticle *gamma = (TParticle *) particles->At(iPart);
 387     if(gamma->GetPdgCode()!=22) continue;
 388     if(gamma->Energy()<0.001022) continue;
 389     TVector3 gammaV3(gamma->Px(),gamma->Py(),gamma->Pz());
 390     double frac=RandomEnergyFraction(1,gamma->Energy());
 391     double Ee1=frac*gamma->Energy();
 392     double Ee2=(1-frac)*gamma->Energy();
 393     double Pe1=sqrt((Ee1+0.000511)*(Ee1-0.000511));
 394     double Pe2=sqrt((Ee2+0.000511)*(Ee2-0.000511));
 395
 396     TVector3 rotAxis(OrthogonalVector(gammaV3));
 397     Float_t az=fRandom->Uniform(TMath::Pi()*2);
 398     rotAxis.Rotate(az,gammaV3);
 399     TVector3 e1V3(gammaV3);
 400     double u=RandomPolarAngle();
 401     e1V3.Rotate(u/Ee1,rotAxis);
 402     e1V3=e1V3.Unit();
 403     e1V3*=Pe1;
 404     TVector3 e2V3(gammaV3);
 405     e2V3.Rotate(-u/Ee2,rotAxis);
 406     e2V3=e2V3.Unit();
 407     e2V3*=Pe2;
 408     // gamma = new TParticle(*gamma);
 409     // particles->RemoveAt(iPart);
 410     gamma->SetFirstDaughter(nPartNew+1);
 411     gamma->SetLastDaughter(nPartNew+2);
 412     // new((*particles)[iPart]) TParticle(*gamma);
 413     // delete gamma;
 414
 415     TLorentzVector vtx;
 416     gamma->ProductionVertex(vtx);
 417     new((*particles)[nPartNew]) TParticle(11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, TLorentzVector(e1V3,Ee1), vtx);
 418     nPartNew++;
 419     new((*particles)[nPartNew]) TParticle(-11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, TLorentzVector(e2V3,Ee2), vtx);
 420     nPartNew++;
 421   }
 422   return nPartNew;
 423 }
 424
 425 //____________________________________________________________
 426 void AliGenParam::Init()
 427 {
 428   // Initialisation
 429
 430     if (TVirtualMC::GetMC()) fDecayer = TVirtualMC::GetMC()->GetDecayer();
 431   //Begin_Html
 432   /*
 433     <img src="picts/AliGenParam.gif">
 434   */
 435   //End_Html
 436     char name[256];
 437     snprintf(name, 256, "pt-parameterisation for %s", GetName());
 438
 439     if (fPtPara) fPtPara->Delete();
 440     fPtPara = new TF1(name, fPtParaFunc, fPtMin, fPtMax,0);
 441     gROOT->GetListOfFunctions()->Remove(fPtPara);
 442   //  Set representation precision to 10 MeV
 443     Int_t npx= Int_t((fPtMax - fPtMin) / fDeltaPt);
 444
 445     fPtPara->SetNpx(npx);
 446
 447     snprintf(name, 256, "y-parameterisation  for %s", GetName());
 448     if (fYPara) fYPara->Delete();
 449     fYPara  = new TF1(name, fYParaFunc, fYMin, fYMax, 0);
 450     gROOT->GetListOfFunctions()->Remove(fYPara);
 451
 452   snprintf(name, 256, "v2-parameterisation for %s", GetName());
 453   if (fV2Para) fV2Para->Delete();
 454   fV2Para  = new TF1(name, fV2ParaFunc, fPtMin, fPtMax, 0);
 455   // fV2Para  = new TF1(name, "2*[0]/(1+TMath::Exp([1]*([2]-x)))-[0]", fPtMin, fPtMax);
 456   // fV2Para->SetParameter(0, 0.236910);
 457   // fV2Para->SetParameter(1, 1.71122);
 458   // fV2Para->SetParameter(2, 0.0827617);
 459   //gROOT->GetListOfFunctions()->Remove(fV2Para);  //TR: necessary?
 460
 461   snprintf(name, 256, "dNdPhi for %s", GetName());
 462   if (fdNdPhi) fdNdPhi->Delete();
 463   fdNdPhi  = new TF1(name, "1+2*[0]*TMath::Cos(2*(x-[1]))", fPhiMin, fPhiMax);
 464   //gROOT->GetListOfFunctions()->Remove(fdNdPhi);  //TR: necessary?
 465
 466     snprintf(name, 256, "pt-for-%s", GetName());
 467     TF1 ptPara(name ,fPtParaFunc, 0, 15, 0);
 468     snprintf(name, 256, "y-for-%s", GetName());
 469     TF1 yPara(name, fYParaFunc, -6, 6, 0);
 470
 471   //
 472   // dN/dy| y=0
 473     Double_t y1=0;
 474     Double_t y2=0;
 475
 476     fdNdy0=fYParaFunc(&y1,&y2);
 477   //
 478   // Integral over generation region
 479 #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
 480     Float_t intYS  = yPara.Integral(fYMin, fYMax,(Double_t*) 0x0,1.e-6);
 481     Float_t intPt0 = ptPara.Integral(0,15,(Double_t *) 0x0,1.e-6);
 482     Float_t intPtS = ptPara.Integral(fPtMin,fPtMax,(Double_t*) 0x0,1.e-6);
 483 #else
 484     Float_t intYS  = yPara.Integral(fYMin, fYMax,1.e-6);
 485     Float_t intPt0 = ptPara.Integral(0,15,1.e-6);
 486     Float_t intPtS = ptPara.Integral(fPtMin,fPtMax,1.e-6);
 487 #endif
 488   Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();    //TR: should probably be done differently in case of anisotropic phi...
 489     if (fAnalog == kAnalog) {
 490         fYWgt  = intYS/fdNdy0;
 491         fPtWgt = intPtS/intPt0;
 492         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
 493     } else {
 494         fYWgt = intYS/fdNdy0;
 495         fPtWgt = (fPtMax-fPtMin)/intPt0;
 496         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
 497     }
 498   //
 499   // particle decay related initialization
 500     fDecayer->SetForceDecay(fForceDecay);
 501     fDecayer->Init();
 502
 503   //
 504     AliGenMC::Init();
 505 }
 506
 507 //____________________________________________________________
 508 void AliGenParam::Generate()
 509 {
 510   //
 511   // Generate 1 event (see Generate(Int_t ntimes) for details
 512   //
 513   GenerateN(1);
 514 }
 515 //____________________________________________________________
 516 void AliGenParam::GenerateN(Int_t ntimes)
 517 {
 518   //
 519   // Generate ntimes*'npart' light and heavy mesons (J/Psi, upsilon or phi, Pion,
 520   // Kaons, Etas, Omegas) and Baryons (proton, antiprotons, neutrons and
 521   // antineutrons in the the desired theta, phi and momentum windows;
 522   // Gaussian smearing on the vertex is done if selected.
 523   // The decay of heavy mesons is done using lujet,
 524   //    and the childern particle are tracked by GEANT
 525   // However, light mesons are directly tracked by GEANT
 526   // setting fForceDecay = nodecay (SetForceDecay(nodecay))
 527   //
 528   //
 529   //  Reinitialize decayer
 530   fDecayer->SetForceDecay(fForceDecay);
 531   fDecayer->Init();
 532
 533   //
 534   Float_t polar[3]= {0,0,0};  // Polarisation of the parent particle (for GEANT tracking)
 535   Float_t origin0[3];         // Origin of the generated parent particle (for GEANT tracking)
 536   Float_t time0;              // Time0 of the generated parent particle
 537   Float_t pt, pl, ptot;       // Transverse, logitudinal and total momenta of the parent particle
 538   Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
 539   Float_t p[3], pc[3],
 540           och[3];             // Momentum, polarisation and origin of the children particles from lujet
 541   Double_t ty, xmt;
 542   Int_t nt, i, j;
 543   Float_t  wgtp, wgtch;
 544   Double_t dummy;
 545   static TClonesArray *particles;
 546   //
 547   if(!particles) particles = new TClonesArray("TParticle",1000);
 548
 549   TDatabasePDG *pDataBase = TDatabasePDG::Instance();
 550   //
 551   Float_t random[6];
 552
 553   // Calculating vertex position per event
 554   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
 555   time0 = fTimeOrigin;
 556   if(fVertexSmear==kPerEvent) {
 557       Vertex();
 558       for (j=0;j<3;j++) origin0[j]=fVertex[j];
 559       time0 = fTime;
 560   }
 561
 562   Int_t ipa=0;
 563
 564   // Generating fNpart particles
 565   fNprimaries = 0;
 566
 567   Int_t nGen = fNpart*ntimes;
 568   while (ipa<nGen) {
 569       while(1) {
 570       //
 571       // particle type
 572           Int_t iPart = fIpParaFunc(fRandom);
 573       Int_t iTemp = iPart;
 574
 575       // custom pdg codes for to destinguish direct photons
 576       if(iPart==220000) iPart=22;
 577
 578           fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;
 579           TParticlePDG *particle = pDataBase->GetParticle(iPart);
 580           Float_t am = particle->Mass();
 581
 582           Rndm(random,2);
 583       //
 584       // y
 585           ty = TMath::TanH(fYPara->GetRandom());
 586
 587       //
 588       // pT
 589           if (fAnalog == kAnalog) {
 590               pt=fPtPara->GetRandom();
 591               wgtp=fParentWeight;
 592               wgtch=fChildWeight;
 593           } else {
 594               pt=fPtMin+random[1]*(fPtMax-fPtMin);
 595               Double_t ptd=pt;
 596               wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
 597               wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
 598           }
 599           xmt=sqrt(pt*pt+am*am);
 600           if (TMath::Abs(ty)==1.) {
 601               ty=0.;
 602               Fatal("AliGenParam",
 603                     "Division by 0: Please check you rapidity range !");
 604           }
 605       //
 606       // phi
 607           //      if(!ipa)
 608           //phi=fEvPlane; //align first particle of each event with event plane
 609           //else{
 610         double v2 = fV2Para->Eval(pt);
 611         fdNdPhi->SetParameter(0,v2);
 612         fdNdPhi->SetParameter(1,fEvPlane);
 613         phi=fdNdPhi->GetRandom();
 614         //     }
 615
 616           pl=xmt*ty/sqrt((1.-ty)*(1.+ty));
 617           theta=TMath::ATan2(pt,pl);
 618       // Cut on theta
 619           if(theta<fThetaMin || theta>fThetaMax) continue;
 620           ptot=TMath::Sqrt(pt*pt+pl*pl);
 621       // Cut on momentum
 622           if(ptot<fPMin || ptot>fPMax) continue;
 623       //
 624           p[0]=pt*TMath::Cos(phi);
 625           p[1]=pt*TMath::Sin(phi);
 626           p[2]=pl;
 627
 628           if(fVertexSmear==kPerTrack) {
 629               Rndm(random,6);
 630               for (j=0;j<3;j++) {
 631                   origin0[j]=
 632                       fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 633                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 634               }
 635               Rndm(random,2);
 636               time0 = fTimeOrigin + fOsigma[2]/TMath::Ccgs()*
 637                 TMath::Cos(2*random[0]*TMath::Pi())*
 638                 TMath::Sqrt(-2*TMath::Log(random[1]));
 639           }
 640
 641       // Looking at fForceDecay :
 642       // if fForceDecay != none Primary particle decays using
 643       // AliPythia and children are tracked by GEANT
 644       //
 645       // if fForceDecay == none Primary particle is tracked by GEANT
 646       // (In the latest, make sure that GEANT actually does all the decays you want)
 647       //
 648           Bool_t decayed = kFALSE;
 649
 650
 651           if (fForceDecay != kNoDecay) {
 652         // Using lujet to decay particle
 653               Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
 654               TLorentzVector pmom(p[0], p[1], p[2], energy);
 655               fDecayer->Decay(iPart,&pmom);
 656         //
 657         // select decay particles
 658               Int_t np=fDecayer->ImportParticles(particles);
 659
 660         iPart=iTemp;
 661         if(iPart==220000){
 662           TParticle *gamma = (TParticle *)particles->At(0);
 663           gamma->SetPdgCode(iPart);
 664           np=VirtualGammaPairProduction(particles,np);
 665         }
 666         if(fForceConv) np=ForceGammaConversion(particles,np);
 667
 668               //  Selecting  GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut;
 669               if (fGeometryAcceptance)
 670                 if (!CheckAcceptanceGeometry(np,particles)) continue;
 671               Int_t ncsel=0;
 672               Int_t* pFlag      = new Int_t[np];
 673               Int_t* pParent    = new Int_t[np];
 674               Int_t* pSelected  = new Int_t[np];
 675               Int_t* trackIt    = new Int_t[np];
 676
 677               for (i=0; i<np; i++) {
 678                   pFlag[i]     =  0;
 679                   pSelected[i] =  0;
 680                   pParent[i]   = -1;
 681               }
 682
 683               if (np >1) {
 684                   decayed = kTRUE;
 685                   TParticle* iparticle =  0;
 686                   Int_t ipF, ipL;
 687                   for (i = 1; i<np ; i++) {
 688                       trackIt[i] = 1;
 689                       iparticle = (TParticle *) particles->At(i);
 690                       Int_t kf = iparticle->GetPdgCode();
 691                       Int_t ks = iparticle->GetStatusCode();
 692             // flagged particle
 693
 694                       if (pFlag[i] == 1) {
 695                           ipF = iparticle->GetFirstDaughter();
 696                           ipL = iparticle->GetLastDaughter();
 697                           if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
 698                           continue;
 699                       }
 700
 701             // flag decay products of particles with long life-time (c tau > .3 mum)
 702
 703                       if (ks != 1) {
 704               //                          TParticlePDG *particle = pDataBase->GetParticle(kf);
 705
 706                           Double_t lifeTime = fDecayer->GetLifetime(kf);
 707               //                          Double_t mass     = particle->Mass();
 708               //                          Double_t width    = particle->Width();
 709                           if (lifeTime > (Double_t) fMaxLifeTime) {
 710                               ipF = iparticle->GetFirstDaughter();
 711                               ipL = iparticle->GetLastDaughter();
 712                               if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
 713                           } else{
 714                               trackIt[i]     = 0;
 715                               pSelected[i]   = 1;
 716                           }
 717                       } // ks==1 ?
 718             //
 719             // children
 720
 721                       if ((ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll || fSelectAll) && trackIt[i])
 722                       {
 723                           if (fCutOnChild) {
 724                               pc[0]=iparticle->Px();
 725                               pc[1]=iparticle->Py();
 726                               pc[2]=iparticle->Pz();
 727                               Bool_t  childok = KinematicSelection(iparticle, 1);
 728                               if(childok) {
 729                                   pSelected[i]  = 1;
 730                                   ncsel++;
 731                               } else {
 732                                   ncsel=-1;
 733                                   break;
 734                               } // child kine cuts
 735                           } else {
 736                               pSelected[i]  = 1;
 737                               ncsel++;
 738                           } // if child selection
 739                       } // select muon
 740                   } // decay particle loop
 741               } // if decay products
 742
 743               Int_t iparent;
 744               if ((fCutOnChild && ncsel >0) || !fCutOnChild){
 745                   ipa++;
 746           //
 747           // Parent
 748
 749
 750                   PushTrack(0, -1, iPart, p, origin0, polar, time0, kPPrimary, nt, wgtp, ((decayed)? 11 : 1));
 751                   pParent[0] = nt;
 752                   KeepTrack(nt);
 753                   fNprimaries++;
 754
 755           //
 756           // Decay Products
 757           //
 758                   for (i = 1; i < np; i++) {
 759                       if (pSelected[i]) {
 760                           TParticle* iparticle = (TParticle *) particles->At(i);
 761                           Int_t kf   = iparticle->GetPdgCode();
 762                           Int_t ksc  = iparticle->GetStatusCode();
 763                           Int_t jpa  = iparticle->GetFirstMother()-1;
 764
 765                           och[0] = origin0[0]+iparticle->Vx();
 766                           och[1] = origin0[1]+iparticle->Vy();
 767                           och[2] = origin0[2]+iparticle->Vz();
 768                           pc[0]  = iparticle->Px();
 769                           pc[1]  = iparticle->Py();
 770                           pc[2]  = iparticle->Pz();
 771
 772                           if (jpa > -1) {
 773                               iparent = pParent[jpa];
 774                           } else {
 775                               iparent = -1;
 776                           }
 777
 778                           PushTrack(fTrackIt * trackIt[i], iparent, kf,
 779                                            pc, och, polar,
 780                                            time0 + iparticle->T(), kPDecay, nt, wgtch, ksc);
 781                           pParent[i] = nt;
 782                           KeepTrack(nt);
 783                           fNprimaries++;
 784                       } // Selected
 785                   } // Particle loop
 786               }  // Decays by Lujet
 787               particles->Clear();
 788               if (pFlag)      delete[] pFlag;
 789               if (pParent)    delete[] pParent;
 790               if (pSelected)  delete[] pSelected;
 791               if (trackIt)    delete[] trackIt;
 792           } // kinematic selection
 793           else  // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
 794           {
 795             gAlice->GetMCApp()->
 796                 PushTrack(fTrackIt,-1,iPart,p,origin0,polar,time0,kPPrimary,nt,wgtp, 1);
 797             ipa++;
 798             fNprimaries++;
 799           }
 800           break;
 801     } // while
 802   } // event loop
 803
 804   SetHighWaterMark(nt);
 805
 806   AliGenEventHeader* header = new AliGenEventHeader("PARAM");
 807   header->SetPrimaryVertex(fVertex);
 808   header->SetInteractionTime(fTime);
 809   header->SetNProduced(fNprimaries);
 810   AddHeader(header);
 811 }
 812 //____________________________________________________________________________________
 813 Float_t AliGenParam::GetRelativeArea(Float_t ptMin, Float_t ptMax, Float_t yMin, Float_t yMax, Float_t phiMin, Float_t phiMax)
 814 {
 815   //
 816   // Normalisation for selected kinematic region
 817   //
 818 #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
 819   Float_t ratio =
 820     fPtPara->Integral(ptMin,ptMax,(Double_t *)0,1.e-6) / fPtPara->Integral( fPtPara->GetXmin(), fPtPara->GetXmax(),(Double_t *)0,1.e-6) *
 821     fYPara->Integral(yMin,yMax,(Double_t *)0,1.e-6)/fYPara->Integral(fYPara->GetXmin(),fYPara->GetXmax(),(Double_t *)0,1.e-6)   *
 822     (phiMax-phiMin)/360.;
 823 #else
 824   Float_t ratio =
 825     fPtPara->Integral(ptMin,ptMax,1.e-6) / fPtPara->Integral( fPtPara->GetXmin(), fPtPara->GetXmax(),1.e-6) *
 826     fYPara->Integral(yMin,yMax,1.e-6)/fYPara->Integral(fYPara->GetXmin(),fYPara->GetXmax(),1.e-6)   *
 827     (phiMax-phiMin)/360.;
 828 #endif
 829   return TMath::Abs(ratio);
 830 }
 831
 832 //____________________________________________________________________________________
 833
 834 void AliGenParam::Draw( const char * /*opt*/)
 835 {
 836     //
 837     // Draw the pT and y Distributions
 838     //
 839      TCanvas *c0 = new TCanvas("c0","Canvas 0",400,10,600,700);
 840      c0->Divide(2,1);
 841      c0->cd(1);
 842      fPtPara->Draw();
 843      fPtPara->GetHistogram()->SetXTitle("p_{T} (GeV)");
 844      c0->cd(2);
 845      fYPara->Draw();
 846      fYPara->GetHistogram()->SetXTitle("y");
 847 }
 848
 849
 850
 851