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