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 Int_t AliGenMUONlib::IpPsiP(TRandom *)
 343 {
 344 // Psi prime composition
 345     return 100443;
 346 }
 347 Int_t AliGenMUONlib::IpJpsiFamily(TRandom *)
 348 {
 349 // J/Psi composition
 350   Int_t ip;
 351   Float_t r = gRandom->Rndm();
 352   if (r < 0.98) {
 353     ip = 443;
 354   } else {
 355     ip = 100443;
 356   }
 357   return ip;
 358 }
 359
 360
 361
 362 //                      Upsilon
 363 //
 364 //
 365 //                  pt-distribution
 366 //____________________________________________________________
 367 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t */*dummy*/ )
 368 {
 369 // Upsilon pT
 370   const Double_t kpt0 = 5.3;
 371   const Double_t kxn  = 2.5;
 372   Double_t x=*px;
 373   //
 374   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 375   return x/TMath::Power(pass1,kxn);
 376 }
 377
 378 Double_t AliGenMUONlib::PtUpsilonPbPb( Double_t *px, Double_t */*dummy*/)
 379 {
 380
 381 //
 382 // Upsilon pT
 383 //
 384 //
 385 // R. Vogt 2002
 386 // PbPb 5.5 TeV
 387 // MRST HO
 388 // mc = 1.4 GeV, pt-kick 1 GeV
 389 //
 390     Float_t x = px[0];
 391     Double_t c[8] = {
 392         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,
 393         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
 394     };
 395     Double_t y;
 396     if (x < 10.) {
 397         Int_t j;
 398         y = c[j = 7];
 399         while (j > 0) y  = y * x +c[--j];
 400         y = x * TMath::Exp(y);
 401     } else {
 402         y = 0.;
 403     }
 404     return y;
 405 }
 406
 407 Double_t AliGenMUONlib::PtUpsilonPP( Double_t *px, Double_t */*dummy*/)
 408 {
 409
 410 //
 411 // Upsilon pT
 412 //
 413 //
 414 // R. Vogt 2002
 415 // pp 14 TeV
 416 // MRST HO
 417 // mc = 1.4 GeV, pt-kick 1 GeV
 418 //
 419     Float_t x = px[0];
 420     Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,
 421                      -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
 422
 423     Double_t y;
 424     if (x < 10.) {
 425         Int_t j;
 426         y = c[j = 7];
 427         while (j > 0) y  = y * x +c[--j];
 428         y = x * TMath::Exp(y);
 429     } else {
 430         y = 0.;
 431     }
 432     return y;
 433 }
 434
 435 //
 436 //                    y-distribution
 437 //
 438 //____________________________________________________________
 439 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t */*dummy*/)
 440 {
 441 // Upsilon y
 442   const Double_t ky0 = 3.;
 443   const Double_t kb=1.;
 444   Double_t yu;
 445   Double_t y=TMath::Abs(*py);
 446   //
 447   if (y < ky0)
 448     yu=kb;
 449   else
 450     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 451   return yu;
 452 }
 453
 454
 455 Double_t AliGenMUONlib::YUpsilonPbPb( Double_t *px, Double_t */*dummy*/)
 456 {
 457
 458 //
 459 // Upsilon y
 460 //
 461 //
 462 // R. Vogt 2002
 463 // PbPb 5.5 TeV
 464 // MRST HO
 465 // mc = 1.4 GeV, pt-kick 1 GeV
 466 //
 467
 468     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
 469                      -2.99753e-09, 1.28895e-05};
 470
 471     Double_t x = px[0];
 472     if (TMath::Abs(x) > 5.55) return 0.;
 473     Int_t j;
 474     Double_t y = c[j = 6];
 475     while (j > 0) y  = y * x +c[--j];
 476     return y;
 477 }
 478
 479 Double_t AliGenMUONlib::YUpsilonPP( Double_t *px, Double_t */*dummy*/)
 480 {
 481
 482 //
 483 // Upsilon y
 484 //
 485 //
 486 // R. Vogt 2002
 487 // p p  14. TeV
 488 // MRST HO
 489 // mc = 1.4 GeV, pt-kick 1 GeV
 490 //
 491     Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04,
 492                      -6.20539e-10, 1.29943e-05};
 493
 494     Double_t x = px[0];
 495     if (TMath::Abs(x) > 6.2) return 0.;
 496     Int_t j;
 497     Double_t y = c[j = 6];
 498     while (j > 0) y  = y * x +c[--j];
 499     return y;
 500 }
 501
 502 //                 particle composition
 503 //
 504 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
 505 {
 506 // y composition
 507     return 553;
 508 }
 509 Int_t AliGenMUONlib::IpUpsilonP(TRandom *)
 510 {
 511 // y composition
 512     return 100553;
 513 }
 514 Int_t AliGenMUONlib::IpUpsilonPP(TRandom *)
 515 {
 516 // y composition
 517     return 200553;
 518 }
 519 Int_t AliGenMUONlib::IpUpsilonFamily(TRandom *)
 520 {
 521 // y composition
 522   Int_t ip;
 523   Float_t r = gRandom->Rndm();
 524
 525   if (r < 0.712) {
 526     ip = 553;
 527   } else if (r < 0.896) {
 528     ip = 100553;
 529   } else {
 530     ip = 200553;
 531   }
 532   return ip;
 533 }
 534
 535
 536 //
 537 //                        Phi
 538 //
 539 //
 540 //    pt-distribution (by scaling of pion distribution)
 541 //____________________________________________________________
 542 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t */*dummy*/)
 543 {
 544 // Phi pT
 545   return PtScal(*px,7);
 546 }
 547 //    y-distribution
 548 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t */*dummy*/)
 549 {
 550 // Phi y
 551     Double_t *dum=0;
 552     return YJpsi(px,dum);
 553 }
 554 //                 particle composition
 555 //
 556 Int_t AliGenMUONlib::IpPhi(TRandom *)
 557 {
 558 // Phi composition
 559     return 333;
 560 }
 561
 562 //
 563 //                        omega
 564 //
 565 //
 566 //    pt-distribution (by scaling of pion distribution)
 567 //____________________________________________________________
 568 Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t */*dummy*/)
 569 {
 570 // Omega pT
 571   return PtScal(*px,5);
 572 }
 573 //    y-distribution
 574 Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t */*dummy*/)
 575 {
 576 // Omega y
 577     Double_t *dum=0;
 578     return YJpsi(px,dum);
 579 }
 580 //                 particle composition
 581 //
 582 Int_t AliGenMUONlib::IpOmega(TRandom *)
 583 {
 584 // Omega composition
 585     return 223;
 586 }
 587
 588
 589 //
 590 //                        Eta
 591 //
 592 //
 593 //    pt-distribution (by scaling of pion distribution)
 594 //____________________________________________________________
 595 Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t */*dummy*/)
 596 {
 597 // Eta pT
 598   return PtScal(*px,3);
 599 }
 600 //    y-distribution
 601 Double_t AliGenMUONlib::YEta( Double_t *px, Double_t */*dummy*/)
 602 {
 603 // Eta y
 604     Double_t *dum=0;
 605     return YJpsi(px,dum);
 606 }
 607 //                 particle composition
 608 //
 609 Int_t AliGenMUONlib::IpEta(TRandom *)
 610 {
 611 // Eta composition
 612     return 221;
 613 }
 614
 615 //
 616 //                        Charm
 617 //
 618 //
 619 //                    pt-distribution
 620 //____________________________________________________________
 621 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t */*dummy*/)
 622 {
 623 // Charm pT
 624   const Double_t kpt0 = 4.08;
 625   const Double_t kxn  = 9.40;
 626
 627   Double_t x=*px;
 628   //
 629   Double_t pass1 = 1.+(x/kpt0);
 630   return x/TMath::Power(pass1,kxn);
 631 }
 632 //                  y-distribution
 633 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t */*dummy*/)
 634 {
 635 // Charm y
 636     Double_t *dum=0;
 637     return YJpsi(px,dum);
 638 }
 639
 640 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
 641 {
 642 // Charm composition
 643     Float_t random;
 644     Int_t ip;
 645 //    411,421,431,4122
 646     random = ran->Rndm();
 647     if (random < 0.5) {
 648         ip=411;
 649     } else if (random < 0.75) {
 650         ip=421;
 651     } else if (random < 0.90) {
 652         ip=431;
 653     } else {
 654         ip=4122;
 655     }
 656     if (ran->Rndm() < 0.5) {ip=-ip;}
 657
 658     return ip;
 659 }
 660
 661
 662 //
 663 //                        Beauty
 664 //
 665 //
 666 //                    pt-distribution
 667 //____________________________________________________________
 668 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t */*dummy*/)
 669 {
 670 // Beauty pT
 671   const Double_t kpt0 = 4.;
 672   const Double_t kxn  = 3.6;
 673   Double_t x=*px;
 674   //
 675   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 676   return x/TMath::Power(pass1,kxn);
 677 }
 678 //                     y-distribution
 679 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t */*dummy*/)
 680 {
 681 // Beauty y
 682     Double_t *dum=0;
 683     return YJpsi(px,dum);
 684 }
 685
 686 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
 687 {
 688 // Beauty Composition
 689     Float_t random;
 690     Int_t ip;
 691     random = ran->Rndm();
 692     if (random < 0.5) {
 693         ip=511;
 694     } else if (random < 0.75) {
 695         ip=521;
 696     } else if (random < 0.90) {
 697         ip=531;
 698     } else {
 699         ip=5122;
 700     }
 701     if (ran->Rndm() < 0.5) {ip=-ip;}
 702
 703     return ip;
 704 }
 705
 706 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
 707 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
 708 {
 709 // Return pointer to pT parameterisation
 710     TString sname = TString(tname);
 711     GenFunc func;
 712     switch (param)
 713     {
 714     case kPhi:
 715         func=PtPhi;
 716         break;
 717     case kOmega:
 718         func=PtOmega;
 719         break;
 720     case kEta:
 721         func=PtEta;
 722         break;
 723     case kJpsiFamily:
 724     case kPsiP:
 725     case kJpsi:
 726         if (sname == "Vogt" || sname == "Vogt PbPb") {
 727             func=PtJpsiPbPb;
 728         } else if (sname == "Vogt pp") {
 729             func=PtJpsiPP;
 730         } else {
 731             func=PtJpsi;
 732         }
 733         break;
 734     case kJpsiFromB:
 735         func = PtJpsiBPbPb;
 736         break;
 737     case kUpsilonFamily:
 738     case kUpsilonP:
 739     case kUpsilonPP:
 740     case kUpsilon:
 741         if (sname == "Vogt" || sname == "Vogt PbPb") {
 742             func=PtUpsilonPbPb;
 743         } else if (sname == "Vogt pp") {
 744             func=PtUpsilonPP;
 745         } else {
 746             func=PtUpsilon;
 747         }
 748         break;
 749     case kCharm:
 750         func=PtCharm;
 751         break;
 752     case kBeauty:
 753         func=PtBeauty;
 754         break;
 755     case kPion:
 756         func=PtPion;
 757         break;
 758     case kKaon:
 759         func=PtKaon;
 760         break;
 761     default:
 762         func=0;
 763         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
 764     }
 765     return func;
 766 }
 767
 768 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
 769 {
 770     TString sname = TString(tname);
 771
 772 // Return pointer to y- parameterisation
 773     GenFunc func;
 774     switch (param)
 775     {
 776     case kPhi:
 777         func=YPhi;
 778         break;
 779     case kEta:
 780         func=YEta;
 781         break;
 782     case kOmega:
 783         func=YOmega;
 784         break;
 785     case kJpsiFamily:
 786     case kPsiP:
 787     case kJpsi:
 788         if (sname == "Vogt" || sname == "Vogt PbPb") {
 789             func=YJpsiPbPb;
 790         } else if (sname == "Vogt pp"){
 791             func=YJpsiPP;
 792         } else {
 793             func=YJpsi;
 794         }
 795         break;
 796     case kJpsiFromB:
 797         func = YJpsiBPbPb;
 798         break;
 799     case kUpsilonFamily:
 800     case kUpsilonP:
 801     case kUpsilonPP:
 802     case kUpsilon:
 803         if (sname == "Vogt" || sname == "Vogt PbPb") {
 804             func=YUpsilonPbPb;
 805         } else if (sname == "Vogt pp") {
 806             func = YUpsilonPP;
 807         } else {
 808             func=YUpsilon;
 809         }
 810         break;
 811     case kCharm:
 812         func=YCharm;
 813         break;
 814     case kBeauty:
 815         func=YBeauty;
 816         break;
 817     case kPion:
 818         func=YPion;
 819         break;
 820     case kKaon:
 821         func=YKaon;
 822         break;
 823     default:
 824         func=0;
 825         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
 826     }
 827     return func;
 828 }
 829 typedef Int_t (*GenFuncIp) (TRandom *);
 830 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* /*tname*/) const
 831 {
 832 // Return pointer to particle type parameterisation
 833     GenFuncIp func;
 834     switch (param)
 835     {
 836     case kPhi:
 837         func=IpPhi;
 838         break;
 839     case kEta:
 840         func=IpEta;
 841         break;
 842     case kOmega:
 843         func=IpOmega;
 844         break;
 845     case kJpsiFamily:
 846         func=IpJpsiFamily;
 847         break;
 848     case kPsiP:
 849         func=IpPsiP;
 850         break;
 851     case kJpsi:
 852     case kJpsiFromB:
 853         func=IpJpsi;
 854         break;
 855     case kUpsilon:
 856         func=IpUpsilon;
 857         break;
 858     case kUpsilonFamily:
 859       func=IpUpsilonFamily;
 860       break;
 861     case kUpsilonP:
 862         func=IpUpsilonP;
 863         break;
 864     case kUpsilonPP:
 865         func=IpUpsilonPP;
 866         break;
 867     case kCharm:
 868         func=IpCharm;
 869         break;
 870     case kBeauty:
 871         func=IpBeauty;
 872         break;
 873     case kPion:
 874         func=IpPion;
 875         break;
 876     case kKaon:
 877         func=IpKaon;
 878         break;
 879     default:
 880         func=0;
 881         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
 882     }
 883     return func;
 884 }
 885
 886
 887
 888 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0,
 889                                    Float_t dx,
 890                                    Int_t n, Int_t no)
 891 {
 892 //
 893 // Neville's alorithm for interpolation
 894 //
 895 // x:  x-value
 896 // y:  Input array
 897 // x0: minimum x
 898 // dx: step size
 899 //  n: number of data points
 900 // no: order of polynom
 901 //
 902     Float_t*  c = new Float_t[n];
 903     Float_t*  d = new Float_t[n];
 904     Int_t m, i;
 905     for (i = 0; i < n; i++) {
 906         c[i] = y[i];
 907         d[i] = y[i];
 908     }
 909
 910     Int_t   ns  = int((x - x0)/dx);
 911
 912     Float_t y1  = y[ns];
 913     ns--;
 914     for (m = 0; m < no; m++) {
 915         for (i = 0; i < n-m; i++) {
 916             Float_t ho = x0 + Float_t(i) * dx - x;
 917             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
 918             Float_t w  = c[i+1] - d[i];
 919             Float_t den = ho-hp;
 920             den = w/den;
 921             d[i] = hp * den;
 922             c[i] = ho * den;
 923         }
 924         Float_t dy;
 925
 926         if (2*ns < (n-m-1)) {
 927             dy  = c[ns+1];
 928         } else {
 929             dy  = d[ns--];
 930         }
 931         y1 += dy;}
 932     delete[] c;
 933     delete[] d;
 934
 935     return y1;
 936 }
 937
 938