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