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.14.6.1  2002/06/10 14:57:41  hristov
  19 Merged with v3-08-02
  20
  21 Revision 1.15  2002/04/17 10:11:51  morsch
  22 Coding Rule violations corrected.
  23
  24 Revision 1.14  2002/02/22 17:26:43  morsch
  25 Eta and omega added.
  26
  27 Revision 1.13  2001/03/27 11:01:04  morsch
  28 Charm pt-distribution corrected. More realistic y-distribution for pi and K.
  29
  30 Revision 1.12  2001/03/09 13:01:41  morsch
  31 - enum constants for paramterisation type (particle family) moved to AliGen*lib.h
  32 - use AliGenGSIlib::kUpsilon, AliGenPHOSlib::kEtaPrime to access the constants
  33
  34 Revision 1.11  2000/11/30 07:12:50  alibrary
  35 Introducing new Rndm and QA classes
  36
  37 Revision 1.10  2000/06/29 21:08:27  morsch
  38 All paramatrisation libraries derive from the pure virtual base class AliGenLib.
  39 This allows to pass a pointer to a library directly to AliGenParam and avoids the
  40 use of function pointers in Config.C.
  41
  42 Revision 1.9  2000/06/14 15:20:56  morsch
  43 Include clean-up (IH)
  44
  45 Revision 1.8  2000/06/09 20:32:11  morsch
  46 All coding rule violations except RS3 corrected
  47
  48 Revision 1.7  2000/05/02 08:12:13  morsch
  49 Coding rule violations corrected.
  50
  51 Revision 1.6  1999/09/29 09:24:14  fca
  52 Introduction of the Copyright and cvs Log
  53
  54 */
  55
  56 // Library class for particle pt and y distributions used for
  57 // muon spectrometer simulations.
  58 // To be used with AliGenParam.
  59 // The following particle typed can be simulated:
  60 // pi, K, phi, omega, eta, J/Psi, Upsilon, charm and beauty mesons.
  61 //
  62 // andreas.morsch@cern.ch
  63 //
  64
  65 #include "TMath.h"
  66 #include "TRandom.h"
  67
  68 #include "AliGenMUONlib.h"
  69
  70 ClassImp(AliGenMUONlib)
  71 //
  72 //  Pions
  73 Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *dummy)
  74 {
  75 //
  76 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
  77 //     POWER LAW FOR PT > 500 MEV
  78 //     MT SCALING BELOW (T=160 MEV)
  79 //
  80   const Double_t kp0 = 1.3;
  81   const Double_t kxn = 8.28;
  82   const Double_t kxlim=0.5;
  83   const Double_t kt=0.160;
  84   const Double_t kxmpi=0.139;
  85   const Double_t kb=1.;
  86   Double_t y, y1, xmpi2, ynorm, a;
  87   Double_t x=*px;
  88   //
  89   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
  90   xmpi2=kxmpi*kxmpi;
  91   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
  92   a=ynorm/y1;
  93   if (x > kxlim)
  94     y=a*TMath::Power(kp0/(kp0+x),kxn);
  95   else
  96     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
  97   return y*x;
  98 }
  99 //
 100 // y-distribution
 101 //
 102 Double_t AliGenMUONlib::YPion( Double_t *py, Double_t *dummy)
 103 {
 104 // Pion y
 105   Double_t y=TMath::Abs(*py);
 106 /*
 107   const Double_t ka    = 7000.;
 108   const Double_t kdy   = 4.;
 109   Double_t ex = y*y/(2*kdy*kdy);
 110   return ka*TMath::Exp(-ex);
 111 */
 112   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
 113
 114 }
 115 //                 particle composition
 116 //
 117 Int_t AliGenMUONlib::IpPion(TRandom *ran)
 118 {
 119 // Pion composition
 120     if (ran->Rndm() < 0.5) {
 121         return  211;
 122     } else {
 123         return -211;
 124     }
 125 }
 126
 127 //____________________________________________________________
 128 //
 129 // Mt-scaling
 130
 131 Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
 132 {
 133   //    SCALING EN MASSE PAR RAPPORT A PTPI
 134   //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
 135   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
 136   //     VALUE MESON/PI AT 5 GEV
 137   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
 138   np--;
 139   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
 140   Double_t fmax2=f5/kfmax[np];
 141   // PIONS
 142   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
 143   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
 144                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
 145   return fmtscal*ptpion;
 146 }
 147 //
 148 // kaon
 149 //
 150 //                pt-distribution
 151 //____________________________________________________________
 152 Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *dummy)
 153 {
 154 // Kaon pT
 155   return PtScal(*px,2);
 156 }
 157
 158 // y-distribution
 159 //____________________________________________________________
 160 Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t *dummy)
 161 {
 162 // Kaon y
 163   Double_t y=TMath::Abs(*py);
 164 /*
 165   const Double_t ka    = 1000.;
 166   const Double_t kdy   = 4.;
 167   //
 168   Double_t ex = y*y/(2*kdy*kdy);
 169   return ka*TMath::Exp(-ex);
 170 */
 171
 172   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
 173 }
 174
 175 //                 particle composition
 176 //
 177 Int_t AliGenMUONlib::IpKaon(TRandom *ran)
 178 {
 179 // Kaon composition
 180     if (ran->Rndm() < 0.5) {
 181         return  321;
 182     } else {
 183         return -321;
 184     }
 185 }
 186
 187 //                    J/Psi
 188 //
 189 //
 190 //                pt-distribution
 191 //____________________________________________________________
 192 Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *dummy)
 193 {
 194 // J/Psi pT
 195   const Double_t kpt0 = 4.;
 196   const Double_t kxn  = 3.6;
 197   Double_t x=*px;
 198   //
 199   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 200   return x/TMath::Power(pass1,kxn);
 201 }
 202 //
 203 //               y-distribution
 204 //____________________________________________________________
 205 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *dummy)
 206 {
 207 // J/psi y
 208   const Double_t ky0 = 4.;
 209   const Double_t kb=1.;
 210   Double_t yj;
 211   Double_t y=TMath::Abs(*py);
 212   //
 213   if (y < ky0)
 214     yj=kb;
 215   else
 216     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 217   return yj;
 218 }
 219 //                 particle composition
 220 //
 221 Int_t AliGenMUONlib::IpJpsi(TRandom *)
 222 {
 223 // J/Psi composition
 224     return 443;
 225 }
 226
 227 //                      Upsilon
 228 //
 229 //
 230 //                  pt-distribution
 231 //____________________________________________________________
 232 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t *dummy )
 233 {
 234 // Upsilon pT
 235   const Double_t kpt0 = 5.3;
 236   const Double_t kxn  = 2.5;
 237   Double_t x=*px;
 238   //
 239   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 240   return x/TMath::Power(pass1,kxn);
 241 }
 242 //
 243 //                    y-distribution
 244 //
 245 //____________________________________________________________
 246 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *dummy)
 247 {
 248 // Upsilon y
 249   const Double_t ky0 = 3.;
 250   const Double_t kb=1.;
 251   Double_t yu;
 252   Double_t y=TMath::Abs(*py);
 253   //
 254   if (y < ky0)
 255     yu=kb;
 256   else
 257     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 258   return yu;
 259 }
 260 //                 particle composition
 261 //
 262 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
 263 {
 264 // y composition
 265     return 553;
 266 }
 267
 268 //
 269 //                        Phi
 270 //
 271 //
 272 //    pt-distribution (by scaling of pion distribution)
 273 //____________________________________________________________
 274 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *dummy)
 275 {
 276 // Phi pT
 277   return PtScal(*px,7);
 278 }
 279 //    y-distribution
 280 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *dummy)
 281 {
 282 // Phi y
 283     Double_t *dum=0;
 284     return YJpsi(px,dum);
 285 }
 286 //                 particle composition
 287 //
 288 Int_t AliGenMUONlib::IpPhi(TRandom *)
 289 {
 290 // Phi composition
 291     return 333;
 292 }
 293
 294 //
 295 //                        omega
 296 //
 297 //
 298 //    pt-distribution (by scaling of pion distribution)
 299 //____________________________________________________________
 300 Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t *dummy)
 301 {
 302 // Omega pT
 303   return PtScal(*px,5);
 304 }
 305 //    y-distribution
 306 Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t *dummy)
 307 {
 308 // Omega y
 309     Double_t *dum=0;
 310     return YJpsi(px,dum);
 311 }
 312 //                 particle composition
 313 //
 314 Int_t AliGenMUONlib::IpOmega(TRandom *)
 315 {
 316 // Omega composition
 317     return 223;
 318 }
 319
 320
 321 //
 322 //                        Eta
 323 //
 324 //
 325 //    pt-distribution (by scaling of pion distribution)
 326 //____________________________________________________________
 327 Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t *dummy)
 328 {
 329 // Eta pT
 330   return PtScal(*px,3);
 331 }
 332 //    y-distribution
 333 Double_t AliGenMUONlib::YEta( Double_t *px, Double_t *dummy)
 334 {
 335 // Eta y
 336     Double_t *dum=0;
 337     return YJpsi(px,dum);
 338 }
 339 //                 particle composition
 340 //
 341 Int_t AliGenMUONlib::IpEta(TRandom *)
 342 {
 343 // Eta composition
 344     return 221;
 345 }
 346
 347 //
 348 //                        Charm
 349 //
 350 //
 351 //                    pt-distribution
 352 //____________________________________________________________
 353 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *dummy)
 354 {
 355 // Charm pT
 356   const Double_t kpt0 = 4.08;
 357   const Double_t kxn  = 9.40;
 358
 359   Double_t x=*px;
 360   //
 361   Double_t pass1 = 1.+(x/kpt0);
 362   return x/TMath::Power(pass1,kxn);
 363 }
 364 //                  y-distribution
 365 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *dummy)
 366 {
 367 // Charm y
 368     Double_t *dum=0;
 369     return YJpsi(px,dum);
 370 }
 371
 372 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
 373 {
 374 // Charm composition
 375     Float_t random;
 376     Int_t ip;
 377 //    411,421,431,4122
 378     random = ran->Rndm();
 379     if (random < 0.5) {
 380         ip=411;
 381     } else if (random < 0.75) {
 382         ip=421;
 383     } else if (random < 0.90) {
 384         ip=431;
 385     } else {
 386         ip=4122;
 387     }
 388     if (ran->Rndm() < 0.5) {ip=-ip;}
 389
 390     return ip;
 391 }
 392
 393
 394 //
 395 //                        Beauty
 396 //
 397 //
 398 //                    pt-distribution
 399 //____________________________________________________________
 400 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *dummy)
 401 {
 402 // Beauty pT
 403   const Double_t kpt0 = 4.;
 404   const Double_t kxn  = 3.6;
 405   Double_t x=*px;
 406   //
 407   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 408   return x/TMath::Power(pass1,kxn);
 409 }
 410 //                     y-distribution
 411 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *dummy)
 412 {
 413 // Beauty y
 414     Double_t *dum=0;
 415     return YJpsi(px,dum);
 416 }
 417
 418 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
 419 {
 420 // Beauty Composition
 421     Float_t random;
 422     Int_t ip;
 423     random = ran->Rndm();
 424     if (random < 0.5) {
 425         ip=511;
 426     } else if (random < 0.75) {
 427         ip=521;
 428     } else if (random < 0.90) {
 429         ip=531;
 430     } else {
 431         ip=5122;
 432     }
 433     if (ran->Rndm() < 0.5) {ip=-ip;}
 434
 435     return ip;
 436 }
 437
 438 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
 439 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
 440 {
 441 // Return pointer to pT parameterisation
 442     GenFunc func;
 443     switch (param)
 444     {
 445     case kPhi:
 446         func=PtPhi;
 447         break;
 448     case kOmega:
 449         func=PtOmega;
 450         break;
 451     case kEta:
 452         func=PtEta;
 453         break;
 454     case kJpsi:
 455         func=PtJpsi;
 456         break;
 457     case kUpsilon:
 458         func=PtUpsilon;
 459         break;
 460     case kCharm:
 461         func=PtCharm;
 462         break;
 463     case kBeauty:
 464         func=PtBeauty;
 465         break;
 466     case kPion:
 467         func=PtPion;
 468         break;
 469     case kKaon:
 470         func=PtKaon;
 471         break;
 472     default:
 473         func=0;
 474         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
 475     }
 476     return func;
 477 }
 478
 479 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
 480 {
 481 // Return pointer to y- parameterisation
 482     GenFunc func;
 483     switch (param)
 484     {
 485     case kPhi:
 486         func=YPhi;
 487         break;
 488     case kEta:
 489         func=YEta;
 490         break;
 491     case kOmega:
 492         func=YOmega;
 493         break;
 494     case kJpsi:
 495         func=YJpsi;
 496         break;
 497     case kUpsilon:
 498         func=YUpsilon;
 499         break;
 500     case kCharm:
 501         func=YCharm;
 502         break;
 503     case kBeauty:
 504         func=YBeauty;
 505         break;
 506     case kPion:
 507         func=YPion;
 508         break;
 509     case kKaon:
 510         func=YKaon;
 511         break;
 512     default:
 513         func=0;
 514         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
 515     }
 516     return func;
 517 }
 518 typedef Int_t (*GenFuncIp) (TRandom *);
 519 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* tname) const
 520 {
 521 // Return pointer to particle type parameterisation
 522     GenFuncIp func;
 523     switch (param)
 524     {
 525     case kPhi:
 526         func=IpPhi;
 527         break;
 528     case kEta:
 529         func=IpEta;
 530         break;
 531     case kOmega:
 532         func=IpOmega;
 533         break;
 534     case kJpsi:
 535         func=IpJpsi;
 536         break;
 537     case kUpsilon:
 538         func=IpUpsilon;
 539         break;
 540     case kCharm:
 541         func=IpCharm;
 542         break;
 543     case kBeauty:
 544         func=IpBeauty;
 545         break;
 546     case kPion:
 547         func=IpPion;
 548         break;
 549     case kKaon:
 550         func=IpKaon;
 551         break;
 552     default:
 553         func=0;
 554         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
 555     }
 556     return func;
 557 }
 558
 559
 560
 561