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