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::PtJpsiCDFscaledPP( Double_t *px, Double_t */*dummy*/)
 177 {
 178 // J/Psi pT
 179 //
 180 // pp 14 TeV
 181 //
 182 // scaled from CDF data at 2 TeV
 183
 184   const Double_t kpt0 = 5.355;
 185   const Double_t kxn  = 3.821;
 186   Double_t x=*px;
 187   //
 188   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 189   return x/TMath::Power(pass1,kxn);
 190 }
 191
 192 Double_t AliGenMUONlib::PtJpsiFlat( Double_t */*px*/, Double_t */*dummy*/ )
 193 {
 194   return 1.;
 195 }
 196
 197 Double_t AliGenMUONlib::PtJpsiPbPb( Double_t *px, Double_t */*dummy*/)
 198 {
 199 // J/Psi pT spectrum
 200 //
 201 // R. Vogt 2002
 202 // PbPb 5.5 TeV
 203 // MRST HO
 204 // mc = 1.4 GeV, pt-kick 1 GeV
 205 //
 206     Float_t x = px[0];
 207     Float_t c[8] = {
 208         -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00,
 209         -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
 210     };
 211
 212     Double_t y;
 213     if (x < 10.) {
 214         Int_t j;
 215         y = c[j = 7];
 216         while (j > 0) y  = y * x +c[--j];
 217         y = x * TMath::Exp(y);
 218     } else {
 219         y = 0.;
 220     }
 221     return y;
 222 }
 223
 224 Double_t AliGenMUONlib::PtJpsiBPbPb( Double_t *px, Double_t */*dummy*/)
 225 {
 226 // J/Psi pT spectrum
 227 // B -> J/Psi X
 228     Double_t x0 =   4.0384;
 229     Double_t  n =   3.0288;
 230
 231     Double_t x = px[0];
 232     Double_t y = x / TMath::Power((1. + (x/x0)*(x/x0)), n);
 233
 234     return y;
 235 }
 236
 237
 238 Double_t AliGenMUONlib::PtJpsiPP( Double_t *px, Double_t */*dummy*/)
 239 {
 240 // J/Psi pT spectrum
 241 //
 242 // R. Vogt 2002
 243 // pp 14 TeV
 244 // MRST HO
 245 // mc = 1.4 GeV, pt-kick 1 GeV
 246 //
 247     Float_t x = px[0];
 248     Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
 249
 250     Double_t y;
 251     if (x < 10.) {
 252         Int_t j;
 253         y = c[j = 3];
 254         while (j > 0) y  = y * x +c[--j];
 255         y = x * TMath::Exp(y);
 256     } else {
 257         y = 0.;
 258     }
 259     return y;
 260 }
 261
 262 //
 263 //               y-distribution
 264 //____________________________________________________________
 265 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t */*dummy*/)
 266 {
 267 // J/psi y
 268   const Double_t ky0 = 4.;
 269   const Double_t kb=1.;
 270   Double_t yj;
 271   Double_t y=TMath::Abs(*py);
 272   //
 273   if (y < ky0)
 274     yj=kb;
 275   else
 276     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 277   return yj;
 278 }
 279
 280 Double_t AliGenMUONlib::YJpsiFlat( Double_t */*py*/, Double_t */*dummy*/ )
 281 {
 282   return 1.;
 283 }
 284
 285
 286 Double_t AliGenMUONlib::YJpsiPbPb( Double_t *px, Double_t */*dummy*/)
 287 {
 288
 289 //
 290 // J/Psi y
 291 //
 292 //
 293 // R. Vogt 2002
 294 // PbPb 5.5 TeV
 295 // MRST HO
 296 // mc = 1.4 GeV, pt-kick 1 GeV
 297 //
 298     Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
 299     Double_t x = TMath::Abs(px[0]);
 300     Double_t y;
 301
 302     if (x < 4.) {
 303         y = 31.754;
 304     } else if (x < 6) {
 305         Int_t j;
 306         y = c[j = 4];
 307         while (j > 0) y  = y * x + c[--j];
 308     } else {
 309         y =0.;
 310     }
 311
 312     return y;
 313 }
 314
 315 Double_t AliGenMUONlib::YJpsiCDFscaled( Double_t *px, Double_t* dummy)
 316 {
 317     // J/Psi y
 318     return AliGenMUONlib::YJpsiPbPb(px, dummy);
 319 }
 320
 321 Double_t AliGenMUONlib::YJpsiCDFscaledPP( Double_t *px, Double_t* dummy)
 322 {
 323     // J/Psi y
 324     return AliGenMUONlib::YJpsiPP(px, dummy);
 325 }
 326
 327 Double_t AliGenMUONlib::YJpsiPP( Double_t *px, Double_t */*dummy*/)
 328 {
 329
 330 //
 331 // J/Psi y
 332 //
 333 //
 334 // R. Vogt 2002
 335 // pp 14  TeV
 336 // MRST HO
 337 // mc = 1.4 GeV, pt-kick 1 GeV
 338 //
 339
 340     Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
 341     Double_t x = TMath::Abs(px[0]);
 342     Double_t y;
 343
 344     if (x < 2.5) {
 345         y = 96.455 - 0.8483 * x * x;
 346     } else if (x < 7.9) {
 347         Int_t j;
 348         y = c[j = 4];
 349         while (j > 0) y  = y * x + c[--j];
 350     } else {
 351         y =0.;
 352     }
 353
 354     return y;
 355 }
 356
 357 Double_t AliGenMUONlib::YJpsiBPbPb( Double_t *px, Double_t */*dummy*/)
 358 {
 359
 360 //
 361 // J/Psi from B->J/Psi X
 362 //
 363 //
 364
 365
 366     Double_t c[7] = {7.37025e-02, 0., -2.94487e-03, 0., 6.07953e-06, 0., 5.39219e-07};
 367
 368     Double_t x = TMath::Abs(px[0]);
 369     Double_t y;
 370
 371     if (x > 6.) {
 372         y = 0.;
 373     } else {
 374         Int_t j;
 375         y = c[j = 6];
 376         while (j > 0) y  = y * x + c[--j];
 377     }
 378
 379     return y;
 380 }
 381
 382
 383
 384 //                 particle composition
 385 //
 386 Int_t AliGenMUONlib::IpJpsi(TRandom *)
 387 {
 388 // J/Psi composition
 389     return 443;
 390 }
 391 Int_t AliGenMUONlib::IpPsiP(TRandom *)
 392 {
 393 // Psi prime composition
 394     return 100443;
 395 }
 396 Int_t AliGenMUONlib::IpJpsiFamily(TRandom *)
 397 {
 398 // J/Psi composition
 399   Int_t ip;
 400   Float_t r = gRandom->Rndm();
 401   if (r < 0.98) {
 402     ip = 443;
 403   } else {
 404     ip = 100443;
 405   }
 406   return ip;
 407 }
 408
 409
 410
 411 //                      Upsilon
 412 //
 413 //
 414 //                  pt-distribution
 415 //____________________________________________________________
 416 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t */*dummy*/ )
 417 {
 418 // Upsilon pT
 419   const Double_t kpt0 = 5.3;
 420   const Double_t kxn  = 2.5;
 421   Double_t x=*px;
 422   //
 423   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 424   return x/TMath::Power(pass1,kxn);
 425 }
 426
 427 Double_t AliGenMUONlib::PtUpsilonCDFscaled( Double_t *px, Double_t */*dummy*/ )
 428 {
 429 // Upsilon pT
 430   const Double_t kpt0 = 7.753;
 431   const Double_t kxn  = 3.042;
 432   Double_t x=*px;
 433   //
 434   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 435   return x/TMath::Power(pass1,kxn);
 436 }
 437
 438 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP( Double_t *px, Double_t */*dummy*/ )
 439 {
 440 // Upsilon pT
 441 //
 442 // pp 14 TeV
 443 //
 444 // scaled from CDF data at 2 TeV
 445
 446   const Double_t kpt0 = 8.610;
 447   const Double_t kxn  = 3.051;
 448   Double_t x=*px;
 449   //
 450   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 451   return x/TMath::Power(pass1,kxn);
 452 }
 453
 454 Double_t AliGenMUONlib::PtUpsilonFlat( Double_t */*px*/, Double_t */*dummy*/ )
 455 {
 456   return 1.;
 457 }
 458
 459 Double_t AliGenMUONlib::PtUpsilonPbPb( Double_t *px, Double_t */*dummy*/)
 460 {
 461
 462 //
 463 // Upsilon pT
 464 //
 465 //
 466 // R. Vogt 2002
 467 // PbPb 5.5 TeV
 468 // MRST HO
 469 // mc = 1.4 GeV, pt-kick 1 GeV
 470 //
 471     Float_t x = px[0];
 472     Double_t c[8] = {
 473         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,
 474         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
 475     };
 476     Double_t y;
 477     if (x < 10.) {
 478         Int_t j;
 479         y = c[j = 7];
 480         while (j > 0) y  = y * x +c[--j];
 481         y = x * TMath::Exp(y);
 482     } else {
 483         y = 0.;
 484     }
 485     return y;
 486 }
 487
 488 Double_t AliGenMUONlib::PtUpsilonPP( Double_t *px, Double_t */*dummy*/)
 489 {
 490
 491 //
 492 // Upsilon pT
 493 //
 494 //
 495 // R. Vogt 2002
 496 // pp 14 TeV
 497 // MRST HO
 498 // mc = 1.4 GeV, pt-kick 1 GeV
 499 //
 500     Float_t x = px[0];
 501     Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,
 502                      -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
 503
 504     Double_t y;
 505     if (x < 10.) {
 506         Int_t j;
 507         y = c[j = 7];
 508         while (j > 0) y  = y * x +c[--j];
 509         y = x * TMath::Exp(y);
 510     } else {
 511         y = 0.;
 512     }
 513     return y;
 514 }
 515
 516 //
 517 //                    y-distribution
 518 //
 519 //____________________________________________________________
 520 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t */*dummy*/)
 521 {
 522 // Upsilon y
 523   const Double_t ky0 = 3.;
 524   const Double_t kb=1.;
 525   Double_t yu;
 526   Double_t y=TMath::Abs(*py);
 527   //
 528   if (y < ky0)
 529     yu=kb;
 530   else
 531     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 532   return yu;
 533 }
 534
 535
 536 Double_t AliGenMUONlib::YUpsilonPbPb( Double_t *px, Double_t */*dummy*/)
 537 {
 538
 539 //
 540 // Upsilon y
 541 //
 542 //
 543 // R. Vogt 2002
 544 // PbPb 5.5 TeV
 545 // MRST HO
 546 // mc = 1.4 GeV, pt-kick 1 GeV
 547 //
 548
 549     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
 550                      -2.99753e-09, 1.28895e-05};
 551
 552     Double_t x = px[0];
 553     if (TMath::Abs(x) > 5.55) return 0.;
 554     Int_t j;
 555     Double_t y = c[j = 6];
 556     while (j > 0) y  = y * x +c[--j];
 557     return y;
 558 }
 559
 560 Double_t AliGenMUONlib::YUpsilonCDFscaled( Double_t *px, Double_t *dummy)
 561 {
 562     // Upsilon y
 563     return AliGenMUONlib::YUpsilonPbPb(px, dummy);
 564
 565 }
 566
 567 Double_t AliGenMUONlib::YUpsilonCDFscaledPP( Double_t *px, Double_t *dummy)
 568 {
 569     // Upsilon y
 570     return AliGenMUONlib::YUpsilonPP(px, dummy);
 571
 572 }
 573
 574 Double_t AliGenMUONlib::YUpsilonFlat( Double_t */*px*/, Double_t */*dummy*/)
 575 {
 576     // Upsilon y
 577     return 1.;
 578
 579 }
 580
 581 Double_t AliGenMUONlib::YUpsilonPP( Double_t *px, Double_t */*dummy*/)
 582 {
 583
 584 //
 585 // Upsilon y
 586 //
 587 //
 588 // R. Vogt 2002
 589 // p p  14. TeV
 590 // MRST HO
 591 // mc = 1.4 GeV, pt-kick 1 GeV
 592 //
 593     Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04,
 594                      -6.20539e-10, 1.29943e-05};
 595
 596     Double_t x = px[0];
 597     if (TMath::Abs(x) > 6.2) return 0.;
 598     Int_t j;
 599     Double_t y = c[j = 6];
 600     while (j > 0) y  = y * x +c[--j];
 601     return y;
 602 }
 603
 604 //                 particle composition
 605 //
 606 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
 607 {
 608 // y composition
 609     return 553;
 610 }
 611 Int_t AliGenMUONlib::IpUpsilonP(TRandom *)
 612 {
 613 // y composition
 614     return 100553;
 615 }
 616 Int_t AliGenMUONlib::IpUpsilonPP(TRandom *)
 617 {
 618 // y composition
 619     return 200553;
 620 }
 621 Int_t AliGenMUONlib::IpUpsilonFamily(TRandom *)
 622 {
 623 // y composition
 624   Int_t ip;
 625   Float_t r = gRandom->Rndm();
 626
 627   if (r < 0.712) {
 628     ip = 553;
 629   } else if (r < 0.896) {
 630     ip = 100553;
 631   } else {
 632     ip = 200553;
 633   }
 634   return ip;
 635 }
 636
 637
 638 //
 639 //                        Phi
 640 //
 641 //
 642 //    pt-distribution (by scaling of pion distribution)
 643 //____________________________________________________________
 644 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t */*dummy*/)
 645 {
 646 // Phi pT
 647   return PtScal(*px,7);
 648 }
 649 //    y-distribution
 650 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t */*dummy*/)
 651 {
 652 // Phi y
 653     Double_t *dum=0;
 654     return YJpsi(px,dum);
 655 }
 656 //                 particle composition
 657 //
 658 Int_t AliGenMUONlib::IpPhi(TRandom *)
 659 {
 660 // Phi composition
 661     return 333;
 662 }
 663
 664 //
 665 //                        omega
 666 //
 667 //
 668 //    pt-distribution (by scaling of pion distribution)
 669 //____________________________________________________________
 670 Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t */*dummy*/)
 671 {
 672 // Omega pT
 673   return PtScal(*px,5);
 674 }
 675 //    y-distribution
 676 Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t */*dummy*/)
 677 {
 678 // Omega y
 679     Double_t *dum=0;
 680     return YJpsi(px,dum);
 681 }
 682 //                 particle composition
 683 //
 684 Int_t AliGenMUONlib::IpOmega(TRandom *)
 685 {
 686 // Omega composition
 687     return 223;
 688 }
 689
 690
 691 //
 692 //                        Eta
 693 //
 694 //
 695 //    pt-distribution (by scaling of pion distribution)
 696 //____________________________________________________________
 697 Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t */*dummy*/)
 698 {
 699 // Eta pT
 700   return PtScal(*px,3);
 701 }
 702 //    y-distribution
 703 Double_t AliGenMUONlib::YEta( Double_t *px, Double_t */*dummy*/)
 704 {
 705 // Eta y
 706     Double_t *dum=0;
 707     return YJpsi(px,dum);
 708 }
 709 //                 particle composition
 710 //
 711 Int_t AliGenMUONlib::IpEta(TRandom *)
 712 {
 713 // Eta composition
 714     return 221;
 715 }
 716
 717 //
 718 //                        Charm
 719 //
 720 //
 721 //                    pt-distribution
 722 //____________________________________________________________
 723 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t */*dummy*/)
 724 {
 725 // Charm pT
 726   const Double_t kpt0 = 2.25;
 727   const Double_t kxn  = 3.17;
 728
 729   Double_t x=*px;
 730   //
 731   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 732   return x/TMath::Power(pass1,kxn);
 733 }
 734
 735 Double_t AliGenMUONlib::PtCharmCentral( Double_t *px, Double_t */*dummy*/)
 736 {
 737 // Charm pT
 738   const Double_t kpt0 = 2.12;
 739   const Double_t kxn  = 2.78;
 740
 741   Double_t x=*px;
 742   //
 743   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 744   return x/TMath::Power(pass1,kxn);
 745 }
 746 //                  y-distribution
 747 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t */*dummy*/)
 748 {
 749 // Charm y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225)
 750 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
 751 // shadowing + kt broadening
 752
 753     Double_t x=px[0];
 754     Double_t c[2]={-2.42985e-03,-2.31001e-04};
 755     Double_t y=1+(c[0]*TMath::Power(x,2))+(c[1]*TMath::Power(x,4));
 756     Double_t ycharm;
 757
 758     if (TMath::Abs(x)>8) {
 759       ycharm=0.;
 760     }
 761     else {
 762       ycharm=TMath::Power(y,3);
 763     }
 764
 765     return ycharm;
 766 }
 767
 768
 769 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
 770 {
 771 // Charm composition
 772     Float_t random;
 773     Int_t ip;
 774 //    411,421,431,4122
 775     random = ran->Rndm();
 776 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3
 777 //  >>>>> cf. tab 4 p 11
 778
 779     if (random < 0.30) {
 780         ip=421;
 781     } else if (random < 0.60) {
 782         ip=-421;
 783     } else if (random < 0.70) {
 784         ip=411;
 785     } else if (random < 0.80) {
 786         ip=-411;
 787     } else if (random < 0.86) {
 788         ip=431;
 789     } else if (random < 0.92) {
 790         ip=-431;
 791     } else if (random < 0.96) {
 792         ip=4122;
 793     } else {
 794         ip=-4122;
 795     }
 796
 797     return ip;
 798 }
 799
 800 //
 801 //                        Beauty
 802 //
 803 //
 804 //                    pt-distribution
 805 //____________________________________________________________
 806 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t */*dummy*/)
 807 {
 808 // Beauty pT
 809   const Double_t kpt0 = 6.53;
 810   const Double_t kxn  = 3.59;
 811   Double_t x=*px;
 812   //
 813   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 814   return x/TMath::Power(pass1,kxn);
 815 }
 816
 817 Double_t AliGenMUONlib::PtBeautyCentral( Double_t *px, Double_t */*dummy*/)
 818 {
 819 // Beauty pT
 820   const Double_t kpt0 = 6.14;
 821   const Double_t kxn  = 2.93;
 822   Double_t x=*px;
 823   //
 824   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 825   return x/TMath::Power(pass1,kxn);
 826 }
 827 //                     y-distribution
 828 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t */*dummy*/)
 829 {
 830 // Beauty y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225)
 831 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
 832 // shadowing + kt broadening
 833
 834     Double_t x=px[0];
 835     Double_t c[2]={-1.27590e-02,-2.42731e-04};
 836     Double_t y=1+c[0]*TMath::Power(x,2)+c[1]*TMath::Power(x,4);
 837     Double_t ybeauty;
 838
 839     if (TMath::Abs(x)>6) {
 840       ybeauty=0.;
 841     }
 842     else {
 843       ybeauty=TMath::Power(y,3);
 844     }
 845
 846     return ybeauty;
 847 }
 848
 849
 850 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
 851 {
 852 // Beauty Composition
 853     Float_t random;
 854     Int_t ip;
 855     random = ran->Rndm();
 856
 857 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3
 858 //  >>>>> cf. tab 4 p 11
 859
 860  if (random < 0.20) {
 861         ip=511;
 862     } else if (random < 0.40) {
 863         ip=-511;
 864     } else if (random < 0.605) {
 865         ip=521;
 866     } else if (random < 0.81) {
 867         ip=-521;
 868     } else if (random < 0.87) {
 869         ip=531;
 870     } else if (random < 0.93) {
 871         ip=-531;
 872     } else if (random < 0.965) {
 873         ip=5122;
 874     } else {
 875         ip=-5122;
 876     }
 877
 878  return ip;
 879 }
 880
 881
 882 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
 883 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
 884 {
 885 // Return pointer to pT parameterisation
 886     TString sname = TString(tname);
 887     GenFunc func;
 888     switch (param)
 889     {
 890     case kPhi:
 891         func=PtPhi;
 892         break;
 893     case kOmega:
 894         func=PtOmega;
 895         break;
 896     case kEta:
 897         func=PtEta;
 898         break;
 899     case kJpsiFamily:
 900     case kPsiP:
 901     case kJpsi:
 902         if (sname == "Vogt" || sname == "Vogt PbPb") {
 903             func=PtJpsiPbPb;
 904         } else if (sname == "Vogt pp") {
 905             func=PtJpsiPP;
 906         } else if (sname == "CDF scaled") {
 907             func=PtJpsiCDFscaled;
 908         } else if (sname == "CDF pp") {
 909             func=PtJpsiCDFscaledPP;
 910         } else if (sname == "Flat") {
 911             func=PtJpsiFlat;
 912         } else {
 913             func=PtJpsi;
 914         }
 915         break;
 916     case kJpsiFromB:
 917         func = PtJpsiBPbPb;
 918         break;
 919     case kUpsilonFamily:
 920     case kUpsilonP:
 921     case kUpsilonPP:
 922     case kUpsilon:
 923         if (sname == "Vogt" || sname == "Vogt PbPb") {
 924             func=PtUpsilonPbPb;
 925         } else if (sname == "Vogt pp") {
 926             func=PtUpsilonPP;
 927         } else if (sname == "CDF scaled") {
 928             func=PtUpsilonCDFscaled;
 929         } else if (sname == "CDF pp") {
 930             func=PtUpsilonCDFscaledPP;
 931         } else if (sname == "Flat") {
 932             func=PtUpsilonFlat;
 933         } else {
 934             func=PtUpsilon;
 935         }
 936         break;
 937     case kCharm:
 938         if (sname == "central") {
 939             func=PtCharmCentral;
 940         } else {
 941             func=PtCharm;
 942         }
 943         break;
 944     case kBeauty:
 945         if (sname == "central") {
 946             func=PtBeautyCentral;
 947         } else {
 948             func=PtBeauty;
 949         }
 950         break;
 951     case kPion:
 952         func=PtPion;
 953         break;
 954     case kKaon:
 955         func=PtKaon;
 956         break;
 957     case kChi_c0:
 958         func=PtChi_c0;
 959         break;
 960     case kChi_c1:
 961         func=PtChi_c1;
 962         break;
 963     case kChi_c2:
 964         func=PtChi_c2;
 965         break;
 966     case kChi_c:
 967         func=PtChi_c;
 968         break;
 969     default:
 970         func=0;
 971         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
 972     }
 973     return func;
 974 }
 975
 976 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
 977 {
 978   //
 979   // Return pointer to y- parameterisation
 980   //
 981     TString sname = TString(tname);
 982     GenFunc func;
 983     switch (param)
 984     {
 985     case kPhi:
 986         func=YPhi;
 987         break;
 988     case kEta:
 989         func=YEta;
 990         break;
 991     case kOmega:
 992         func=YOmega;
 993         break;
 994     case kJpsiFamily:
 995     case kPsiP:
 996     case kJpsi:
 997         if (sname == "Vogt" || sname == "Vogt PbPb") {
 998             func=YJpsiPbPb;
 999         } else if (sname == "Vogt pp"){
1000             func=YJpsiPP;
1001         } else if (sname == "CDF scaled") {
1002             func=YJpsiCDFscaled;
1003         } else if (sname == "CDF pp") {
1004             func=YJpsiCDFscaledPP;
1005         } else if (sname == "Flat") {
1006             func=YJpsiFlat;
1007         } else {
1008             func=YJpsi;
1009         }
1010         break;
1011     case kJpsiFromB:
1012         func = YJpsiBPbPb;
1013         break;
1014     case kUpsilonFamily:
1015     case kUpsilonP:
1016     case kUpsilonPP:
1017     case kUpsilon:
1018         if (sname == "Vogt" || sname == "Vogt PbPb") {
1019             func=YUpsilonPbPb;
1020         } else if (sname == "Vogt pp") {
1021             func = YUpsilonPP;
1022         } else if (sname == "CDF scaled") {
1023             func=YUpsilonCDFscaled;
1024         } else if (sname == "CDF pp") {
1025             func=YUpsilonCDFscaledPP;
1026         } else if (sname == "Flat") {
1027             func=YUpsilonFlat;
1028         } else {
1029             func=YUpsilon;
1030         }
1031         break;
1032     case kCharm:
1033         func=YCharm;
1034         break;
1035     case kBeauty:
1036         func=YBeauty;
1037         break;
1038     case kPion:
1039         func=YPion;
1040         break;
1041     case kKaon:
1042         func=YKaon;
1043         break;
1044     case kChi_c0:
1045         func=YChi_c0;
1046         break;
1047     case kChi_c1:
1048         func=YChi_c1;
1049         break;
1050     case kChi_c2:
1051         func=YChi_c2;
1052         break;
1053     case kChi_c:
1054         func=YChi_c;
1055         break;
1056     default:
1057         func=0;
1058         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
1059     }
1060     return func;
1061 }
1062
1063 //
1064 //                    Chi
1065 //
1066 //
1067 //                pt-distribution
1068 //____________________________________________________________
1069 Double_t AliGenMUONlib::PtChi_c0( Double_t *px, Double_t */*dummy*/)
1070 {
1071 // Chi_c1 pT
1072   const Double_t kpt0 = 4.;
1073   const Double_t kxn  = 3.6;
1074   Double_t x=*px;
1075   //
1076   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1077   return x/TMath::Power(pass1,kxn);
1078 }
1079 Double_t AliGenMUONlib::PtChi_c1( Double_t *px, Double_t */*dummy*/)
1080 {
1081 // Chi_c1 pT
1082   const Double_t kpt0 = 4.;
1083   const Double_t kxn  = 3.6;
1084   Double_t x=*px;
1085   //
1086   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1087   return x/TMath::Power(pass1,kxn);
1088 }
1089 Double_t AliGenMUONlib::PtChi_c2( Double_t *px, Double_t */*dummy*/)
1090 {
1091 // Chi_c2 pT
1092   const Double_t kpt0 = 4.;
1093   const Double_t kxn  = 3.6;
1094   Double_t x=*px;
1095   //
1096   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1097   return x/TMath::Power(pass1,kxn);
1098 }
1099 Double_t AliGenMUONlib::PtChi_c( Double_t *px, Double_t */*dummy*/)
1100 {
1101 // Chi_c family pT
1102   const Double_t kpt0 = 4.;
1103   const Double_t kxn  = 3.6;
1104   Double_t x=*px;
1105   //
1106   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1107   return x/TMath::Power(pass1,kxn);
1108 }
1109
1110 //
1111 //               y-distribution
1112 //____________________________________________________________
1113 Double_t AliGenMUONlib::YChi_c0(Double_t *py, Double_t */*dummy*/)
1114 {
1115 // Chi-1c y
1116   const Double_t ky0 = 4.;
1117  const Double_t kb=1.;
1118   Double_t yj;
1119   Double_t y=TMath::Abs(*py);
1120   //
1121   if (y < ky0)
1122     yj=kb;
1123   else
1124     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
1125   return yj;
1126 }
1127
1128 Double_t AliGenMUONlib::YChi_c1(Double_t *py, Double_t */*dummy*/)
1129 {
1130 // Chi-1c y
1131   const Double_t ky0 = 4.;
1132   const Double_t kb=1.;
1133   Double_t yj;
1134   Double_t y=TMath::Abs(*py);
1135   //
1136   if (y < ky0)
1137     yj=kb;
1138   else
1139     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
1140   return yj;
1141 }
1142
1143 Double_t AliGenMUONlib::YChi_c2(Double_t *py, Double_t */*dummy*/)
1144 {
1145 // Chi-2c y
1146   const Double_t ky0 = 4.;
1147   const Double_t kb=1.;
1148   Double_t yj;
1149   Double_t y=TMath::Abs(*py);
1150   //
1151   if (y < ky0)
1152     yj=kb;
1153   else
1154     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
1155   return yj;
1156 }
1157
1158 Double_t AliGenMUONlib::YChi_c(Double_t *py, Double_t */*dummy*/)
1159 {
1160 // Chi_c family y
1161   const Double_t ky0 = 4.;
1162   const Double_t kb=1.;
1163   Double_t yj;
1164   Double_t y=TMath::Abs(*py);
1165   //
1166   if (y < ky0)
1167     yj=kb;
1168   else
1169     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
1170   return yj;
1171 }
1172
1173 //                 particle composition
1174 //
1175 Int_t AliGenMUONlib::IpChi_c0(TRandom *)
1176 {
1177 // Chi composition
1178     return 10441;
1179 }
1180 //
1181 Int_t AliGenMUONlib::IpChi_c1(TRandom *)
1182 {
1183 // Chi composition
1184     return 20443;
1185 }
1186 Int_t AliGenMUONlib::IpChi_c2(TRandom *)
1187 {
1188 // Chi_c2 prime composition
1189     return 445;
1190 }
1191 Int_t AliGenMUONlib::IpChi_c(TRandom *)
1192 {
1193 // Chi composition
1194   Int_t ip;
1195   Float_t r = gRandom->Rndm();
1196   if (r < 0.001) {
1197     ip = 10441;
1198   } else if( r < 0.377 ) {
1199     ip = 20443;
1200   } else {
1201     ip = 445;
1202   }
1203   return ip;
1204 }
1205
1206
1207 //_____________________________________________________________
1208
1209 typedef Int_t (*GenFuncIp) (TRandom *);
1210 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* /*tname*/) const
1211 {
1212 // Return pointer to particle type parameterisation
1213     GenFuncIp func;
1214     switch (param)
1215     {
1216     case kPhi:
1217         func=IpPhi;
1218         break;
1219     case kEta:
1220         func=IpEta;
1221         break;
1222     case kOmega:
1223         func=IpOmega;
1224         break;
1225     case kJpsiFamily:
1226         func=IpJpsiFamily;
1227         break;
1228     case kPsiP:
1229         func=IpPsiP;
1230         break;
1231     case kJpsi:
1232     case kJpsiFromB:
1233         func=IpJpsi;
1234         break;
1235     case kUpsilon:
1236         func=IpUpsilon;
1237         break;
1238     case kUpsilonFamily:
1239       func=IpUpsilonFamily;
1240       break;
1241     case kUpsilonP:
1242         func=IpUpsilonP;
1243         break;
1244     case kUpsilonPP:
1245         func=IpUpsilonPP;
1246         break;
1247     case kCharm:
1248         func=IpCharm;
1249         break;
1250     case kBeauty:
1251         func=IpBeauty;
1252         break;
1253     case kPion:
1254         func=IpPion;
1255         break;
1256     case kKaon:
1257         func=IpKaon;
1258         break;
1259     case kChi_c0:
1260         func=IpChi_c0;
1261         break;
1262     case kChi_c1:
1263         func=IpChi_c1;
1264         break;
1265     case kChi_c2:
1266         func=IpChi_c2;
1267         break;
1268     case kChi_c:
1269         func=IpChi_c;
1270         break;
1271     default:
1272         func=0;
1273         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
1274     }
1275     return func;
1276 }
1277
1278
1279
1280 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0,
1281                                    Float_t dx,
1282                                    Int_t n, Int_t no)
1283 {
1284 //
1285 // Neville's alorithm for interpolation
1286 //
1287 // x:  x-value
1288 // y:  Input array
1289 // x0: minimum x
1290 // dx: step size
1291 //  n: number of data points
1292 // no: order of polynom
1293 //
1294     Float_t*  c = new Float_t[n];
1295     Float_t*  d = new Float_t[n];
1296     Int_t m, i;
1297     for (i = 0; i < n; i++) {
1298         c[i] = y[i];
1299         d[i] = y[i];
1300     }
1301
1302     Int_t   ns  = int((x - x0)/dx);
1303
1304     Float_t y1  = y[ns];
1305     ns--;
1306     for (m = 0; m < no; m++) {
1307         for (i = 0; i < n-m; i++) {
1308             Float_t ho = x0 + Float_t(i) * dx - x;
1309             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
1310             Float_t w  = c[i+1] - d[i];
1311             Float_t den = ho-hp;
1312             den = w/den;
1313             d[i] = hp * den;
1314             c[i] = ho * den;
1315         }
1316         Float_t dy;
1317
1318         if (2*ns < (n-m-1)) {
1319             dy  = c[ns+1];
1320         } else {
1321             dy  = d[ns--];
1322         }
1323         y1 += dy;}
1324     delete[] c;
1325     delete[] d;
1326
1327     return y1;
1328 }
1329
1330