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