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