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