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