1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.21 2003/03/25 12:01:11 morsch
19 Quarkonia for pp @ 14 TeV added.
21 Revision 1.20 2003/03/13 11:54:39 morsch
22 Limited pT range for parameterized Upsilon and J/Psi pT distributions.
24 Revision 1.19 2003/02/24 16:46:11 morsch
25 New parameterisation for Psi and Upsilon (PbPb)
27 Revision 1.18 2002/11/07 09:13:42 morsch
28 Use "Vogt" to label new distributions.
30 Revision 1.17 2002/11/07 09:06:10 morsch
31 J/Psi and Upsilon pt and y distributions from R. Vogt 2002 added.
33 Revision 1.16 2002/10/14 14:55:35 hristov
34 Merging the VirtualMC branch to the main development branch (HEAD)
36 Revision 1.14.6.1 2002/06/10 14:57:41 hristov
39 Revision 1.15 2002/04/17 10:11:51 morsch
40 Coding Rule violations corrected.
42 Revision 1.14 2002/02/22 17:26:43 morsch
45 Revision 1.13 2001/03/27 11:01:04 morsch
46 Charm pt-distribution corrected. More realistic y-distribution for pi and K.
48 Revision 1.12 2001/03/09 13:01:41 morsch
49 - enum constants for paramterisation type (particle family) moved to AliGen*lib.h
50 - use AliGenGSIlib::kUpsilon, AliGenPHOSlib::kEtaPrime to access the constants
52 Revision 1.11 2000/11/30 07:12:50 alibrary
53 Introducing new Rndm and QA classes
55 Revision 1.10 2000/06/29 21:08:27 morsch
56 All paramatrisation libraries derive from the pure virtual base class AliGenLib.
57 This allows to pass a pointer to a library directly to AliGenParam and avoids the
58 use of function pointers in Config.C.
60 Revision 1.9 2000/06/14 15:20:56 morsch
63 Revision 1.8 2000/06/09 20:32:11 morsch
64 All coding rule violations except RS3 corrected
66 Revision 1.7 2000/05/02 08:12:13 morsch
67 Coding rule violations corrected.
69 Revision 1.6 1999/09/29 09:24:14 fca
70 Introduction of the Copyright and cvs Log
74 // Library class for particle pt and y distributions used for
75 // muon spectrometer simulations.
76 // To be used with AliGenParam.
77 // The following particle typed can be simulated:
78 // pi, K, phi, omega, eta, J/Psi, Upsilon, charm and beauty mesons.
80 // andreas.morsch@cern.ch
86 #include "AliGenMUONlib.h"
88 ClassImp(AliGenMUONlib)
91 Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *dummy)
94 // PT-PARAMETERIZATION CDF, PRL 61(88) 1819
95 // POWER LAW FOR PT > 500 MEV
96 // MT SCALING BELOW (T=160 MEV)
98 const Double_t kp0 = 1.3;
99 const Double_t kxn = 8.28;
100 const Double_t kxlim=0.5;
101 const Double_t kt=0.160;
102 const Double_t kxmpi=0.139;
103 const Double_t kb=1.;
104 Double_t y, y1, xmpi2, ynorm, a;
107 y1=TMath::Power(kp0/(kp0+kxlim),kxn);
109 ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
112 y=a*TMath::Power(kp0/(kp0+x),kxn);
114 y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
120 Double_t AliGenMUONlib::YPion( Double_t *py, Double_t *dummy)
123 Double_t y=TMath::Abs(*py);
125 const Double_t ka = 7000.;
126 const Double_t kdy = 4.;
127 Double_t ex = y*y/(2*kdy*kdy);
128 return ka*TMath::Exp(-ex);
130 return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
133 // particle composition
135 Int_t AliGenMUONlib::IpPion(TRandom *ran)
138 if (ran->Rndm() < 0.5) {
145 //____________________________________________________________
149 Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
151 // SCALING EN MASSE PAR RAPPORT A PTPI
152 // MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
153 const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
154 // VALUE MESON/PI AT 5 GEV
155 const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
157 Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
158 Double_t fmax2=f5/kfmax[np];
160 Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
161 Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
162 (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
163 return fmtscal*ptpion;
169 //____________________________________________________________
170 Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *dummy)
173 return PtScal(*px,2);
177 //____________________________________________________________
178 Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t *dummy)
181 Double_t y=TMath::Abs(*py);
183 const Double_t ka = 1000.;
184 const Double_t kdy = 4.;
186 Double_t ex = y*y/(2*kdy*kdy);
187 return ka*TMath::Exp(-ex);
190 return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
193 // particle composition
195 Int_t AliGenMUONlib::IpKaon(TRandom *ran)
198 if (ran->Rndm() < 0.5) {
209 //____________________________________________________________
210 Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *dummy)
213 const Double_t kpt0 = 4.;
214 const Double_t kxn = 3.6;
217 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
218 return x/TMath::Power(pass1,kxn);
221 Double_t AliGenMUONlib::PtJpsiPbPb( Double_t *px, Double_t *dummy)
228 // mc = 1.4 GeV, pt-kick 1 GeV
232 -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00,
233 -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
240 while (j > 0) y = y * x +c[--j];
241 y = x * TMath::Exp(y);
247 Double_t AliGenMUONlib::PtJpsiPP( Double_t *px, Double_t *dummy)
254 // mc = 1.4 GeV, pt-kick 1 GeV
257 Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
263 while (j > 0) y = y * x +c[--j];
264 y = x * TMath::Exp(y);
273 //____________________________________________________________
274 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *dummy)
277 const Double_t ky0 = 4.;
278 const Double_t kb=1.;
280 Double_t y=TMath::Abs(*py);
285 yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
290 Double_t AliGenMUONlib::YJpsiPbPb( Double_t *px, Double_t *dummy)
300 // mc = 1.4 GeV, pt-kick 1 GeV
302 Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
303 Double_t x = TMath::Abs(px[0]);
311 while (j > 0) y = y * x + c[--j];
319 Double_t AliGenMUONlib::YJpsiPP( Double_t *px, Double_t *dummy)
329 // mc = 1.4 GeV, pt-kick 1 GeV
332 Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
333 Double_t x = TMath::Abs(px[0]);
337 y = 96.455 - 0.8483 * x * x;
338 } else if (x < 7.9) {
341 while (j > 0) y = y * x + c[--j];
349 // particle composition
351 Int_t AliGenMUONlib::IpJpsi(TRandom *)
355 Float_t r = gRandom->Rndm();
368 //____________________________________________________________
369 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t *dummy )
372 const Double_t kpt0 = 5.3;
373 const Double_t kxn = 2.5;
376 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
377 return x/TMath::Power(pass1,kxn);
380 Double_t AliGenMUONlib::PtUpsilonPbPb( Double_t *px, Double_t *dummy)
390 // mc = 1.4 GeV, pt-kick 1 GeV
394 -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,
395 -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
401 while (j > 0) y = y * x +c[--j];
402 y = x * TMath::Exp(y);
409 Double_t AliGenMUONlib::PtUpsilonPP( Double_t *px, Double_t *dummy)
419 // mc = 1.4 GeV, pt-kick 1 GeV
422 Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,
423 -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
429 while (j > 0) y = y * x +c[--j];
430 y = x * TMath::Exp(y);
440 //____________________________________________________________
441 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *dummy)
444 const Double_t ky0 = 3.;
445 const Double_t kb=1.;
447 Double_t y=TMath::Abs(*py);
452 yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
457 Double_t AliGenMUONlib::YUpsilonPbPb( Double_t *px, Double_t *dummy)
467 // mc = 1.4 GeV, pt-kick 1 GeV
470 Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
471 -2.99753e-09, 1.28895e-05};
474 if (TMath::Abs(x) > 5.55) return 0.;
476 Double_t y = c[j = 6];
477 while (j > 0) y = y * x +c[--j];
481 Double_t AliGenMUONlib::YUpsilonPP( Double_t *px, Double_t *dummy)
491 // mc = 1.4 GeV, pt-kick 1 GeV
493 Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04,
494 -6.20539e-10, 1.29943e-05};
497 if (TMath::Abs(x) > 6.2) return 0.;
499 Double_t y = c[j = 6];
500 while (j > 0) y = y * x +c[--j];
504 // particle composition
506 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
510 Float_t r = gRandom->Rndm();
514 } else if (r < 0.896) {
527 // pt-distribution (by scaling of pion distribution)
528 //____________________________________________________________
529 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *dummy)
532 return PtScal(*px,7);
535 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *dummy)
539 return YJpsi(px,dum);
541 // particle composition
543 Int_t AliGenMUONlib::IpPhi(TRandom *)
553 // pt-distribution (by scaling of pion distribution)
554 //____________________________________________________________
555 Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t *dummy)
558 return PtScal(*px,5);
561 Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t *dummy)
565 return YJpsi(px,dum);
567 // particle composition
569 Int_t AliGenMUONlib::IpOmega(TRandom *)
580 // pt-distribution (by scaling of pion distribution)
581 //____________________________________________________________
582 Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t *dummy)
585 return PtScal(*px,3);
588 Double_t AliGenMUONlib::YEta( Double_t *px, Double_t *dummy)
592 return YJpsi(px,dum);
594 // particle composition
596 Int_t AliGenMUONlib::IpEta(TRandom *)
607 //____________________________________________________________
608 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *dummy)
611 const Double_t kpt0 = 4.08;
612 const Double_t kxn = 9.40;
616 Double_t pass1 = 1.+(x/kpt0);
617 return x/TMath::Power(pass1,kxn);
620 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *dummy)
624 return YJpsi(px,dum);
627 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
633 random = ran->Rndm();
636 } else if (random < 0.75) {
638 } else if (random < 0.90) {
643 if (ran->Rndm() < 0.5) {ip=-ip;}
654 //____________________________________________________________
655 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *dummy)
658 const Double_t kpt0 = 4.;
659 const Double_t kxn = 3.6;
662 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
663 return x/TMath::Power(pass1,kxn);
666 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *dummy)
670 return YJpsi(px,dum);
673 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
675 // Beauty Composition
678 random = ran->Rndm();
681 } else if (random < 0.75) {
683 } else if (random < 0.90) {
688 if (ran->Rndm() < 0.5) {ip=-ip;}
693 typedef Double_t (*GenFunc) (Double_t*, Double_t*);
694 GenFunc AliGenMUONlib::GetPt(Int_t param, const char* tname) const
696 // Return pointer to pT parameterisation
697 TString sname = TString(tname);
711 if (sname == "Vogt" || sname == "Vogt PbPb") {
713 } else if (sname == "Vogt pp") {
720 if (sname == "Vogt" || sname == "Vogt PbPb") {
722 } else if (sname == "Vogt pp") {
742 printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
747 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
749 TString sname = TString(tname);
751 // Return pointer to y- parameterisation
765 if (sname == "Vogt" || sname == "Vogt PbPb") {
767 } else if (sname == "Vogt pp"){
775 if (sname == "Vogt" || sname == "Vogt PbPb") {
777 } else if (sname == "Vogt pp") {
797 printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
801 typedef Int_t (*GenFuncIp) (TRandom *);
802 GenFuncIp AliGenMUONlib::GetIp(Int_t param, const char* tname) const
804 // Return pointer to particle type parameterisation
837 printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
844 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0,
849 // Neville's alorithm for interpolation
855 // n: number of data points
856 // no: order of polynom
858 Float_t* c = new Float_t[n];
859 Float_t* d = new Float_t[n];
861 for (i = 0; i < n; i++) {
866 Int_t ns = int((x - x0)/dx);
870 for (m = 0; m < no; m++) {
871 for (i = 0; i < n-m; i++) {
872 Float_t ho = x0 + Float_t(i) * dx - x;
873 Float_t hp = x0 + Float_t(i+m+1) * dx - x;
874 Float_t w = c[i+1] - d[i];
882 if (2*ns < (n-m-1)) {