EVGEN/AliGenEMlib.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 /* $Id: AliGenEMlib.cxx 30052 2008-11-25 14:54:18Z morsch $ */
  17
  18 /////////////////////////////////////////////////////////////////////////////
  19 //                                                                         //
  20 // Implementation of AliGenEMlib for electron, di-electron, and photon     //
  21 // cocktail calculations.                                                  //
  22 // It is based on AliGenGSIlib.                                            //
  23 //                                                                         //
  24 // Responsible: R.Averbeck@gsi.de                                          //
  25 //                                                                         //
  26 /////////////////////////////////////////////////////////////////////////////
  27
  28
  29 #include "TMath.h"
  30 #include "TRandom.h"
  31 #include "TString.h"
  32 #include "AliGenEMlib.h"
  33
  34
  35 ClassImp(AliGenEMlib)
  36
  37 AliGenEMlib::Centbins_t AliGenEMlib::fCentbin=k2030;
  38
  39 Double_t AliGenEMlib::CrossOverLc(const double a, const double b, const double x){
  40 if(x<b-a/2) return 1.0;
  41   else if(x>b+a/2) return 0.0;
  42   else return cos(((x-b)/a+0.5)*TMath::Pi())/2+0.5;
  43 }
  44 Double_t AliGenEMlib::CrossOverRc(const double a, const double b, const double x){
  45   return 1-CrossOverLc(a,b,x);
  46 }
  47
  48 const Double_t AliGenEMlib::fpTparam[8][14] = {
  49   // charged pion
  50   { 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 0.0000000000, 1.0000000000, 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 9999.0000000, 1.0000000000, 0.0000000000, 0.0000000000 },   //
  51   { 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 0.0000000000, 1.0000000000, 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 9999.0000000, 1.0000000000, 0.0000000000, 0.0000000000 },   //
  52   { 3.217746e+03, 1.632570e+00, 9.340162e+00, 1.0000, 2.650000e+00, 4.200635e+00, 9.039589e+08, 1.191548e-01, 6.907293e+00, 1.0000, 6.750000e+00, 4.099970e+00, 6.036772e+01, 5.928279e+00 },   //10-20
  53   { 2.125480e+03, 1.711882e+00, 9.585665e+00, 1.0000, 3.100000e+00, 5.041511e+00, 2.431808e+08, 1.155071e-01, 6.574663e+00, 1.0000, 6.250000e+00, 1.842070e+00, 3.928902e+01, 5.860970e+00 },   //20-30
  54   { 1.577897e+03, 1.411948e+00, 8.638815e+00, 1.0000, 2.550000e+00, 4.066432e+00, 3.774439e+08, 1.000330e-01, 6.535971e+00, 1.0000, 6.750000e+00, 2.482514e+00, 3.495494e+01, 5.954321e+00 },   //30-40
  55   { 7.061859e+02, 1.223810e+00, 8.532113e+00, 1.0000, 1.350000e+00, 1.956213e+00, 1.318169e+04, 5.658401e-01, 7.157575e+00, 1.0000, 5.250000e+00, 3.900000e+00, 1.958841e+01, 6.114787e+00 },   //40-50
  56   { 7.061859e+02, 1.223810e+00, 8.532113e+00, 1.0000, 1.350000e+00, 1.956213e+00, 1.318169e+04, 5.658401e-01, 7.157575e+00, 1.0000, 5.250000e+00, 3.900000e+00, 1.958841e+01, 6.114787e+00 },   //
  57   { 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 0.0000000000, 1.0000000000, 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 9999.0000000, 1.0000000000, 0.0000000000, 0.0000000000 } }; //pp
  58
  59 const Double_t AliGenEMlib::fv2param[2][8][9] = { {
  60     // charged pion
  61     {  3.971545e-02, 2.111261e+00, 6.499885e-02, 3.659836e+00, 7.812234e+00, 1.479471e-02, 1.241708e+00, 2.817950e+00, 1.910663e-01 }   //0-5
  62     ,{ 7.424082e-02, 2.417588e+00, 7.891084e-02, 2.232841e+00, 3.398147e+00, 3.410206e-02, 4.138276e-01,-1.152430e+00, 5.729093e-01 }   //5-10
  63     ,{ 1.094608e-01, 2.357420e+00, 7.614904e-02, 2.294245e+00, 3.538364e+00, 4.932739e-02, 4.557926e-01, 9.218702e-03, 6.428540e-01 }   //10-20
  64     ,{ 1.456377e-01, 2.322408e+00, 7.747166e-02, 2.148424e+00, 3.238113e+00, 5.414722e-02, 4.042938e-01, 1.903040e-01, 6.816790e-01 }   //20-30
  65     ,{ 1.745154e-01, 2.234489e+00, 7.962225e-02, 2.007075e+00, 2.925536e+00, 5.499138e-02, 3.958957e-01, 1.780793e+00, 4.852930e-01 }   //30-40
  66     ,{ 2.384302e-01, 1.578935e+00, 9.234643e-02, 4.328926e+00, 1.020669e+01, 1.001340e-01, 2.433905e+00, 2.966673e+00, 2.646807e-01 }   //40-50
  67     ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }   //
  68     ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }   //pp no v2
  69   },{
  70     // charged kaon
  71     {  0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }  //
  72     ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }  //
  73     ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }  //
  74     ,{ 1.432314e-01, 1.607297e+00, 2.296868e-01, 1.821156e+00, 2.951652e+00,-2.201385e-01, 1.110419e-01, 8.228701e+00, 4.469488e-01 }  //20-30
  75     ,{ 1.691586e-01, 1.617165e+00, 2.159350e-01, 1.649338e+00, 2.607635e+00,-9.536279e+00, 3.860086e-01, 1.802996e+01, 2.509780e-01 }  //30-40
  76     ,{ 1.733831e-01, 1.712705e+00, 1.935862e-01, 1.745523e+00, 2.845436e+00,-5.772356e-01, 6.861616e-02, 5.292878e+00, 9.058815e-01 }  //40-50
  77     ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }  //
  78     ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }  //pp no v2
  79   } };
  80
  81 //==========================================================================
  82 //
  83 //              Definition of Particle Distributions
  84 //
  85 //==========================================================================
  86
  87 //--------------------------------------------------------------------------
  88 //
  89 //                              Pizero
  90 //
  91 //--------------------------------------------------------------------------
  92 Int_t AliGenEMlib::IpPizero(TRandom *)
  93 {
  94 // Return pizero pdg code
  95   return 111;
  96 }
  97
  98 Double_t AliGenEMlib::PtPizero( const Double_t *px, const Double_t */*dummy*/ )
  99 {
 100 // Generate pizero pT distribution from modified Hagedorn parameterization
 101 // taken from fit to unidentified hadrons in pp at 7 TeV
 102   // const Double_t kc=0.000565;
 103   // const Double_t kp0=0.2472;
 104   // const Double_t kp1=4.354;
 105   // const Double_t kn=7.007;
 106   Double_t invYield;
 107   Double_t x=*px;
 108
 109   // invYield = kc/TMath::Power(kp0+x/kp1,kn);
 110
 111   if(!x)return 0;
 112   const Double_t *c=fpTparam[fCentbin];
 113
 114   invYield = CrossOverLc(c[5],c[4],x)*c[0]/TMath::Power(c[3]+x/c[1],c[2]);
 115   invYield +=CrossOverRc(c[5],c[4],x)*c[6]/TMath::Power(c[9]+x/c[7],c[8])*CrossOverLc(c[11],c[10],x);
 116   invYield +=CrossOverRc(c[11],c[10],x)*c[12]/TMath::Power(x,c[13]);
 117
 118   return invYield*(2*TMath::Pi()*x);
 119 }
 120
 121 Double_t AliGenEMlib::YPizero( const Double_t *py, const Double_t */*dummy*/ )
 122 {
 123   return YFlat(*py);
 124
 125 }
 126
 127 Double_t AliGenEMlib::V2Pizero( const Double_t *px, const Double_t */*dummy*/ )
 128 {
 129   return V2Param(px,fv2param[0][fCentbin]);
 130 }
 131
 132 //--------------------------------------------------------------------------
 133 //
 134 //                              Eta
 135 //
 136 //--------------------------------------------------------------------------
 137 Int_t AliGenEMlib::IpEta(TRandom *)
 138 {
 139 // Return eta pdg code
 140   return 221;
 141 }
 142
 143 Double_t AliGenEMlib::PtEta( const Double_t *px, const Double_t */*dummy*/ )
 144 {
 145 // Eta pT
 146   return MtScal(*px,1);
 147 }
 148
 149 Double_t AliGenEMlib::YEta( const Double_t *py, const Double_t */*dummy*/ )
 150 {
 151   return YFlat(*py);
 152 }
 153
 154 Double_t AliGenEMlib::V2Eta( const Double_t *px, const Double_t */*dummy*/ )
 155 {
 156   return V2Param(px,fv2param[1][fCentbin]);
 157 }
 158
 159 //--------------------------------------------------------------------------
 160 //
 161 //                              Rho
 162 //
 163 //--------------------------------------------------------------------------
 164 Int_t AliGenEMlib::IpRho(TRandom *)
 165 {
 166 // Return rho pdg code
 167   return 113;
 168 }
 169
 170 Double_t AliGenEMlib::PtRho( const Double_t *px, const Double_t */*dummy*/ )
 171 {
 172 // Rho pT
 173   return MtScal(*px,2);
 174 }
 175
 176 Double_t AliGenEMlib::YRho( const Double_t *py, const Double_t */*dummy*/ )
 177 {
 178   return YFlat(*py);
 179
 180 }
 181
 182 //--------------------------------------------------------------------------
 183 //
 184 //                              Omega
 185 //
 186 //--------------------------------------------------------------------------
 187 Int_t AliGenEMlib::IpOmega(TRandom *)
 188 {
 189 // Return omega pdg code
 190   return 223;
 191 }
 192
 193 Double_t AliGenEMlib::PtOmega( const Double_t *px, const Double_t */*dummy*/ )
 194 {
 195 // Omega pT
 196   return MtScal(*px,3);
 197 }
 198
 199 Double_t AliGenEMlib::YOmega( const Double_t *py, const Double_t */*dummy*/ )
 200 {
 201   return YFlat(*py);
 202
 203 }
 204
 205 //--------------------------------------------------------------------------
 206 //
 207 //                              Etaprime
 208 //
 209 //--------------------------------------------------------------------------
 210 Int_t AliGenEMlib::IpEtaprime(TRandom *)
 211 {
 212 // Return etaprime pdg code
 213   return 331;
 214 }
 215
 216 Double_t AliGenEMlib::PtEtaprime( const Double_t *px, const Double_t */*dummy*/ )
 217 {
 218 // Eta pT
 219   return MtScal(*px,4);
 220 }
 221
 222 Double_t AliGenEMlib::YEtaprime( const Double_t *py, const Double_t */*dummy*/ )
 223 {
 224   return YFlat(*py);
 225
 226 }
 227
 228 //--------------------------------------------------------------------------
 229 //
 230 //                              Phi
 231 //
 232 //--------------------------------------------------------------------------
 233 Int_t AliGenEMlib::IpPhi(TRandom *)
 234 {
 235 // Return phi pdg code
 236   return 333;
 237 }
 238
 239 Double_t AliGenEMlib::PtPhi( const Double_t *px, const Double_t */*dummy*/ )
 240 {
 241 // Phi pT
 242   return MtScal(*px,5);
 243 }
 244
 245 Double_t AliGenEMlib::YPhi( const Double_t *py, const Double_t */*dummy*/ )
 246 {
 247   return YFlat(*py);
 248
 249 }
 250
 251 Double_t AliGenEMlib::YFlat(Double_t /*y*/)
 252 {
 253 //--------------------------------------------------------------------------
 254 //
 255 //                    flat rapidity distribution
 256 //
 257 //--------------------------------------------------------------------------
 258
 259   Double_t dNdy = 1.;
 260
 261   return dNdy;
 262
 263 }
 264
 265 //=============================================================
 266 //
 267 //                    Mt-scaling
 268 //
 269 //=============================================================
 270 //
 271  Double_t AliGenEMlib::MtScal(Double_t pt, Int_t np)
 272 {
 273 // Function for the calculation of the Pt distribution for a
 274 // given particle np, from the pizero Pt distribution using
 275 // mt scaling.
 276
 277 // MASS   0=>PIZERO, 1=>ETA, 2=>RHO, 3=>OMEGA, 4=>ETAPRIME, 5=>PHI
 278
 279   const Double_t khm[6] = {0.13498, 0.54751, 0.7755, 0.78265, 0.95778, 1.01946};
 280
 281   Double_t scaledPt = sqrt(pt*pt + khm[np]*khm[np] - khm[0]*khm[0]);
 282   Double_t scaledYield = PtPizero(&scaledPt, (Double_t*) 0);
 283
 284   //     VALUE MESON/PI AT 5 GEV
 285
 286   Double_t normPt = 5.;
 287   Double_t scaledNormPt = sqrt(normPt*normPt + khm[np]*khm[np] - khm[0]*khm[0]);
 288   const Double_t kfmax[6]={1., 0.48, 1.0, 0.9, 0.25, 0.4};
 289
 290   Double_t norm = kfmax[np] * (PtPizero(&normPt, (Double_t*) 0) / PtPizero(&scaledNormPt, (Double_t*) 0));
 291
 292   return norm*(pt/scaledPt)*scaledYield;
 293 }
 294
 295 Double_t AliGenEMlib::V2Param(const Double_t *px, const Double_t *par)
 296 {
 297   // Very general parametrization of the v2
 298
 299   return TMath::Max(CrossOverLc(par[4],par[3],px[0])*(2*par[0]/(1+TMath::Exp(par[1]*(par[2]-px[0])))-par[0])+CrossOverRc(par[4],par[3],px[0])*((par[8]-par[5])/(1+TMath::Exp(par[6]*(px[0]-par[7])))+par[5]),0.0);
 300 }
 301
 302 Double_t AliGenEMlib::V2Flat(const Double_t */*px*/, const Double_t */*param*/)
 303 {
 304   // Flat v2
 305
 306   return 1.0;
 307 }
 308
 309 //==========================================================================
 310 //
 311 //                     Set Getters
 312 //
 313 //==========================================================================
 314
 315 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
 316
 317 typedef Int_t (*GenFuncIp) (TRandom *);
 318
 319 GenFunc AliGenEMlib::GetPt(Int_t param, const char * tname) const
 320 {
 321 // Return pointer to pT parameterisation
 322    GenFunc func=0;
 323    TString sname(tname);
 324
 325    switch (param)
 326     {
 327     case kPizero:
 328       func=PtPizero;
 329       break;
 330     case kEta:
 331       func=PtEta;
 332       break;
 333     case kRho:
 334       func=PtRho;
 335       break;
 336     case kOmega:
 337       func=PtOmega;
 338       break;
 339     case kEtaprime:
 340       func=PtEtaprime;
 341       break;
 342     case kPhi:
 343       func=PtPhi;
 344       break;
 345
 346     default:
 347       func=0;
 348       printf("<AliGenEMlib::GetPt> unknown parametrisation\n");
 349     }
 350    return func;
 351 }
 352
 353 GenFunc AliGenEMlib::GetY(Int_t param, const char * tname) const
 354 {
 355 // Return pointer to y- parameterisation
 356    GenFunc func=0;
 357    TString sname(tname);
 358
 359    switch (param)
 360     {
 361     case kPizero:
 362          func=YPizero;
 363          break;
 364     case kEta:
 365          func=YEta;
 366          break;
 367     case kRho:
 368          func=YRho;
 369          break;
 370     case kOmega:
 371          func=YOmega;
 372          break;
 373     case kEtaprime:
 374          func=YEtaprime;
 375          break;
 376     case kPhi:
 377          func=YPhi;
 378          break;
 379
 380     default:
 381         func=0;
 382         printf("<AliGenEMlib::GetY> unknown parametrisation\n");
 383     }
 384     return func;
 385 }
 386
 387 GenFuncIp AliGenEMlib::GetIp(Int_t param, const char * tname) const
 388 {
 389 // Return pointer to particle type parameterisation
 390    GenFuncIp func=0;
 391    TString sname(tname);
 392
 393    switch (param)
 394     {
 395     case kPizero:
 396          func=IpPizero;
 397          break;
 398     case kEta:
 399          func=IpEta;
 400          break;
 401     case kRho:
 402          func=IpRho;
 403          break;
 404     case kOmega:
 405          func=IpOmega;
 406          break;
 407     case kEtaprime:
 408          func=IpEtaprime;
 409          break;
 410     case kPhi:
 411          func=IpPhi;
 412          break;
 413
 414     default:
 415         func=0;
 416         printf("<AliGenEMlib::GetIp> unknown parametrisation\n");
 417     }
 418     return func;
 419 }
 420
 421 GenFunc AliGenEMlib::GetV2(Int_t param, const char * tname) const
 422 {
 423   // Return pointer to v2-parameterisation
 424   GenFunc func=0;
 425   TString sname(tname);
 426
 427   switch (param)
 428     {
 429     case kPizero:
 430       func=V2Pizero;
 431       break;
 432     case kEta:
 433       func=V2Eta;
 434       break;
 435     case kRho:
 436       func=V2Pizero;
 437       break;
 438     case kOmega:
 439       func=V2Pizero;
 440       break;
 441     case kEtaprime:
 442       func=V2Pizero;
 443       break;
 444     case kPhi:
 445       func=V2Pizero;
 446       break;
 447
 448     default:
 449       func=0;
 450       printf("<AliGenEMlib::GetV2> unknown parametrisation\n");
 451     }
 452   return func;
 453 }