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