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