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