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