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