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