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