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