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