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 'npart' of light and heavy mesons (J/Psi, upsilon or phi, Pion,
 278   // Kaons, Etas, Omegas) and Baryons (proton, antiprotons, neutrons and
 279   // antineutrons in the the desired theta, phi and momentum windows;
 280   // Gaussian smearing on the vertex is done if selected.
 281   // The decay of heavy mesons is done using lujet,
 282   //    and the childern particle are tracked by GEANT
 283   // However, light mesons are directly tracked by GEANT
 284   // setting fForceDecay = nodecay (SetForceDecay(nodecay))
 285   //
 286   //
 287   //  Reinitialize decayer
 288   fDecayer->SetForceDecay(fForceDecay);
 289   fDecayer->Init();
 290
 291   //
 292   Float_t polar[3]= {0,0,0};  // Polarisation of the parent particle (for GEANT tracking)
 293   Float_t origin0[3];         // Origin of the generated parent particle (for GEANT tracking)
 294   Float_t time0;              // Time0 of the generated parent particle
 295   Float_t pt, pl, ptot;       // Transverse, logitudinal and total momenta of the parent particle
 296   Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
 297   Float_t p[3], pc[3],
 298           och[3];             // Momentum, polarisation and origin of the children particles from lujet
 299   Double_t ty, xmt;
 300   Int_t nt, i, j;
 301   Float_t  wgtp, wgtch;
 302   Double_t dummy;
 303   static TClonesArray *particles;
 304   //
 305   if(!particles) particles = new TClonesArray("TParticle",1000);
 306
 307   TDatabasePDG *pDataBase = TDatabasePDG::Instance();
 308   //
 309   Float_t random[6];
 310
 311   // Calculating vertex position per event
 312   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
 313   time0 = fTimeOrigin;
 314   if(fVertexSmear==kPerEvent) {
 315       Vertex();
 316       for (j=0;j<3;j++) origin0[j]=fVertex[j];
 317       time0 = fTime;
 318   }
 319
 320   Int_t ipa=0;
 321
 322   // Generating fNpart particles
 323   fNprimaries = 0;
 324
 325   while (ipa<fNpart) {
 326       while(1) {
 327       //
 328       // particle type
 329           Int_t iPart = fIpParaFunc(fRandom);
 330           fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;
 331           TParticlePDG *particle = pDataBase->GetParticle(iPart);
 332           Float_t am = particle->Mass();
 333
 334           Rndm(random,2);
 335       //
 336       // y
 337           ty = TMath::TanH(fYPara->GetRandom());
 338       //
 339       // pT
 340           if (fAnalog == kAnalog) {
 341               pt=fPtPara->GetRandom();
 342               wgtp=fParentWeight;
 343               wgtch=fChildWeight;
 344           } else {
 345               pt=fPtMin+random[1]*(fPtMax-fPtMin);
 346               Double_t ptd=pt;
 347               wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
 348               wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
 349           }
 350           xmt=sqrt(pt*pt+am*am);
 351           if (TMath::Abs(ty)==1.) {
 352               ty=0.;
 353               Fatal("AliGenParam",
 354                     "Division by 0: Please check you rapidity range !");
 355           }
 356       //
 357       // phi
 358       // if(!ipa)
 359       //        phi=fEvPlane; //align first particle of each event with event plane
 360       // else{
 361         double v2 = fV2Para->Eval(pt);
 362         fdNdPhi->SetParameter(0,v2);
 363         fdNdPhi->SetParameter(1,fEvPlane);
 364         phi=fdNdPhi->GetRandom();
 365       // }
 366
 367           pl=xmt*ty/sqrt((1.-ty)*(1.+ty));
 368           theta=TMath::ATan2(pt,pl);
 369       // Cut on theta
 370           if(theta<fThetaMin || theta>fThetaMax) continue;
 371           ptot=TMath::Sqrt(pt*pt+pl*pl);
 372       // Cut on momentum
 373           if(ptot<fPMin || ptot>fPMax) continue;
 374       //
 375           p[0]=pt*TMath::Cos(phi);
 376           p[1]=pt*TMath::Sin(phi);
 377           p[2]=pl;
 378           if(fVertexSmear==kPerTrack) {
 379               Rndm(random,6);
 380               for (j=0;j<3;j++) {
 381                   origin0[j]=
 382                       fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 383                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 384               }
 385               Rndm(random,2);
 386               time0 = fTimeOrigin + fOsigma[2]/TMath::Ccgs()*
 387                 TMath::Cos(2*random[0]*TMath::Pi())*
 388                 TMath::Sqrt(-2*TMath::Log(random[1]));
 389           }
 390
 391       // Looking at fForceDecay :
 392       // if fForceDecay != none Primary particle decays using
 393       // AliPythia and children are tracked by GEANT
 394       //
 395       // if fForceDecay == none Primary particle is tracked by GEANT
 396       // (In the latest, make sure that GEANT actually does all the decays you want)
 397       //
 398           Bool_t decayed = kFALSE;
 399
 400
 401           if (fForceDecay != kNoDecay) {
 402         // Using lujet to decay particle
 403               Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
 404               TLorentzVector pmom(p[0], p[1], p[2], energy);
 405               fDecayer->Decay(iPart,&pmom);
 406         //
 407         // select decay particles
 408               Int_t np=fDecayer->ImportParticles(particles);
 409
 410               //  Selecting  GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut;
 411               if (fGeometryAcceptance)
 412                 if (!CheckAcceptanceGeometry(np,particles)) continue;
 413               Int_t ncsel=0;
 414               Int_t* pFlag      = new Int_t[np];
 415               Int_t* pParent    = new Int_t[np];
 416               Int_t* pSelected  = new Int_t[np];
 417               Int_t* trackIt    = new Int_t[np];
 418
 419               for (i=0; i<np; i++) {
 420                   pFlag[i]     =  0;
 421                   pSelected[i] =  0;
 422                   pParent[i]   = -1;
 423               }
 424
 425               if (np >1) {
 426                   decayed = kTRUE;
 427                   TParticle* iparticle =  0;
 428                   Int_t ipF, ipL;
 429                   for (i = 1; i<np ; i++) {
 430                       trackIt[i] = 1;
 431                       iparticle = (TParticle *) particles->At(i);
 432                       Int_t kf = iparticle->GetPdgCode();
 433                       Int_t ks = iparticle->GetStatusCode();
 434             // flagged particle
 435
 436                       if (pFlag[i] == 1) {
 437                           ipF = iparticle->GetFirstDaughter();
 438                           ipL = iparticle->GetLastDaughter();
 439                           if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
 440                           continue;
 441                       }
 442
 443             // flag decay products of particles with long life-time (c tau > .3 mum)
 444
 445                       if (ks != 1) {
 446               //                          TParticlePDG *particle = pDataBase->GetParticle(kf);
 447
 448                           Double_t lifeTime = fDecayer->GetLifetime(kf);
 449               //                          Double_t mass     = particle->Mass();
 450               //                          Double_t width    = particle->Width();
 451                           if (lifeTime > (Double_t) fMaxLifeTime) {
 452                               ipF = iparticle->GetFirstDaughter();
 453                               ipL = iparticle->GetLastDaughter();
 454                               if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
 455                           } else{
 456                               trackIt[i]     = 0;
 457                               pSelected[i]   = 1;
 458                           }
 459                       } // ks==1 ?
 460             //
 461             // children
 462
 463                       if ((ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll || fSelectAll) && trackIt[i])
 464                       {
 465                           if (fCutOnChild) {
 466                               pc[0]=iparticle->Px();
 467                               pc[1]=iparticle->Py();
 468                               pc[2]=iparticle->Pz();
 469                               Bool_t  childok = KinematicSelection(iparticle, 1);
 470                               if(childok) {
 471                                   pSelected[i]  = 1;
 472                                   ncsel++;
 473                               } else {
 474                                   ncsel=-1;
 475                                   break;
 476                               } // child kine cuts
 477                           } else {
 478                               pSelected[i]  = 1;
 479                               ncsel++;
 480                           } // if child selection
 481                       } // select muon
 482                   } // decay particle loop
 483               } // if decay products
 484
 485               Int_t iparent;
 486               if ((fCutOnChild && ncsel >0) || !fCutOnChild){
 487                   ipa++;
 488           //
 489           // Parent
 490
 491
 492                   PushTrack(0, -1, iPart, p, origin0, polar, time0, kPPrimary, nt, wgtp, ((decayed)? 11 : 1));
 493                   pParent[0] = nt;
 494                   KeepTrack(nt);
 495                   fNprimaries++;
 496
 497           //
 498           // Decay Products
 499           //
 500                   for (i = 1; i < np; i++) {
 501                       if (pSelected[i]) {
 502                           TParticle* iparticle = (TParticle *) particles->At(i);
 503                           Int_t kf   = iparticle->GetPdgCode();
 504                           Int_t ksc  = iparticle->GetStatusCode();
 505                           Int_t jpa  = iparticle->GetFirstMother()-1;
 506
 507                           och[0] = origin0[0]+iparticle->Vx();
 508                           och[1] = origin0[1]+iparticle->Vy();
 509                           och[2] = origin0[2]+iparticle->Vz();
 510                           pc[0]  = iparticle->Px();
 511                           pc[1]  = iparticle->Py();
 512                           pc[2]  = iparticle->Pz();
 513
 514                           if (jpa > -1) {
 515                               iparent = pParent[jpa];
 516                           } else {
 517                               iparent = -1;
 518                           }
 519
 520                           PushTrack(fTrackIt * trackIt[i], iparent, kf,
 521                                            pc, och, polar,
 522                                            time0 + iparticle->T(), kPDecay, nt, wgtch, ksc);
 523                           pParent[i] = nt;
 524                           KeepTrack(nt);
 525                           fNprimaries++;
 526                       } // Selected
 527                   } // Particle loop
 528               }  // Decays by Lujet
 529               particles->Clear();
 530               if (pFlag)      delete[] pFlag;
 531               if (pParent)    delete[] pParent;
 532               if (pSelected)  delete[] pSelected;
 533               if (trackIt)    delete[] trackIt;
 534           } // kinematic selection
 535           else  // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
 536           {
 537             gAlice->GetMCApp()->
 538                 PushTrack(fTrackIt,-1,iPart,p,origin0,polar,time0,kPPrimary,nt,wgtp, 1);
 539             ipa++;
 540             fNprimaries++;
 541           }
 542           break;
 543     } // while
 544   } // event loop
 545
 546   SetHighWaterMark(nt);
 547
 548   AliGenEventHeader* header = new AliGenEventHeader("PARAM");
 549   header->SetPrimaryVertex(fVertex);
 550   header->SetInteractionTime(fTime);
 551   header->SetNProduced(fNprimaries);
 552   AddHeader(header);
 553 }
 554 //____________________________________________________________________________________
 555 Float_t AliGenParam::GetRelativeArea(Float_t ptMin, Float_t ptMax, Float_t yMin, Float_t yMax, Float_t phiMin, Float_t phiMax)
 556 {
 557   //
 558   // Normalisation for selected kinematic region
 559   //
 560   Float_t ratio =
 561     fPtPara->Integral(ptMin,ptMax,(Double_t *)0,1.e-6) / fPtPara->Integral( fPtPara->GetXmin(), fPtPara->GetXmax(),(Double_t *)0,1.e-6) *
 562     fYPara->Integral(yMin,yMax,(Double_t *)0,1.e-6)/fYPara->Integral(fYPara->GetXmin(),fYPara->GetXmax(),(Double_t *)0,1.e-6)   *
 563     (phiMax-phiMin)/360.;
 564   return TMath::Abs(ratio);
 565 }
 566
 567 //____________________________________________________________________________________
 568
 569 void AliGenParam::Draw( const char * /*opt*/)
 570 {
 571     //
 572     // Draw the pT and y Distributions
 573     //
 574      TCanvas *c0 = new TCanvas("c0","Canvas 0",400,10,600,700);
 575      c0->Divide(2,1);
 576      c0->cd(1);
 577      fPtPara->Draw();
 578      fPtPara->GetHistogram()->SetXTitle("p_{T} (GeV)");
 579      c0->cd(2);
 580      fYPara->Draw();
 581      fYPara->GetHistogram()->SetXTitle("y");
 582 }
 583
 584
 585
 586