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 /* $Id$ */
  17
  18 // Class to generate particles from using paramtrized pT and y distributions.
  19 // Distributions are obtained from pointer to object of type AliGenLib.
  20 // (For example AliGenMUONlib)
  21 // Decays are performed using Pythia.
  22 // andreas.morsch@cern.ch
  23
  24 #include <TCanvas.h>
  25 #include <TClonesArray.h>
  26 #include <TDatabasePDG.h>
  27 #include <TF1.h>
  28 #include <TH1F.h>
  29 #include <TLorentzVector.h>
  30 #include <TMath.h>
  31 #include <TParticle.h>
  32 #include <TParticlePDG.h>
  33 #include <TROOT.h>
  34 #include <TVirtualMC.h>
  35
  36 #include "AliDecayer.h"
  37 #include "AliGenMUONlib.h"
  38 #include "AliGenParam.h"
  39 #include "AliMC.h"
  40 #include "AliRun.h"
  41
  42 ClassImp(AliGenParam)
  43
  44 //------------------------------------------------------------
  45
  46   //Begin_Html
  47   /*
  48     <img src="picts/AliGenParam.gif">
  49   */
  50   //End_Html
  51
  52 //____________________________________________________________
  53     AliGenParam::AliGenParam():
  54         fPtParaFunc(0),
  55         fYParaFunc(0),
  56         fIpParaFunc(0),
  57         fPtPara(0),
  58         fYPara(0),
  59         fParam(0),
  60         fdNdy0(0.),
  61         fYWgt(0.),
  62         fPtWgt(0.),
  63         fBias(0.),
  64         fTrials(0),
  65         fDeltaPt(0.01),
  66         fDecayer(0)
  67 {
  68 // Default constructor
  69 }
  70 //____________________________________________________________
  71 AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library,  Int_t param, char* tname)
  72     :AliGenMC(npart),
  73      fPtParaFunc(Library->GetPt(param, tname)),
  74      fYParaFunc (Library->GetY (param, tname)),
  75      fIpParaFunc(Library->GetIp(param, tname)),
  76      fPtPara(0),
  77      fYPara(0),
  78      fParam(param),
  79      fdNdy0(0.),
  80      fYWgt(0.),
  81      fPtWgt(0.),
  82      fBias(0.),
  83      fTrials(0),
  84      fDeltaPt(0.01),
  85      fDecayer(0)
  86 {
  87 // Constructor using number of particles parameterisation id and library
  88     fName = "Param";
  89     fTitle= "Particle Generator using pT and y parameterisation";
  90     fAnalog = kAnalog;
  91     SetForceDecay();
  92 }
  93 //____________________________________________________________
  94 AliGenParam::AliGenParam(Int_t npart, Int_t param, const char* tname, const char* name):
  95     AliGenMC(npart),
  96     fPtParaFunc(0),
  97     fYParaFunc (0),
  98     fIpParaFunc(0),
  99     fPtPara(0),
 100     fYPara(0),
 101     fParam(param),
 102     fdNdy0(0.),
 103     fYWgt(0.),
 104     fPtWgt(0.),
 105     fBias(0.),
 106     fTrials(0),
 107     fDeltaPt(0.01),
 108     fDecayer(0)
 109 {
 110 // Constructor using parameterisation id and number of particles
 111 //
 112     fName = name;
 113     fTitle= "Particle Generator using pT and y parameterisation";
 114
 115     AliGenLib* pLibrary = new AliGenMUONlib();
 116     fPtParaFunc = pLibrary->GetPt(param, tname);
 117     fYParaFunc  = pLibrary->GetY (param, tname);
 118     fIpParaFunc = pLibrary->GetIp(param, tname);
 119
 120     fAnalog = kAnalog;
 121     fChildSelect.Set(5);
 122     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
 123     SetForceDecay();
 124     SetCutOnChild();
 125     SetChildMomentumRange();
 126     SetChildPtRange();
 127     SetChildPhiRange();
 128     SetChildThetaRange();
 129 }
 130 //____________________________________________________________
 131
 132 AliGenParam::AliGenParam(Int_t npart, Int_t param,
 133                          Double_t (*PtPara) (Double_t*, Double_t*),
 134                          Double_t (*YPara ) (Double_t* ,Double_t*),
 135                          Int_t    (*IpPara) (TRandom *))
 136     :AliGenMC(npart),
 137
 138      fPtParaFunc(PtPara),
 139      fYParaFunc(YPara),
 140      fIpParaFunc(IpPara),
 141      fPtPara(0),
 142      fYPara(0),
 143      fParam(param),
 144      fdNdy0(0.),
 145      fYWgt(0.),
 146      fPtWgt(0.),
 147      fBias(0.),
 148      fTrials(0),
 149      fDeltaPt(0.01),
 150      fDecayer(0)
 151 {
 152 // Constructor
 153 // Gines Martinez 1/10/99
 154     fName = "Param";
 155     fTitle= "Particle Generator using pT and y parameterisation";
 156
 157     fAnalog = kAnalog;
 158     fChildSelect.Set(5);
 159     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
 160     SetForceDecay();
 161     SetCutOnChild();
 162     SetChildMomentumRange();
 163     SetChildPtRange();
 164     SetChildPhiRange();
 165     SetChildThetaRange();
 166 }
 167
 168 //____________________________________________________________
 169 AliGenParam::~AliGenParam()
 170 {
 171 // Destructor
 172     delete  fPtPara;
 173     delete  fYPara;
 174 }
 175
 176 //____________________________________________________________
 177 void AliGenParam::Init()
 178 {
 179 // Initialisation
 180
 181     if (gMC) fDecayer = gMC->GetDecayer();
 182   //Begin_Html
 183   /*
 184     <img src="picts/AliGenParam.gif">
 185   */
 186   //End_Html
 187     char name[256];
 188     sprintf(name, "pt-parameterisation for %s", GetName());
 189
 190     if (fPtPara) fPtPara->Delete();
 191     fPtPara = new TF1(name, fPtParaFunc, fPtMin, fPtMax,0);
 192     gROOT->GetListOfFunctions()->Remove(fPtPara);
 193 //  Set representation precision to 10 MeV
 194     Int_t npx= Int_t((fPtMax - fPtMin) / fDeltaPt);
 195
 196     fPtPara->SetNpx(npx);
 197
 198     sprintf(name, "y-parameterisation  for %s", GetName());
 199     if (fYPara) fYPara->Delete();
 200     fYPara  = new TF1(name, fYParaFunc, fYMin, fYMax, 0);
 201     gROOT->GetListOfFunctions()->Remove(fYPara);
 202
 203
 204     sprintf(name, "pt-for-%s", GetName());
 205     TF1 ptPara(name ,fPtParaFunc, 0, 15, 0);
 206     sprintf(name, "y-for-%s", GetName());
 207     TF1 yPara(name, fYParaFunc, -6, 6, 0);
 208
 209 //
 210 // dN/dy| y=0
 211     Double_t y1=0;
 212     Double_t y2=0;
 213
 214     fdNdy0=fYParaFunc(&y1,&y2);
 215 //
 216 // Integral over generation region
 217     Float_t intYS  = yPara.Integral(fYMin, fYMax,(Double_t*) 0x0,1.e-6);
 218     Float_t intPt0 = ptPara.Integral(0,15,(Double_t *) 0x0,1.e-6);
 219     Float_t intPtS = ptPara.Integral(fPtMin,fPtMax,(Double_t*) 0x0,1.e-6);
 220     Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();
 221     if (fAnalog == kAnalog) {
 222         fYWgt  = intYS/fdNdy0;
 223         fPtWgt = intPtS/intPt0;
 224         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
 225     } else {
 226         fYWgt = intYS/fdNdy0;
 227         fPtWgt = (fPtMax-fPtMin)/intPt0;
 228         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
 229     }
 230 //
 231 // particle decay related initialization
 232     fDecayer->SetForceDecay(fForceDecay);
 233     fDecayer->Init();
 234
 235 //
 236     AliGenMC::Init();
 237 }
 238
 239 //____________________________________________________________
 240 void AliGenParam::Generate()
 241 {
 242 //
 243 // Generate 'npart' of light and heavy mesons (J/Psi, upsilon or phi, Pion,
 244 // Kaons, Etas, Omegas) and Baryons (proton, antiprotons, neutrons and
 245 // antineutrons in the the desired theta, phi and momentum windows;
 246 // Gaussian smearing on the vertex is done if selected.
 247 // The decay of heavy mesons is done using lujet,
 248 //    and the childern particle are tracked by GEANT
 249 // However, light mesons are directly tracked by GEANT
 250 // setting fForceDecay = nodecay (SetForceDecay(nodecay))
 251 //
 252 //
 253 //  Reinitialize decayer
 254   fDecayer->SetForceDecay(fForceDecay);
 255   fDecayer->Init();
 256
 257 //
 258   Float_t polar[3]= {0,0,0};  // Polarisation of the parent particle (for GEANT tracking)
 259   Float_t origin0[3];         // Origin of the generated parent particle (for GEANT tracking)
 260   Float_t pt, pl, ptot;       // Transverse, logitudinal and total momenta of the parent particle
 261   Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
 262   Float_t p[3], pc[3],
 263           och[3];             // Momentum, polarisation and origin of the children particles from lujet
 264   Double_t ty, xmt;
 265   Int_t nt, i, j;
 266   Float_t  wgtp, wgtch;
 267   Double_t dummy;
 268   static TClonesArray *particles;
 269   //
 270   if(!particles) particles = new TClonesArray("TParticle",1000);
 271
 272   TDatabasePDG *pDataBase = TDatabasePDG::Instance();
 273   //
 274   Float_t random[6];
 275
 276 // Calculating vertex position per event
 277   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
 278   if(fVertexSmear==kPerEvent) {
 279       Vertex();
 280       for (j=0;j<3;j++) origin0[j]=fVertex[j];
 281   }
 282
 283   Int_t ipa=0;
 284
 285 // Generating fNpart particles
 286   while (ipa<fNpart) {
 287       while(1) {
 288 //
 289 // particle type
 290           Int_t iPart = fIpParaFunc(fRandom);
 291           fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;
 292           TParticlePDG *particle = pDataBase->GetParticle(iPart);
 293           Float_t am = particle->Mass();
 294
 295           Rndm(random,2);
 296 //
 297 // phi
 298           phi=fPhiMin+random[0]*(fPhiMax-fPhiMin);
 299 //
 300 // y
 301           ty = TMath::TanH(fYPara->GetRandom());
 302 //
 303 // pT
 304           if (fAnalog == kAnalog) {
 305               pt=fPtPara->GetRandom();
 306               wgtp=fParentWeight;
 307               wgtch=fChildWeight;
 308           } else {
 309               pt=fPtMin+random[1]*(fPtMax-fPtMin);
 310               Double_t ptd=pt;
 311               wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
 312               wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
 313           }
 314           xmt=sqrt(pt*pt+am*am);
 315           if (TMath::Abs(ty)==1.) {
 316               ty=0.;
 317               Fatal("AliGenParam",
 318                     "Division by 0: Please check you rapidity range !");
 319           }
 320
 321           pl=xmt*ty/sqrt(1.-ty*ty);
 322           theta=TMath::ATan2(pt,pl);
 323 // Cut on theta
 324           if(theta<fThetaMin || theta>fThetaMax) continue;
 325           ptot=TMath::Sqrt(pt*pt+pl*pl);
 326 // Cut on momentum
 327           if(ptot<fPMin || ptot>fPMax) continue;
 328 //
 329           p[0]=pt*TMath::Cos(phi);
 330           p[1]=pt*TMath::Sin(phi);
 331           p[2]=pl;
 332           if(fVertexSmear==kPerTrack) {
 333               Rndm(random,6);
 334               for (j=0;j<3;j++) {
 335                   origin0[j]=
 336                       fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 337                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 338               }
 339           }
 340
 341 // Looking at fForceDecay :
 342 // if fForceDecay != none Primary particle decays using
 343 // AliPythia and children are tracked by GEANT
 344 //
 345 // if fForceDecay == none Primary particle is tracked by GEANT
 346 // (In the latest, make sure that GEANT actually does all the decays you want)
 347 //
 348
 349           if (fForceDecay != kNoDecay) {
 350 // Using lujet to decay particle
 351               Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
 352               TLorentzVector pmom(p[0], p[1], p[2], energy);
 353               fDecayer->Decay(iPart,&pmom);
 354 //
 355 // select decay particles
 356               Int_t np=fDecayer->ImportParticles(particles);
 357
 358               //  Selecting  GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut;
 359               if (fGeometryAcceptance)
 360                 if (!CheckAcceptanceGeometry(np,particles)) continue;
 361               Int_t ncsel=0;
 362               Int_t* pFlag      = new Int_t[np];
 363               Int_t* pParent    = new Int_t[np];
 364               Int_t* pSelected  = new Int_t[np];
 365               Int_t* trackIt    = new Int_t[np];
 366
 367               for (i=0; i<np; i++) {
 368                   pFlag[i]     =  0;
 369                   pSelected[i] =  0;
 370                   pParent[i]   = -1;
 371               }
 372
 373               if (np >1) {
 374                   TParticle* iparticle =  (TParticle *) particles->At(0);
 375                   Int_t ipF, ipL;
 376                   for (i = 1; i<np ; i++) {
 377                       trackIt[i] = 1;
 378                       iparticle = (TParticle *) particles->At(i);
 379                       Int_t kf = iparticle->GetPdgCode();
 380                       Int_t ks = iparticle->GetStatusCode();
 381 // flagged particle
 382
 383                       if (pFlag[i] == 1) {
 384                           ipF = iparticle->GetFirstDaughter();
 385                           ipL = iparticle->GetLastDaughter();
 386                           if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
 387                           continue;
 388                       }
 389
 390 // flag decay products of particles with long life-time (c tau > .3 mum)
 391
 392                       if (ks != 1) {
 393 //                        TParticlePDG *particle = pDataBase->GetParticle(kf);
 394
 395                           Double_t lifeTime = fDecayer->GetLifetime(kf);
 396 //                        Double_t mass     = particle->Mass();
 397 //                        Double_t width    = particle->Width();
 398                           if (lifeTime > (Double_t) fMaxLifeTime) {
 399                               ipF = iparticle->GetFirstDaughter();
 400                               ipL = iparticle->GetLastDaughter();
 401                               if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
 402                           } else{
 403                               trackIt[i]     = 0;
 404                               pSelected[i]   = 1;
 405                           }
 406                       } // ks==1 ?
 407 //
 408 // children
 409
 410                       if (ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll && trackIt[i])
 411                       {
 412                           if (fCutOnChild) {
 413                               pc[0]=iparticle->Px();
 414                               pc[1]=iparticle->Py();
 415                               pc[2]=iparticle->Pz();
 416                               Bool_t  childok = KinematicSelection(iparticle, 1);
 417                               if(childok) {
 418                                   pSelected[i]  = 1;
 419                                   ncsel++;
 420                               } else {
 421                                   ncsel=-1;
 422                                   break;
 423                               } // child kine cuts
 424                           } else {
 425                               pSelected[i]  = 1;
 426                               ncsel++;
 427                           } // if child selection
 428                       } // select muon
 429                   } // decay particle loop
 430               } // if decay products
 431
 432               Int_t iparent;
 433               if ((fCutOnChild && ncsel >0) || !fCutOnChild){
 434                   ipa++;
 435 //
 436 // Parent
 437                   PushTrack(0, -1, iPart, p, origin0, polar, 0, kPPrimary, nt, wgtp);
 438                   pParent[0] = nt;
 439                   KeepTrack(nt);
 440 //
 441 // Decay Products
 442 //
 443                   for (i = 1; i < np; i++) {
 444                       if (pSelected[i]) {
 445                           TParticle* iparticle = (TParticle *) particles->At(i);
 446                           Int_t kf  = iparticle->GetPdgCode();
 447                           Int_t ipa = iparticle->GetFirstMother()-1;
 448
 449                           och[0] = origin0[0]+iparticle->Vx()/10;
 450                           och[1] = origin0[1]+iparticle->Vy()/10;
 451                           och[2] = origin0[2]+iparticle->Vz()/10;
 452                           pc[0]  = iparticle->Px();
 453                           pc[1]  = iparticle->Py();
 454                           pc[2]  = iparticle->Pz();
 455
 456                           if (ipa > -1) {
 457                               iparent = pParent[ipa];
 458                           } else {
 459                               iparent = -1;
 460                           }
 461
 462                           PushTrack(fTrackIt*trackIt[i], iparent, kf,
 463                                            pc, och, polar,
 464                                            0, kPDecay, nt, wgtch);
 465                           pParent[i] = nt;
 466                           KeepTrack(nt);
 467                       } // Selected
 468                   } // Particle loop
 469               }  // Decays by Lujet
 470               particles->Clear();
 471               if (pFlag)      delete[] pFlag;
 472               if (pParent)    delete[] pParent;
 473               if (pSelected)  delete[] pSelected;
 474               if (trackIt)    delete[] trackIt;
 475           } // kinematic selection
 476           else  // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
 477           {
 478             gAlice->GetMCApp()->
 479                 PushTrack(fTrackIt,-1,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp);
 480             ipa++;
 481           }
 482           break;
 483     } // while
 484   } // event loop
 485   SetHighWaterMark(nt);
 486 }
 487 //____________________________________________________________________________________
 488 Float_t AliGenParam::GetRelativeArea(Float_t ptMin, Float_t ptMax, Float_t yMin, Float_t yMax, Float_t phiMin, Float_t phiMax)
 489 {
 490 //
 491 // Normalisation for selected kinematic region
 492 //
 493   Float_t ratio =
 494     fPtPara->Integral(ptMin,ptMax,(Double_t *)0,1.e-6) / fPtPara->Integral( fPtPara->GetXmin(), fPtPara->GetXmax(),(Double_t *)0,1.e-6) *
 495     fYPara->Integral(yMin,yMax,(Double_t *)0,1.e-6)/fYPara->Integral(fYPara->GetXmin(),fYPara->GetXmax(),(Double_t *)0,1.e-6)   *
 496     (phiMax-phiMin)/360.;
 497   return TMath::Abs(ratio);
 498 }
 499
 500 //____________________________________________________________________________________
 501
 502 void AliGenParam::Draw( const char * /*opt*/)
 503 {
 504     //
 505     // Draw the pT and y Distributions
 506     //
 507      TCanvas *c0 = new TCanvas("c0","Canvas 0",400,10,600,700);
 508      c0->Divide(2,1);
 509      c0->cd(1);
 510      fPtPara->Draw();
 511      fPtPara->GetHistogram()->SetXTitle("p_{T} (GeV)");
 512      c0->cd(2);
 513      fYPara->Draw();
 514      fYPara->GetHistogram()->SetXTitle("y");
 515 }
 516
 517
 518
 519