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 /*
  17 $Log$
  18 Revision 1.27  2001/02/02 15:21:10  morsch
  19 Set high water mark after last particle.
  20 Use Vertex() method for Vertex.
  21
  22 Revision 1.26  2000/12/21 16:24:06  morsch
  23 Coding convention clean-up
  24
  25 Revision 1.25  2000/11/30 07:12:50  alibrary
  26 Introducing new Rndm and QA classes
  27
  28 Revision 1.24  2000/10/18 19:11:27  hristov
  29 Division by zero fixed
  30
  31 Revision 1.23  2000/10/02 21:28:06  fca
  32 Removal of useless dependecies via forward declarations
  33
  34 Revision 1.22  2000/09/12 14:14:55  morsch
  35 Call fDecayer->ForceDecay() at the beginning of Generate().
  36
  37 Revision 1.21  2000/09/08 15:39:01  morsch
  38 Handle the case fForceDecay=all during the generation, i.e. select all secondaries.
  39
  40 Revision 1.20  2000/09/06 14:35:44  morsch
  41 Use AliDecayerPythia for particle decays.
  42
  43 Revision 1.19  2000/07/11 18:24:56  fca
  44 Coding convention corrections + few minor bug fixes
  45
  46 Revision 1.18  2000/06/29 21:08:27  morsch
  47 All paramatrisation libraries derive from the pure virtual base class AliGenLib.
  48 This allows to pass a pointer to a library directly to AliGenParam and avoids the
  49 use of function pointers in Config.C.
  50
  51 Revision 1.17  2000/06/09 20:33:30  morsch
  52 All coding rule violations except RS3 corrected
  53
  54 Revision 1.16  2000/05/02 07:51:31  morsch
  55 - Control precision of pT sampling TF1::SetNpx(..)
  56 - Correct initialisation of child-cuts in all constructors.
  57 - Most coding rule violations corrected.
  58
  59 Revision 1.15  2000/04/03 15:42:12  morsch
  60 Cuts on primary particles are separated from those on the decay products. Methods
  61 SetChildMomentumRange, SetChildPtRange, SetChildPhiRange, SetChildThetaRange added.
  62
  63 Revision 1.14  1999/11/09 07:38:48  fca
  64 Changes for compatibility with version 2.23 of ROOT
  65
  66 Revision 1.13  1999/11/04 11:30:31  fca
  67 Correct the logics for SetForceDecay
  68
  69 Revision 1.12  1999/11/03 17:43:20  fca
  70 New version from G.Martinez & A.Morsch
  71
  72 Revision 1.11  1999/09/29 09:24:14  fca
  73 Introduction of the Copyright and cvs Log
  74
  75 */
  76
  77 #include "AliGenParam.h"
  78 #include "AliDecayerPythia.h"
  79 #include "AliGenMUONlib.h"
  80 #include "AliRun.h"
  81 #include <TParticle.h>
  82 #include <TParticlePDG.h>
  83 #include <TDatabasePDG.h>
  84 #include <TLorentzVector.h>
  85
  86 #include <TF1.h>
  87
  88 ClassImp(AliGenParam)
  89
  90 //------------------------------------------------------------
  91
  92   //Begin_Html
  93   /*
  94     <img src="picts/AliGenParam.gif">
  95   */
  96   //End_Html
  97
  98 //____________________________________________________________
  99 //____________________________________________________________
 100 AliGenParam::AliGenParam()
 101 {
 102 // Deafault constructor
 103     fPtPara = 0;
 104     fYPara  = 0;
 105     fParam  = 0;
 106     fAnalog = kAnalog;
 107     SetCutOnChild();
 108     SetChildMomentumRange();
 109     SetChildPtRange();
 110     SetChildPhiRange();
 111     SetChildThetaRange();
 112     SetDeltaPt();
 113 //
 114 // Set random number generator
 115     sRandom = fRandom;
 116 }
 117
 118 AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library,  Int_t param, char* tname):AliGenerator(npart)
 119 {
 120 // Constructor using number of particles parameterisation id and library
 121
 122     fPtParaFunc = Library->GetPt(param, tname);
 123     fYParaFunc  = Library->GetY (param, tname);
 124     fIpParaFunc = Library->GetIp(param, tname);
 125
 126     fPtPara = 0;
 127     fYPara  = 0;
 128     fParam  = param;
 129     fAnalog = kAnalog;
 130     fChildSelect.Set(5);
 131     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
 132     SetForceDecay();
 133     SetCutOnChild();
 134     SetChildMomentumRange();
 135     SetChildPtRange();
 136     SetChildPhiRange();
 137     SetChildThetaRange();
 138     SetDeltaPt();
 139 //
 140 // Set random number generator
 141     sRandom = fRandom;
 142 }
 143
 144 //____________________________________________________________
 145
 146 AliGenParam::AliGenParam(Int_t npart, Int_t param, char* tname):AliGenerator(npart)
 147 {
 148 // Constructor using parameterisation id and number of particles
 149 //
 150     AliGenLib* pLibrary = new AliGenMUONlib();
 151
 152     fPtParaFunc = pLibrary->GetPt(param, tname);
 153     fYParaFunc  = pLibrary->GetY (param, tname);
 154     fIpParaFunc = pLibrary->GetIp(param, tname);
 155
 156     fPtPara = 0;
 157     fYPara  = 0;
 158     fParam  = param;
 159     fAnalog = kAnalog;
 160     fChildSelect.Set(5);
 161     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
 162     SetForceDecay();
 163     SetCutOnChild();
 164     SetChildMomentumRange();
 165     SetChildPtRange();
 166     SetChildPhiRange();
 167     SetChildThetaRange();
 168     SetDeltaPt();
 169 }
 170
 171 AliGenParam::AliGenParam(Int_t npart, Int_t param,
 172                          Double_t (*PtPara) (Double_t*, Double_t*),
 173                          Double_t (*YPara ) (Double_t* ,Double_t*),
 174                          Int_t    (*IpPara) (TRandom *))
 175     :AliGenerator(npart)
 176 {
 177 // Constructor
 178 // Gines Martinez 1/10/99
 179     fPtParaFunc = PtPara;
 180     fYParaFunc  = YPara;
 181     fIpParaFunc = IpPara;
 182 //
 183     fPtPara = 0;
 184     fYPara  = 0;
 185     fParam  = param;
 186     fAnalog = kAnalog;
 187     fChildSelect.Set(5);
 188     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
 189     SetForceDecay();
 190     SetCutOnChild();
 191     SetChildMomentumRange();
 192     SetChildPtRange();
 193     SetChildPhiRange();
 194     SetChildThetaRange();
 195     SetDeltaPt();
 196 }
 197
 198
 199 AliGenParam::AliGenParam(const AliGenParam & Paramd)
 200 {
 201 // copy constructor
 202 }
 203
 204 //____________________________________________________________
 205 AliGenParam::~AliGenParam()
 206 {
 207 // Destructor
 208     delete  fPtPara;
 209     delete  fYPara;
 210 }
 211
 212 //____________________________________________________________
 213 void AliGenParam::Init()
 214 {
 215 // Initialisation
 216
 217     fDecayer = new AliDecayerPythia();
 218   //Begin_Html
 219   /*
 220     <img src="picts/AliGenParam.gif">
 221   */
 222   //End_Html
 223
 224     fPtPara = new TF1("Pt-Parametrization",fPtParaFunc,fPtMin,fPtMax,0);
 225 //  Set representation precision to 10 MeV
 226     Int_t npx= Int_t((fPtMax-fPtMin)/fDeltaPt);
 227
 228     fPtPara->SetNpx(npx);
 229
 230     fYPara  = new TF1("Y -Parametrization",fYParaFunc,fYMin,fYMax,0);
 231     TF1* ptPara = new TF1("Pt-Parametrization",fPtParaFunc,0,15,0);
 232     TF1* yPara  = new TF1("Y -Parametrization",fYParaFunc,-6,6,0);
 233
 234 //
 235 // dN/dy| y=0
 236     Double_t y1=0;
 237     Double_t y2=0;
 238
 239     fdNdy0=fYParaFunc(&y1,&y2);
 240 //
 241 // Integral over generation region
 242     Float_t intYS  = yPara ->Integral(fYMin, fYMax);
 243     Float_t intPt0 = ptPara->Integral(0,15);
 244     Float_t intPtS = ptPara->Integral(fPtMin,fPtMax);
 245     Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();
 246     if (fAnalog == kAnalog) {
 247         fYWgt  = intYS/fdNdy0;
 248         fPtWgt = intPtS/intPt0;
 249         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
 250     } else {
 251         fYWgt = intYS/fdNdy0;
 252         fPtWgt = (fPtMax-fPtMin)/intPt0;
 253         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
 254     }
 255 //
 256 // particle decay related initialization
 257     fDecayer->SetForceDecay(fForceDecay);
 258     fDecayer->Init();
 259
 260 //
 261     switch (fForceDecay)
 262     {
 263     case kSemiElectronic:
 264     case kDiElectron:
 265     case kBJpsiDiElectron:
 266     case kBPsiPrimeDiElectron:
 267         fChildSelect[0]=11;
 268         break;
 269     case kSemiMuonic:
 270     case kDiMuon:
 271     case kBJpsiDiMuon:
 272     case kBPsiPrimeDiMuon:
 273         fChildSelect[0]=13;
 274         break;
 275     case kPiToMu:
 276         fChildSelect[0]=13;
 277         break;
 278     case kKaToMu:
 279         fChildSelect[0]=13;
 280         break;
 281     case kHadronicD:
 282 // Implement me !!
 283         break;
 284     case kNoDecay:
 285         break;
 286     case kAll:
 287         break;
 288     }
 289 }
 290
 291 //____________________________________________________________
 292 void AliGenParam::Generate()
 293 {
 294 //
 295 // Generate 'npart' of light and heavy mesons (J/Psi, upsilon or phi, Pion,
 296 // Kaons, Etas, Omegas) and Baryons (proton, antiprotons, neutrons and
 297 // antineutrons in the the desired theta, phi and momentum windows;
 298 // Gaussian smearing on the vertex is done if selected.
 299 // The decay of heavy mesons is done using lujet,
 300 //    and the childern particle are tracked by GEANT
 301 // However, light mesons are directly tracked by GEANT
 302 // setting fForceDecay = nodecay (SetForceDecay(nodecay))
 303 //
 304
 305     fDecayer->ForceDecay();
 306   Float_t polar[3]= {0,0,0};  // Polarisation of the parent particle (for GEANT tracking)
 307   Float_t origin0[3];         // Origin of the generated parent particle (for GEANT tracking)
 308   Float_t pt, pl, ptot;       // Transverse, logitudinal and total momenta of the parent particle
 309   Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
 310   Float_t p[3], pc[3],
 311           och[3], pch[10][3]; // Momentum, polarisation and origin of the children particles from lujet
 312   Float_t ty, xmt;
 313   Int_t nt, i, j, kfch[10];
 314   Float_t  wgtp, wgtch;
 315   Double_t dummy;
 316   static TClonesArray *particles;
 317   //
 318   if(!particles) particles=new TClonesArray("TParticle",1000);
 319
 320   static TDatabasePDG *pDataBase = new TDatabasePDG();
 321   if(!pDataBase) pDataBase = new TDatabasePDG();
 322   //
 323   Float_t random[6];
 324
 325 // Calculating vertex position per event
 326   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
 327   if(fVertexSmear==kPerEvent) {
 328 //      Rndm(random,6);
 329 //      for (j=0;j<3;j++) {
 330 //        origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 331 //            TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 332 //            TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 333 //      }
 334 //    }
 335       Vertex();
 336       for (j=0;j<3;j++) origin0[j]=fVertex[j];
 337   }
 338
 339   Int_t ipa=0;
 340 // Generating fNpart particles
 341   while (ipa<fNpart) {
 342       while(1) {
 343 //
 344 // particle type
 345           Int_t iPart = fIpParaFunc(fRandom);
 346           fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;
 347           TParticlePDG *particle = pDataBase->GetParticle(iPart);
 348           Float_t am = particle->Mass();
 349
 350           Rndm(random,2);
 351 //
 352 // phi
 353           phi=fPhiMin+random[0]*(fPhiMax-fPhiMin);
 354 //
 355 // y
 356           ty=Float_t(TMath::TanH(fYPara->GetRandom()));
 357 //
 358 // pT
 359           if (fAnalog == kAnalog) {
 360               pt=fPtPara->GetRandom();
 361               wgtp=fParentWeight;
 362               wgtch=fChildWeight;
 363           } else {
 364               pt=fPtMin+random[1]*(fPtMax-fPtMin);
 365               Double_t ptd=pt;
 366               wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
 367               wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
 368           }
 369           xmt=sqrt(pt*pt+am*am);
 370       if (TMath::Abs(ty)==1) ty=0;
 371           pl=xmt*ty/sqrt(1.-ty*ty);
 372           theta=TMath::ATan2(pt,pl);
 373 // Cut on theta
 374           if(theta<fThetaMin || theta>fThetaMax) continue;
 375           ptot=TMath::Sqrt(pt*pt+pl*pl);
 376 // Cut on momentum
 377           if(ptot<fPMin || ptot>fPMax) continue;
 378           p[0]=pt*TMath::Cos(phi);
 379           p[1]=pt*TMath::Sin(phi);
 380           p[2]=pl;
 381           if(fVertexSmear==kPerTrack) {
 382               Rndm(random,6);
 383               for (j=0;j<3;j++) {
 384                   origin0[j]=
 385                       fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 386                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 387               }
 388           }
 389
 390 // Looking at fForceDecay :
 391 // if fForceDecay != none Primary particle decays using
 392 // AliPythia and children are tracked by GEANT
 393 //
 394 // if fForceDecay == none Primary particle is tracked by GEANT
 395 // (In the latest, make sure that GEANT actually does all the decays you want)
 396 //
 397           if (fForceDecay != kNoDecay) {
 398 // Using lujet to decay particle
 399               Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
 400               TLorentzVector pmom(p[0], p[1], p[2], energy);
 401               fDecayer->Decay(iPart,&pmom);
 402 //
 403 // select decay particles
 404               Int_t np=fDecayer->ImportParticles(particles);
 405               Int_t ncsel=0;
 406               for (i = 1; i<np; i++) {
 407                   TParticle *  iparticle = (TParticle *) particles->At(i);
 408                   Int_t kf = iparticle->GetPdgCode();
 409 //
 410 // children
 411                   if (ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll)
 412                   {
 413                       pc[0]=iparticle->Px();
 414                       pc[1]=iparticle->Py();
 415                       pc[2]=iparticle->Pz();
 416                       och[0]=origin0[0]+iparticle->Vx()/10;
 417                       och[1]=origin0[1]+iparticle->Vy()/10;
 418                       och[2]=origin0[2]+iparticle->Vz()/10;
 419                       if (fCutOnChild) {
 420                           Float_t ptChild=TMath::Sqrt(pc[0]*pc[0]+pc[1]*pc[1]);
 421                           Float_t pChild=TMath::Sqrt(ptChild*ptChild+pc[2]*pc[2]);
 422                           Float_t thetaChild=TMath::ATan2(ptChild,pc[2]);
 423                           Float_t phiChild=TMath::ATan2(pc[1],pc[0]);
 424                           Bool_t childok =
 425                               ((ptChild    > fChildPtMin    && ptChild    <fChildPtMax)      &&
 426                                (pChild     > fChildPMin     && pChild     <fChildPMax)       &&
 427                                (thetaChild > fChildThetaMin && thetaChild <fChildThetaMax)   &&
 428                                (phiChild   > fChildPhiMin   && phiChild   <fChildPhiMax));
 429                           if(childok)
 430                           {
 431                               pch[ncsel][0]=pc[0];
 432                               pch[ncsel][1]=pc[1];
 433                               pch[ncsel][2]=pc[2];
 434                               kfch[ncsel]=kf;
 435                               ncsel++;
 436                           } else {
 437                               ncsel=-1;
 438                               break;
 439                           } // child kine cuts
 440                       } else {
 441                           pch[ncsel][0]=pc[0];
 442                           pch[ncsel][1]=pc[1];
 443                           pch[ncsel][2]=pc[2];
 444                           kfch[ncsel]=kf;
 445                           ncsel++;
 446                       } // if child selection
 447                   } // select muon
 448               } // decay particle loop
 449               Int_t iparent;
 450               if ((fCutOnChild && ncsel >0) || !fCutOnChild){
 451                   ipa++;
 452 //
 453 // parent
 454                   gAlice->
 455                       SetTrack(0,-1,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp);
 456                   iparent=nt;
 457                   gAlice->KeepTrack(nt);
 458                   for (i=0; i< ncsel; i++) {
 459                       gAlice->SetTrack(fTrackIt,iparent,kfch[i],
 460                                        &pch[i][0],och,polar,
 461                                        0,kPDecay,nt,wgtch);
 462                       gAlice->KeepTrack(nt);
 463                   }
 464               }  // Decays by Lujet
 465           } // kinematic selection
 466           else  // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
 467           {
 468             gAlice->
 469                 SetTrack(fTrackIt,-1,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp);
 470             ipa++;
 471           }
 472           break;
 473     } // while
 474   } // event loop
 475   gAlice->SetHighWaterMark(nt);
 476
 477 }
 478
 479 Bool_t AliGenParam::ChildSelected(Int_t ip)
 480 {
 481 // True if particle is in list of selected children
 482     for (Int_t i=0; i<5; i++)
 483     {
 484         if (fChildSelect[i]==ip) return kTRUE;
 485     }
 486     return kFALSE;
 487 }
 488
 489 Bool_t AliGenParam::KinematicSelection(TParticle *particle)
 490 {
 491 // Perform kinematic cuts
 492     Float_t px=particle->Px();
 493     Float_t py=particle->Py();
 494     Float_t pz=particle->Pz();
 495 //
 496 // momentum cut
 497     Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
 498     if (p > fPMax || p < fPMin)
 499     {
 500 //      printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
 501         return kFALSE;
 502     }
 503     Float_t pt=TMath::Sqrt(px*px+py*py);
 504
 505 //
 506 // theta cut
 507     Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(p)));
 508     if (theta > fThetaMax || theta < fThetaMin)
 509     {
 510 //      printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
 511         return kFALSE;
 512     }
 513
 514     return kTRUE;
 515 }
 516
 517
 518 AliGenParam& AliGenParam::operator=(const  AliGenParam& rhs)
 519 {
 520 // Assignment operator
 521     return *this;
 522 }
 523
 524
 525