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