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