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