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