EVGEN/AliGenHIJINGpara.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 /*
  17 $Log$
  18 Revision 1.18  2002/12/11 11:58:11  morsch
  19 Bug in formula for pi0 energy for decay corrected.
  20
  21 Revision 1.17  2002/12/10 17:44:57  morsch
  22 Correct mother child relation for pi0.
  23
  24 Revision 1.16  2002/11/28 11:46:15  morsch
  25 Don't track pi0 if already decayed.
  26
  27 Revision 1.15  2002/11/28 11:38:53  morsch
  28 Typo corrected.
  29
  30 Revision 1.14  2002/11/26 17:12:36  morsch
  31 Decay pi0 if requested.
  32
  33 Revision 1.13  2002/10/14 14:55:35  hristov
  34 Merging the VirtualMC branch to the main development branch (HEAD)
  35
  36 Revision 1.11.4.1  2002/07/24 08:56:28  alibrary
  37 Updating EVGEN on TVirtulaMC
  38
  39 Revision 1.12  2002/06/19 06:56:34  hristov
  40 Memory leak corrected
  41
  42 Revision 1.11  2002/03/20 10:21:13  hristov
  43 Set fPtMax to 15 GeV in order to avoid some numerical problems
  44
  45 Revision 1.10  2001/10/15 16:44:46  morsch
  46 - Possibility for vertex distribution truncation.
  47 - Write mc header with vertex position.
  48
  49 Revision 1.9  2001/07/27 17:09:36  morsch
  50 Use local SetTrack, KeepTrack and SetHighWaterMark methods
  51 to delegate either to local stack or to stack owned by AliRun.
  52 (Piotr Skowronski, A.M.)
  53
  54 Revision 1.8  2001/07/20 11:03:58  morsch
  55 Issue warning message if used outside allowed eta range (-8 to 8).
  56
  57 Revision 1.7  2001/07/17 12:41:01  morsch
  58 - Calculation of fraction of event corresponding to selected pt-range corrected
  59 (R. Turrisi)
  60 - Parent weight corrected.
  61
  62 Revision 1.6  2001/05/16 14:57:10  alibrary
  63 New files for folders and Stack
  64
  65 Revision 1.5  2000/12/21 16:24:06  morsch
  66 Coding convention clean-up
  67
  68 Revision 1.4  2000/11/30 07:12:50  alibrary
  69 Introducing new Rndm and QA classes
  70
  71 Revision 1.3  2000/10/02 21:28:06  fca
  72 Removal of useless dependecies via forward declarations
  73
  74 Revision 1.2  2000/07/11 18:24:55  fca
  75 Coding convention corrections + few minor bug fixes
  76
  77 Revision 1.1  2000/06/09 20:20:30  morsch
  78 Same class as previously in AliSimpleGen.cxx
  79 All coding rule violations except RS3 corrected (AM)
  80
  81 */
  82
  83 // Parameterisation of pi and K, eta and pt distributions
  84 // used for the ALICE TDRs.
  85 // eta: according to HIJING (shadowing + quenching)
  86 // pT : according to CDF measurement at 1.8 TeV
  87 // Author: andreas.morsch@cern.ch
  88
  89
  90 //Begin_Html
  91 /*
  92 <img src="picts/AliGeneratorClass.gif">
  93 </pre>
  94 <br clear=left>
  95 <font size=+2 color=red>
  96 <p>The responsible person for this module is
  97 <a href="mailto:andreas.morsch@cern.ch">Andreas Morsch</a>.
  98 </font>
  99 <pre>
 100 */
 101 //End_Html
 102 //                                                               //
 103 ///////////////////////////////////////////////////////////////////
 104
 105 #include <TArrayF.h>
 106 #include <TClonesArray.h>
 107 #include <TDatabasePDG.h>
 108 #include <TF1.h>
 109 #include <TParticle.h>
 110 #include <TPDGCode.h>
 111
 112 #include "AliConst.h"
 113 #include "AliDecayer.h"
 114 #include "AliDecayerPythia.h"
 115 #include "AliGenEventHeader.h"
 116 #include "AliGenHIJINGpara.h"
 117 #include "AliRun.h"
 118
 119 ClassImp(AliGenHIJINGpara)
 120
 121 AliGenHIJINGpara::AliGenHIJINGpara(const AliGenHIJINGpara & para)
 122 {
 123 // copy constructor
 124 }
 125
 126 //_____________________________________________________________________________
 127 static Double_t ptpi(Double_t *px, Double_t *)
 128 {
 129   //
 130   //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
 131   //     POWER LAW FOR PT > 500 MEV
 132   //     MT SCALING BELOW (T=160 MEV)
 133   //
 134   const Double_t kp0 = 1.3;
 135   const Double_t kxn = 8.28;
 136   const Double_t kxlim=0.5;
 137   const Double_t kt=0.160;
 138   const Double_t kxmpi=0.139;
 139   const Double_t kb=1.;
 140   Double_t y, y1, xmpi2, ynorm, a;
 141   Double_t x=*px;
 142   //
 143   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
 144   xmpi2=kxmpi*kxmpi;
 145   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
 146   a=ynorm/y1;
 147   if (x > kxlim)
 148     y=a*TMath::Power(kp0/(kp0+x),kxn);
 149   else
 150     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
 151   return y*x;
 152 }
 153
 154 //_____________________________________________________________________________
 155 static Double_t ptscal(Double_t pt, Int_t np)
 156 {
 157     //    SCALING EN MASSE PAR RAPPORT A PTPI
 158     //     MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
 159     const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
 160     //     VALUE MESON/PI AT 5 GEV
 161     const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
 162     np--;
 163     Double_t f5=TMath::Power(((
 164         sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
 165     Double_t fmax2=f5/kfmax[np];
 166     // PIONS
 167     Double_t ptpion=100.*ptpi(&pt, (Double_t*) 0);
 168     Double_t fmtscal=TMath::Power(((
 169         sqrt(pt*pt+0.018215)+2.)/ (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/
 170         fmax2;
 171     return fmtscal*ptpion;
 172 }
 173
 174 //_____________________________________________________________________________
 175 static Double_t ptka( Double_t *px, Double_t *)
 176 {
 177     //
 178     // pt parametrisation for k
 179     //
 180     return ptscal(*px,2);
 181 }
 182
 183
 184 //_____________________________________________________________________________
 185 static Double_t etapic( Double_t *py, Double_t *)
 186 {
 187   //
 188   // eta parametrisation for pi
 189   //
 190     const Double_t ka1    = 4913.;
 191     const Double_t ka2    = 1819.;
 192     const Double_t keta1  = 0.22;
 193     const Double_t keta2  = 3.66;
 194     const Double_t kdeta1 = 1.47;
 195     const Double_t kdeta2 = 1.51;
 196     Double_t y=TMath::Abs(*py);
 197     //
 198     Double_t ex1 = (y-keta1)*(y-keta1)/(2*kdeta1*kdeta1);
 199     Double_t ex2 = (y-keta2)*(y-keta2)/(2*kdeta2*kdeta2);
 200     return ka1*TMath::Exp(-ex1)+ka2*TMath::Exp(-ex2);
 201 }
 202
 203 //_____________________________________________________________________________
 204 static Double_t etakac( Double_t *py, Double_t *)
 205 {
 206     //
 207     // eta parametrisation for ka
 208     //
 209     const Double_t ka1    = 497.6;
 210     const Double_t ka2    = 215.6;
 211     const Double_t keta1  = 0.79;
 212     const Double_t keta2  = 4.09;
 213     const Double_t kdeta1 = 1.54;
 214     const Double_t kdeta2 = 1.40;
 215     Double_t y=TMath::Abs(*py);
 216     //
 217     Double_t ex1 = (y-keta1)*(y-keta1)/(2*kdeta1*kdeta1);
 218     Double_t ex2 = (y-keta2)*(y-keta2)/(2*kdeta2*kdeta2);
 219     return ka1*TMath::Exp(-ex1)+ka2*TMath::Exp(-ex2);
 220 }
 221
 222 //_____________________________________________________________________________
 223 AliGenHIJINGpara::AliGenHIJINGpara()
 224   :AliGenerator()
 225 {
 226     //
 227     // Default constructor
 228     //
 229     fPtpi    =  0;
 230     fPtka    =  0;
 231     fETApic  =  0;
 232     fETAkac  =  0;
 233     fDecayer =  0;
 234     fNt      = -1;
 235     SetCutVertexZ();
 236     SetPtRange();
 237     SetPi0Decays();
 238 }
 239
 240 //_____________________________________________________________________________
 241 AliGenHIJINGpara::AliGenHIJINGpara(Int_t npart)
 242   :AliGenerator(npart)
 243 {
 244   //
 245   // Standard constructor
 246   //
 247     fName="HIGINGpara";
 248     fTitle="HIJING Parametrisation Particle Generator";
 249     fPtpi    =  0;
 250     fPtka    =  0;
 251     fETApic  =  0;
 252     fETAkac  =  0;
 253     fDecayer =  0;
 254     fNt      = -1;
 255     SetCutVertexZ();
 256     SetPtRange();
 257     SetPi0Decays();
 258 }
 259
 260 //_____________________________________________________________________________
 261 AliGenHIJINGpara::~AliGenHIJINGpara()
 262 {
 263   //
 264   // Standard destructor
 265   //
 266     delete fPtpi;
 267     delete fPtka;
 268     delete fETApic;
 269     delete fETAkac;
 270 }
 271
 272 //_____________________________________________________________________________
 273 void AliGenHIJINGpara::Init()
 274 {
 275   //
 276   // Initialise the HIJING parametrisation
 277   //
 278     Float_t etaMin =-TMath::Log(TMath::Tan(
 279         TMath::Min((Double_t)fThetaMax/2,TMath::Pi()/2-1.e-10)));
 280     Float_t etaMax = -TMath::Log(TMath::Tan(
 281         TMath::Max((Double_t)fThetaMin/2,1.e-10)));
 282     fPtpi   = new TF1("ptpi",&ptpi,0,20,0);
 283     fPtka   = new TF1("ptka",&ptka,0,20,0);
 284     fETApic = new TF1("etapic",&etapic,etaMin,etaMax,0);
 285     fETAkac = new TF1("etakac",&etakac,etaMin,etaMax,0);
 286
 287     TF1 etaPic0("etapic",&etapic,-7,7,0);
 288     TF1 etaKac0("etakac",&etakac,-7,7,0);
 289
 290     TF1 ptPic0("ptpi",&ptpi,0.,15.,0);
 291     TF1 ptKac0("ptka",&ptka,0.,15.,0);
 292
 293     Float_t intETApi  = etaPic0.Integral(-0.5, 0.5);
 294     Float_t intETAka  = etaKac0.Integral(-0.5, 0.5);
 295     Float_t scalePi   = 7316/(intETApi/1.5);
 296     Float_t scaleKa   =  684/(intETAka/2.0);
 297
 298 //  Fraction of events corresponding to the selected pt-range
 299     Float_t intPt    = (0.877*ptPic0.Integral(0, 15)+
 300                         0.123*ptKac0.Integral(0, 15));
 301     Float_t intPtSel = (0.877*ptPic0.Integral(fPtMin, fPtMax)+
 302                         0.123*ptKac0.Integral(fPtMin, fPtMax));
 303     Float_t ptFrac   = intPtSel/intPt;
 304
 305 //  Fraction of events corresponding to the selected eta-range
 306     Float_t intETASel  = (scalePi*etaPic0.Integral(etaMin, etaMax)+
 307                           scaleKa*etaKac0.Integral(etaMin, etaMax));
 308 //  Fraction of events corresponding to the selected phi-range
 309     Float_t phiFrac    = (fPhiMax-fPhiMin)/2/TMath::Pi();
 310
 311     fParentWeight = Float_t(fNpart)/(intETASel*ptFrac*phiFrac);
 312
 313     printf("%s: The number of particles in the selected kinematic region corresponds to %f percent of a full event\n ",
 314            ClassName(),100.*fParentWeight);
 315
 316 // Issue warning message if etaMin or etaMax are outside the alowed range
 317 // of the parametrization
 318     if (etaMin < -8.001 || etaMax > 8.001) {
 319         printf("\n \n WARNING FROM AliGenHIJINGPara !");
 320         printf("\n YOU ARE USING THE PARAMETERISATION OUTSIDE ");
 321         printf("\n THE ALLOWED PSEUDORAPIDITY RANGE (-8. - 8.)");
 322         printf("\n YOUR LIMITS: %f %f \n \n ", etaMin, etaMax);
 323     }
 324 //
 325 //
 326     if (fPi0Decays) fDecayer = new AliDecayerPythia();
 327 }
 328
 329
 330 //_____________________________________________________________________________
 331 void AliGenHIJINGpara::Generate()
 332 {
 333   //
 334   // Generate one trigger
 335   //
 336
 337
 338     const Float_t kRaKpic=0.14;
 339     const Float_t kBorne=1/(1+kRaKpic);
 340     Float_t polar[3]= {0,0,0};
 341     //
 342     const Int_t kPions[3] = {kPi0, kPiPlus, kPiMinus};
 343     const Int_t kKaons[4] = {kK0Long, kK0Short, kKPlus, kKMinus};
 344     //
 345     Float_t origin[3];
 346     Float_t pt, pl, ptot;
 347     Float_t phi, theta;
 348     Float_t p[3];
 349     Int_t i, part, j;
 350     //
 351     TF1 *ptf;
 352     TF1 *etaf;
 353     //
 354     Float_t random[6];
 355     //
 356     for (j=0;j<3;j++) origin[j]=fOrigin[j];
 357
 358     if(fVertexSmear == kPerEvent) {
 359         Float_t dv[3];
 360         dv[2] = 1.e10;
 361         while(TMath::Abs(dv[2]) > fCutVertexZ*fOsigma[2]) {
 362             Rndm(random,6);
 363             for (j=0; j < 3; j++) {
 364                 dv[j] = fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 365                     TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 366             }
 367         }
 368         for (j=0; j < 3; j++) origin[j] += dv[j];
 369     } // if kPerEvent
 370     TArrayF eventVertex;
 371     eventVertex.Set(3);
 372     eventVertex[0] = origin[0];
 373     eventVertex[1] = origin[1];
 374     eventVertex[2] = origin[2];
 375
 376     for(i=0;i<fNpart;i++) {
 377         while(1) {
 378             Rndm(random,3);
 379             if(random[0]<kBorne) {
 380                 part=kPions[Int_t (random[1]*3)];
 381                 ptf=fPtpi;
 382                 etaf=fETApic;
 383             } else {
 384                 part=kKaons[Int_t (random[1]*4)];
 385                 ptf=fPtka;
 386                 etaf=fETAkac;
 387             }
 388             phi=fPhiMin+random[2]*(fPhiMax-fPhiMin);
 389             theta=2*TMath::ATan(TMath::Exp(-etaf->GetRandom()));
 390             if(theta<fThetaMin || theta>fThetaMax) continue;
 391             pt=ptf->GetRandom();
 392             pl=pt/TMath::Tan(theta);
 393             ptot=TMath::Sqrt(pt*pt+pl*pl);
 394             if(ptot<fPMin || ptot>fPMax) continue;
 395             p[0]=pt*TMath::Cos(phi);
 396             p[1]=pt*TMath::Sin(phi);
 397             p[2]=pl;
 398             if(fVertexSmear==kPerTrack) {
 399                 Rndm(random,6);
 400                 for (j=0;j<3;j++) {
 401                     origin[j]=fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 402                         TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 403                 }
 404             }
 405             if (part == kPi0 && fPi0Decays){
 406 //
 407 //          Decay pi0 if requested
 408                 SetTrack(0,-1,part,p,origin,polar,0,kPPrimary,fNt,fParentWeight);
 409                 KeepTrack(fNt);
 410                 DecayPi0(origin, p);
 411             } else {
 412                 SetTrack(fTrackIt,-1,part,p,origin,polar,0,kPPrimary,fNt,fParentWeight);
 413                 KeepTrack(fNt);
 414             }
 415
 416             break;
 417         }
 418         SetHighWaterMark(fNt);
 419     }
 420 //
 421
 422 // Header
 423     AliGenEventHeader* header = new AliGenEventHeader("HIJINGparam");
 424 // Event Vertex
 425     header->SetPrimaryVertex(eventVertex);
 426     gAlice->SetGenEventHeader(header);
 427 }
 428
 429 AliGenHIJINGpara& AliGenHIJINGpara::operator=(const  AliGenHIJINGpara& rhs)
 430 {
 431 // Assignment operator
 432     return *this;
 433 }
 434
 435 void AliGenHIJINGpara::SetPtRange(Float_t ptmin, Float_t ptmax) {
 436     AliGenerator::SetPtRange(ptmin, ptmax);
 437 }
 438
 439 void AliGenHIJINGpara::DecayPi0(Float_t* orig, Float_t * p)
 440 {
 441 //
 442 //    Decay the pi0
 443 //    and put decay products on the stack
 444 //
 445     static TClonesArray *particles;
 446     if(!particles) particles = new TClonesArray("TParticle",1000);
 447 //
 448     const Float_t kMass = TDatabasePDG::Instance()->GetParticle(kPi0)->Mass();
 449     Float_t       e     = TMath::Sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]+ kMass * kMass);
 450 //
 451 //  Decay the pi0
 452     TLorentzVector pmom(p[0], p[1], p[2], e);
 453     fDecayer->Decay(kPi0, &pmom);
 454
 455 //
 456 // Put decay particles on the stack
 457 //
 458     Float_t polar[3] = {0., 0., 0.};
 459     Int_t np = fDecayer->ImportParticles(particles);
 460     Int_t nt;
 461     for (Int_t i = 1; i < np; i++)
 462     {
 463         TParticle* iParticle =  (TParticle *) particles->At(i);
 464         p[0] = iParticle->Px();
 465         p[1] = iParticle->Py();
 466         p[2] = iParticle->Pz();
 467         Int_t part = iParticle->GetPdgCode();
 468
 469         SetTrack(fTrackIt, fNt, part, p, orig, polar, 0, kPDecay, nt, fParentWeight);
 470         KeepTrack(nt);
 471     }
 472     fNt = nt;
 473 }