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