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