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::PtJpsiPbPb( Double_t *px, Double_t */*dummy*/)
 166 {
 167 // J/Psi pT spectrum
 168 //
 169 // R. Vogt 2002
 170 // PbPb 5.5 TeV
 171 // MRST HO
 172 // mc = 1.4 GeV, pt-kick 1 GeV
 173 //
 174     Float_t x = px[0];
 175     Float_t c[8] = {
 176         -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00,
 177         -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
 178     };
 179
 180     Double_t y;
 181     if (x < 10.) {
 182         Int_t j;
 183         y = c[j = 7];
 184         while (j > 0) y  = y * x +c[--j];
 185         y = x * TMath::Exp(y);
 186     } else {
 187         y = 0.;
 188     }
 189     return y;
 190 }
 191
 192 Double_t AliGenMUONlib::PtJpsiBPbPb( Double_t *px, Double_t */*dummy*/)
 193 {
 194 // J/Psi pT spectrum
 195 // B -> J/Psi X
 196     Double_t x0 =   4.0384;
 197     Double_t  n =   3.0288;
 198
 199     Double_t x = px[0];
 200     Double_t y = x / TMath::Power((1. + (x/x0)*(x/x0)), n);
 201
 202     return y;
 203 }
 204
 205
 206 Double_t AliGenMUONlib::PtJpsiPP( Double_t *px, Double_t */*dummy*/)
 207 {
 208 // J/Psi pT spectrum
 209 //
 210 // R. Vogt 2002
 211 // pp 14 TeV
 212 // MRST HO
 213 // mc = 1.4 GeV, pt-kick 1 GeV
 214 //
 215     Float_t x = px[0];
 216     Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
 217
 218     Double_t y;
 219     if (x < 10.) {
 220         Int_t j;
 221         y = c[j = 3];
 222         while (j > 0) y  = y * x +c[--j];
 223         y = x * TMath::Exp(y);
 224     } else {
 225         y = 0.;
 226     }
 227     return y;
 228 }
 229
 230 //
 231 //               y-distribution
 232 //____________________________________________________________
 233 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t */*dummy*/)
 234 {
 235 // J/psi y
 236   const Double_t ky0 = 4.;
 237   const Double_t kb=1.;
 238   Double_t yj;
 239   Double_t y=TMath::Abs(*py);
 240   //
 241   if (y < ky0)
 242     yj=kb;
 243   else
 244     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 245   return yj;
 246 }
 247
 248
 249 Double_t AliGenMUONlib::YJpsiPbPb( Double_t *px, Double_t */*dummy*/)
 250 {
 251
 252 //
 253 // J/Psi y
 254 //
 255 //
 256 // R. Vogt 2002
 257 // PbPb 5.5 TeV
 258 // MRST HO
 259 // mc = 1.4 GeV, pt-kick 1 GeV
 260 //
 261     Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
 262     Double_t x = TMath::Abs(px[0]);
 263     Double_t y;
 264
 265     if (x < 4.) {
 266         y = 31.754;
 267     } else if (x < 6) {
 268         Int_t j;
 269         y = c[j = 4];
 270         while (j > 0) y  = y * x + c[--j];
 271     } else {
 272         y =0.;
 273     }
 274
 275     return y;
 276 }
 277
 278 Double_t AliGenMUONlib::YJpsiPP( Double_t *px, Double_t */*dummy*/)
 279 {
 280
 281 //
 282 // J/Psi y
 283 //
 284 //
 285 // R. Vogt 2002
 286 // pp 14  TeV
 287 // MRST HO
 288 // mc = 1.4 GeV, pt-kick 1 GeV
 289 //
 290
 291     Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
 292     Double_t x = TMath::Abs(px[0]);
 293     Double_t y;
 294
 295     if (x < 2.5) {
 296         y = 96.455 - 0.8483 * x * x;
 297     } else if (x < 7.9) {
 298         Int_t j;
 299         y = c[j = 4];
 300         while (j > 0) y  = y * x + c[--j];
 301     } else {
 302         y =0.;
 303     }
 304
 305     return y;
 306 }
 307
 308 Double_t AliGenMUONlib::YJpsiBPbPb( Double_t *px, Double_t */*dummy*/)
 309 {
 310
 311 //
 312 // J/Psi from B->J/Psi X
 313 //
 314 //
 315
 316
 317     Double_t c[7] = {7.37025e-02, 0., -2.94487e-03, 0., 6.07953e-06, 0., 5.39219e-07};
 318
 319     Double_t x = TMath::Abs(px[0]);
 320     Double_t y;
 321
 322     if (x > 6.) {
 323         y = 0.;
 324     } else {
 325         Int_t j;
 326         y = c[j = 6];
 327         while (j > 0) y  = y * x + c[--j];
 328     }
 329
 330     return y;
 331 }
 332
 333
 334
 335 //                 particle composition
 336 //
 337 Int_t AliGenMUONlib::IpJpsi(TRandom *)
 338 {
 339 // J/Psi composition
 340     return 443;
 341 }
 342
 343 //                      Upsilon
 344 //
 345 //
 346 //                  pt-distribution
 347 //____________________________________________________________
 348 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t */*dummy*/ )
 349 {
 350 // Upsilon pT
 351   const Double_t kpt0 = 5.3;
 352   const Double_t kxn  = 2.5;
 353   Double_t x=*px;
 354   //
 355   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 356   return x/TMath::Power(pass1,kxn);
 357 }
 358
 359 Double_t AliGenMUONlib::PtUpsilonPbPb( Double_t *px, Double_t */*dummy*/)
 360 {
 361
 362 //
 363 // Upsilon pT
 364 //
 365 //
 366 // R. Vogt 2002
 367 // PbPb 5.5 TeV
 368 // MRST HO
 369 // mc = 1.4 GeV, pt-kick 1 GeV
 370 //
 371     Float_t x = px[0];
 372     Double_t c[8] = {
 373         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,
 374         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
 375     };
 376     Double_t y;
 377     if (x < 10.) {
 378         Int_t j;
 379         y = c[j = 7];
 380         while (j > 0) y  = y * x +c[--j];
 381         y = x * TMath::Exp(y);
 382     } else {
 383         y = 0.;
 384     }
 385     return y;
 386 }
 387
 388 Double_t AliGenMUONlib::PtUpsilonPP( Double_t *px, Double_t */*dummy*/)
 389 {
 390
 391 //
 392 // Upsilon pT
 393 //
 394 //
 395 // R. Vogt 2002
 396 // pp 14 TeV
 397 // MRST HO
 398 // mc = 1.4 GeV, pt-kick 1 GeV
 399 //
 400     Float_t x = px[0];
 401     Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,
 402                      -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
 403
 404     Double_t y;
 405     if (x < 10.) {
 406         Int_t j;
 407         y = c[j = 7];
 408         while (j > 0) y  = y * x +c[--j];
 409         y = x * TMath::Exp(y);
 410     } else {
 411         y = 0.;
 412     }
 413     return y;
 414 }
 415
 416 //
 417 //                    y-distribution
 418 //
 419 //____________________________________________________________
 420 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t */*dummy*/)
 421 {
 422 // Upsilon y
 423   const Double_t ky0 = 3.;
 424   const Double_t kb=1.;
 425   Double_t yu;
 426   Double_t y=TMath::Abs(*py);
 427   //
 428   if (y < ky0)
 429     yu=kb;
 430   else
 431     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 432   return yu;
 433 }
 434
 435
 436 Double_t AliGenMUONlib::YUpsilonPbPb( Double_t *px, Double_t */*dummy*/)
 437 {
 438
 439 //
 440 // Upsilon y
 441 //
 442 //
 443 // R. Vogt 2002
 444 // PbPb 5.5 TeV
 445 // MRST HO
 446 // mc = 1.4 GeV, pt-kick 1 GeV
 447 //
 448
 449     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
 450                      -2.99753e-09, 1.28895e-05};
 451
 452     Double_t x = px[0];
 453     if (TMath::Abs(x) > 5.55) return 0.;
 454     Int_t j;
 455     Double_t y = c[j = 6];
 456     while (j > 0) y  = y * x +c[--j];
 457     return y;
 458 }
 459
 460 Double_t AliGenMUONlib::YUpsilonPP( Double_t *px, Double_t */*dummy*/)
 461 {
 462
 463 //
 464 // Upsilon y
 465 //
 466 //
 467 // R. Vogt 2002
 468 // p p  14. TeV
 469 // MRST HO
 470 // mc = 1.4 GeV, pt-kick 1 GeV
 471 //
 472     Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04,
 473                      -6.20539e-10, 1.29943e-05};
 474
 475     Double_t x = px[0];
 476     if (TMath::Abs(x) > 6.2) return 0.;
 477     Int_t j;
 478     Double_t y = c[j = 6];
 479     while (j > 0) y  = y * x +c[--j];
 480     return y;
 481 }
 482
 483 //                 particle composition
 484 //
 485 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
 486 {
 487 // y composition
 488     return 553;
 489 }
 490
 491 //
 492 //                        Phi
 493 //
 494 //
 495 //    pt-distribution (by scaling of pion distribution)
 496 //____________________________________________________________
 497 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t */*dummy*/)
 498 {
 499 // Phi pT
 500   return PtScal(*px,7);
 501 }
 502 //    y-distribution
 503 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t */*dummy*/)
 504 {
 505 // Phi y
 506     Double_t *dum=0;
 507     return YJpsi(px,dum);
 508 }
 509 //                 particle composition
 510 //
 511 Int_t AliGenMUONlib::IpPhi(TRandom *)
 512 {
 513 // Phi composition
 514     return 333;
 515 }
 516
 517 //
 518 //                        omega
 519 //
 520 //
 521 //    pt-distribution (by scaling of pion distribution)
 522 //____________________________________________________________
 523 Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t */*dummy*/)
 524 {
 525 // Omega pT
 526   return PtScal(*px,5);
 527 }
 528 //    y-distribution
 529 Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t */*dummy*/)
 530 {
 531 // Omega y
 532     Double_t *dum=0;
 533     return YJpsi(px,dum);
 534 }
 535 //                 particle composition
 536 //
 537 Int_t AliGenMUONlib::IpOmega(TRandom *)
 538 {
 539 // Omega composition
 540     return 223;
 541 }
 542
 543
 544 //
 545 //                        Eta
 546 //
 547 //
 548 //    pt-distribution (by scaling of pion distribution)
 549 //____________________________________________________________
 550 Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t */*dummy*/)
 551 {
 552 // Eta pT
 553   return PtScal(*px,3);
 554 }
 555 //    y-distribution
 556 Double_t AliGenMUONlib::YEta( Double_t *px, Double_t */*dummy*/)
 557 {
 558 // Eta y
 559     Double_t *dum=0;
 560     return YJpsi(px,dum);
 561 }
 562 //                 particle composition
 563 //
 564 Int_t AliGenMUONlib::IpEta(TRandom *)
 565 {
 566 // Eta composition
 567     return 221;
 568 }
 569
 570 //
 571 //                        Charm
 572 //
 573 //
 574 //                    pt-distribution
 575 //____________________________________________________________
 576 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t */*dummy*/)
 577 {
 578 // Charm pT
 579   const Double_t kpt0 = 4.08;
 580   const Double_t kxn  = 9.40;
 581
 582   Double_t x=*px;
 583   //
 584   Double_t pass1 = 1.+(x/kpt0);
 585   return x/TMath::Power(pass1,kxn);
 586 }
 587 //                  y-distribution
 588 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t */*dummy*/)
 589 {
 590 // Charm y
 591     Double_t *dum=0;
 592     return YJpsi(px,dum);
 593 }
 594
 595 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
 596 {
 597 // Charm composition
 598     Float_t random;
 599     Int_t ip;
 600 //    411,421,431,4122
 601     random = ran->Rndm();
 602     if (random < 0.5) {
 603         ip=411;
 604     } else if (random < 0.75) {
 605         ip=421;
 606     } else if (random < 0.90) {
 607         ip=431;
 608     } else {
 609         ip=4122;
 610     }
 611     if (ran->Rndm() < 0.5) {ip=-ip;}
 612
 613     return ip;
 614 }
 615
 616
 617 //
 618 //                        Beauty
 619 //
 620 //
 621 //                    pt-distribution
 622 //____________________________________________________________
 623 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t */*dummy*/)
 624 {
 625 // Beauty pT
 626   const Double_t kpt0 = 4.;
 627   const Double_t kxn  = 3.6;
 628   Double_t x=*px;
 629   //
 630   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 631   return x/TMath::Power(pass1,kxn);
 632 }
 633 //                     y-distribution
 634 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t */*dummy*/)
 635 {
 636 // Beauty y
 637     Double_t *dum=0;
 638     return YJpsi(px,dum);
 639 }
 640
 641 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
 642 {
 643 // Beauty Composition
 644     Float_t random;
 645     Int_t ip;
 646     random = ran->Rndm();
 647     if (random < 0.5) {
 648         ip=511;
 649     } else if (random < 0.75) {
 650         ip=521;
 651     } else if (random < 0.90) {
 652         ip=531;
 653     } else {
 654         ip=5122;
 655     }
 656     if (ran->Rndm() < 0.5) {ip=-ip;}
 657
 658     return ip;
 659 }
 660
 661 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
 662 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
 663 {
 664 // Return pointer to pT parameterisation
 665     TString sname = TString(tname);
 666     GenFunc func;
 667     switch (param)
 668     {
 669     case kPhi:
 670         func=PtPhi;
 671         break;
 672     case kOmega:
 673         func=PtOmega;
 674         break;
 675     case kEta:
 676         func=PtEta;
 677         break;
 678     case kJpsi:
 679         if (sname == "Vogt" || sname == "Vogt PbPb") {
 680             func=PtJpsiPbPb;
 681         } else if (sname == "Vogt pp") {
 682             func=PtJpsiPP;
 683         } else {
 684             func=PtJpsi;
 685         }
 686         break;
 687     case kJpsiFromB:
 688         func = PtJpsiBPbPb;
 689         break;
 690     case kUpsilon:
 691         if (sname == "Vogt" || sname == "Vogt PbPb") {
 692             func=PtUpsilonPbPb;
 693         } else if (sname == "Vogt pp") {
 694             func=PtUpsilonPP;
 695         } else {
 696             func=PtUpsilon;
 697         }
 698         break;
 699     case kCharm:
 700         func=PtCharm;
 701         break;
 702     case kBeauty:
 703         func=PtBeauty;
 704         break;
 705     case kPion:
 706         func=PtPion;
 707         break;
 708     case kKaon:
 709         func=PtKaon;
 710         break;
 711     default:
 712         func=0;
 713         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
 714     }
 715     return func;
 716 }
 717
 718 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
 719 {
 720     TString sname = TString(tname);
 721
 722 // Return pointer to y- parameterisation
 723     GenFunc func;
 724     switch (param)
 725     {
 726     case kPhi:
 727         func=YPhi;
 728         break;
 729     case kEta:
 730         func=YEta;
 731         break;
 732     case kOmega:
 733         func=YOmega;
 734         break;
 735     case kJpsi:
 736         if (sname == "Vogt" || sname == "Vogt PbPb") {
 737             func=YJpsiPbPb;
 738         } else if (sname == "Vogt pp"){
 739             func=YJpsiPP;
 740         } else {
 741             func=YJpsi;
 742         }
 743
 744         break;
 745     case kJpsiFromB:
 746         func = YJpsiBPbPb;
 747         break;
 748     case kUpsilon:
 749         if (sname == "Vogt" || sname == "Vogt PbPb") {
 750             func=YUpsilonPbPb;
 751         } else if (sname == "Vogt pp") {
 752             func = YUpsilonPP;
 753         } else {
 754             func=YUpsilon;
 755         }
 756         break;
 757     case kCharm:
 758         func=YCharm;
 759         break;
 760     case kBeauty:
 761         func=YBeauty;
 762         break;
 763     case kPion:
 764         func=YPion;
 765         break;
 766     case kKaon:
 767         func=YKaon;
 768         break;
 769     default:
 770         func=0;
 771         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
 772     }
 773     return func;
 774 }
 775 typedef Int_t (*GenFuncIp) (TRandom *);
 776 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* /*tname*/) const
 777 {
 778 // Return pointer to particle type parameterisation
 779     GenFuncIp func;
 780     switch (param)
 781     {
 782     case kPhi:
 783         func=IpPhi;
 784         break;
 785     case kEta:
 786         func=IpEta;
 787         break;
 788     case kOmega:
 789         func=IpOmega;
 790         break;
 791     case kJpsi:
 792     case kJpsiFromB:
 793         func=IpJpsi;
 794         break;
 795     case kUpsilon:
 796         func=IpUpsilon;
 797         break;
 798     case kCharm:
 799         func=IpCharm;
 800         break;
 801     case kBeauty:
 802         func=IpBeauty;
 803         break;
 804     case kPion:
 805         func=IpPion;
 806         break;
 807     case kKaon:
 808         func=IpKaon;
 809         break;
 810     default:
 811         func=0;
 812         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
 813     }
 814     return func;
 815 }
 816
 817
 818
 819 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0,
 820                                    Float_t dx,
 821                                    Int_t n, Int_t no)
 822 {
 823 //
 824 // Neville's alorithm for interpolation
 825 //
 826 // x:  x-value
 827 // y:  Input array
 828 // x0: minimum x
 829 // dx: step size
 830 //  n: number of data points
 831 // no: order of polynom
 832 //
 833     Float_t*  c = new Float_t[n];
 834     Float_t*  d = new Float_t[n];
 835     Int_t m, i;
 836     for (i = 0; i < n; i++) {
 837         c[i] = y[i];
 838         d[i] = y[i];
 839     }
 840
 841     Int_t   ns  = int((x - x0)/dx);
 842
 843     Float_t y1  = y[ns];
 844     ns--;
 845     for (m = 0; m < no; m++) {
 846         for (i = 0; i < n-m; i++) {
 847             Float_t ho = x0 + Float_t(i) * dx - x;
 848             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
 849             Float_t w  = c[i+1] - d[i];
 850             Float_t den = ho-hp;
 851             den = w/den;
 852             d[i] = hp * den;
 853             c[i] = ho * den;
 854         }
 855         Float_t dy;
 856
 857         if (2*ns < (n-m-1)) {
 858             dy  = c[ns+1];
 859         } else {
 860             dy  = d[ns--];
 861         }
 862         y1 += dy;}
 863     delete[] c;
 864     delete[] d;
 865
 866     return y1;
 867 }
 868
 869