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 /*
  17 $Log$
  18 Revision 1.20  2003/03/13 11:54:39  morsch
  19 Limited pT range for parameterized Upsilon and J/Psi pT distributions.
  20
  21 Revision 1.19  2003/02/24 16:46:11  morsch
  22 New parameterisation for Psi and Upsilon (PbPb)
  23
  24 Revision 1.18  2002/11/07 09:13:42  morsch
  25 Use "Vogt" to label new distributions.
  26
  27 Revision 1.17  2002/11/07 09:06:10  morsch
  28 J/Psi and Upsilon pt and y distributions from R. Vogt 2002 added.
  29
  30 Revision 1.16  2002/10/14 14:55:35  hristov
  31 Merging the VirtualMC branch to the main development branch (HEAD)
  32
  33 Revision 1.14.6.1  2002/06/10 14:57:41  hristov
  34 Merged with v3-08-02
  35
  36 Revision 1.15  2002/04/17 10:11:51  morsch
  37 Coding Rule violations corrected.
  38
  39 Revision 1.14  2002/02/22 17:26:43  morsch
  40 Eta and omega added.
  41
  42 Revision 1.13  2001/03/27 11:01:04  morsch
  43 Charm pt-distribution corrected. More realistic y-distribution for pi and K.
  44
  45 Revision 1.12  2001/03/09 13:01:41  morsch
  46 - enum constants for paramterisation type (particle family) moved to AliGen*lib.h
  47 - use AliGenGSIlib::kUpsilon, AliGenPHOSlib::kEtaPrime to access the constants
  48
  49 Revision 1.11  2000/11/30 07:12:50  alibrary
  50 Introducing new Rndm and QA classes
  51
  52 Revision 1.10  2000/06/29 21:08:27  morsch
  53 All paramatrisation libraries derive from the pure virtual base class AliGenLib.
  54 This allows to pass a pointer to a library directly to AliGenParam and avoids the
  55 use of function pointers in Config.C.
  56
  57 Revision 1.9  2000/06/14 15:20:56  morsch
  58 Include clean-up (IH)
  59
  60 Revision 1.8  2000/06/09 20:32:11  morsch
  61 All coding rule violations except RS3 corrected
  62
  63 Revision 1.7  2000/05/02 08:12:13  morsch
  64 Coding rule violations corrected.
  65
  66 Revision 1.6  1999/09/29 09:24:14  fca
  67 Introduction of the Copyright and cvs Log
  68
  69 */
  70
  71 // Library class for particle pt and y distributions used for
  72 // muon spectrometer simulations.
  73 // To be used with AliGenParam.
  74 // The following particle typed can be simulated:
  75 // pi, K, phi, omega, eta, J/Psi, Upsilon, charm and beauty mesons.
  76 //
  77 // andreas.morsch@cern.ch
  78 //
  79
  80 #include "TMath.h"
  81 #include "TRandom.h"
  82
  83 #include "AliGenMUONlib.h"
  84
  85 ClassImp(AliGenMUONlib)
  86 //
  87 //  Pions
  88 Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *dummy)
  89 {
  90 //
  91 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
  92 //     POWER LAW FOR PT > 500 MEV
  93 //     MT SCALING BELOW (T=160 MEV)
  94 //
  95   const Double_t kp0 = 1.3;
  96   const Double_t kxn = 8.28;
  97   const Double_t kxlim=0.5;
  98   const Double_t kt=0.160;
  99   const Double_t kxmpi=0.139;
 100   const Double_t kb=1.;
 101   Double_t y, y1, xmpi2, ynorm, a;
 102   Double_t x=*px;
 103   //
 104   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
 105   xmpi2=kxmpi*kxmpi;
 106   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
 107   a=ynorm/y1;
 108   if (x > kxlim)
 109     y=a*TMath::Power(kp0/(kp0+x),kxn);
 110   else
 111     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
 112   return y*x;
 113 }
 114 //
 115 // y-distribution
 116 //
 117 Double_t AliGenMUONlib::YPion( Double_t *py, Double_t *dummy)
 118 {
 119 // Pion y
 120   Double_t y=TMath::Abs(*py);
 121 /*
 122   const Double_t ka    = 7000.;
 123   const Double_t kdy   = 4.;
 124   Double_t ex = y*y/(2*kdy*kdy);
 125   return ka*TMath::Exp(-ex);
 126 */
 127   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
 128
 129 }
 130 //                 particle composition
 131 //
 132 Int_t AliGenMUONlib::IpPion(TRandom *ran)
 133 {
 134 // Pion composition
 135     if (ran->Rndm() < 0.5) {
 136         return  211;
 137     } else {
 138         return -211;
 139     }
 140 }
 141
 142 //____________________________________________________________
 143 //
 144 // Mt-scaling
 145
 146 Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
 147 {
 148   //    SCALING EN MASSE PAR RAPPORT A PTPI
 149   //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
 150   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
 151   //     VALUE MESON/PI AT 5 GEV
 152   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
 153   np--;
 154   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
 155   Double_t fmax2=f5/kfmax[np];
 156   // PIONS
 157   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
 158   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
 159                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
 160   return fmtscal*ptpion;
 161 }
 162 //
 163 // kaon
 164 //
 165 //                pt-distribution
 166 //____________________________________________________________
 167 Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *dummy)
 168 {
 169 // Kaon pT
 170   return PtScal(*px,2);
 171 }
 172
 173 // y-distribution
 174 //____________________________________________________________
 175 Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t *dummy)
 176 {
 177 // Kaon y
 178   Double_t y=TMath::Abs(*py);
 179 /*
 180   const Double_t ka    = 1000.;
 181   const Double_t kdy   = 4.;
 182   //
 183   Double_t ex = y*y/(2*kdy*kdy);
 184   return ka*TMath::Exp(-ex);
 185 */
 186
 187   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
 188 }
 189
 190 //                 particle composition
 191 //
 192 Int_t AliGenMUONlib::IpKaon(TRandom *ran)
 193 {
 194 // Kaon composition
 195     if (ran->Rndm() < 0.5) {
 196         return  321;
 197     } else {
 198         return -321;
 199     }
 200 }
 201
 202 //                    J/Psi
 203 //
 204 //
 205 //                pt-distribution
 206 //____________________________________________________________
 207 Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *dummy)
 208 {
 209 // J/Psi pT
 210   const Double_t kpt0 = 4.;
 211   const Double_t kxn  = 3.6;
 212   Double_t x=*px;
 213   //
 214   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 215   return x/TMath::Power(pass1,kxn);
 216 }
 217
 218 Double_t AliGenMUONlib::PtJpsiPbPb( Double_t *px, Double_t *dummy)
 219 {
 220 // J/Psi pT spectrum
 221 //
 222 // R. Vogt 2002
 223 // PbPb 5.5 TeV
 224 // MRST HO
 225 // mc = 1.4 GeV, pt-kick 1 GeV
 226 //
 227     Float_t x = px[0];
 228     Float_t c[8] = {
 229         -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00,
 230         -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
 231     };
 232
 233     Double_t y;
 234     if (x < 10.) {
 235         Int_t j;
 236         y = c[j = 7];
 237         while (j > 0) y  = y * x +c[--j];
 238         y = x * TMath::Exp(y);
 239     } else {
 240         y = 0.;
 241     }
 242     return y;
 243 }
 244 Double_t AliGenMUONlib::PtJpsiPP( Double_t *px, Double_t *dummy)
 245 {
 246 // J/Psi pT spectrum
 247 //
 248 // R. Vogt 2002
 249 // pp 14 TeV
 250 // MRST HO
 251 // mc = 1.4 GeV, pt-kick 1 GeV
 252 //
 253     Float_t x = px[0];
 254     Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
 255
 256     Double_t y;
 257     if (x < 10.) {
 258         Int_t j;
 259         y = c[j = 3];
 260         while (j > 0) y  = y * x +c[--j];
 261         y = x * TMath::Exp(y);
 262     } else {
 263         y = 0.;
 264     }
 265     return y;
 266 }
 267
 268 //
 269 //               y-distribution
 270 //____________________________________________________________
 271 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *dummy)
 272 {
 273 // J/psi y
 274   const Double_t ky0 = 4.;
 275   const Double_t kb=1.;
 276   Double_t yj;
 277   Double_t y=TMath::Abs(*py);
 278   //
 279   if (y < ky0)
 280     yj=kb;
 281   else
 282     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 283   return yj;
 284 }
 285
 286
 287 Double_t AliGenMUONlib::YJpsiPbPb( Double_t *px, Double_t *dummy)
 288 {
 289
 290 //
 291 // J/Psi y
 292 //
 293 //
 294 // R. Vogt 2002
 295 // PbPb 5.5 TeV
 296 // MRST HO
 297 // mc = 1.4 GeV, pt-kick 1 GeV
 298 //
 299     Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
 300     Double_t x = TMath::Abs(px[0]);
 301     Double_t y;
 302
 303     if (x < 4.) {
 304         y = 31.754;
 305     } else if (x < 6) {
 306         Int_t j;
 307         y = c[j = 4];
 308         while (j > 0) y  = y * x + c[--j];
 309     } else {
 310         y =0.;
 311     }
 312
 313     return y;
 314 }
 315
 316 Double_t AliGenMUONlib::YJpsiPP( Double_t *px, Double_t *dummy)
 317 {
 318
 319 //
 320 // J/Psi y
 321 //
 322 //
 323 // R. Vogt 2002
 324 // pp 14  TeV
 325 // MRST HO
 326 // mc = 1.4 GeV, pt-kick 1 GeV
 327 //
 328
 329     Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
 330     Double_t x = TMath::Abs(px[0]);
 331     Double_t y;
 332
 333     if (x < 2.5) {
 334         y = 96.455 - 0.8483 * x * x;
 335     } else if (x < 7.9) {
 336         Int_t j;
 337         y = c[j = 4];
 338         while (j > 0) y  = y * x + c[--j];
 339     } else {
 340         y =0.;
 341     }
 342
 343     return y;
 344 }
 345
 346 //                 particle composition
 347 //
 348 Int_t AliGenMUONlib::IpJpsi(TRandom *)
 349 {
 350 // J/Psi composition
 351     return 443;
 352 }
 353
 354 //                      Upsilon
 355 //
 356 //
 357 //                  pt-distribution
 358 //____________________________________________________________
 359 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t *dummy )
 360 {
 361 // Upsilon pT
 362   const Double_t kpt0 = 5.3;
 363   const Double_t kxn  = 2.5;
 364   Double_t x=*px;
 365   //
 366   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 367   return x/TMath::Power(pass1,kxn);
 368 }
 369
 370 Double_t AliGenMUONlib::PtUpsilonPbPb( Double_t *px, Double_t *dummy)
 371 {
 372
 373 //
 374 // Upsilon pT
 375 //
 376 //
 377 // R. Vogt 2002
 378 // PbPb 5.5 TeV
 379 // MRST HO
 380 // mc = 1.4 GeV, pt-kick 1 GeV
 381 //
 382     Float_t x = px[0];
 383     Double_t c[8] = {
 384         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,
 385         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
 386     };
 387     Double_t y;
 388     if (x < 10.) {
 389         Int_t j;
 390         y = c[j = 7];
 391         while (j > 0) y  = y * x +c[--j];
 392         y = x * TMath::Exp(y);
 393     } else {
 394         y = 0.;
 395     }
 396     return y;
 397 }
 398
 399 Double_t AliGenMUONlib::PtUpsilonPP( Double_t *px, Double_t *dummy)
 400 {
 401
 402 //
 403 // Upsilon pT
 404 //
 405 //
 406 // R. Vogt 2002
 407 // pp 14 TeV
 408 // MRST HO
 409 // mc = 1.4 GeV, pt-kick 1 GeV
 410 //
 411     Float_t x = px[0];
 412     Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,
 413                      -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
 414
 415     Double_t y;
 416     if (x < 10.) {
 417         Int_t j;
 418         y = c[j = 7];
 419         while (j > 0) y  = y * x +c[--j];
 420         y = x * TMath::Exp(y);
 421     } else {
 422         y = 0.;
 423     }
 424     return y;
 425 }
 426
 427 //
 428 //                    y-distribution
 429 //
 430 //____________________________________________________________
 431 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *dummy)
 432 {
 433 // Upsilon y
 434   const Double_t ky0 = 3.;
 435   const Double_t kb=1.;
 436   Double_t yu;
 437   Double_t y=TMath::Abs(*py);
 438   //
 439   if (y < ky0)
 440     yu=kb;
 441   else
 442     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 443   return yu;
 444 }
 445
 446
 447 Double_t AliGenMUONlib::YUpsilonPbPb( Double_t *px, Double_t *dummy)
 448 {
 449
 450 //
 451 // Upsilon y
 452 //
 453 //
 454 // R. Vogt 2002
 455 // PbPb 5.5 TeV
 456 // MRST HO
 457 // mc = 1.4 GeV, pt-kick 1 GeV
 458 //
 459
 460     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
 461                      -2.99753e-09, 1.28895e-05};
 462
 463     Double_t x = px[0];
 464     if (TMath::Abs(x) > 5.55) return 0.;
 465     Int_t j;
 466     Double_t y = c[j = 6];
 467     while (j > 0) y  = y * x +c[--j];
 468     return y;
 469 }
 470
 471 Double_t AliGenMUONlib::YUpsilonPP( Double_t *px, Double_t *dummy)
 472 {
 473
 474 //
 475 // Upsilon y
 476 //
 477 //
 478 // R. Vogt 2002
 479 // p p  14. TeV
 480 // MRST HO
 481 // mc = 1.4 GeV, pt-kick 1 GeV
 482 //
 483     Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04,
 484                      -6.20539e-10, 1.29943e-05};
 485
 486     Double_t x = px[0];
 487     if (TMath::Abs(x) > 6.2) return 0.;
 488     Int_t j;
 489     Double_t y = c[j = 6];
 490     while (j > 0) y  = y * x +c[--j];
 491     return y;
 492 }
 493
 494 //                 particle composition
 495 //
 496 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
 497 {
 498 // y composition
 499     return 553;
 500 }
 501
 502 //
 503 //                        Phi
 504 //
 505 //
 506 //    pt-distribution (by scaling of pion distribution)
 507 //____________________________________________________________
 508 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *dummy)
 509 {
 510 // Phi pT
 511   return PtScal(*px,7);
 512 }
 513 //    y-distribution
 514 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *dummy)
 515 {
 516 // Phi y
 517     Double_t *dum=0;
 518     return YJpsi(px,dum);
 519 }
 520 //                 particle composition
 521 //
 522 Int_t AliGenMUONlib::IpPhi(TRandom *)
 523 {
 524 // Phi composition
 525     return 333;
 526 }
 527
 528 //
 529 //                        omega
 530 //
 531 //
 532 //    pt-distribution (by scaling of pion distribution)
 533 //____________________________________________________________
 534 Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t *dummy)
 535 {
 536 // Omega pT
 537   return PtScal(*px,5);
 538 }
 539 //    y-distribution
 540 Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t *dummy)
 541 {
 542 // Omega y
 543     Double_t *dum=0;
 544     return YJpsi(px,dum);
 545 }
 546 //                 particle composition
 547 //
 548 Int_t AliGenMUONlib::IpOmega(TRandom *)
 549 {
 550 // Omega composition
 551     return 223;
 552 }
 553
 554
 555 //
 556 //                        Eta
 557 //
 558 //
 559 //    pt-distribution (by scaling of pion distribution)
 560 //____________________________________________________________
 561 Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t *dummy)
 562 {
 563 // Eta pT
 564   return PtScal(*px,3);
 565 }
 566 //    y-distribution
 567 Double_t AliGenMUONlib::YEta( Double_t *px, Double_t *dummy)
 568 {
 569 // Eta y
 570     Double_t *dum=0;
 571     return YJpsi(px,dum);
 572 }
 573 //                 particle composition
 574 //
 575 Int_t AliGenMUONlib::IpEta(TRandom *)
 576 {
 577 // Eta composition
 578     return 221;
 579 }
 580
 581 //
 582 //                        Charm
 583 //
 584 //
 585 //                    pt-distribution
 586 //____________________________________________________________
 587 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *dummy)
 588 {
 589 // Charm pT
 590   const Double_t kpt0 = 4.08;
 591   const Double_t kxn  = 9.40;
 592
 593   Double_t x=*px;
 594   //
 595   Double_t pass1 = 1.+(x/kpt0);
 596   return x/TMath::Power(pass1,kxn);
 597 }
 598 //                  y-distribution
 599 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *dummy)
 600 {
 601 // Charm y
 602     Double_t *dum=0;
 603     return YJpsi(px,dum);
 604 }
 605
 606 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
 607 {
 608 // Charm composition
 609     Float_t random;
 610     Int_t ip;
 611 //    411,421,431,4122
 612     random = ran->Rndm();
 613     if (random < 0.5) {
 614         ip=411;
 615     } else if (random < 0.75) {
 616         ip=421;
 617     } else if (random < 0.90) {
 618         ip=431;
 619     } else {
 620         ip=4122;
 621     }
 622     if (ran->Rndm() < 0.5) {ip=-ip;}
 623
 624     return ip;
 625 }
 626
 627
 628 //
 629 //                        Beauty
 630 //
 631 //
 632 //                    pt-distribution
 633 //____________________________________________________________
 634 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *dummy)
 635 {
 636 // Beauty pT
 637   const Double_t kpt0 = 4.;
 638   const Double_t kxn  = 3.6;
 639   Double_t x=*px;
 640   //
 641   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 642   return x/TMath::Power(pass1,kxn);
 643 }
 644 //                     y-distribution
 645 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *dummy)
 646 {
 647 // Beauty y
 648     Double_t *dum=0;
 649     return YJpsi(px,dum);
 650 }
 651
 652 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
 653 {
 654 // Beauty Composition
 655     Float_t random;
 656     Int_t ip;
 657     random = ran->Rndm();
 658     if (random < 0.5) {
 659         ip=511;
 660     } else if (random < 0.75) {
 661         ip=521;
 662     } else if (random < 0.90) {
 663         ip=531;
 664     } else {
 665         ip=5122;
 666     }
 667     if (ran->Rndm() < 0.5) {ip=-ip;}
 668
 669     return ip;
 670 }
 671
 672 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
 673 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
 674 {
 675 // Return pointer to pT parameterisation
 676     TString sname = TString(tname);
 677     GenFunc func;
 678     switch (param)
 679     {
 680     case kPhi:
 681         func=PtPhi;
 682         break;
 683     case kOmega:
 684         func=PtOmega;
 685         break;
 686     case kEta:
 687         func=PtEta;
 688         break;
 689     case kJpsi:
 690         if (sname == "Vogt" || sname == "Vogt PbPb") {
 691             func=PtJpsiPbPb;
 692         } else if (sname == "Vogt pp") {
 693             func=PtJpsiPP;
 694         } else {
 695             func=PtJpsi;
 696         }
 697         break;
 698     case kUpsilon:
 699         if (sname == "Vogt" || sname == "Vogt PbPb") {
 700             func=PtUpsilonPbPb;
 701         } else if (sname == "Vogt pp") {
 702             func=PtUpsilonPP;
 703         } else {
 704             func=PtUpsilon;
 705         }
 706         break;
 707     case kCharm:
 708         func=PtCharm;
 709         break;
 710     case kBeauty:
 711         func=PtBeauty;
 712         break;
 713     case kPion:
 714         func=PtPion;
 715         break;
 716     case kKaon:
 717         func=PtKaon;
 718         break;
 719     default:
 720         func=0;
 721         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
 722     }
 723     return func;
 724 }
 725
 726 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
 727 {
 728     TString sname = TString(tname);
 729
 730 // Return pointer to y- parameterisation
 731     GenFunc func;
 732     switch (param)
 733     {
 734     case kPhi:
 735         func=YPhi;
 736         break;
 737     case kEta:
 738         func=YEta;
 739         break;
 740     case kOmega:
 741         func=YOmega;
 742         break;
 743     case kJpsi:
 744         if (sname == "Vogt" || sname == "Vogt PbPb") {
 745             func=YJpsiPbPb;
 746         } else if (sname == "Vogt pp"){
 747             func=YJpsiPP;
 748         } else {
 749             func=YJpsi;
 750         }
 751
 752         break;
 753     case kUpsilon:
 754         if (sname == "Vogt" || sname == "Vogt PbPb") {
 755             func=YUpsilonPbPb;
 756         } else if (sname == "Vogt pp") {
 757             func = YUpsilonPP;
 758         } else {
 759             func=YUpsilon;
 760         }
 761         break;
 762     case kCharm:
 763         func=YCharm;
 764         break;
 765     case kBeauty:
 766         func=YBeauty;
 767         break;
 768     case kPion:
 769         func=YPion;
 770         break;
 771     case kKaon:
 772         func=YKaon;
 773         break;
 774     default:
 775         func=0;
 776         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
 777     }
 778     return func;
 779 }
 780 typedef Int_t (*GenFuncIp) (TRandom *);
 781 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* tname) const
 782 {
 783 // Return pointer to particle type parameterisation
 784     GenFuncIp func;
 785     switch (param)
 786     {
 787     case kPhi:
 788         func=IpPhi;
 789         break;
 790     case kEta:
 791         func=IpEta;
 792         break;
 793     case kOmega:
 794         func=IpOmega;
 795         break;
 796     case kJpsi:
 797         func=IpJpsi;
 798         break;
 799     case kUpsilon:
 800         func=IpUpsilon;
 801         break;
 802     case kCharm:
 803         func=IpCharm;
 804         break;
 805     case kBeauty:
 806         func=IpBeauty;
 807         break;
 808     case kPion:
 809         func=IpPion;
 810         break;
 811     case kKaon:
 812         func=IpKaon;
 813         break;
 814     default:
 815         func=0;
 816         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
 817     }
 818     return func;
 819 }
 820
 821
 822
 823 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0,
 824                                    Float_t dx,
 825                                    Int_t n, Int_t no)
 826 {
 827 //
 828 // Neville's alorithm for interpolation
 829 //
 830 // x:  x-value
 831 // y:  Input array
 832 // x0: minimum x
 833 // dx: step size
 834 //  n: number of data points
 835 // no: order of polynom
 836 //
 837     Float_t*  c = new Float_t[n];
 838     Float_t*  d = new Float_t[n];
 839     Int_t m, i;
 840     for (i = 0; i < n; i++) {
 841         c[i] = y[i];
 842         d[i] = y[i];
 843     }
 844
 845     Int_t   ns  = int((x - x0)/dx);
 846
 847     Float_t y1  = y[ns];
 848     ns--;
 849     for (m = 0; m < no; m++) {
 850         for (i = 0; i < n-m; i++) {
 851             Float_t ho = x0 + Float_t(i) * dx - x;
 852             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
 853             Float_t w  = c[i+1] - d[i];
 854             Float_t den = ho-hp;
 855             den = w/den;
 856             d[i] = hp * den;
 857             c[i] = ho * den;
 858         }
 859         Float_t dy;
 860
 861         if (2*ns < (n-m-1)) {
 862             dy  = c[ns+1];
 863         } else {
 864             dy  = d[ns--];
 865         }
 866         y1 += dy;}
 867     delete[] c;
 868     delete[] d;
 869
 870     return y1;
 871 }
 872
 873