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