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