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