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