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.6  1999/09/29 09:24:14  fca
  19 Introduction of the Copyright and cvs Log
  20
  21 */
  22
  23 #include "AliGenMUONlib.h"
  24 #include "AliRun.h"
  25 ClassImp(AliGenMUONlib)
  26 //
  27 //  Pions
  28 Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *dummy)
  29 {
  30 //
  31 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
  32 //     POWER LAW FOR PT > 500 MEV
  33 //     MT SCALING BELOW (T=160 MEV)
  34 //
  35   const Double_t kp0 = 1.3;
  36   const Double_t kxn = 8.28;
  37   const Double_t kxlim=0.5;
  38   const Double_t kt=0.160;
  39   const Double_t kxmpi=0.139;
  40   const Double_t kb=1.;
  41   Double_t y, y1, xmpi2, ynorm, a;
  42   Double_t x=*px;
  43   //
  44   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
  45   xmpi2=kxmpi*kxmpi;
  46   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
  47   a=ynorm/y1;
  48   if (x > kxlim)
  49     y=a*TMath::Power(kp0/(kp0+x),kxn);
  50   else
  51     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
  52   return y*x;
  53 }
  54 //
  55 // y-distribution
  56 //
  57 Double_t AliGenMUONlib::YPion( Double_t *py, Double_t *dummy)
  58 {
  59 // Pion y
  60   const Double_t ka    = 7000.;
  61   const Double_t kdy   = 4.;
  62
  63   Double_t y=TMath::Abs(*py);
  64   //
  65   Double_t ex = y*y/(2*kdy*kdy);
  66   return ka*TMath::Exp(-ex);
  67 }
  68 //                 particle composition
  69 //
  70 Int_t AliGenMUONlib::IpPion()
  71 {
  72 // Pion composition
  73     Float_t random[1];
  74     gMC->Rndm(random,1);
  75     if (random[0] < 0.5) {
  76         return  211;
  77     } else {
  78         return -211;
  79     }
  80 }
  81
  82 //____________________________________________________________
  83 //
  84 // Mt-scaling
  85
  86 Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
  87 {
  88   //    SCALING EN MASSE PAR RAPPORT A PTPI
  89   //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
  90   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
  91   //     VALUE MESON/PI AT 5 GEV
  92   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
  93   np--;
  94   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
  95   Double_t fmax2=f5/kfmax[np];
  96   // PIONS
  97   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
  98   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
  99                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
 100   return fmtscal*ptpion;
 101 }
 102 //
 103 // kaon
 104 //
 105 //                pt-distribution
 106 //____________________________________________________________
 107 Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *dummy)
 108 {
 109 // Kaon pT
 110   return PtScal(*px,2);
 111 }
 112
 113 // y-distribution
 114 //____________________________________________________________
 115 Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t *dummy)
 116 {
 117 // Kaon y
 118   const Double_t ka    = 1000.;
 119   const Double_t kdy   = 4.;
 120
 121
 122   Double_t y=TMath::Abs(*py);
 123   //
 124   Double_t ex = y*y/(2*kdy*kdy);
 125   return ka*TMath::Exp(-ex);
 126 }
 127
 128 //                 particle composition
 129 //
 130 Int_t AliGenMUONlib::IpKaon()
 131 {
 132 // Kaon composition
 133     Float_t random[1];
 134     gMC->Rndm(random,1);
 135     if (random[0] < 0.5) {
 136         return  321;
 137     } else {
 138         return -321;
 139     }
 140 }
 141
 142 //                    J/Psi
 143 //
 144 //
 145 //                pt-distribution
 146 //____________________________________________________________
 147 Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *dummy)
 148 {
 149 // J/Psi pT
 150   const Double_t kpt0 = 4.;
 151   const Double_t kxn  = 3.6;
 152   Double_t x=*px;
 153   //
 154   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 155   return x/TMath::Power(pass1,kxn);
 156 }
 157 //
 158 //               y-distribution
 159 //____________________________________________________________
 160 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *dummy)
 161 {
 162 // J/psi y
 163   const Double_t ky0 = 4.;
 164   const Double_t kb=1.;
 165   Double_t yj;
 166   Double_t y=TMath::Abs(*py);
 167   //
 168   if (y < ky0)
 169     yj=kb;
 170   else
 171     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 172   return yj;
 173 }
 174 //                 particle composition
 175 //
 176 Int_t AliGenMUONlib::IpJpsi()
 177 {
 178 // J/Psi composition
 179     return 443;
 180 }
 181
 182 //                      Upsilon
 183 //
 184 //
 185 //                  pt-distribution
 186 //____________________________________________________________
 187 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t *dummy )
 188 {
 189 // Upsilon pT
 190   const Double_t kpt0 = 5.3;
 191   const Double_t kxn  = 2.5;
 192   Double_t x=*px;
 193   //
 194   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 195   return x/TMath::Power(pass1,kxn);
 196 }
 197 //
 198 //                    y-distribution
 199 //
 200 //____________________________________________________________
 201 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *dummy)
 202 {
 203 // Upsilon y
 204   const Double_t ky0 = 3.;
 205   const Double_t kb=1.;
 206   Double_t yu;
 207   Double_t y=TMath::Abs(*py);
 208   //
 209   if (y < ky0)
 210     yu=kb;
 211   else
 212     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 213   return yu;
 214 }
 215 //                 particle composition
 216 //
 217 Int_t AliGenMUONlib::IpUpsilon()
 218 {
 219 // y composition
 220     return 553;
 221 }
 222
 223 //
 224 //                        Phi
 225 //
 226 //
 227 //    pt-distribution (by scaling of pion distribution)
 228 //____________________________________________________________
 229 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *dummy)
 230 {
 231 // Phi pT
 232   return PtScal(*px,7);
 233 }
 234 //    y-distribution
 235 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *dummy)
 236 {
 237 // Phi y
 238     Double_t *dum=0;
 239     return YJpsi(px,dum);
 240 }
 241 //                 particle composition
 242 //
 243 Int_t AliGenMUONlib::IpPhi()
 244 {
 245 // Phi composition
 246     return 41;
 247 }
 248
 249 //
 250 //                        Charm
 251 //
 252 //
 253 //                    pt-distribution
 254 //____________________________________________________________
 255 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *dummy)
 256 {
 257 // Charm pT
 258   const Double_t kpt0 = 4.08;
 259   const Double_t kxn  = 9.40;
 260   Double_t x=*px;
 261   //
 262   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 263   return x/TMath::Power(pass1,kxn);
 264 }
 265 //                  y-distribution
 266 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *dummy)
 267 {
 268 // Charm y
 269     Double_t *dum=0;
 270     return YJpsi(px,dum);
 271 }
 272
 273 Int_t AliGenMUONlib::IpCharm()
 274 {
 275 // Charm composition
 276     Float_t random[2];
 277     Int_t ip;
 278 //    411,421,431,4122
 279     gMC->Rndm(random,2);
 280     if (random[0] < 0.5) {
 281         ip=411;
 282     } else if (random[0] < 0.75) {
 283         ip=421;
 284     } else if (random[0] < 0.90) {
 285         ip=431;
 286     } else {
 287         ip=4122;
 288     }
 289     if (random[1] < 0.5) {ip=-ip;}
 290
 291     return ip;
 292 }
 293
 294
 295 //
 296 //                        Beauty
 297 //
 298 //
 299 //                    pt-distribution
 300 //____________________________________________________________
 301 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *dummy)
 302 {
 303 // Beauty pT
 304   const Double_t kpt0 = 4.;
 305   const Double_t kxn  = 3.6;
 306   Double_t x=*px;
 307   //
 308   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 309   return x/TMath::Power(pass1,kxn);
 310 }
 311 //                     y-distribution
 312 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *dummy)
 313 {
 314 // Beauty y
 315     Double_t *dum=0;
 316     return YJpsi(px,dum);
 317 }
 318
 319 Int_t AliGenMUONlib::IpBeauty()
 320 {
 321 // Beauty Composition
 322     Float_t random[2];
 323     Int_t ip;
 324     gMC->Rndm(random,2);
 325     if (random[0] < 0.5) {
 326         ip=511;
 327     } else if (random[0] < 0.75) {
 328         ip=521;
 329     } else if (random[0] < 0.90) {
 330         ip=531;
 331     } else {
 332         ip=5122;
 333     }
 334     if (random[1] < 0.5) {ip=-ip;}
 335
 336     return ip;
 337 }
 338
 339 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
 340 GenFunc AliGenMUONlib::GetPt(Param_t param)
 341 {
 342 // Return pointer to pT parameterisation
 343     GenFunc func;
 344     switch (param)
 345     {
 346     case phi_p:
 347         func=PtPhi;
 348         break;
 349     case jpsi_p:
 350         func=PtJpsi;
 351         break;
 352     case upsilon_p:
 353         func=PtUpsilon;
 354         break;
 355     case charm_p:
 356         func=PtCharm;
 357         break;
 358     case beauty_p:
 359         func=PtBeauty;
 360         break;
 361     case pion_p:
 362         func=PtPion;
 363         break;
 364     case kaon_p:
 365         func=PtKaon;
 366         break;
 367     default:
 368         func=0;
 369         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
 370     }
 371     return func;
 372 }
 373
 374 GenFunc AliGenMUONlib::GetY(Param_t param)
 375 {
 376 // Return pointer to y- parameterisation
 377     GenFunc func;
 378     switch (param)
 379     {
 380     case phi_p:
 381         func=YPhi;
 382         break;
 383     case jpsi_p:
 384         func=YJpsi;
 385         break;
 386     case upsilon_p:
 387         func=YUpsilon;
 388         break;
 389     case charm_p:
 390         func=YCharm;
 391         break;
 392     case beauty_p:
 393         func=YBeauty;
 394         break;
 395     case pion_p:
 396         func=YPion;
 397         break;
 398     case kaon_p:
 399         func=YKaon;
 400         break;
 401     default:
 402         func=0;
 403         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
 404     }
 405     return func;
 406 }
 407 typedef Int_t (*GenFuncIp) ();
 408 GenFuncIp AliGenMUONlib::GetIp(Param_t param)
 409 {
 410 // Return pointer to particle type parameterisation
 411     GenFuncIp func;
 412     switch (param)
 413     {
 414     case phi_p:
 415         func=IpPhi;
 416         break;
 417     case jpsi_p:
 418         func=IpJpsi;
 419         break;
 420     case upsilon_p:
 421         func=IpUpsilon;
 422         break;
 423     case charm_p:
 424         func=IpCharm;
 425         break;
 426     case beauty_p:
 427         func=IpBeauty;
 428         break;
 429     case pion_p:
 430         func=IpPion;
 431         break;
 432     case kaon_p:
 433         func=IpKaon;
 434         break;
 435     default:
 436         func=0;
 437         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
 438     }
 439     return func;
 440 }
 441
 442
 443
 444