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