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