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