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