EVGEN/AliGenMUONlib.cxx

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