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