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