EVGEN/AliGenMUONlib.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$ */
  17
  18 // Library class for particle pt and y distributions used for
  19 // muon spectrometer simulations.
  20 // To be used with AliGenParam.
  21 // The following particle typed can be simulated:
  22 // pi, K, phi, omega, eta, J/Psi, Upsilon, charm and beauty mesons.
  23 //
  24 // andreas.morsch@cern.ch
  25 //
  26
  27 #include "TMath.h"
  28 #include "TRandom.h"
  29
  30 #include "AliGenMUONlib.h"
  31
  32 ClassImp(AliGenMUONlib)
  33 //
  34 //  Pions
  35 Double_t AliGenMUONlib::PtPion(const Double_t *px, const Double_t* /*dummy*/)
  36 {
  37 //
  38 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
  39 //     POWER LAW FOR PT > 500 MEV
  40 //     MT SCALING BELOW (T=160 MEV)
  41 //
  42   const Double_t kp0 = 1.3;
  43   const Double_t kxn = 8.28;
  44   const Double_t kxlim=0.5;
  45   const Double_t kt=0.160;
  46   const Double_t kxmpi=0.139;
  47   const Double_t kb=1.;
  48   Double_t y, y1, xmpi2, ynorm, a;
  49   Double_t x=*px;
  50   //
  51   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
  52   xmpi2=kxmpi*kxmpi;
  53   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
  54   a=ynorm/y1;
  55   if (x > kxlim)
  56     y=a*TMath::Power(kp0/(kp0+x),kxn);
  57   else
  58     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
  59   return y*x;
  60 }
  61 //
  62 // y-distribution
  63 //
  64 Double_t AliGenMUONlib::YPion( const Double_t *py, const Double_t */*dummy*/)
  65 {
  66 // Pion y
  67   Double_t y=TMath::Abs(*py);
  68 /*
  69   const Double_t ka    = 7000.;
  70   const Double_t kdy   = 4.;
  71   Double_t ex = y*y/(2*kdy*kdy);
  72   return ka*TMath::Exp(-ex);
  73 */
  74   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
  75
  76 }
  77 //                 particle composition
  78 //
  79 Int_t AliGenMUONlib::IpPion(TRandom *ran)
  80 {
  81 // Pion composition
  82     if (ran->Rndm() < 0.5) {
  83         return  211;
  84     } else {
  85         return -211;
  86     }
  87 }
  88
  89 //____________________________________________________________
  90 //
  91 // Mt-scaling
  92
  93 Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
  94 {
  95   //    SCALING EN MASSE PAR RAPPORT A PTPI
  96   //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
  97   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
  98   //     VALUE MESON/PI AT 5 GEV
  99   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
 100   np--;
 101   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
 102   Double_t fmax2=f5/kfmax[np];
 103   // PIONS
 104   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
 105   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
 106                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
 107   return fmtscal*ptpion;
 108 }
 109 //
 110 // kaon
 111 //
 112 //                pt-distribution
 113 //____________________________________________________________
 114 Double_t AliGenMUONlib::PtKaon( const Double_t *px, const Double_t */*dummy*/)
 115 {
 116 // Kaon pT
 117   return PtScal(*px,2);
 118 }
 119
 120 // y-distribution
 121 //____________________________________________________________
 122 Double_t AliGenMUONlib::YKaon( const Double_t *py, const Double_t */*dummy*/)
 123 {
 124 // Kaon y
 125   Double_t y=TMath::Abs(*py);
 126 /*
 127   const Double_t ka    = 1000.;
 128   const Double_t kdy   = 4.;
 129   //
 130   Double_t ex = y*y/(2*kdy*kdy);
 131   return ka*TMath::Exp(-ex);
 132 */
 133
 134   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
 135 }
 136
 137 //                 particle composition
 138 //
 139 Int_t AliGenMUONlib::IpKaon(TRandom *ran)
 140 {
 141 // Kaon composition
 142     if (ran->Rndm() < 0.5) {
 143         return  321;
 144     } else {
 145         return -321;
 146     }
 147 }
 148
 149 //                    J/Psi
 150 //
 151 //
 152 //                pt-distribution
 153 //____________________________________________________________
 154 Double_t AliGenMUONlib::PtJpsiPP7000( const Double_t *px, const Double_t */*dummy*/)
 155 {
 156 // J/Psi pT
 157 //
 158 // pp 7 TeV
 159 // using ALICE data at 2.5<y<4, see arXiv:1103.2394
 160
 161   const Double_t kpt0 = 2.44;
 162   const Double_t kxn  = 3.9;
 163   Double_t x=*px;
 164   //
 165   Double_t pass1 = 1.+0.36*(x/kpt0)*(x/kpt0);
 166   return x/TMath::Power(pass1,kxn);
 167 }
 168
 169 Double_t AliGenMUONlib::PtJpsiPP2760( const Double_t *px, const Double_t */*dummy*/)
 170 {
 171 // J/Psi pT
 172 //
 173 // pp 2.76 TeV
 174 // from the fit of RHIC + LHC data, see arXiv:1103.2394
 175
 176   const Double_t kpt0 = 2.31;
 177   const Double_t kxn  = 3.9;
 178   Double_t x=*px;
 179   //
 180   Double_t pass1 = 1.+0.36*(x/kpt0)*(x/kpt0);
 181   return x/TMath::Power(pass1,kxn);
 182 }
 183
 184 Double_t AliGenMUONlib::PtJpsiPbPb2760( const Double_t *px, const Double_t *dummy)
 185 {
 186 // J/Psi pT
 187 //
 188 // PbPb 2.76 TeV, for EKS98 with minimum bias shadowing factor 0.66
 189 //
 190   Double_t c[5] = {6.01022e-01, 4.70988e-02, -2.27917e-03, 3.09885e-05, 1.31955e-06};
 191   Double_t x=*px;
 192   Double_t y;
 193   Int_t j;
 194   y = c[j = 4];
 195   while (j > 0) y  = y * x + c[--j];
 196   //
 197   Double_t d = 1.+c[4]*TMath::Power(x,4);
 198   return y/d * AliGenMUONlib::PtJpsiPP2760(px,dummy);
 199 }
 200
 201 Double_t AliGenMUONlib::PtJpsi( const Double_t *px, const Double_t */*dummy*/)
 202 {
 203 // J/Psi pT
 204   const Double_t kpt0 = 4.;
 205   const Double_t kxn  = 3.6;
 206   Double_t x=*px;
 207   //
 208   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 209   return x/TMath::Power(pass1,kxn);
 210 }
 211
 212 Double_t AliGenMUONlib::PtJpsiCDFscaled( const Double_t *px, const Double_t */*dummy*/)
 213 {
 214 // J/Psi pT
 215 //
 216 // PbPb 5.5 TeV
 217 // scaled from CDF data at 2 TeV
 218 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
 219
 220   const Double_t kpt0 = 5.100;
 221   const Double_t kxn  = 4.102;
 222   Double_t x=*px;
 223   //
 224   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 225   return x/TMath::Power(pass1,kxn);
 226 }
 227
 228 Double_t AliGenMUONlib::PtJpsiCDFscaledPP( const Double_t *px, const Double_t */*dummy*/)
 229 {
 230 // J/Psi pT
 231 //
 232 // pp 14 TeV
 233 // scaled from CDF data at 2 TeV
 234
 235   const Double_t kpt0 = 5.630;
 236   const Double_t kxn  = 4.071;
 237   Double_t x=*px;
 238   //
 239   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 240   return x/TMath::Power(pass1,kxn);
 241 }
 242
 243 Double_t AliGenMUONlib::PtJpsiCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
 244 {
 245 // J/Psi pT
 246 //
 247 // pp 10 TeV
 248 // scaled from CDF data at 2 TeV
 249
 250   const Double_t kpt0 = 5.334;
 251   const Double_t kxn  = 4.071;
 252   Double_t x=*px;
 253   //
 254   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 255   return x/TMath::Power(pass1,kxn);
 256 }
 257
 258 Double_t AliGenMUONlib::PtJpsiCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
 259 {
 260 // J/Psi pT
 261 //
 262 // pp 8.8 TeV
 263 // scaled from CDF data at 2 TeV
 264 //
 265   const Double_t kpt0 = 5.245;
 266   const Double_t kxn  = 4.071;
 267   Double_t x=*px;
 268   //
 269   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 270   return x/TMath::Power(pass1,kxn);
 271 }
 272
 273 Double_t AliGenMUONlib::PtJpsiCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
 274 {
 275 // J/Psi pT
 276 //
 277 // pp 7 TeV
 278 // scaled from CDF data at 2 TeV
 279
 280   const Double_t kpt0 = 5.072;
 281   const Double_t kxn  = 4.071;
 282   Double_t x=*px;
 283   //
 284   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 285   return x/TMath::Power(pass1,kxn);
 286 }
 287
 288 Double_t AliGenMUONlib::PtJpsiCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
 289 {
 290 // J/Psi pT
 291 //
 292 // pp 3.94 TeV
 293 // scaled from CDF data at 2 TeV
 294 //
 295   const Double_t kpt0 = 4.647;
 296   const Double_t kxn  = 4.071;
 297   Double_t x=*px;
 298   //
 299   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 300   return x/TMath::Power(pass1,kxn);
 301 }
 302
 303 Double_t AliGenMUONlib::PtJpsiCDFscaledPP3( const Double_t *px, const Double_t */*dummy*/)
 304 {
 305 // J/Psi pT
 306 //
 307 // pp 2.76 TeV
 308 // scaled from CDF data at 1.9 TeV
 309 //
 310   const Double_t kpt0 = 4.435;
 311   const Double_t kxn  = 4.071;
 312   Double_t x=*px;
 313   //
 314   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 315   return x/TMath::Power(pass1,kxn);
 316 }
 317
 318 Double_t AliGenMUONlib::PtJpsiCDFscaledPP2( const Double_t *px, const Double_t */*dummy*/)
 319 {
 320 // J/Psi pT
 321 //
 322 // pp 1.9 TeV
 323 // fit of the CDF data at 1.9 TeV
 324 //
 325   const Double_t kpt0 = 4.233;
 326   const Double_t kxn  = 4.071;
 327   Double_t x=*px;
 328   //
 329   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 330   return x/TMath::Power(pass1,kxn);
 331 }
 332
 333 Double_t AliGenMUONlib::PtJpsiCDFscaledPPb9( const Double_t *px, const Double_t *dummy)
 334 {
 335 // J/Psi pT
 336 //
 337 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.80
 338 //
 339   Double_t c[5] = {6.42774e-01, 1.86168e-02, -6.77296e-04, 8.93512e-06, 1.31586e-07};
 340   Double_t x=*px;
 341   Double_t y;
 342   Int_t j;
 343   y = c[j = 4];
 344   while (j > 0) y  = y * x + c[--j];
 345   //
 346   Double_t d = 1.+c[4]*TMath::Power(x,4);
 347   return y/d * AliGenMUONlib::PtJpsiCDFscaledPP9(px,dummy);
 348 }
 349
 350 Double_t AliGenMUONlib::PtJpsiCDFscaledPbP9( const Double_t *px, const Double_t *dummy)
 351 {
 352 // J/Psi pT
 353 //
 354 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.80
 355 //
 356   Double_t c[5] = {8.58557e-01, 5.39791e-02, -4.75180e-03, 2.49463e-04, 5.52396e-05};
 357   Double_t x=*px;
 358   Double_t y;
 359   Int_t j;
 360   y = c[j = 4];
 361   while (j > 0) y  = y * x + c[--j];
 362   //
 363   Double_t d = 1.+c[4]*TMath::Power(x,4);
 364   return y/d * AliGenMUONlib::PtJpsiCDFscaledPP9(px,dummy);
 365 }
 366
 367 Double_t AliGenMUONlib::PtJpsiCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
 368 {
 369 // J/Psi pT
 370 //
 371 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.66
 372 //
 373   Double_t c[5] = {6.01022e-01, 4.70988e-02, -2.27917e-03, 3.09885e-05, 1.31955e-06};
 374   Double_t x=*px;
 375   Double_t y;
 376   Int_t j;
 377   y = c[j = 4];
 378   while (j > 0) y  = y * x + c[--j];
 379   //
 380   Double_t d = 1.+c[4]*TMath::Power(x,4);
 381   return y/d * AliGenMUONlib::PtJpsiCDFscaledPP4(px,dummy);
 382 }
 383
 384 Double_t AliGenMUONlib::PtJpsiFlat( const Double_t */*px*/, const Double_t */*dummy*/ )
 385 {
 386   return 1.;
 387 }
 388
 389 Double_t AliGenMUONlib::PtJpsiPbPb( const Double_t *px, const Double_t */*dummy*/)
 390 {
 391 // J/Psi pT spectrum
 392 //
 393 // R. Vogt 2002
 394 // PbPb 5.5 TeV
 395 // MRST HO
 396 // mc = 1.4 GeV, pt-kick 1 GeV
 397 //
 398     Float_t x = px[0];
 399     Float_t c[8] = {
 400         -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00,
 401         -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
 402     };
 403
 404     Double_t y;
 405     if (x < 10.) {
 406         Int_t j;
 407         y = c[j = 7];
 408         while (j > 0) y  = y * x +c[--j];
 409         y = x * TMath::Exp(y);
 410     } else {
 411         y = 0.;
 412     }
 413     return y;
 414 }
 415
 416 Double_t AliGenMUONlib::PtJpsiBPbPb( const Double_t *px, const Double_t */*dummy*/)
 417 {
 418 // J/Psi pT spectrum
 419 // B -> J/Psi X
 420     Double_t x0 =   4.0384;
 421     Double_t  n =   3.0288;
 422
 423     Double_t x = px[0];
 424     Double_t y = x / TMath::Power((1. + (x/x0)*(x/x0)), n);
 425
 426     return y;
 427 }
 428
 429
 430 Double_t AliGenMUONlib::PtJpsiPP( const Double_t *px, const Double_t */*dummy*/)
 431 {
 432 // J/Psi pT spectrum
 433 //
 434 // R. Vogt 2002
 435 // pp 14 TeV
 436 // MRST HO
 437 // mc = 1.4 GeV, pt-kick 1 GeV
 438 //
 439     Float_t x = px[0];
 440     Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
 441
 442     Double_t y;
 443     if (x < 10.) {
 444         Int_t j;
 445         y = c[j = 3];
 446         while (j > 0) y  = y * x +c[--j];
 447         y = x * TMath::Exp(y);
 448     } else {
 449         y = 0.;
 450     }
 451     return y;
 452 }
 453
 454 //
 455 //               y-distribution
 456 //____________________________________________________________
 457 Double_t AliGenMUONlib::YJpsiPP7000( const Double_t *px, const Double_t */*dummy*/)
 458 {
 459 // J/Psi y
 460 //
 461 // pp 7 TeV
 462 // from the fit of RHIC + LHC data, see arXiv:1103.2394
 463 //
 464     Double_t x = px[0]/7.72;
 465     x = x*x;
 466     Double_t y = TMath::Exp(-x/0.383/0.383/2);
 467     if(x > 1) y=0;
 468     return y;
 469 }
 470
 471 Double_t AliGenMUONlib::YJpsiPP2760( const Double_t *px, const Double_t */*dummy*/)
 472 {
 473 // J/Psi y
 474 //
 475 // pp 2.76 TeV
 476 // from the fit of RHIC + LHC data, see arXiv:1103.2394
 477 //
 478     Double_t x = px[0]/6.79;
 479     x = x*x;
 480     Double_t y = TMath::Exp(-x/0.383/0.383/2);
 481     if(x > 1) y=0;
 482     return y;
 483 }
 484
 485 Double_t AliGenMUONlib::YJpsiPbPb2760( const Double_t *px, const Double_t *dummy)
 486 {
 487 // J/Psi y
 488 //
 489 // PbPb 2.76 TeV, for EKS98 with minimum bias shadowing factor 0.66
 490 //
 491     Double_t c[4] = {5.95228e-01, 9.45069e-03, 2.44710e-04, -1.32894e-05};
 492     Double_t x = px[0]*px[0];
 493     Double_t y;
 494     Int_t j;
 495     y = c[j = 3];
 496     while (j > 0) y  = y * x + c[--j];
 497     if(y<0) y=0;
 498
 499     return y * AliGenMUONlib::YJpsiPP2760(px,dummy);
 500 }
 501
 502 Double_t AliGenMUONlib::YJpsi(const Double_t *py, const Double_t */*dummy*/)
 503 {
 504 // J/psi y
 505   const Double_t ky0 = 4.;
 506   const Double_t kb=1.;
 507   Double_t yj;
 508   Double_t y=TMath::Abs(*py);
 509   //
 510   if (y < ky0)
 511     yj=kb;
 512   else
 513     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 514   return yj;
 515 }
 516
 517 Double_t AliGenMUONlib::YJpsiFlat( const Double_t */*py*/, const Double_t */*dummy*/ )
 518 {
 519   return 1.;
 520 }
 521
 522
 523 Double_t AliGenMUONlib::YJpsiPbPb( const Double_t *px, const Double_t */*dummy*/)
 524 {
 525
 526 //
 527 // J/Psi y
 528 //
 529 //
 530 // R. Vogt 2002
 531 // PbPb 5.5 TeV
 532 // MRST HO
 533 // mc = 1.4 GeV, pt-kick 1 GeV
 534 //
 535     Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
 536     Double_t x = TMath::Abs(px[0]);
 537     Double_t y;
 538
 539     if (x < 4.) {
 540         y = 31.754;
 541     } else if (x < 6) {
 542         Int_t j;
 543         y = c[j = 4];
 544         while (j > 0) y  = y * x + c[--j];
 545     } else {
 546         y =0.;
 547     }
 548
 549     return y;
 550 }
 551
 552 Double_t AliGenMUONlib::YJpsiCDFscaled( const Double_t *px, const Double_t* dummy)
 553 {
 554     // J/Psi y
 555     return AliGenMUONlib::YJpsiPbPb(px, dummy);
 556 }
 557
 558 Double_t AliGenMUONlib::YJpsiCDFscaledPP( const Double_t *px, const Double_t* dummy)
 559 {
 560     // J/Psi y
 561     return AliGenMUONlib::YJpsiPP(px, dummy);
 562 }
 563
 564 Double_t AliGenMUONlib::YJpsiCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
 565 {
 566 // J/Psi y
 567 //
 568 // pp 10 TeV
 569 // scaled from YJpsiPP(14 TeV) using 10 TeV / 14 TeV ratio of y-spectra in LO pQCD.
 570 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
 571 //
 572
 573     Double_t c[5] = {2.46681e+01, 8.91486e+01, -3.21227e+01, 3.63075e+00, -1.32047e-01};
 574
 575     Double_t x = TMath::Abs(px[0]);
 576     Double_t y;
 577
 578     if (x < 3.2) {
 579         y = 98.523 - 1.3664 * x * x;
 580     } else if (x < 7.5) {
 581         Int_t j;
 582         y = c[j = 4];
 583         while (j > 0) y  = y * x + c[--j];
 584     } else {
 585         y =0.;
 586     }
 587
 588     if(y<0) y=0;
 589
 590     return y;
 591 }
 592
 593 Double_t AliGenMUONlib::YJpsiCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
 594 {
 595 // J/Psi y
 596 //
 597 // pp 8.8 TeV
 598 // rescaling of YJpsiPP(14 TeV) using 8.8 TeV / 14 TeV ratio of y-spectra in LO QCD
 599 //
 600     Double_t c[5] = {3.33882e+02, -1.30980e+02, 2.59082e+01, -3.08935e+00, 1.56375e-01};
 601     Double_t x = TMath::Abs(px[0]);
 602     Double_t y;
 603
 604     if (x < 3.7) {
 605         y = 99.236 - 1.5498 * x * x;
 606     } else if (x < 7.4) {
 607         Int_t j;
 608         y = c[j = 4];
 609         while (j > 0) y  = y * x + c[--j];
 610     } else {
 611         y =0.;
 612     }
 613
 614     if(y<0) y=0;
 615
 616     return y;
 617 }
 618
 619 Double_t AliGenMUONlib::YJpsiCDFscaledPP9dummy(Double_t px)
 620 {
 621     return AliGenMUONlib::YJpsiCDFscaledPP9(&px, (Double_t*) 0);
 622 }
 623
 624 Double_t AliGenMUONlib::YJpsiCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
 625 {
 626 // J/Psi y
 627 //
 628 // pp 7 TeV
 629 // scaled from YJpsiPP(14 TeV) using 7 TeV / 14 TeV ratio of y-spectra in LO pQCD.
 630 //
 631
 632     Double_t c[5] = {6.71181e+02, -3.69240e+02, 8.89644e+01, -1.04937e+01, 4.80959e-01};
 633
 634     Double_t x = TMath::Abs(px[0]);
 635     Double_t y;
 636
 637     if (x < 4.0) {
 638         y = 100.78 - 1.8353 * x * x;
 639     } else if (x < 7.3) {
 640         Int_t j;
 641         y = c[j = 4];
 642         while (j > 0) y  = y * x + c[--j];
 643     } else {
 644         y =0.;
 645     }
 646
 647     if(y<0) y=0;
 648
 649     return y;
 650 }
 651
 652 Double_t AliGenMUONlib::YJpsiCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
 653 {
 654 // J/Psi y
 655 //
 656 // pp 3.94 TeV
 657 // rescaling of YJpsiPP(14 TeV) using 3.94 TeV / 14 TeV ratio of y-spectra in LO QCD
 658 //
 659     Double_t c[5] = {4.00785e+02, -1.41159e+01, -3.28599e+01, 5.53048e+00, -2.45151e-01};
 660     Double_t x = TMath::Abs(px[0]);
 661     Double_t y;
 662
 663     if (x < 5.5) {
 664         y = 107.389 - 2.7454 * x * x;
 665     } else if (x < 7.0) {
 666         Int_t j;
 667         y = c[j = 4];
 668         while (j > 0) y  = y * x + c[--j];
 669     } else {
 670         y =0.;
 671     }
 672
 673     if(y<0) y=0;
 674
 675     return y;
 676 }
 677
 678 Double_t AliGenMUONlib::YJpsiCDFscaledPP3( const Double_t *px, const Double_t *dummy)
 679 {
 680 // J/Psi y
 681     return AliGenMUONlib::YJpsiPP2760(px, dummy);
 682 }
 683
 684 Double_t AliGenMUONlib::YJpsiCDFscaledPP2( const Double_t *px, const Double_t */*dummy*/)
 685 {
 686 // J/Psi y
 687 //
 688 // pp 1.9 TeV
 689 // from the fit of RHIC + LHC data, see arXiv:1103.2394
 690 //
 691     Double_t x = px[0]/6.42;
 692     x = x*x;
 693     Double_t y = TMath::Exp(-x/0.383/0.383/2);
 694     if(x > 1) y=0;
 695     return y;
 696 }
 697
 698 Double_t AliGenMUONlib::YJpsiPP( const Double_t *px, const Double_t */*dummy*/)
 699 {
 700
 701 //
 702 // J/Psi y
 703 //
 704 //
 705 // R. Vogt 2002
 706 // pp 14  TeV
 707 // MRST HO
 708 // mc = 1.4 GeV, pt-kick 1 GeV
 709 //
 710
 711     Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
 712     Double_t x = TMath::Abs(px[0]);
 713     Double_t y;
 714
 715     if (x < 2.5) {
 716         y = 96.455 - 0.8483 * x * x;
 717     } else if (x < 7.9) {
 718         Int_t j;
 719         y = c[j = 4];
 720         while (j > 0) y  = y * x + c[--j];
 721     } else {
 722         y =0.;
 723     }
 724
 725     return y;
 726 }
 727
 728 Double_t AliGenMUONlib::YJpsiCDFscaledPPb9( const Double_t *px, const Double_t */*dummy*/)
 729 {
 730 // J/Psi y
 731 //
 732 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.80
 733 //
 734     Double_t c[7] = {7.52296e-01, 2.49917e-02, 3.36500e-03, 1.91187e-03, 2.92154e-04,
 735                      -4.16509e-05,-7.62709e-06};
 736     Double_t y;
 737     Double_t x = px[0] + 0.47;              // rapidity shift
 738     Int_t j;
 739     y = c[j = 6];
 740     while (j > 0) y = y * x + c[--j];
 741     if(y<0) y=0;
 742
 743     return y * AliGenMUONlib::YJpsiCDFscaledPP9dummy(x);
 744 }
 745
 746 Double_t AliGenMUONlib::YJpsiCDFscaledPbP9( const Double_t *px, const Double_t */*dummy*/)
 747 {
 748 // J/Psi y
 749 //
 750 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.80
 751 //
 752     Double_t c[7] = {7.52296e-01, 2.49917e-02, 3.36500e-03, 1.91187e-03, 2.92154e-04,
 753                      -4.16509e-05,-7.62709e-06};
 754     Double_t y;
 755     Double_t x = -px[0] + 0.47;              // rapidity shift
 756     Int_t j;
 757     y = c[j = 6];
 758     while (j > 0) y = y * x + c[--j];
 759     if(y<0) y=0;
 760
 761     return y * AliGenMUONlib::YJpsiCDFscaledPP9dummy(x);
 762 }
 763
 764 Double_t AliGenMUONlib::YJpsiCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
 765 {
 766 // J/Psi y
 767 //
 768 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.66
 769 //
 770     Double_t c[4] = {5.95228e-01, 9.45069e-03, 2.44710e-04, -1.32894e-05};
 771     Double_t x = px[0]*px[0];
 772     Double_t y;
 773     Int_t j;
 774     y = c[j = 3];
 775     while (j > 0) y  = y * x + c[--j];
 776     if(y<0) y=0;
 777
 778     return y * AliGenMUONlib::YJpsiCDFscaledPP4(px,dummy);
 779 }
 780
 781 Double_t AliGenMUONlib::YJpsiBPbPb( const Double_t *px, const Double_t */*dummy*/)
 782 {
 783
 784 //
 785 // J/Psi from B->J/Psi X
 786 //
 787 //
 788
 789
 790     Double_t c[7] = {7.37025e-02, 0., -2.94487e-03, 0., 6.07953e-06, 0., 5.39219e-07};
 791
 792     Double_t x = TMath::Abs(px[0]);
 793     Double_t y;
 794
 795     if (x > 6.) {
 796         y = 0.;
 797     } else {
 798         Int_t j;
 799         y = c[j = 6];
 800         while (j > 0) y  = y * x + c[--j];
 801     }
 802
 803     return y;
 804 }
 805
 806
 807
 808 //                 particle composition
 809 //
 810 Int_t AliGenMUONlib::IpJpsi(TRandom *)
 811 {
 812 // J/Psi composition
 813     return 443;
 814 }
 815 Int_t AliGenMUONlib::IpPsiP(TRandom *)
 816 {
 817 // Psi prime composition
 818     return 100443;
 819 }
 820 Int_t AliGenMUONlib::IpJpsiFamily(TRandom *)
 821 {
 822 // J/Psi composition
 823   Int_t ip;
 824   Float_t r = gRandom->Rndm();
 825   if (r < 0.98) {
 826     ip = 443;
 827   } else {
 828     ip = 100443;
 829   }
 830   return ip;
 831 }
 832
 833
 834
 835 //                      Upsilon
 836 //
 837 //
 838 //                  pt-distribution
 839 //____________________________________________________________
 840 Double_t AliGenMUONlib::PtUpsilon( const Double_t *px, const Double_t */*dummy*/ )
 841 {
 842 // Upsilon pT
 843   const Double_t kpt0 = 5.3;
 844   const Double_t kxn  = 2.5;
 845   Double_t x=*px;
 846   //
 847   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 848   return x/TMath::Power(pass1,kxn);
 849 }
 850
 851 Double_t AliGenMUONlib::PtUpsilonCDFscaled( const Double_t *px, const Double_t */*dummy*/ )
 852 {
 853 // Upsilon pT
 854   const Double_t kpt0 = 7.753;
 855   const Double_t kxn  = 3.042;
 856   Double_t x=*px;
 857   //
 858   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 859   return x/TMath::Power(pass1,kxn);
 860 }
 861
 862 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP( const Double_t *px, const Double_t */*dummy*/ )
 863 {
 864 // Upsilon pT
 865 //
 866 // pp 14 TeV
 867 //
 868 // scaled from CDF data at 2 TeV
 869
 870   const Double_t kpt0 = 8.610;
 871   const Double_t kxn  = 3.051;
 872   Double_t x=*px;
 873   //
 874   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 875   return x/TMath::Power(pass1,kxn);
 876 }
 877
 878 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
 879 {
 880 // Upsilon pT
 881 //
 882 // pp 10 TeV
 883 //
 884 // scaled from CDF data at 2 TeV
 885
 886   const Double_t kpt0 = 8.235;
 887   const Double_t kxn  = 3.051;
 888   Double_t x=*px;
 889   //
 890   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 891   return x/TMath::Power(pass1,kxn);
 892 }
 893
 894 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
 895 {
 896 // Upsilon pT
 897 //
 898 // pp 8.8 TeV
 899 // scaled from CDF data at 2 TeV
 900 //
 901   const Double_t kpt0 = 8.048;
 902   const Double_t kxn  = 3.051;
 903   Double_t x=*px;
 904   //
 905   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 906   return x/TMath::Power(pass1,kxn);
 907 }
 908
 909 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
 910 {
 911 // Upsilon pT
 912 //
 913 // pp 7 TeV
 914 //
 915 // scaled from CDF data at 2 TeV
 916
 917   const Double_t kpt0 = 7.817;
 918   const Double_t kxn  = 3.051;
 919   Double_t x=*px;
 920   //
 921   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 922   return x/TMath::Power(pass1,kxn);
 923 }
 924
 925 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
 926 {
 927 // Upsilon pT
 928 //
 929 // pp 3.94 TeV
 930 // scaled from CDF data at 2 TeV
 931 //
 932   const Double_t kpt0 = 7.189;
 933   const Double_t kxn  = 3.051;
 934   Double_t x=*px;
 935   //
 936   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 937   return x/TMath::Power(pass1,kxn);
 938 }
 939
 940 Double_t AliGenMUONlib::PtUpsilonCDFscaledPPb9( const Double_t *px, const Double_t *dummy)
 941 {
 942 // Upsilon pT
 943 //
 944 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.90
 945 //
 946   Double_t c[5] = {7.64952e-01, 1.12501e-04, 4.96038e-04, -3.03198e-05, 3.74035e-06};
 947   Double_t x=*px;
 948   Double_t y;
 949   Int_t j;
 950   y = c[j = 4];
 951   while (j > 0) y  = y * x + c[--j];
 952   //
 953   Double_t d = 1.+c[4]*TMath::Power(x,4);
 954   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP9(px,dummy);
 955 }
 956
 957 Double_t AliGenMUONlib::PtUpsilonCDFscaledPbP9( const Double_t *px, const Double_t *dummy)
 958 {
 959 // Upsilon pT
 960 //
 961 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.90
 962 //
 963   Double_t c[5] = {1.09881e+00, 3.08329e-03, -2.00356e-04, 8.28991e-06, 2.52576e-06};
 964   Double_t x=*px;
 965   Double_t y;
 966   Int_t j;
 967   y = c[j = 4];
 968   while (j > 0) y  = y * x + c[--j];
 969   //
 970   Double_t d = 1.+c[4]*TMath::Power(x,4);
 971   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP9(px,dummy);
 972 }
 973
 974 Double_t AliGenMUONlib::PtUpsilonCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
 975 {
 976 // Upsilon pT
 977 //
 978 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.85
 979 //
 980   Double_t c[5] = {8.65872e-01, 2.05465e-03, 2.56063e-04, -1.65598e-05, 2.29209e-06};
 981   Double_t x=*px;
 982   Double_t y;
 983   Int_t j;
 984   y = c[j = 4];
 985   while (j > 0) y  = y * x + c[--j];
 986   //
 987   Double_t d = 1.+c[4]*TMath::Power(x,4);
 988   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP4(px,dummy);
 989 }
 990
 991 Double_t AliGenMUONlib::PtUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/ )
 992 {
 993   return 1.;
 994 }
 995
 996 Double_t AliGenMUONlib::PtUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
 997 {
 998 //
 999 // Upsilon pT
1000 //
1001 //
1002 // R. Vogt 2002
1003 // PbPb 5.5 TeV
1004 // MRST HO
1005 // mc = 1.4 GeV, pt-kick 1 GeV
1006 //
1007     Float_t x = px[0];
1008     Double_t c[8] = {
1009         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,
1010         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
1011     };
1012     Double_t y;
1013     if (x < 10.) {
1014         Int_t j;
1015         y = c[j = 7];
1016         while (j > 0) y  = y * x +c[--j];
1017         y = x * TMath::Exp(y);
1018     } else {
1019         y = 0.;
1020     }
1021     return y;
1022 }
1023
1024 Double_t AliGenMUONlib::PtUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
1025 {
1026 //
1027 // Upsilon pT
1028 //
1029 //
1030 // R. Vogt 2002
1031 // pp 14 TeV
1032 // MRST HO
1033 // mc = 1.4 GeV, pt-kick 1 GeV
1034 //
1035     Float_t x = px[0];
1036     Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,
1037                      -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
1038
1039     Double_t y;
1040     if (x < 10.) {
1041         Int_t j;
1042         y = c[j = 7];
1043         while (j > 0) y  = y * x +c[--j];
1044         y = x * TMath::Exp(y);
1045     } else {
1046         y = 0.;
1047     }
1048     return y;
1049 }
1050
1051 //
1052 //                    y-distribution
1053 //
1054 //____________________________________________________________
1055 Double_t AliGenMUONlib::YUpsilon(const Double_t *py, const Double_t */*dummy*/)
1056 {
1057 // Upsilon y
1058   const Double_t ky0 = 3.;
1059   const Double_t kb=1.;
1060   Double_t yu;
1061   Double_t y=TMath::Abs(*py);
1062   //
1063   if (y < ky0)
1064     yu=kb;
1065   else
1066     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
1067   return yu;
1068 }
1069
1070
1071 Double_t AliGenMUONlib::YUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
1072 {
1073
1074 //
1075 // Upsilon y
1076 //
1077 //
1078 // R. Vogt 2002
1079 // PbPb 5.5 TeV
1080 // MRST HO
1081 // mc = 1.4 GeV, pt-kick 1 GeV
1082 //
1083
1084     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
1085                      -2.99753e-09, 1.28895e-05};
1086     Double_t x = TMath::Abs(px[0]);
1087     if (x > 5.55) return 0.;
1088     Int_t j;
1089     Double_t y = c[j = 6];
1090     while (j > 0) y  = y * x +c[--j];
1091     return y;
1092 }
1093
1094 Double_t AliGenMUONlib::YUpsilonCDFscaled( const Double_t *px, const Double_t *dummy)
1095 {
1096     // Upsilon y
1097     return AliGenMUONlib::YUpsilonPbPb(px, dummy);
1098
1099 }
1100
1101 Double_t AliGenMUONlib::YUpsilonCDFscaledPP( const Double_t *px, const Double_t *dummy)
1102 {
1103     // Upsilon y
1104     return AliGenMUONlib::YUpsilonPP(px, dummy);
1105
1106 }
1107
1108 Double_t AliGenMUONlib::YUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/)
1109 {
1110     // Upsilon y
1111     return 1.;
1112
1113 }
1114
1115 Double_t AliGenMUONlib::YUpsilonCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
1116 {
1117 // Upsilon y
1118 //
1119 // pp 10 TeV
1120 // scaled from YUpsilonPP(14 TeV) using 10 TeV / 14 TeV ratio of y-spectra in LO pQCD.
1121 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
1122 //
1123     Double_t c[4] = {1., -2.17877e-02, -6.52830e-04, 1.40578e-05};
1124     Double_t x = TMath::Abs(px[0]);
1125     if (x > 6.1) return 0.;
1126     Int_t j;
1127     Double_t y = c[j = 3];
1128     while (j > 0) y  = y * x*x +c[--j];
1129     return y;
1130 }
1131
1132 Double_t AliGenMUONlib::YUpsilonCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
1133 {
1134 // Upsilon y
1135 //
1136 // pp 8.8 TeV
1137 // rescaling of YUpsilonPP(14 TeV) using 8.8 TeV / 14 TeV ratio of y-spectra in LO QCD
1138 //
1139     Double_t c[4] = {1., -2.37621e-02, -6.29610e-04, 1.47976e-05};
1140     Double_t x = TMath::Abs(px[0]);
1141     if (x > 6.1) return 0.;
1142     Int_t j;
1143     Double_t y = c[j = 3];
1144     while (j > 0) y  = y * x*x +c[--j];
1145     return y;
1146 }
1147
1148 Double_t AliGenMUONlib::YUpsilonCDFscaledPP9dummy(Double_t px)
1149 {
1150     return AliGenMUONlib::YUpsilonCDFscaledPP9(&px, (Double_t*) 0);
1151 }
1152
1153 Double_t AliGenMUONlib::YUpsilonCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
1154 {
1155 // Upsilon y
1156 //
1157 // pp 7 TeV
1158 // scaled from YUpsilonPP(14 TeV) using 7 TeV / 14 TeV ratio of y-spectra in LO pQCD.
1159 //
1160     Double_t c[4] = {1., -2.61009e-02, -6.83937e-04, 1.78451e-05};
1161     Double_t x = TMath::Abs(px[0]);
1162     if (x > 6.0) return 0.;
1163     Int_t j;
1164     Double_t y = c[j = 3];
1165     while (j > 0) y  = y * x*x +c[--j];
1166     return y;
1167 }
1168
1169 Double_t AliGenMUONlib::YUpsilonCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
1170 {
1171 // Upsilon y
1172 //
1173 // pp 3.94 TeV
1174 // rescaling of YUpsilonPP(14 TeV) using 3.94 TeV / 14 TeV ratio of y-spectra in LO QCD
1175 //
1176     Double_t c[4] = {1., -3.91924e-02, -4.26184e-04, 2.10914e-05};
1177     Double_t x = TMath::Abs(px[0]);
1178     if (x > 5.7) return 0.;
1179     Int_t j;
1180     Double_t y = c[j = 3];
1181     while (j > 0) y  = y * x*x +c[--j];
1182
1183     return y;
1184 }
1185
1186 Double_t AliGenMUONlib::YUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
1187 {
1188
1189 //
1190 // Upsilon y
1191 //
1192 //
1193 // R. Vogt 2002
1194 // p p  14. TeV
1195 // MRST HO
1196 // mc = 1.4 GeV, pt-kick 1 GeV
1197 //
1198     Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04,
1199                      -6.20539e-10, 1.29943e-05};
1200     Double_t x = TMath::Abs(px[0]);
1201     if (x > 6.2) return 0.;
1202     Int_t j;
1203     Double_t y = c[j = 6];
1204     while (j > 0) y  = y * x +c[--j];
1205     return y;
1206 }
1207
1208 Double_t AliGenMUONlib::YUpsilonCDFscaledPPb9( const Double_t *px, const Double_t */*dummy*/)
1209 {
1210 // Upsilon y
1211 //
1212 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.90
1213 //
1214     Double_t c[7] = {8.71829e-01, 4.77467e-02, 8.09671e-03, 6.45294e-04, -2.15730e-04,
1215                      -4.67538e-05,-2.11683e-06};
1216     Double_t y;
1217     Double_t x = px[0] + 0.47;              // rapidity shift
1218     Int_t j;
1219     y = c[j = 6];
1220     while (j > 0) y = y * x + c[--j];
1221     if(y<0) y=0;
1222
1223     return y * AliGenMUONlib::YUpsilonCDFscaledPP9dummy(x);
1224 }
1225
1226 Double_t AliGenMUONlib::YUpsilonCDFscaledPbP9( const Double_t *px, const Double_t */*dummy*/)
1227 {
1228 // Upsilon y
1229 //
1230 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.90
1231 //
1232     Double_t c[7] = {8.71829e-01, 4.77467e-02, 8.09671e-03, 6.45294e-04, -2.15730e-04,
1233                      -4.67538e-05,-2.11683e-06};
1234     Double_t y;
1235     Double_t x = -px[0] + 0.47;              // rapidity shift
1236     Int_t j;
1237     y = c[j = 6];
1238     while (j > 0) y = y * x + c[--j];
1239     if(y<0) y=0;
1240
1241     return y * AliGenMUONlib::YUpsilonCDFscaledPP9dummy(x);
1242 }
1243
1244 Double_t AliGenMUONlib::YUpsilonCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
1245 {
1246 // Upsilon y
1247 //
1248 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.85
1249 //
1250     Double_t c[4] = {8.27837e-01, 1.70115e-02, -1.26046e-03, 1.52091e-05};
1251     Double_t x = px[0]*px[0];
1252     Double_t y;
1253     Int_t j;
1254     y = c[j = 3];
1255     while (j > 0) y  = y * x + c[--j];
1256     if(y<0) y=0;
1257
1258     return y * AliGenMUONlib::YUpsilonCDFscaledPP4(px,dummy);
1259 }
1260
1261
1262 //                 particle composition
1263 //
1264 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
1265 {
1266 // y composition
1267     return 553;
1268 }
1269 Int_t AliGenMUONlib::IpUpsilonP(TRandom *)
1270 {
1271 // y composition
1272     return 100553;
1273 }
1274 Int_t AliGenMUONlib::IpUpsilonPP(TRandom *)
1275 {
1276 // y composition
1277     return 200553;
1278 }
1279 Int_t AliGenMUONlib::IpUpsilonFamily(TRandom *)
1280 {
1281 // y composition
1282   Int_t ip;
1283   Float_t r = gRandom->Rndm();
1284
1285   if (r < 0.712) {
1286     ip = 553;
1287   } else if (r < 0.896) {
1288     ip = 100553;
1289   } else {
1290     ip = 200553;
1291   }
1292   return ip;
1293 }
1294
1295
1296 //
1297 //                        Phi
1298 //
1299 //
1300 //    pt-distribution (by scaling of pion distribution)
1301 //____________________________________________________________
1302 Double_t AliGenMUONlib::PtPhi( const Double_t *px, const Double_t */*dummy*/)
1303 {
1304 // Phi pT
1305   return PtScal(*px,7);
1306 }
1307 //    y-distribution
1308 Double_t AliGenMUONlib::YPhi( const Double_t *px, const Double_t */*dummy*/)
1309 {
1310 // Phi y
1311     Double_t *dum=0;
1312     return YJpsi(px,dum);
1313 }
1314 //                 particle composition
1315 //
1316 Int_t AliGenMUONlib::IpPhi(TRandom *)
1317 {
1318 // Phi composition
1319     return 333;
1320 }
1321
1322 //
1323 //                        omega
1324 //
1325 //
1326 //    pt-distribution (by scaling of pion distribution)
1327 //____________________________________________________________
1328 Double_t AliGenMUONlib::PtOmega( const Double_t *px, const Double_t */*dummy*/)
1329 {
1330 // Omega pT
1331   return PtScal(*px,5);
1332 }
1333 //    y-distribution
1334 Double_t AliGenMUONlib::YOmega( const Double_t *px, const Double_t */*dummy*/)
1335 {
1336 // Omega y
1337     Double_t *dum=0;
1338     return YJpsi(px,dum);
1339 }
1340 //                 particle composition
1341 //
1342 Int_t AliGenMUONlib::IpOmega(TRandom *)
1343 {
1344 // Omega composition
1345     return 223;
1346 }
1347
1348
1349 //
1350 //                        Eta
1351 //
1352 //
1353 //    pt-distribution (by scaling of pion distribution)
1354 //____________________________________________________________
1355 Double_t AliGenMUONlib::PtEta( const Double_t *px, const Double_t */*dummy*/)
1356 {
1357 // Eta pT
1358   return PtScal(*px,3);
1359 }
1360 //    y-distribution
1361 Double_t AliGenMUONlib::YEta( const Double_t *px, const Double_t */*dummy*/)
1362 {
1363 // Eta y
1364     Double_t *dum=0;
1365     return YJpsi(px,dum);
1366 }
1367 //                 particle composition
1368 //
1369 Int_t AliGenMUONlib::IpEta(TRandom *)
1370 {
1371 // Eta composition
1372     return 221;
1373 }
1374
1375 //
1376 //                        Charm
1377 //
1378 //
1379 //                    pt-distribution
1380 //____________________________________________________________
1381 Double_t AliGenMUONlib::PtCharm( const Double_t *px, const Double_t */*dummy*/)
1382 {
1383 // Charm pT
1384   const Double_t kpt0 = 2.25;
1385   const Double_t kxn  = 3.17;
1386   Double_t x=*px;
1387   //
1388   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1389   return x/TMath::Power(pass1,kxn);
1390 }
1391
1392 Double_t AliGenMUONlib::PtCharmCentral( const Double_t *px, const Double_t */*dummy*/)
1393 {
1394 // Charm pT
1395   const Double_t kpt0 = 2.12;
1396   const Double_t kxn  = 2.78;
1397   Double_t x=*px;
1398   //
1399   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1400   return x/TMath::Power(pass1,kxn);
1401 }
1402 Double_t AliGenMUONlib::PtCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1403 {
1404 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
1405 // PtCharmFiMjSkPP = PtCharmF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
1406 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
1407 // dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
1408 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
1409 // calculations for the following inputs:
1410 // Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11
1411 // for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV
1412 // for j=0,1 & 2 respectively;
1413 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
1414 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2
1415 // for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
1416 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
1417 // June 2008, Smbat.Grigoryan@cern.ch
1418
1419 // Charm pT
1420 // Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
1421 // for pp collisions at 14 TeV with one c-cbar pair per event.
1422 // Corresponding NLO total cross section is 5.68 mb
1423
1424
1425   const Double_t kpt0 = 2.2930;
1426   const Double_t kxn  = 3.1196;
1427   Double_t c[3]={-5.2180e-01,1.8753e-01,2.8669e-02};
1428   Double_t x=*px;
1429   //
1430   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1431   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1432 }
1433 Double_t AliGenMUONlib::PtCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1434 {
1435 // Charm pT
1436 // Corresponding NLO total cross section is 6.06 mb
1437   const Double_t kpt0 = 2.8669;
1438   const Double_t kxn  = 3.1044;
1439   Double_t c[3]={-4.6714e-01,1.5005e-01,4.5003e-02};
1440   Double_t x=*px;
1441   //
1442   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1443   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1444 }
1445 Double_t AliGenMUONlib::PtCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1446 {
1447 // Charm pT
1448 // Corresponding NLO total cross section is 6.06 mb
1449   const Double_t kpt0 = 1.8361;
1450   const Double_t kxn  = 3.2966;
1451   Double_t c[3]={-6.1550e-01,2.6498e-01,1.0728e-02};
1452   Double_t x=*px;
1453   //
1454   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1455   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1456 }
1457 Double_t AliGenMUONlib::PtCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
1458 {
1459 // Charm pT
1460 // Corresponding NLO total cross section is 7.69 mb
1461   const Double_t kpt0 = 2.1280;
1462   const Double_t kxn  = 3.1397;
1463   Double_t c[3]={-5.4021e-01,2.0944e-01,2.5211e-02};
1464   Double_t x=*px;
1465   //
1466   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1467   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1468 }
1469 Double_t AliGenMUONlib::PtCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
1470 {
1471 // Charm pT
1472 // Corresponding NLO total cross section is 4.81 mb
1473   const Double_t kpt0 = 2.4579;
1474   const Double_t kxn  = 3.1095;
1475   Double_t c[3]={-5.1497e-01,1.7532e-01,3.2429e-02};
1476   Double_t x=*px;
1477   //
1478   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1479   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1480 }
1481 Double_t AliGenMUONlib::PtCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
1482 {
1483 // Charm pT
1484 // Corresponding NLO total cross section is 14.09 mb
1485   const Double_t kpt0 = 2.1272;
1486   const Double_t kxn  = 3.1904;
1487   Double_t c[3]={-4.6088e-01,2.1918e-01,2.3055e-02};
1488   Double_t x=*px;
1489   //
1490   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1491   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1492 }
1493 Double_t AliGenMUONlib::PtCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
1494 {
1495 // Charm pT
1496 // Corresponding NLO total cross section is 1.52 mb
1497   const Double_t kpt0 = 2.8159;
1498   const Double_t kxn  = 3.0857;
1499   Double_t c[3]={-6.4691e-01,2.0289e-01,2.4922e-02};
1500   Double_t x=*px;
1501   //
1502   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1503   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1504 }
1505 Double_t AliGenMUONlib::PtCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
1506 {
1507 // Charm pT
1508 // Corresponding NLO total cross section is 3.67 mb
1509   const Double_t kpt0 = 2.7297;
1510   const Double_t kxn  = 3.3019;
1511   Double_t c[3]={-6.2216e-01,1.9031e-01,1.5341e-02};
1512   Double_t x=*px;
1513   //
1514   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1515   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1516 }
1517 Double_t AliGenMUONlib::PtCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
1518 {
1519 // Charm pT
1520 // Corresponding NLO total cross section is 3.38 mb
1521   const Double_t kpt0 = 2.3894;
1522   const Double_t kxn  = 3.1075;
1523   Double_t c[3]={-4.9742e-01,1.7032e-01,2.5994e-02};
1524   Double_t x=*px;
1525   //
1526   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1527   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1528 }
1529 Double_t AliGenMUONlib::PtCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
1530 {
1531 // Charm pT
1532 // Corresponding NLO total cross section is 10.37 mb
1533   const Double_t kpt0 = 2.0187;
1534   const Double_t kxn  = 3.3011;
1535   Double_t c[3]={-3.9869e-01,2.9248e-01,1.1763e-02};
1536   Double_t x=*px;
1537   //
1538   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1539   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1540 }
1541 Double_t AliGenMUONlib::PtCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
1542 {
1543 // Charm pT
1544 // Corresponding NLO total cross section is 7.22 mb
1545   const Double_t kpt0 = 2.1089;
1546   const Double_t kxn  = 3.1848;
1547   Double_t c[3]={-4.6275e-01,1.8114e-01,2.1363e-02};
1548   Double_t x=*px;
1549   //
1550   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1551   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1552 }
1553
1554 //                  y-distribution
1555 Double_t AliGenMUONlib::YCharm( const Double_t *px, const Double_t */*dummy*/)
1556 {
1557 // Charm y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225)
1558 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
1559 // shadowing + kt broadening
1560
1561     Double_t x=px[0];
1562     Double_t c[2]={-2.42985e-03,-2.31001e-04};
1563     Double_t y=1+(c[0]*TMath::Power(x,2))+(c[1]*TMath::Power(x,4));
1564     Double_t ycharm;
1565
1566     if (TMath::Abs(x)>8) {
1567       ycharm=0.;
1568     }
1569     else {
1570       ycharm=TMath::Power(y,3);
1571     }
1572
1573     return ycharm;
1574 }
1575 Double_t AliGenMUONlib::YCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1576 {
1577 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
1578 // YCharmFiMjSkPP = YCharmF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
1579 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
1580 // dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
1581 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
1582 // calculations for the following inputs:
1583 // Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11
1584 // for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV
1585 // for j=0,1 & 2 respectively;
1586 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
1587 // with a/b = 1/1,1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for
1588 // k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
1589 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
1590 // June 2008, Smbat.Grigoryan@cern.ch
1591
1592 // Charm y
1593 // Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
1594 // for pp collisions at 14 TeV with one c-cbar pair per event.
1595 // Corresponding NLO total cross section is 5.68 mb
1596
1597     Double_t x=px[0];
1598     Double_t c[2]={7.0909e-03,6.1967e-05};
1599     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1600     Double_t ycharm;
1601
1602     if (TMath::Abs(x)>9) {
1603       ycharm=0.;
1604     }
1605     else {
1606       ycharm=TMath::Power(y,3);
1607     }
1608
1609     return ycharm;
1610 }
1611 Double_t AliGenMUONlib::YCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1612 {
1613 // Charm y
1614 // Corresponding NLO total cross section is 6.06 mb
1615     Double_t x=px[0];
1616     Double_t c[2]={6.9707e-03,6.0971e-05};
1617     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1618     Double_t ycharm;
1619
1620     if (TMath::Abs(x)>9) {
1621       ycharm=0.;
1622     }
1623     else {
1624       ycharm=TMath::Power(y,3);
1625     }
1626
1627     return ycharm;
1628 }
1629 Double_t AliGenMUONlib::YCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1630 {
1631 // Charm y
1632 // Corresponding NLO total cross section is 6.06 mb
1633     Double_t x=px[0];
1634     Double_t c[2]={7.1687e-03,6.5303e-05};
1635     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1636     Double_t ycharm;
1637
1638     if (TMath::Abs(x)>9) {
1639       ycharm=0.;
1640     }
1641     else {
1642       ycharm=TMath::Power(y,3);
1643     }
1644
1645     return ycharm;
1646 }
1647 Double_t AliGenMUONlib::YCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
1648 {
1649 // Charm y
1650 // Corresponding NLO total cross section is 7.69 mb
1651     Double_t x=px[0];
1652     Double_t c[2]={5.9090e-03,7.1854e-05};
1653     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1654     Double_t ycharm;
1655
1656     if (TMath::Abs(x)>9) {
1657       ycharm=0.;
1658     }
1659     else {
1660       ycharm=TMath::Power(y,3);
1661     }
1662
1663     return ycharm;
1664 }
1665 Double_t AliGenMUONlib::YCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
1666 {
1667 // Charm y
1668 // Corresponding NLO total cross section is 4.81 mb
1669     Double_t x=px[0];
1670     Double_t c[2]={8.0882e-03,5.5872e-05};
1671     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1672     Double_t ycharm;
1673
1674     if (TMath::Abs(x)>9) {
1675       ycharm=0.;
1676     }
1677     else {
1678       ycharm=TMath::Power(y,3);
1679     }
1680
1681     return ycharm;
1682 }
1683 Double_t AliGenMUONlib::YCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
1684 {
1685 // Charm y
1686 // Corresponding NLO total cross section is 14.09 mb
1687     Double_t x=px[0];
1688     Double_t c[2]={7.2520e-03,6.2691e-05};
1689     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1690     Double_t ycharm;
1691
1692     if (TMath::Abs(x)>9) {
1693       ycharm=0.;
1694     }
1695     else {
1696       ycharm=TMath::Power(y,3);
1697     }
1698
1699     return ycharm;
1700 }
1701 Double_t AliGenMUONlib::YCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
1702 {
1703 // Charm y
1704 // Corresponding NLO total cross section is 1.52 mb
1705     Double_t x=px[0];
1706     Double_t c[2]={1.1040e-04,1.4498e-04};
1707     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1708     Double_t ycharm;
1709
1710     if (TMath::Abs(x)>9) {
1711       ycharm=0.;
1712     }
1713     else {
1714       ycharm=TMath::Power(y,3);
1715     }
1716
1717     return ycharm;
1718 }
1719 Double_t AliGenMUONlib::YCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
1720 {
1721 // Charm y
1722 // Corresponding NLO total cross section is 3.67 mb
1723     Double_t x=px[0];
1724     Double_t c[2]={-3.1328e-03,1.8270e-04};
1725     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1726     Double_t ycharm;
1727
1728     if (TMath::Abs(x)>9) {
1729       ycharm=0.;
1730     }
1731     else {
1732       ycharm=TMath::Power(y,3);
1733     }
1734
1735     return ycharm;
1736 }
1737 Double_t AliGenMUONlib::YCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
1738 {
1739 // Charm y
1740 // Corresponding NLO total cross section is 3.38 mb
1741     Double_t x=px[0];
1742     Double_t c[2]={7.0865e-03,6.2532e-05};
1743     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1744     Double_t ycharm;
1745
1746     if (TMath::Abs(x)>9) {
1747       ycharm=0.;
1748     }
1749     else {
1750       ycharm=TMath::Power(y,3);
1751     }
1752
1753     return ycharm;
1754 }
1755 Double_t AliGenMUONlib::YCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
1756 {
1757 // Charm y
1758 // Corresponding NLO total cross section is 10.37 mb
1759     Double_t x=px[0];
1760     Double_t c[2]={7.7070e-03,5.3533e-05};
1761     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1762     Double_t ycharm;
1763
1764     if (TMath::Abs(x)>9) {
1765       ycharm=0.;
1766     }
1767     else {
1768       ycharm=TMath::Power(y,3);
1769     }
1770
1771     return ycharm;
1772 }
1773 Double_t AliGenMUONlib::YCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
1774 {
1775 // Charm y
1776 // Corresponding NLO total cross section is 7.22 mb
1777     Double_t x=px[0];
1778     Double_t c[2]={7.9195e-03,5.3823e-05};
1779     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1780     Double_t ycharm;
1781
1782     if (TMath::Abs(x)>9) {
1783       ycharm=0.;
1784     }
1785     else {
1786       ycharm=TMath::Power(y,3);
1787     }
1788
1789     return ycharm;
1790 }
1791
1792
1793 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
1794 {
1795 // Charm composition
1796     Float_t random;
1797     Int_t ip;
1798 //    411,421,431,4122
1799     random = ran->Rndm();
1800 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3
1801 //  >>>>> cf. tab 4 p 11
1802
1803     if (random < 0.30) {
1804         ip=421;
1805     } else if (random < 0.60) {
1806         ip=-421;
1807     } else if (random < 0.70) {
1808         ip=411;
1809     } else if (random < 0.80) {
1810         ip=-411;
1811     } else if (random < 0.86) {
1812         ip=431;
1813     } else if (random < 0.92) {
1814         ip=-431;
1815     } else if (random < 0.96) {
1816         ip=4122;
1817     } else {
1818         ip=-4122;
1819     }
1820
1821     return ip;
1822 }
1823
1824 //
1825 //                        Beauty
1826 //
1827 //
1828 //                    pt-distribution
1829 //____________________________________________________________
1830 Double_t AliGenMUONlib::PtBeauty( const Double_t *px, const Double_t */*dummy*/)
1831 {
1832 // Beauty pT
1833   const Double_t kpt0 = 6.53;
1834   const Double_t kxn  = 3.59;
1835   Double_t x=*px;
1836   //
1837   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1838   return x/TMath::Power(pass1,kxn);
1839 }
1840
1841 Double_t AliGenMUONlib::PtBeautyCentral( const Double_t *px, const Double_t */*dummy*/)
1842 {
1843 // Beauty pT
1844   const Double_t kpt0 = 6.14;
1845   const Double_t kxn  = 2.93;
1846   Double_t x=*px;
1847   //
1848   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1849   return x/TMath::Power(pass1,kxn);
1850 }
1851 Double_t AliGenMUONlib::PtBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1852 {
1853 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
1854 // PtBeautyFiMjSkPP = PtBeautyF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
1855 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
1856 // dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
1857 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
1858 // calculations for the following inputs:
1859 // Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004
1860 // for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV
1861 // for j=0,1 & 2 respectively;
1862 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
1863 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for
1864 // k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
1865 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
1866 // June 2008, Smbat.Grigoryan@cern.ch
1867
1868 // Beauty pT
1869 // Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
1870 // for pp collisions at 14 TeV with one b-bbar pair per event.
1871 // Corresponding NLO total cross section is 0.494 mb
1872
1873   const Double_t kpt0 = 8.0575;
1874   const Double_t kxn  = 3.1921;
1875   Double_t x=*px;
1876   //
1877   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1878   return x/TMath::Power(pass1,kxn);
1879 }
1880 Double_t AliGenMUONlib::PtBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1881 {
1882 // Beauty pT
1883 // Corresponding NLO total cross section is 0.445 mb
1884   const Double_t kpt0 = 8.6239;
1885   const Double_t kxn  = 3.2911;
1886   Double_t x=*px;
1887   //
1888   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1889   return x/TMath::Power(pass1,kxn);
1890 }
1891 Double_t AliGenMUONlib::PtBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1892 {
1893 // Beauty pT
1894 // Corresponding NLO total cross section is 0.445 mb
1895   const Double_t kpt0 = 7.3367;
1896   const Double_t kxn  = 3.0692;
1897   Double_t x=*px;
1898   //
1899   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1900   return x/TMath::Power(pass1,kxn);
1901 }
1902 Double_t AliGenMUONlib::PtBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
1903 {
1904 // Beauty pT
1905 // Corresponding NLO total cross section is 0.518 mb
1906   const Double_t kpt0 = 7.6409;
1907   const Double_t kxn  = 3.1364;
1908   Double_t x=*px;
1909   //
1910   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1911   return x/TMath::Power(pass1,kxn);
1912 }
1913 Double_t AliGenMUONlib::PtBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
1914 {
1915 // Beauty pT
1916 // Corresponding NLO total cross section is 0.384 mb
1917   const Double_t kpt0 = 8.4948;
1918   const Double_t kxn  = 3.2546;
1919   Double_t x=*px;
1920   //
1921   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1922   return x/TMath::Power(pass1,kxn);
1923 }
1924 Double_t AliGenMUONlib::PtBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
1925 {
1926 // Beauty pT
1927 // Corresponding NLO total cross section is 0.648 mb
1928   const Double_t kpt0 = 7.6631;
1929   const Double_t kxn  = 3.1621;
1930   Double_t x=*px;
1931   //
1932   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1933   return x/TMath::Power(pass1,kxn);
1934 }
1935 Double_t AliGenMUONlib::PtBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
1936 {
1937 // Beauty pT
1938 // Corresponding NLO total cross section is 0.294 mb
1939   const Double_t kpt0 = 8.7245;
1940   const Double_t kxn  = 3.2213;
1941   Double_t x=*px;
1942   //
1943   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1944   return x/TMath::Power(pass1,kxn);
1945 }
1946 Double_t AliGenMUONlib::PtBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
1947 {
1948 // Beauty pT
1949 // Corresponding NLO total cross section is 0.475 mb
1950   const Double_t kpt0 = 8.5296;
1951   const Double_t kxn  = 3.2187;
1952   Double_t x=*px;
1953   //
1954   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1955   return x/TMath::Power(pass1,kxn);
1956 }
1957 Double_t AliGenMUONlib::PtBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
1958 {
1959 // Beauty pT
1960 // Corresponding NLO total cross section is 0.324 mb
1961   const Double_t kpt0 = 7.9440;
1962   const Double_t kxn  = 3.1614;
1963   Double_t x=*px;
1964   //
1965   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1966   return x/TMath::Power(pass1,kxn);
1967 }
1968 Double_t AliGenMUONlib::PtBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
1969 {
1970 // Beauty pT
1971 // Corresponding NLO total cross section is 0.536 mb
1972   const Double_t kpt0 = 8.2408;
1973   const Double_t kxn  = 3.3029;
1974   Double_t x=*px;
1975   //
1976   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1977   return x/TMath::Power(pass1,kxn);
1978 }
1979 Double_t AliGenMUONlib::PtBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
1980 {
1981 // Beauty pT
1982 // Corresponding NLO total cross section is 0.420 mb
1983   const Double_t kpt0 = 7.8041;
1984   const Double_t kxn  = 3.2094;
1985   Double_t x=*px;
1986   //
1987   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1988   return x/TMath::Power(pass1,kxn);
1989 }
1990
1991 //                     y-distribution
1992 Double_t AliGenMUONlib::YBeauty( const Double_t *px, const Double_t */*dummy*/)
1993 {
1994 // Beauty y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225)
1995 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
1996 // shadowing + kt broadening
1997
1998     Double_t x=px[0];
1999     Double_t c[2]={-1.27590e-02,-2.42731e-04};
2000     Double_t y=1+c[0]*TMath::Power(x,2)+c[1]*TMath::Power(x,4);
2001     Double_t ybeauty;
2002
2003     if (TMath::Abs(x)>6) {
2004       ybeauty=0.;
2005     }
2006     else {
2007       ybeauty=TMath::Power(y,3);
2008     }
2009
2010     return ybeauty;
2011 }
2012 Double_t AliGenMUONlib::YBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
2013 {
2014 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
2015 // YBeautyFiMjSkPP = YBeautyF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
2016 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
2017 // dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
2018 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
2019 // calculations for the following inputs:
2020 // Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004
2021 // for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV
2022 // for j=0,1 & 2 respectively;
2023 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
2024 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2
2025 // for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
2026 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
2027 // June 2008, Smbat.Grigoryan@cern.ch
2028
2029 // Beauty y
2030 // Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
2031 // for pp collisions at 14 TeV with one b-bbar pair per event.
2032 // Corresponding NLO total cross section is 0.494 mb
2033
2034
2035     Double_t x=px[0];
2036     Double_t c[2]={1.2350e-02,9.2667e-05};
2037     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
2038     Double_t ybeauty;
2039
2040     if (TMath::Abs(x)>7.6) {
2041       ybeauty=0.;
2042     }
2043     else {
2044       ybeauty=TMath::Power(y,3);
2045     }
2046
2047     return ybeauty;
2048 }
2049 Double_t AliGenMUONlib::YBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
2050 {
2051 // Beauty y
2052 // Corresponding NLO total cross section is 0.445 mb
2053     Double_t x=px[0];
2054     Double_t c[2]={1.2292e-02,9.1847e-05};
2055     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
2056     Double_t ybeauty;
2057
2058     if (TMath::Abs(x)>7.6) {
2059       ybeauty=0.;
2060     }
2061     else {
2062       ybeauty=TMath::Power(y,3);
2063     }
2064
2065     return ybeauty;
2066 }
2067 Double_t AliGenMUONlib::YBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
2068 {
2069 // Beauty y
2070 // Corresponding NLO total cross section is 0.445 mb
2071     Double_t x=px[0];
2072     Double_t c[2]={1.2436e-02,9.3709e-05};
2073     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
2074     Double_t ybeauty;
2075
2076     if (TMath::Abs(x)>7.6) {
2077       ybeauty=0.;
2078     }
2079     else {
2080       ybeauty=TMath::Power(y,3);
2081     }
2082
2083     return ybeauty;
2084 }
2085 Double_t AliGenMUONlib::YBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
2086 {
2087 // Beauty y
2088 // Corresponding NLO total cross section is 0.518 mb
2089     Double_t x=px[0];
2090     Double_t c[2]={1.1714e-02,1.0068e-04};
2091     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
2092     Double_t ybeauty;
2093
2094     if (TMath::Abs(x)>7.6) {
2095       ybeauty=0.;
2096     }
2097     else {
2098       ybeauty=TMath::Power(y,3);
2099     }
2100
2101     return ybeauty;
2102 }
2103 Double_t AliGenMUONlib::YBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
2104 {
2105 // Beauty y
2106 // Corresponding NLO total cross section is 0.384 mb
2107     Double_t x=px[0];
2108     Double_t c[2]={1.2944e-02,8.5500e-05};
2109     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
2110     Double_t ybeauty;
2111
2112     if (TMath::Abs(x)>7.6) {
2113       ybeauty=0.;
2114     }
2115     else {
2116       ybeauty=TMath::Power(y,3);
2117     }
2118
2119     return ybeauty;
2120 }
2121 Double_t AliGenMUONlib::YBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
2122 {
2123 // Beauty y
2124 // Corresponding NLO total cross section is 0.648 mb
2125     Double_t x=px[0];
2126     Double_t c[2]={1.2455e-02,9.2713e-05};
2127     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
2128     Double_t ybeauty;
2129
2130     if (TMath::Abs(x)>7.6) {
2131       ybeauty=0.;
2132     }
2133     else {
2134       ybeauty=TMath::Power(y,3);
2135     }
2136
2137     return ybeauty;
2138 }
2139 Double_t AliGenMUONlib::YBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
2140 {
2141 // Beauty y
2142 // Corresponding NLO total cross section is 0.294 mb
2143     Double_t x=px[0];
2144     Double_t c[2]={1.0897e-02,1.1878e-04};
2145     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
2146     Double_t ybeauty;
2147
2148     if (TMath::Abs(x)>7.6) {
2149       ybeauty=0.;
2150     }
2151     else {
2152       ybeauty=TMath::Power(y,3);
2153     }
2154
2155     return ybeauty;
2156 }
2157 Double_t AliGenMUONlib::YBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
2158 {
2159 // Beauty y
2160 // Corresponding NLO total cross section is 0.475 mb
2161     Double_t x=px[0];
2162     Double_t c[2]={1.0912e-02,1.1858e-04};
2163     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
2164     Double_t ybeauty;
2165
2166     if (TMath::Abs(x)>7.6) {
2167       ybeauty=0.;
2168     }
2169     else {
2170       ybeauty=TMath::Power(y,3);
2171     }
2172
2173     return ybeauty;
2174 }
2175 Double_t AliGenMUONlib::YBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
2176 {
2177 // Beauty y
2178 // Corresponding NLO total cross section is 0.324 mb
2179     Double_t x=px[0];
2180     Double_t c[2]={1.2378e-02,9.2490e-05};
2181     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
2182     Double_t ybeauty;
2183
2184     if (TMath::Abs(x)>7.6) {
2185       ybeauty=0.;
2186     }
2187     else {
2188       ybeauty=TMath::Power(y,3);
2189     }
2190
2191     return ybeauty;
2192 }
2193 Double_t AliGenMUONlib::YBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
2194 {
2195 // Beauty y
2196 // Corresponding NLO total cross section is 0.536 mb
2197     Double_t x=px[0];
2198     Double_t c[2]={1.2886e-02,8.2912e-05};
2199     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
2200     Double_t ybeauty;
2201
2202     if (TMath::Abs(x)>7.6) {
2203       ybeauty=0.;
2204     }
2205     else {
2206       ybeauty=TMath::Power(y,3);
2207     }
2208
2209     return ybeauty;
2210 }
2211 Double_t AliGenMUONlib::YBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
2212 {
2213 // Beauty y
2214 // Corresponding NLO total cross section is 0.420 mb
2215     Double_t x=px[0];
2216     Double_t c[2]={1.3106e-02,8.0115e-05};
2217     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
2218     Double_t ybeauty;
2219
2220     if (TMath::Abs(x)>7.6) {
2221       ybeauty=0.;
2222     }
2223     else {
2224       ybeauty=TMath::Power(y,3);
2225     }
2226
2227     return ybeauty;
2228 }
2229
2230 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
2231 {
2232 // Beauty Composition
2233     Float_t random;
2234     Int_t ip;
2235     random = ran->Rndm();
2236
2237 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3
2238 //  >>>>> cf. tab 4 p 11
2239
2240  if (random < 0.20) {
2241         ip=511;
2242     } else if (random < 0.40) {
2243         ip=-511;
2244     } else if (random < 0.605) {
2245         ip=521;
2246     } else if (random < 0.81) {
2247         ip=-521;
2248     } else if (random < 0.87) {
2249         ip=531;
2250     } else if (random < 0.93) {
2251         ip=-531;
2252     } else if (random < 0.965) {
2253         ip=5122;
2254     } else {
2255         ip=-5122;
2256     }
2257
2258  return ip;
2259 }
2260
2261
2262 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
2263 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
2264 {
2265 // Return pointer to pT parameterisation
2266     TString sname = TString(tname);
2267     GenFunc func;
2268     switch (param)
2269     {
2270     case kPhi:
2271         func=PtPhi;
2272         break;
2273     case kOmega:
2274         func=PtOmega;
2275         break;
2276     case kEta:
2277         func=PtEta;
2278         break;
2279     case kJpsiFamily:
2280     case kPsiP:
2281     case kChic1:
2282     case kChic2:
2283     case kJpsi:
2284         if (sname == "Vogt" || sname == "Vogt PbPb") {
2285             func=PtJpsiPbPb;
2286         } else if (sname == "Vogt pp") {
2287             func=PtJpsiPP;
2288         } else if (sname == "pp 7") {
2289             func=PtJpsiPP7000;
2290         } else if (sname == "pp 2.76") {
2291             func=PtJpsiPP2760;
2292         } else if (sname == "PbPb 2.76") {
2293             func=PtJpsiPbPb2760;
2294         } else if (sname == "CDF scaled") {
2295             func=PtJpsiCDFscaled;
2296         } else if (sname == "CDF pp") {
2297             func=PtJpsiCDFscaledPP;
2298         } else if (sname == "CDF pp 10") {
2299             func=PtJpsiCDFscaledPP10;
2300         } else if (sname == "CDF pp 8.8") {
2301             func=PtJpsiCDFscaledPP9;
2302         } else if (sname == "CDF pp 7" || sname == "CDF pp 7 flat y") {
2303             func=PtJpsiCDFscaledPP7;
2304         } else if (sname == "CDF pp 3.94") {
2305             func=PtJpsiCDFscaledPP4;
2306         } else if (sname == "CDF pp 2.76") {
2307             func=PtJpsiCDFscaledPP3;
2308         } else if (sname == "CDF pp 1.9") {
2309             func=PtJpsiCDFscaledPP2;
2310         } else if (sname == "CDF pPb 8.8") {
2311             func=PtJpsiCDFscaledPPb9;
2312         } else if (sname == "CDF Pbp 8.8") {
2313             func=PtJpsiCDFscaledPbP9;
2314         } else if (sname == "CDF PbPb 3.94") {
2315             func=PtJpsiCDFscaledPbPb4;
2316         } else if (sname == "Flat" || sname == "CDF pp 7 flat pt") {
2317             func=PtJpsiFlat;
2318         } else {
2319             func=PtJpsi;
2320         }
2321         break;
2322     case kJpsiFromB:
2323         func = PtJpsiBPbPb;
2324         break;
2325     case kUpsilonFamily:
2326     case kUpsilonP:
2327     case kUpsilonPP:
2328     case kUpsilon:
2329         if (sname == "Vogt" || sname == "Vogt PbPb") {
2330             func=PtUpsilonPbPb;
2331         } else if (sname == "Vogt pp") {
2332             func=PtUpsilonPP;
2333         } else if (sname == "CDF scaled") {
2334             func=PtUpsilonCDFscaled;
2335         } else if (sname == "CDF pp") {
2336             func=PtUpsilonCDFscaledPP;
2337         } else if (sname == "CDF pp 10") {
2338             func=PtUpsilonCDFscaledPP10;
2339         } else if (sname == "CDF pp 8.8") {
2340             func=PtUpsilonCDFscaledPP9;
2341         } else if (sname == "CDF pp 7") {
2342             func=PtUpsilonCDFscaledPP7;
2343         } else if (sname == "CDF pp 3.94") {
2344             func=PtUpsilonCDFscaledPP4;
2345         } else if (sname == "CDF pPb 8.8") {
2346             func=PtUpsilonCDFscaledPPb9;
2347         } else if (sname == "CDF Pbp 8.8") {
2348             func=PtUpsilonCDFscaledPbP9;
2349         } else if (sname == "CDF PbPb 3.94") {
2350             func=PtUpsilonCDFscaledPbPb4;
2351         } else if (sname == "Flat") {
2352             func=PtUpsilonFlat;
2353         } else {
2354             func=PtUpsilon;
2355         }
2356         break;
2357     case kCharm:
2358         if (sname == "F0M0S0 pp") {
2359             func=PtCharmF0M0S0PP;
2360         } else if (sname == "F1M0S0 pp") {
2361             func=PtCharmF1M0S0PP;
2362         } else if (sname == "F2M0S0 pp") {
2363             func=PtCharmF2M0S0PP;
2364         } else if (sname == "F0M1S0 pp") {
2365             func=PtCharmF0M1S0PP;
2366         } else if (sname == "F0M2S0 pp") {
2367             func=PtCharmF0M2S0PP;
2368         } else if (sname == "F0M0S1 pp") {
2369             func=PtCharmF0M0S1PP;
2370         } else if (sname == "F0M0S2 pp") {
2371             func=PtCharmF0M0S2PP;
2372         } else if (sname == "F0M0S3 pp") {
2373             func=PtCharmF0M0S3PP;
2374         } else if (sname == "F0M0S4 pp") {
2375             func=PtCharmF0M0S4PP;
2376         } else if (sname == "F0M0S5 pp") {
2377             func=PtCharmF0M0S5PP;
2378         } else if (sname == "F0M0S6 pp") {
2379             func=PtCharmF0M0S6PP;
2380         } else if (sname == "central") {
2381             func=PtCharmCentral;
2382         } else {
2383             func=PtCharm;
2384         }
2385         break;
2386     case kBeauty:
2387         if (sname == "F0M0S0 pp") {
2388             func=PtBeautyF0M0S0PP;
2389         } else if (sname == "F1M0S0 pp") {
2390             func=PtBeautyF1M0S0PP;
2391         } else if (sname == "F2M0S0 pp") {
2392             func=PtBeautyF2M0S0PP;
2393         } else if (sname == "F0M1S0 pp") {
2394             func=PtBeautyF0M1S0PP;
2395         } else if (sname == "F0M2S0 pp") {
2396             func=PtBeautyF0M2S0PP;
2397         } else if (sname == "F0M0S1 pp") {
2398             func=PtBeautyF0M0S1PP;
2399         } else if (sname == "F0M0S2 pp") {
2400             func=PtBeautyF0M0S2PP;
2401         } else if (sname == "F0M0S3 pp") {
2402             func=PtBeautyF0M0S3PP;
2403         } else if (sname == "F0M0S4 pp") {
2404             func=PtBeautyF0M0S4PP;
2405         } else if (sname == "F0M0S5 pp") {
2406             func=PtBeautyF0M0S5PP;
2407         } else if (sname == "F0M0S6 pp") {
2408             func=PtBeautyF0M0S6PP;
2409         } else if (sname == "central") {
2410             func=PtBeautyCentral;
2411         } else {
2412             func=PtBeauty;
2413         }
2414         break;
2415     case kPion:
2416         func=PtPion;
2417         break;
2418     case kKaon:
2419         func=PtKaon;
2420         break;
2421     case kChic0:
2422         func=PtChic0;
2423         break;
2424     case kChic:
2425         func=PtChic;
2426         break;
2427     default:
2428         func=0;
2429         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
2430     }
2431     return func;
2432 }
2433
2434 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
2435 {
2436   //
2437   // Return pointer to y- parameterisation
2438   //
2439     TString sname = TString(tname);
2440     GenFunc func;
2441     switch (param)
2442     {
2443     case kPhi:
2444         func=YPhi;
2445         break;
2446     case kEta:
2447         func=YEta;
2448         break;
2449     case kOmega:
2450         func=YOmega;
2451         break;
2452     case kJpsiFamily:
2453     case kPsiP:
2454     case kChic1:
2455     case kChic2:
2456     case kJpsi:
2457         if (sname == "Vogt" || sname == "Vogt PbPb") {
2458             func=YJpsiPbPb;
2459         } else if (sname == "Vogt pp"){
2460             func=YJpsiPP;
2461         } else if (sname == "pp 7") {
2462             func=YJpsiPP7000;
2463         } else if (sname == "pp 2.76") {
2464             func=YJpsiPP2760;
2465         } else if (sname == "PbPb 2.76") {
2466             func=YJpsiPbPb2760;
2467         } else if (sname == "CDF scaled") {
2468             func=YJpsiCDFscaled;
2469         } else if (sname == "CDF pp") {
2470             func=YJpsiCDFscaledPP;
2471         } else if (sname == "CDF pp 10") {
2472             func=YJpsiCDFscaledPP10;
2473         } else if (sname == "CDF pp 8.8") {
2474             func=YJpsiCDFscaledPP9;
2475         } else if (sname == "CDF pp 7" || sname == "CDF pp 7 flat pt") {
2476             func=YJpsiCDFscaledPP7;
2477         } else if (sname == "CDF pp 3.94") {
2478             func=YJpsiCDFscaledPP4;
2479         } else if (sname == "CDF pp 2.76") {
2480             func=YJpsiCDFscaledPP3;
2481         } else if (sname == "CDF pp 1.9") {
2482             func=YJpsiCDFscaledPP2;
2483         } else if (sname == "CDF pPb 8.8") {
2484             func=YJpsiCDFscaledPPb9;
2485         } else if (sname == "CDF Pbp 8.8") {
2486             func=YJpsiCDFscaledPbP9;
2487         } else if (sname == "CDF PbPb 3.94") {
2488             func=YJpsiCDFscaledPbPb4;
2489         } else if (sname == "Flat" || sname == "CDF pp 7 flat y") {
2490             func=YJpsiFlat;
2491         } else {
2492             func=YJpsi;
2493         }
2494         break;
2495     case kJpsiFromB:
2496         func = YJpsiBPbPb;
2497         break;
2498     case kUpsilonFamily:
2499     case kUpsilonP:
2500     case kUpsilonPP:
2501     case kUpsilon:
2502         if (sname == "Vogt" || sname == "Vogt PbPb") {
2503             func=YUpsilonPbPb;
2504         } else if (sname == "Vogt pp") {
2505             func = YUpsilonPP;
2506         } else if (sname == "CDF scaled") {
2507             func=YUpsilonCDFscaled;
2508         } else if (sname == "CDF pp") {
2509             func=YUpsilonCDFscaledPP;
2510         } else if (sname == "CDF pp 10") {
2511             func=YUpsilonCDFscaledPP10;
2512         } else if (sname == "CDF pp 8.8") {
2513             func=YUpsilonCDFscaledPP9;
2514         } else if (sname == "CDF pp 7") {
2515             func=YUpsilonCDFscaledPP7;
2516         } else if (sname == "CDF pp 3.94") {
2517             func=YUpsilonCDFscaledPP4;
2518         } else if (sname == "CDF pPb 8.8") {
2519             func=YUpsilonCDFscaledPPb9;
2520         } else if (sname == "CDF Pbp 8.8") {
2521             func=YUpsilonCDFscaledPbP9;
2522         } else if (sname == "CDF PbPb 3.94") {
2523             func=YUpsilonCDFscaledPbPb4;
2524         } else if (sname == "Flat") {
2525             func=YUpsilonFlat;
2526         } else {
2527             func=YUpsilon;
2528         }
2529         break;
2530     case kCharm:
2531         if (sname == "F0M0S0 pp") {
2532             func=YCharmF0M0S0PP;
2533         } else if (sname == "F1M0S0 pp") {
2534             func=YCharmF1M0S0PP;
2535         } else if (sname == "F2M0S0 pp") {
2536             func=YCharmF2M0S0PP;
2537         } else if (sname == "F0M1S0 pp") {
2538             func=YCharmF0M1S0PP;
2539         } else if (sname == "F0M2S0 pp") {
2540             func=YCharmF0M2S0PP;
2541         } else if (sname == "F0M0S1 pp") {
2542             func=YCharmF0M0S1PP;
2543         } else if (sname == "F0M0S2 pp") {
2544             func=YCharmF0M0S2PP;
2545         } else if (sname == "F0M0S3 pp") {
2546             func=YCharmF0M0S3PP;
2547         } else if (sname == "F0M0S4 pp") {
2548             func=YCharmF0M0S4PP;
2549         } else if (sname == "F0M0S5 pp") {
2550             func=YCharmF0M0S5PP;
2551         } else if (sname == "F0M0S6 pp") {
2552             func=YCharmF0M0S6PP;
2553         } else {
2554             func=YCharm;
2555         }
2556         break;
2557     case kBeauty:
2558         if (sname == "F0M0S0 pp") {
2559             func=YBeautyF0M0S0PP;
2560         } else if (sname == "F1M0S0 pp") {
2561             func=YBeautyF1M0S0PP;
2562         } else if (sname == "F2M0S0 pp") {
2563             func=YBeautyF2M0S0PP;
2564         } else if (sname == "F0M1S0 pp") {
2565             func=YBeautyF0M1S0PP;
2566         } else if (sname == "F0M2S0 pp") {
2567             func=YBeautyF0M2S0PP;
2568         } else if (sname == "F0M0S1 pp") {
2569             func=YBeautyF0M0S1PP;
2570         } else if (sname == "F0M0S2 pp") {
2571             func=YBeautyF0M0S2PP;
2572         } else if (sname == "F0M0S3 pp") {
2573             func=YBeautyF0M0S3PP;
2574         } else if (sname == "F0M0S4 pp") {
2575             func=YBeautyF0M0S4PP;
2576         } else if (sname == "F0M0S5 pp") {
2577             func=YBeautyF0M0S5PP;
2578         } else if (sname == "F0M0S6 pp") {
2579             func=YBeautyF0M0S6PP;
2580         } else {
2581             func=YBeauty;
2582         }
2583         break;
2584     case kPion:
2585         func=YPion;
2586         break;
2587     case kKaon:
2588         func=YKaon;
2589         break;
2590     case kChic0:
2591         func=YChic0;
2592         break;
2593     case kChic:
2594         func=YChic;
2595         break;
2596     default:
2597         func=0;
2598         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
2599     }
2600     return func;
2601 }
2602
2603 //
2604 //                    Chi
2605 //
2606 //
2607 //                pt-distribution
2608 //____________________________________________________________
2609 Double_t AliGenMUONlib::PtChic0( const Double_t *px, const Double_t */*dummy*/)
2610 {
2611 // Chi_c1 pT
2612   const Double_t kpt0 = 4.;
2613   const Double_t kxn  = 3.6;
2614   Double_t x=*px;
2615   //
2616   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2617   return x/TMath::Power(pass1,kxn);
2618 }
2619 Double_t AliGenMUONlib::PtChic1( const Double_t *px, const Double_t */*dummy*/)
2620 {
2621 // Chi_c1 pT
2622   const Double_t kpt0 = 4.;
2623   const Double_t kxn  = 3.6;
2624   Double_t x=*px;
2625   //
2626   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2627   return x/TMath::Power(pass1,kxn);
2628 }
2629 Double_t AliGenMUONlib::PtChic2( const Double_t *px, const Double_t */*dummy*/)
2630 {
2631 // Chi_c2 pT
2632   const Double_t kpt0 = 4.;
2633   const Double_t kxn  = 3.6;
2634   Double_t x=*px;
2635   //
2636   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2637   return x/TMath::Power(pass1,kxn);
2638 }
2639 Double_t AliGenMUONlib::PtChic( const Double_t *px, const Double_t */*dummy*/)
2640 {
2641 // Chi_c family pT
2642   const Double_t kpt0 = 4.;
2643   const Double_t kxn  = 3.6;
2644   Double_t x=*px;
2645   //
2646   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2647   return x/TMath::Power(pass1,kxn);
2648 }
2649
2650 //
2651 //               y-distribution
2652 //____________________________________________________________
2653 Double_t AliGenMUONlib::YChic0(const Double_t *py, const Double_t */*dummy*/)
2654 {
2655 // Chi-1c y
2656   const Double_t ky0 = 4.;
2657  const Double_t kb=1.;
2658   Double_t yj;
2659   Double_t y=TMath::Abs(*py);
2660   //
2661   if (y < ky0)
2662     yj=kb;
2663   else
2664     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2665   return yj;
2666 }
2667
2668 Double_t AliGenMUONlib::YChic1(const Double_t *py, const Double_t */*dummy*/)
2669 {
2670 // Chi-1c y
2671   const Double_t ky0 = 4.;
2672   const Double_t kb=1.;
2673   Double_t yj;
2674   Double_t y=TMath::Abs(*py);
2675   //
2676   if (y < ky0)
2677     yj=kb;
2678   else
2679     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2680   return yj;
2681 }
2682
2683 Double_t AliGenMUONlib::YChic2(const Double_t *py, const Double_t */*dummy*/)
2684 {
2685 // Chi-2c y
2686   const Double_t ky0 = 4.;
2687   const Double_t kb=1.;
2688   Double_t yj;
2689   Double_t y=TMath::Abs(*py);
2690   //
2691   if (y < ky0)
2692     yj=kb;
2693   else
2694     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2695   return yj;
2696 }
2697
2698 Double_t AliGenMUONlib::YChic(const Double_t *py, const Double_t */*dummy*/)
2699 {
2700 // Chi_c family y
2701   const Double_t ky0 = 4.;
2702   const Double_t kb=1.;
2703   Double_t yj;
2704   Double_t y=TMath::Abs(*py);
2705   //
2706   if (y < ky0)
2707     yj=kb;
2708   else
2709     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2710   return yj;
2711 }
2712
2713 //                 particle composition
2714 //
2715 Int_t AliGenMUONlib::IpChic0(TRandom *)
2716 {
2717 // Chi composition
2718     return 10441;
2719 }
2720 //
2721 Int_t AliGenMUONlib::IpChic1(TRandom *)
2722 {
2723 // Chi composition
2724     return 20443;
2725 }
2726 Int_t AliGenMUONlib::IpChic2(TRandom *)
2727 {
2728 // Chi_c2 prime composition
2729     return 445;
2730 }
2731 Int_t AliGenMUONlib::IpChic(TRandom *)
2732 {
2733 // Chi composition
2734   Int_t ip;
2735   Float_t r = gRandom->Rndm();
2736   if (r < 0.001) {
2737     ip = 10441;
2738   } else if( r < 0.377 ) {
2739     ip = 20443;
2740   } else {
2741     ip = 445;
2742   }
2743   return ip;
2744 }
2745
2746
2747 //_____________________________________________________________
2748
2749 typedef Int_t (*GenFuncIp) (TRandom *);
2750 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* /*tname*/) const
2751 {
2752 // Return pointer to particle type parameterisation
2753     GenFuncIp func;
2754     switch (param)
2755     {
2756     case kPhi:
2757         func=IpPhi;
2758         break;
2759     case kEta:
2760         func=IpEta;
2761         break;
2762     case kOmega:
2763         func=IpOmega;
2764         break;
2765     case kJpsiFamily:
2766         func=IpJpsiFamily;
2767         break;
2768     case kPsiP:
2769         func=IpPsiP;
2770         break;
2771     case kJpsi:
2772     case kJpsiFromB:
2773         func=IpJpsi;
2774         break;
2775     case kUpsilon:
2776         func=IpUpsilon;
2777         break;
2778     case kUpsilonFamily:
2779       func=IpUpsilonFamily;
2780       break;
2781     case kUpsilonP:
2782         func=IpUpsilonP;
2783         break;
2784     case kUpsilonPP:
2785         func=IpUpsilonPP;
2786         break;
2787     case kCharm:
2788         func=IpCharm;
2789         break;
2790     case kBeauty:
2791         func=IpBeauty;
2792         break;
2793     case kPion:
2794         func=IpPion;
2795         break;
2796     case kKaon:
2797         func=IpKaon;
2798         break;
2799     case kChic0:
2800         func=IpChic0;
2801         break;
2802     case kChic1:
2803         func=IpChic1;
2804         break;
2805     case kChic2:
2806         func=IpChic2;
2807         break;
2808     case kChic:
2809         func=IpChic;
2810         break;
2811     default:
2812         func=0;
2813         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
2814     }
2815     return func;
2816 }
2817
2818
2819
2820 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0,
2821                                    Float_t dx,
2822                                    Int_t n, Int_t no)
2823 {
2824 //
2825 // Neville's alorithm for interpolation
2826 //
2827 // x:  x-value
2828 // y:  Input array
2829 // x0: minimum x
2830 // dx: step size
2831 //  n: number of data points
2832 // no: order of polynom
2833 //
2834     Float_t*  c = new Float_t[n];
2835     Float_t*  d = new Float_t[n];
2836     Int_t m, i;
2837     for (i = 0; i < n; i++) {
2838         c[i] = y[i];
2839         d[i] = y[i];
2840     }
2841
2842     Int_t   ns  = int((x - x0)/dx);
2843
2844     Float_t y1  = y[ns];
2845     ns--;
2846     for (m = 0; m < no; m++) {
2847         for (i = 0; i < n-m; i++) {
2848             Float_t ho = x0 + Float_t(i) * dx - x;
2849             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
2850             Float_t w  = c[i+1] - d[i];
2851             Float_t den = ho-hp;
2852             den = w/den;
2853             d[i] = hp * den;
2854             c[i] = ho * den;
2855         }
2856         Float_t dy;
2857
2858         if (2*ns < (n-m-1)) {
2859             dy  = c[ns+1];
2860         } else {
2861             dy  = d[ns--];
2862         }
2863         y1 += dy;}
2864     delete[] c;
2865     delete[] d;
2866
2867     return y1;
2868 }
2869
2870