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