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