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.17  2002/11/07 09:06:10  morsch
  19 J/Psi and Upsilon pt and y distributions from R. Vogt 2002 added.
  20
  21 Revision 1.16  2002/10/14 14:55:35  hristov
  22 Merging the VirtualMC branch to the main development branch (HEAD)
  23
  24 Revision 1.14.6.1  2002/06/10 14:57:41  hristov
  25 Merged with v3-08-02
  26
  27 Revision 1.15  2002/04/17 10:11:51  morsch
  28 Coding Rule violations corrected.
  29
  30 Revision 1.14  2002/02/22 17:26:43  morsch
  31 Eta and omega added.
  32
  33 Revision 1.13  2001/03/27 11:01:04  morsch
  34 Charm pt-distribution corrected. More realistic y-distribution for pi and K.
  35
  36 Revision 1.12  2001/03/09 13:01:41  morsch
  37 - enum constants for paramterisation type (particle family) moved to AliGen*lib.h
  38 - use AliGenGSIlib::kUpsilon, AliGenPHOSlib::kEtaPrime to access the constants
  39
  40 Revision 1.11  2000/11/30 07:12:50  alibrary
  41 Introducing new Rndm and QA classes
  42
  43 Revision 1.10  2000/06/29 21:08:27  morsch
  44 All paramatrisation libraries derive from the pure virtual base class AliGenLib.
  45 This allows to pass a pointer to a library directly to AliGenParam and avoids the
  46 use of function pointers in Config.C.
  47
  48 Revision 1.9  2000/06/14 15:20:56  morsch
  49 Include clean-up (IH)
  50
  51 Revision 1.8  2000/06/09 20:32:11  morsch
  52 All coding rule violations except RS3 corrected
  53
  54 Revision 1.7  2000/05/02 08:12:13  morsch
  55 Coding rule violations corrected.
  56
  57 Revision 1.6  1999/09/29 09:24:14  fca
  58 Introduction of the Copyright and cvs Log
  59
  60 */
  61
  62 // Library class for particle pt and y distributions used for
  63 // muon spectrometer simulations.
  64 // To be used with AliGenParam.
  65 // The following particle typed can be simulated:
  66 // pi, K, phi, omega, eta, J/Psi, Upsilon, charm and beauty mesons.
  67 //
  68 // andreas.morsch@cern.ch
  69 //
  70
  71 #include "TMath.h"
  72 #include "TRandom.h"
  73
  74 #include "AliGenMUONlib.h"
  75
  76 ClassImp(AliGenMUONlib)
  77 //
  78 //  Pions
  79 Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *dummy)
  80 {
  81 //
  82 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
  83 //     POWER LAW FOR PT > 500 MEV
  84 //     MT SCALING BELOW (T=160 MEV)
  85 //
  86   const Double_t kp0 = 1.3;
  87   const Double_t kxn = 8.28;
  88   const Double_t kxlim=0.5;
  89   const Double_t kt=0.160;
  90   const Double_t kxmpi=0.139;
  91   const Double_t kb=1.;
  92   Double_t y, y1, xmpi2, ynorm, a;
  93   Double_t x=*px;
  94   //
  95   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
  96   xmpi2=kxmpi*kxmpi;
  97   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
  98   a=ynorm/y1;
  99   if (x > kxlim)
 100     y=a*TMath::Power(kp0/(kp0+x),kxn);
 101   else
 102     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
 103   return y*x;
 104 }
 105 //
 106 // y-distribution
 107 //
 108 Double_t AliGenMUONlib::YPion( Double_t *py, Double_t *dummy)
 109 {
 110 // Pion y
 111   Double_t y=TMath::Abs(*py);
 112 /*
 113   const Double_t ka    = 7000.;
 114   const Double_t kdy   = 4.;
 115   Double_t ex = y*y/(2*kdy*kdy);
 116   return ka*TMath::Exp(-ex);
 117 */
 118   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
 119
 120 }
 121 //                 particle composition
 122 //
 123 Int_t AliGenMUONlib::IpPion(TRandom *ran)
 124 {
 125 // Pion composition
 126     if (ran->Rndm() < 0.5) {
 127         return  211;
 128     } else {
 129         return -211;
 130     }
 131 }
 132
 133 //____________________________________________________________
 134 //
 135 // Mt-scaling
 136
 137 Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
 138 {
 139   //    SCALING EN MASSE PAR RAPPORT A PTPI
 140   //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
 141   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
 142   //     VALUE MESON/PI AT 5 GEV
 143   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
 144   np--;
 145   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
 146   Double_t fmax2=f5/kfmax[np];
 147   // PIONS
 148   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
 149   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
 150                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
 151   return fmtscal*ptpion;
 152 }
 153 //
 154 // kaon
 155 //
 156 //                pt-distribution
 157 //____________________________________________________________
 158 Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *dummy)
 159 {
 160 // Kaon pT
 161   return PtScal(*px,2);
 162 }
 163
 164 // y-distribution
 165 //____________________________________________________________
 166 Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t *dummy)
 167 {
 168 // Kaon y
 169   Double_t y=TMath::Abs(*py);
 170 /*
 171   const Double_t ka    = 1000.;
 172   const Double_t kdy   = 4.;
 173   //
 174   Double_t ex = y*y/(2*kdy*kdy);
 175   return ka*TMath::Exp(-ex);
 176 */
 177
 178   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
 179 }
 180
 181 //                 particle composition
 182 //
 183 Int_t AliGenMUONlib::IpKaon(TRandom *ran)
 184 {
 185 // Kaon composition
 186     if (ran->Rndm() < 0.5) {
 187         return  321;
 188     } else {
 189         return -321;
 190     }
 191 }
 192
 193 //                    J/Psi
 194 //
 195 //
 196 //                pt-distribution
 197 //____________________________________________________________
 198 Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *dummy)
 199 {
 200 // J/Psi pT
 201   const Double_t kpt0 = 4.;
 202   const Double_t kxn  = 3.6;
 203   Double_t x=*px;
 204   //
 205   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 206   return x/TMath::Power(pass1,kxn);
 207 }
 208
 209 Double_t AliGenMUONlib::PtJpsiPbPb( Double_t *px, Double_t *dummy)
 210 {
 211 // J/Psi pT
 212 //
 213 // R. Vogt 2002
 214 // PbPb 5.5 TeV
 215 // MRST HO
 216 // mc = 1.4 GeV, pt-kick 1 GeV
 217 //
 218
 219     Float_t ptJpsi[100] = {
 220         0.0000e-01,  4.5870e+01,  6.5200e+01,  7.1740e+01,  6.5090e+01,
 221         5.5070e+01,  4.9420e+01,  3.9780e+01,  3.2390e+01,  2.8120e+01,
 222         2.3870e+01,  1.9540e+01,  1.6510e+01,  1.4180e+01,  1.2050e+01,
 223         1.0390e+01,  8.7970e+00,  7.8680e+00,  6.7710e+00,  5.9360e+00,
 224         5.3460e+00,  4.5670e+00,  4.6500e+00,  3.9360e+00,  3.5070e+00,
 225         3.2070e+00,  2.8310e+00,  2.6340e+00,  2.4900e+00,  2.2410e+00,
 226         2.1090e+00,  1.9070e+00,  1.7360e+00,  1.6120e+00,  1.5450e+00,
 227         1.4350e+00,  1.3890e+00,  1.2610e+00,  1.0880e+00,  1.0930e+00,
 228         1.0680e+00,  9.2500e-01,  8.6790e-01,  8.1790e-01,  7.9770e-01,
 229         7.4660e-01,  7.3110e-01,  6.5120e-01,  6.8140e-01,  5.7960e-01,
 230         5.8210e-01,  5.4640e-01,  5.1700e-01,  5.0760e-01,  4.8280e-01,
 231         4.5360e-01,  4.4910e-01,  4.2410e-01,  4.2100e-01,  3.9530e-01,
 232         3.7220e-01,  3.4840e-01,  3.4550e-01,  3.3000e-01,  3.1670e-01,
 233         3.1470e-01,  2.8920e-01,  2.7650e-01,  2.6860e-01,  2.5390e-01,
 234         2.4190e-01,  2.5200e-01,  2.2960e-01,  2.2540e-01,  2.0950e-01,
 235         2.0250e-01,  1.8720e-01,  1.8200e-01,  1.7860e-01,  1.8290e-01,
 236         1.6970e-01,  1.7130e-01,  1.6310e-01,  1.5500e-01,  1.5100e-01,
 237         1.5770e-01,  1.4240e-01,  1.4560e-01,  1.3330e-01,  1.4190e-01,
 238         1.2010e-01,  1.2430e-01,  1.2430e-01,  1.1340e-01,  1.1840e-01,
 239         1.1380e-01,  1.0330e-01,  1.0130e-01,  1.0390e-01,  9.5810e-02
 240     };
 241     Float_t x = px[0] * px[0];
 242     if (x  < 1.5 || x > 100) {
 243         return 0.0;
 244     } else {
 245         Float_t y =  Interpolate(x, ptJpsi, 0.5, 1., 100, 4);
 246         return px[0] * y;
 247     }
 248 }
 249 //
 250 //               y-distribution
 251 //____________________________________________________________
 252 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *dummy)
 253 {
 254 // J/psi y
 255   const Double_t ky0 = 4.;
 256   const Double_t kb=1.;
 257   Double_t yj;
 258   Double_t y=TMath::Abs(*py);
 259   //
 260   if (y < ky0)
 261     yj=kb;
 262   else
 263     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 264   return yj;
 265 }
 266
 267
 268 Double_t AliGenMUONlib::YJpsiPbPb( Double_t *px, Double_t *dummy)
 269 {
 270
 271 //
 272 // J/Psi y
 273 //
 274 //
 275 // R. Vogt 2002
 276 // PbPb 5.5 TeV
 277 // MRST HO
 278 // mc = 1.4 GeV, pt-kick 1 GeV
 279 //
 280
 281     Float_t yJpsi[62] =
 282         {
 283           0.3981E-03, 0.1169E-01, 0.6143E-01, 0.3554E+00, 0.1249E+01
 284         , 0.1677E+01, 0.3634E+01, 0.5414E+01, 0.8242E+01, 0.1102E+02
 285         , 0.1353E+02, 0.1964E+02, 0.2357E+02, 0.2662E+02, 0.3023E+02
 286         , 0.3250E+02, 0.3137E+02, 0.3243E+02, 0.3120E+02, 0.3249E+02
 287         , 0.3166E+02, 0.3104E+02, 0.3203E+02, 0.3149E+02, 0.3117E+02
 288         , 0.3210E+02, 0.3170E+02, 0.3279E+02, 0.3079E+02, 0.3208E+02
 289         , 0.3218E+02, 0.3218E+02, 0.3208E+02, 0.3079E+02, 0.3279E+02
 290         , 0.3170E+02, 0.3210E+02, 0.3118E+02, 0.3149E+02, 0.3203E+02
 291         , 0.3104E+02, 0.3167E+02, 0.3250E+02, 0.3120E+02, 0.3243E+02
 292         , 0.3137E+02, 0.3250E+02, 0.3022E+02, 0.2662E+02, 0.2357E+02
 293         , 0.1964E+02, 0.1353E+02, 0.1102E+02, 0.8242E+01, 0.5414E+01
 294         , 0.3634E+01, 0.1677E+01, 0.1249E+01, 0.3554E+00, 0.6142E-01
 295         , 0.1169E-01, 0.3981E-03};
 296
 297     return Interpolate(px[0], yJpsi, -7.625, 0.25, 62, 2);
 298 }
 299
 300 //                 particle composition
 301 //
 302 Int_t AliGenMUONlib::IpJpsi(TRandom *)
 303 {
 304 // J/Psi composition
 305     return 443;
 306 }
 307
 308 //                      Upsilon
 309 //
 310 //
 311 //                  pt-distribution
 312 //____________________________________________________________
 313 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t *dummy )
 314 {
 315 // Upsilon pT
 316   const Double_t kpt0 = 5.3;
 317   const Double_t kxn  = 2.5;
 318   Double_t x=*px;
 319   //
 320   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 321   return x/TMath::Power(pass1,kxn);
 322 }
 323
 324 Double_t AliGenMUONlib::PtUpsilonPbPb( Double_t *px, Double_t *dummy)
 325 {
 326
 327 //
 328 // Upsilon pT
 329 //
 330 //
 331 // R. Vogt 2002
 332 // PbPb 5.5 TeV
 333 // MRST HO
 334 // mc = 1.4 GeV, pt-kick 1 GeV
 335 //
 336
 337     Float_t ptUps[100] =
 338         {
 339             0.0000e-01, -1.5290e-02,  2.6020e-01,  3.4220e-01,  3.3710e-01,
 340             3.1880e-01,  3.3420e-01,  2.7740e-01,  2.3730e-01,  2.0640e-01,
 341             1.7690e-01,  1.6190e-01,  1.4500e-01,  1.3310e-01,  1.1440e-01,
 342             1.0800e-01,  1.0210e-01,  8.4690e-02,  8.0050e-02,  7.0710e-02,
 343             6.4160e-02,  6.5200e-02,  6.6890e-02,  6.0600e-02,  5.4030e-02,
 344             5.1140e-02,  4.6120e-02,  4.4800e-02,  4.2490e-02,  4.1440e-02,
 345             4.0310e-02,  3.7110e-02,  3.5890e-02,  3.5420e-02,  3.0370e-02,
 346             2.9970e-02,  3.0770e-02,  2.6380e-02,  2.7740e-02,  2.6690e-02,
 347             2.4210e-02,  2.5200e-02,  2.3760e-02,  2.1370e-02,  2.2290e-02,
 348             2.2700e-02,  2.0110e-02,  1.9320e-02,  1.8830e-02,  1.9910e-02,
 349             1.9740e-02,  1.8460e-02,  1.8240e-02,  1.6740e-02,  1.6140e-02,
 350             1.7340e-02,  1.5950e-02,  1.5430e-02,  1.4780e-02,  1.2750e-02,
 351             1.4370e-02,  1.2810e-02,  1.2900e-02,  1.1070e-02,  1.1830e-02,
 352             1.1150e-02,  1.1260e-02,  1.1610e-02,  1.0700e-02,  1.1600e-02,
 353             1.0390e-02,  1.0280e-02,  1.0180e-02,  1.0030e-02,  9.6050e-03,
 354             8.8050e-03,  8.9680e-03,  9.0120e-03,  8.4110e-03,  8.6660e-03,
 355             8.3060e-03,  8.5850e-03,  8.2600e-03,  8.3800e-03,  8.4200e-03,
 356             7.5690e-03,  7.2100e-03,  7.1230e-03,  7.3350e-03,  7.1980e-03,
 357             6.7500e-03,  6.6190e-03,  6.3370e-03,  6.6270e-03,  6.8290e-03,
 358             6.0880e-03,  6.6310e-03,  6.0490e-03,  5.8900e-03,  5.6100e-03
 359         };
 360     Float_t x = px[0] * px[0];
 361     if (x  < 1.5 || x > 100) {
 362         return 0.0;
 363     } else {
 364         Float_t y =  Interpolate(x, ptUps, 0.5, 1., 100, 4);
 365         return px[0] * y;
 366     }
 367 }
 368
 369 //
 370 //                    y-distribution
 371 //
 372 //____________________________________________________________
 373 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *dummy)
 374 {
 375 // Upsilon y
 376   const Double_t ky0 = 3.;
 377   const Double_t kb=1.;
 378   Double_t yu;
 379   Double_t y=TMath::Abs(*py);
 380   //
 381   if (y < ky0)
 382     yu=kb;
 383   else
 384     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 385   return yu;
 386 }
 387
 388
 389 Double_t AliGenMUONlib::YUpsilonPbPb( Double_t *px, Double_t *dummy)
 390 {
 391
 392 //
 393 // Upsilon y
 394 //
 395 //
 396 // R. Vogt 2002
 397 // PbPb 5.5 TeV
 398 // MRST HO
 399 // mc = 1.4 GeV, pt-kick 1 GeV
 400 //
 401
 402     Float_t yUps[52] =
 403         {
 404             0.000000, 0.000065, 0.000767, 0.004332, 0.014790,
 405             0.029370, 0.052060, 0.077930, 0.122500, 0.135800,
 406             0.184000, 0.207900, 0.228300, 0.258500, 0.269500,
 407             0.288500, 0.316200, 0.304100, 0.315800, 0.323300,
 408             0.322400, 0.322600, 0.345500, 0.338100, 0.331900,
 409             0.343700, 0.343700, 0.331900, 0.338100, 0.345500,
 410             0.322600, 0.322400, 0.323300, 0.315800, 0.304100,
 411             0.316200, 0.288500, 0.269500, 0.258500, 0.228300,
 412             0.207900, 0.184000, 0.135800, 0.122500, 0.077930,
 413             0.052060, 0.029380, 0.014780, 0.004332, 0.000767,
 414             0.6479E-04, 0.1013E-06
 415         };
 416
 417
 418
 419     return Interpolate(px[0], yUps, -6.5, 0.25, 52, 2);
 420 }
 421
 422 //                 particle composition
 423 //
 424 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
 425 {
 426 // y composition
 427     return 553;
 428 }
 429
 430 //
 431 //                        Phi
 432 //
 433 //
 434 //    pt-distribution (by scaling of pion distribution)
 435 //____________________________________________________________
 436 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *dummy)
 437 {
 438 // Phi pT
 439   return PtScal(*px,7);
 440 }
 441 //    y-distribution
 442 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *dummy)
 443 {
 444 // Phi y
 445     Double_t *dum=0;
 446     return YJpsi(px,dum);
 447 }
 448 //                 particle composition
 449 //
 450 Int_t AliGenMUONlib::IpPhi(TRandom *)
 451 {
 452 // Phi composition
 453     return 333;
 454 }
 455
 456 //
 457 //                        omega
 458 //
 459 //
 460 //    pt-distribution (by scaling of pion distribution)
 461 //____________________________________________________________
 462 Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t *dummy)
 463 {
 464 // Omega pT
 465   return PtScal(*px,5);
 466 }
 467 //    y-distribution
 468 Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t *dummy)
 469 {
 470 // Omega y
 471     Double_t *dum=0;
 472     return YJpsi(px,dum);
 473 }
 474 //                 particle composition
 475 //
 476 Int_t AliGenMUONlib::IpOmega(TRandom *)
 477 {
 478 // Omega composition
 479     return 223;
 480 }
 481
 482
 483 //
 484 //                        Eta
 485 //
 486 //
 487 //    pt-distribution (by scaling of pion distribution)
 488 //____________________________________________________________
 489 Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t *dummy)
 490 {
 491 // Eta pT
 492   return PtScal(*px,3);
 493 }
 494 //    y-distribution
 495 Double_t AliGenMUONlib::YEta( Double_t *px, Double_t *dummy)
 496 {
 497 // Eta y
 498     Double_t *dum=0;
 499     return YJpsi(px,dum);
 500 }
 501 //                 particle composition
 502 //
 503 Int_t AliGenMUONlib::IpEta(TRandom *)
 504 {
 505 // Eta composition
 506     return 221;
 507 }
 508
 509 //
 510 //                        Charm
 511 //
 512 //
 513 //                    pt-distribution
 514 //____________________________________________________________
 515 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *dummy)
 516 {
 517 // Charm pT
 518   const Double_t kpt0 = 4.08;
 519   const Double_t kxn  = 9.40;
 520
 521   Double_t x=*px;
 522   //
 523   Double_t pass1 = 1.+(x/kpt0);
 524   return x/TMath::Power(pass1,kxn);
 525 }
 526 //                  y-distribution
 527 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *dummy)
 528 {
 529 // Charm y
 530     Double_t *dum=0;
 531     return YJpsi(px,dum);
 532 }
 533
 534 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
 535 {
 536 // Charm composition
 537     Float_t random;
 538     Int_t ip;
 539 //    411,421,431,4122
 540     random = ran->Rndm();
 541     if (random < 0.5) {
 542         ip=411;
 543     } else if (random < 0.75) {
 544         ip=421;
 545     } else if (random < 0.90) {
 546         ip=431;
 547     } else {
 548         ip=4122;
 549     }
 550     if (ran->Rndm() < 0.5) {ip=-ip;}
 551
 552     return ip;
 553 }
 554
 555
 556 //
 557 //                        Beauty
 558 //
 559 //
 560 //                    pt-distribution
 561 //____________________________________________________________
 562 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *dummy)
 563 {
 564 // Beauty pT
 565   const Double_t kpt0 = 4.;
 566   const Double_t kxn  = 3.6;
 567   Double_t x=*px;
 568   //
 569   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 570   return x/TMath::Power(pass1,kxn);
 571 }
 572 //                     y-distribution
 573 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *dummy)
 574 {
 575 // Beauty y
 576     Double_t *dum=0;
 577     return YJpsi(px,dum);
 578 }
 579
 580 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
 581 {
 582 // Beauty Composition
 583     Float_t random;
 584     Int_t ip;
 585     random = ran->Rndm();
 586     if (random < 0.5) {
 587         ip=511;
 588     } else if (random < 0.75) {
 589         ip=521;
 590     } else if (random < 0.90) {
 591         ip=531;
 592     } else {
 593         ip=5122;
 594     }
 595     if (ran->Rndm() < 0.5) {ip=-ip;}
 596
 597     return ip;
 598 }
 599
 600 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
 601 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
 602 {
 603 // Return pointer to pT parameterisation
 604     TString sname = TString(tname);
 605     GenFunc func;
 606     switch (param)
 607     {
 608     case kPhi:
 609         func=PtPhi;
 610         break;
 611     case kOmega:
 612         func=PtOmega;
 613         break;
 614     case kEta:
 615         func=PtEta;
 616         break;
 617     case kJpsi:
 618         if (sname == "Vogt") {
 619             func=PtJpsiPbPb;
 620         } else {
 621             func=PtJpsi;
 622         }
 623         break;
 624     case kUpsilon:
 625         if (sname == "Vogt") {
 626             func=PtUpsilonPbPb;
 627         } else {
 628             func=PtUpsilon;
 629         }
 630         break;
 631     case kCharm:
 632         func=PtCharm;
 633         break;
 634     case kBeauty:
 635         func=PtBeauty;
 636         break;
 637     case kPion:
 638         func=PtPion;
 639         break;
 640     case kKaon:
 641         func=PtKaon;
 642         break;
 643     default:
 644         func=0;
 645         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
 646     }
 647     return func;
 648 }
 649
 650 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
 651 {
 652     TString sname = TString(tname);
 653
 654 // Return pointer to y- parameterisation
 655     GenFunc func;
 656     switch (param)
 657     {
 658     case kPhi:
 659         func=YPhi;
 660         break;
 661     case kEta:
 662         func=YEta;
 663         break;
 664     case kOmega:
 665         func=YOmega;
 666         break;
 667     case kJpsi:
 668         if (sname == "Vogt") {
 669             func=YJpsiPbPb;
 670         } else {
 671             func=YJpsi;
 672         }
 673         break;
 674     case kUpsilon:
 675         if (sname == "Vogt") {
 676             func=YUpsilonPbPb;
 677         } else {
 678             func=YUpsilon;
 679         }
 680         break;
 681     case kCharm:
 682         func=YCharm;
 683         break;
 684     case kBeauty:
 685         func=YBeauty;
 686         break;
 687     case kPion:
 688         func=YPion;
 689         break;
 690     case kKaon:
 691         func=YKaon;
 692         break;
 693     default:
 694         func=0;
 695         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
 696     }
 697     return func;
 698 }
 699 typedef Int_t (*GenFuncIp) (TRandom *);
 700 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* tname) const
 701 {
 702 // Return pointer to particle type parameterisation
 703     GenFuncIp func;
 704     switch (param)
 705     {
 706     case kPhi:
 707         func=IpPhi;
 708         break;
 709     case kEta:
 710         func=IpEta;
 711         break;
 712     case kOmega:
 713         func=IpOmega;
 714         break;
 715     case kJpsi:
 716         func=IpJpsi;
 717         break;
 718     case kUpsilon:
 719         func=IpUpsilon;
 720         break;
 721     case kCharm:
 722         func=IpCharm;
 723         break;
 724     case kBeauty:
 725         func=IpBeauty;
 726         break;
 727     case kPion:
 728         func=IpPion;
 729         break;
 730     case kKaon:
 731         func=IpKaon;
 732         break;
 733     default:
 734         func=0;
 735         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
 736     }
 737     return func;
 738 }
 739
 740
 741
 742 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0,
 743                                    Float_t dx,
 744                                    Int_t n, Int_t no)
 745 {
 746 //
 747 // Neville's alorithm for interpolation
 748 //
 749 // x:  x-value
 750 // y:  Input array
 751 // x0: minimum x
 752 // dx: step size
 753 //  n: number of data points
 754 // no: order of polynom
 755 //
 756     Float_t*  c = new Float_t[n];
 757     Float_t*  d = new Float_t[n];
 758     Int_t m, i;
 759     for (i = 0; i < n; i++) {
 760         c[i] = y[i];
 761         d[i] = y[i];
 762     }
 763
 764     Int_t   ns  = int((x - x0)/dx);
 765
 766     Float_t y1  = y[ns];
 767     ns--;
 768     for (m = 0; m < no; m++) {
 769         for (i = 0; i < n-m; i++) {
 770             Float_t ho = x0 + Float_t(i) * dx - x;
 771             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
 772             Float_t w  = c[i+1] - d[i];
 773             Float_t den = ho-hp;
 774             den = w/den;
 775             d[i] = hp * den;
 776             c[i] = ho * den;
 777         }
 778         Float_t dy;
 779
 780         if (2*ns < (n-m-1)) {
 781             dy  = c[ns+1];
 782         } else {
 783             dy  = d[ns--];
 784         }
 785         y1 += dy;}
 786     delete[] c;
 787     delete[] d;
 788
 789     return y1;
 790 }
 791
 792