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