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