EVGEN/AliGenMUONlib.cxx

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