X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=EVGEN%2FAliGenMUONlib.cxx;h=6207544baf69a2b0598c5f53ec3fa697e5b3ccc4;hp=2ec6a0563ff635cb48d30ea2fa953057d3ede44d;hb=88e5db438d136a52020afa73aed7ee0fa5d70417;hpb=34f60c01e27a5f194bd3152deb556a4bc3d933e9 diff --git a/EVGEN/AliGenMUONlib.cxx b/EVGEN/AliGenMUONlib.cxx index 2ec6a0563ff..6207544baf6 100644 --- a/EVGEN/AliGenMUONlib.cxx +++ b/EVGEN/AliGenMUONlib.cxx @@ -13,29 +13,16 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.11 2000/11/30 07:12:50 alibrary -Introducing new Rndm and QA classes - -Revision 1.10 2000/06/29 21:08:27 morsch -All paramatrisation libraries derive from the pure virtual base class AliGenLib. -This allows to pass a pointer to a library directly to AliGenParam and avoids the -use of function pointers in Config.C. - -Revision 1.9 2000/06/14 15:20:56 morsch -Include clean-up (IH) - -Revision 1.8 2000/06/09 20:32:11 morsch -All coding rule violations except RS3 corrected +/* $Id$ */ -Revision 1.7 2000/05/02 08:12:13 morsch -Coding rule violations corrected. - -Revision 1.6 1999/09/29 09:24:14 fca -Introduction of the Copyright and cvs Log - -*/ +// Library class for particle pt and y distributions used for +// muon spectrometer simulations. +// To be used with AliGenParam. +// The following particle typed can be simulated: +// pi, K, phi, omega, eta, J/Psi, Upsilon, charm and beauty mesons. +// +// andreas.morsch@cern.ch +// #include "TMath.h" #include "TRandom.h" @@ -45,7 +32,7 @@ Introduction of the Copyright and cvs Log ClassImp(AliGenMUONlib) // // Pions -Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *dummy) +Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t* /*dummy*/) { // // PT-PARAMETERIZATION CDF, PRL 61(88) 1819 @@ -74,16 +61,18 @@ Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *dummy) // // y-distribution // -Double_t AliGenMUONlib::YPion( Double_t *py, Double_t *dummy) +Double_t AliGenMUONlib::YPion( Double_t *py, Double_t */*dummy*/) { // Pion y + Double_t y=TMath::Abs(*py); +/* const Double_t ka = 7000.; const Double_t kdy = 4.; - - Double_t y=TMath::Abs(*py); - // Double_t ex = y*y/(2*kdy*kdy); return ka*TMath::Exp(-ex); +*/ + return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02; + } // particle composition // @@ -122,7 +111,7 @@ Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np) // // pt-distribution //____________________________________________________________ -Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *dummy) +Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t */*dummy*/) { // Kaon pT return PtScal(*px,2); @@ -130,17 +119,19 @@ Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *dummy) // y-distribution //____________________________________________________________ -Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t *dummy) +Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t */*dummy*/) { // Kaon y + Double_t y=TMath::Abs(*py); +/* const Double_t ka = 1000.; const Double_t kdy = 4.; - - - Double_t y=TMath::Abs(*py); // Double_t ex = y*y/(2*kdy*kdy); return ka*TMath::Exp(-ex); +*/ + + return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02; } // particle composition @@ -160,7 +151,7 @@ Int_t AliGenMUONlib::IpKaon(TRandom *ran) // // pt-distribution //____________________________________________________________ -Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *dummy) +Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t */*dummy*/) { // J/Psi pT const Double_t kpt0 = 4.; @@ -170,10 +161,76 @@ Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *dummy) Double_t pass1 = 1.+(x/kpt0)*(x/kpt0); return x/TMath::Power(pass1,kxn); } + +Double_t AliGenMUONlib::PtJpsiPbPb( Double_t *px, Double_t */*dummy*/) +{ +// J/Psi pT spectrum +// +// R. Vogt 2002 +// PbPb 5.5 TeV +// MRST HO +// mc = 1.4 GeV, pt-kick 1 GeV +// + Float_t x = px[0]; + Float_t c[8] = { + -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00, + -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05 + }; + + Double_t y; + if (x < 10.) { + Int_t j; + y = c[j = 7]; + while (j > 0) y = y * x +c[--j]; + y = x * TMath::Exp(y); + } else { + y = 0.; + } + return y; +} + +Double_t AliGenMUONlib::PtJpsiBPbPb( Double_t *px, Double_t */*dummy*/) +{ +// J/Psi pT spectrum +// B -> J/Psi X + Double_t x0 = 4.0384; + Double_t n = 3.0288; + + Double_t x = px[0]; + Double_t y = x / TMath::Power((1. + (x/x0)*(x/x0)), n); + + return y; +} + + +Double_t AliGenMUONlib::PtJpsiPP( Double_t *px, Double_t */*dummy*/) +{ +// J/Psi pT spectrum +// +// R. Vogt 2002 +// pp 14 TeV +// MRST HO +// mc = 1.4 GeV, pt-kick 1 GeV +// + Float_t x = px[0]; + Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03}; + + Double_t y; + if (x < 10.) { + Int_t j; + y = c[j = 3]; + while (j > 0) y = y * x +c[--j]; + y = x * TMath::Exp(y); + } else { + y = 0.; + } + return y; +} + // // y-distribution //____________________________________________________________ -Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *dummy) +Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t */*dummy*/) { // J/psi y const Double_t ky0 = 4.; @@ -187,6 +244,94 @@ Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *dummy) yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2); return yj; } + + +Double_t AliGenMUONlib::YJpsiPbPb( Double_t *px, Double_t */*dummy*/) +{ + +// +// J/Psi y +// +// +// R. Vogt 2002 +// PbPb 5.5 TeV +// MRST HO +// mc = 1.4 GeV, pt-kick 1 GeV +// + Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01}; + Double_t x = TMath::Abs(px[0]); + Double_t y; + + if (x < 4.) { + y = 31.754; + } else if (x < 6) { + Int_t j; + y = c[j = 4]; + while (j > 0) y = y * x + c[--j]; + } else { + y =0.; + } + + return y; +} + +Double_t AliGenMUONlib::YJpsiPP( Double_t *px, Double_t */*dummy*/) +{ + +// +// J/Psi y +// +// +// R. Vogt 2002 +// pp 14 TeV +// MRST HO +// mc = 1.4 GeV, pt-kick 1 GeV +// + + Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01}; + Double_t x = TMath::Abs(px[0]); + Double_t y; + + if (x < 2.5) { + y = 96.455 - 0.8483 * x * x; + } else if (x < 7.9) { + Int_t j; + y = c[j = 4]; + while (j > 0) y = y * x + c[--j]; + } else { + y =0.; + } + + return y; +} + +Double_t AliGenMUONlib::YJpsiBPbPb( Double_t *px, Double_t */*dummy*/) +{ + +// +// J/Psi from B->J/Psi X +// +// + + + Double_t c[7] = {7.37025e-02, 0., -2.94487e-03, 0., 6.07953e-06, 0., 5.39219e-07}; + + Double_t x = TMath::Abs(px[0]); + Double_t y; + + if (x > 6.) { + y = 0.; + } else { + Int_t j; + y = c[j = 6]; + while (j > 0) y = y * x + c[--j]; + } + + return y; +} + + + // particle composition // Int_t AliGenMUONlib::IpJpsi(TRandom *) @@ -194,13 +339,32 @@ Int_t AliGenMUONlib::IpJpsi(TRandom *) // J/Psi composition return 443; } +Int_t AliGenMUONlib::IpPsiP(TRandom *) +{ +// Psi prime composition + return 100443; +} +Int_t AliGenMUONlib::IpJpsiFamily(TRandom *) +{ +// J/Psi composition + Int_t ip; + Float_t r = gRandom->Rndm(); + if (r < 0.98) { + ip = 443; + } else { + ip = 100443; + } + return ip; +} + + // Upsilon // // // pt-distribution //____________________________________________________________ -Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t *dummy ) +Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t */*dummy*/ ) { // Upsilon pT const Double_t kpt0 = 5.3; @@ -210,11 +374,69 @@ Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t *dummy ) Double_t pass1 = 1.+(x/kpt0)*(x/kpt0); return x/TMath::Power(pass1,kxn); } + +Double_t AliGenMUONlib::PtUpsilonPbPb( Double_t *px, Double_t */*dummy*/) +{ + +// +// Upsilon pT +// +// +// R. Vogt 2002 +// PbPb 5.5 TeV +// MRST HO +// mc = 1.4 GeV, pt-kick 1 GeV +// + Float_t x = px[0]; + Double_t c[8] = { + -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00, + -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05 + }; + Double_t y; + if (x < 10.) { + Int_t j; + y = c[j = 7]; + while (j > 0) y = y * x +c[--j]; + y = x * TMath::Exp(y); + } else { + y = 0.; + } + return y; +} + +Double_t AliGenMUONlib::PtUpsilonPP( Double_t *px, Double_t */*dummy*/) +{ + +// +// Upsilon pT +// +// +// R. Vogt 2002 +// pp 14 TeV +// MRST HO +// mc = 1.4 GeV, pt-kick 1 GeV +// + Float_t x = px[0]; + Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00, + -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06}; + + Double_t y; + if (x < 10.) { + Int_t j; + y = c[j = 7]; + while (j > 0) y = y * x +c[--j]; + y = x * TMath::Exp(y); + } else { + y = 0.; + } + return y; +} + // // y-distribution // //____________________________________________________________ -Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *dummy) +Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t */*dummy*/) { // Upsilon y const Double_t ky0 = 3.; @@ -228,6 +450,55 @@ Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *dummy) yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2); return yu; } + + +Double_t AliGenMUONlib::YUpsilonPbPb( Double_t *px, Double_t */*dummy*/) +{ + +// +// Upsilon y +// +// +// R. Vogt 2002 +// PbPb 5.5 TeV +// MRST HO +// mc = 1.4 GeV, pt-kick 1 GeV +// + + Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04, + -2.99753e-09, 1.28895e-05}; + + Double_t x = px[0]; + if (TMath::Abs(x) > 5.55) return 0.; + Int_t j; + Double_t y = c[j = 6]; + while (j > 0) y = y * x +c[--j]; + return y; +} + +Double_t AliGenMUONlib::YUpsilonPP( Double_t *px, Double_t */*dummy*/) +{ + +// +// Upsilon y +// +// +// R. Vogt 2002 +// p p 14. TeV +// MRST HO +// mc = 1.4 GeV, pt-kick 1 GeV +// + Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04, + -6.20539e-10, 1.29943e-05}; + + Double_t x = px[0]; + if (TMath::Abs(x) > 6.2) return 0.; + Int_t j; + Double_t y = c[j = 6]; + while (j > 0) y = y * x +c[--j]; + return y; +} + // particle composition // Int_t AliGenMUONlib::IpUpsilon(TRandom *) @@ -235,6 +506,32 @@ Int_t AliGenMUONlib::IpUpsilon(TRandom *) // y composition return 553; } +Int_t AliGenMUONlib::IpUpsilonP(TRandom *) +{ +// y composition + return 100553; +} +Int_t AliGenMUONlib::IpUpsilonPP(TRandom *) +{ +// y composition + return 200553; +} +Int_t AliGenMUONlib::IpUpsilonFamily(TRandom *) +{ +// y composition + Int_t ip; + Float_t r = gRandom->Rndm(); + + if (r < 0.712) { + ip = 553; + } else if (r < 0.896) { + ip = 100553; + } else { + ip = 200553; + } + return ip; +} + // // Phi @@ -242,13 +539,13 @@ Int_t AliGenMUONlib::IpUpsilon(TRandom *) // // pt-distribution (by scaling of pion distribution) //____________________________________________________________ -Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *dummy) +Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t */*dummy*/) { // Phi pT return PtScal(*px,7); } // y-distribution -Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *dummy) +Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t */*dummy*/) { // Phi y Double_t *dum=0; @@ -259,7 +556,60 @@ Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *dummy) Int_t AliGenMUONlib::IpPhi(TRandom *) { // Phi composition - return 41; + return 333; +} + +// +// omega +// +// +// pt-distribution (by scaling of pion distribution) +//____________________________________________________________ +Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t */*dummy*/) +{ +// Omega pT + return PtScal(*px,5); +} +// y-distribution +Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t */*dummy*/) +{ +// Omega y + Double_t *dum=0; + return YJpsi(px,dum); +} +// particle composition +// +Int_t AliGenMUONlib::IpOmega(TRandom *) +{ +// Omega composition + return 223; +} + + +// +// Eta +// +// +// pt-distribution (by scaling of pion distribution) +//____________________________________________________________ +Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t */*dummy*/) +{ +// Eta pT + return PtScal(*px,3); +} +// y-distribution +Double_t AliGenMUONlib::YEta( Double_t *px, Double_t */*dummy*/) +{ +// Eta y + Double_t *dum=0; + return YJpsi(px,dum); +} +// particle composition +// +Int_t AliGenMUONlib::IpEta(TRandom *) +{ +// Eta composition + return 221; } // @@ -268,18 +618,19 @@ Int_t AliGenMUONlib::IpPhi(TRandom *) // // pt-distribution //____________________________________________________________ -Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *dummy) +Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t */*dummy*/) { // Charm pT const Double_t kpt0 = 4.08; const Double_t kxn = 9.40; + Double_t x=*px; // - Double_t pass1 = 1.+(x/kpt0)*(x/kpt0); + Double_t pass1 = 1.+(x/kpt0); return x/TMath::Power(pass1,kxn); } // y-distribution -Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *dummy) +Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t */*dummy*/) { // Charm y Double_t *dum=0; @@ -314,7 +665,7 @@ Int_t AliGenMUONlib::IpCharm(TRandom *ran) // // pt-distribution //____________________________________________________________ -Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *dummy) +Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t */*dummy*/) { // Beauty pT const Double_t kpt0 = 4.; @@ -325,7 +676,7 @@ Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *dummy) return x/TMath::Power(pass1,kxn); } // y-distribution -Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *dummy) +Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t */*dummy*/) { // Beauty y Double_t *dum=0; @@ -353,21 +704,48 @@ Int_t AliGenMUONlib::IpBeauty(TRandom *ran) } typedef Double_t (*GenFunc) (Double_t*, Double_t*); -GenFunc AliGenMUONlib::GetPt(Int_t param, const char* tname) +GenFunc AliGenMUONlib::GetPt(Int_t param, const char* tname) const { // Return pointer to pT parameterisation + TString sname = TString(tname); GenFunc func; switch (param) { case kPhi: func=PtPhi; break; + case kOmega: + func=PtOmega; + break; + case kEta: + func=PtEta; + break; + case kJpsiFamily: + case kPsiP: case kJpsi: - func=PtJpsi; + if (sname == "Vogt" || sname == "Vogt PbPb") { + func=PtJpsiPbPb; + } else if (sname == "Vogt pp") { + func=PtJpsiPP; + } else { + func=PtJpsi; + } break; - case kUpsilon: - func=PtUpsilon; + case kJpsiFromB: + func = PtJpsiBPbPb; break; + case kUpsilonFamily: + case kUpsilonP: + case kUpsilonPP: + case kUpsilon: + if (sname == "Vogt" || sname == "Vogt PbPb") { + func=PtUpsilonPbPb; + } else if (sname == "Vogt pp") { + func=PtUpsilonPP; + } else { + func=PtUpsilon; + } + break; case kCharm: func=PtCharm; break; @@ -387,8 +765,10 @@ GenFunc AliGenMUONlib::GetPt(Int_t param, const char* tname) return func; } -GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) +GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const { + TString sname = TString(tname); + // Return pointer to y- parameterisation GenFunc func; switch (param) @@ -396,11 +776,37 @@ GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) case kPhi: func=YPhi; break; + case kEta: + func=YEta; + break; + case kOmega: + func=YOmega; + break; + case kJpsiFamily: + case kPsiP: case kJpsi: - func=YJpsi; + if (sname == "Vogt" || sname == "Vogt PbPb") { + func=YJpsiPbPb; + } else if (sname == "Vogt pp"){ + func=YJpsiPP; + } else { + func=YJpsi; + } + break; + case kJpsiFromB: + func = YJpsiBPbPb; break; + case kUpsilonFamily: + case kUpsilonP: + case kUpsilonPP: case kUpsilon: - func=YUpsilon; + if (sname == "Vogt" || sname == "Vogt PbPb") { + func=YUpsilonPbPb; + } else if (sname == "Vogt pp") { + func = YUpsilonPP; + } else { + func=YUpsilon; + } break; case kCharm: func=YCharm; @@ -421,7 +827,7 @@ GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) return func; } typedef Int_t (*GenFuncIp) (TRandom *); -GenFuncIp AliGenMUONlib::GetIp(Int_t param, const char* tname) +GenFuncIp AliGenMUONlib::GetIp(Int_t param, const char* /*tname*/) const { // Return pointer to particle type parameterisation GenFuncIp func; @@ -430,12 +836,34 @@ GenFuncIp AliGenMUONlib::GetIp(Int_t param, const char* tname) case kPhi: func=IpPhi; break; + case kEta: + func=IpEta; + break; + case kOmega: + func=IpOmega; + break; + case kJpsiFamily: + func=IpJpsiFamily; + break; + case kPsiP: + func=IpPsiP; + break; case kJpsi: + case kJpsiFromB: func=IpJpsi; break; case kUpsilon: func=IpUpsilon; break; + case kUpsilonFamily: + func=IpUpsilonFamily; + break; + case kUpsilonP: + func=IpUpsilonP; + break; + case kUpsilonPP: + func=IpUpsilonPP; + break; case kCharm: func=IpCharm; break; @@ -457,4 +885,54 @@ GenFuncIp AliGenMUONlib::GetIp(Int_t param, const char* tname) +Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0, + Float_t dx, + Int_t n, Int_t no) +{ +// +// Neville's alorithm for interpolation +// +// x: x-value +// y: Input array +// x0: minimum x +// dx: step size +// n: number of data points +// no: order of polynom +// + Float_t* c = new Float_t[n]; + Float_t* d = new Float_t[n]; + Int_t m, i; + for (i = 0; i < n; i++) { + c[i] = y[i]; + d[i] = y[i]; + } + + Int_t ns = int((x - x0)/dx); + + Float_t y1 = y[ns]; + ns--; + for (m = 0; m < no; m++) { + for (i = 0; i < n-m; i++) { + Float_t ho = x0 + Float_t(i) * dx - x; + Float_t hp = x0 + Float_t(i+m+1) * dx - x; + Float_t w = c[i+1] - d[i]; + Float_t den = ho-hp; + den = w/den; + d[i] = hp * den; + c[i] = ho * den; + } + Float_t dy; + + if (2*ns < (n-m-1)) { + dy = c[ns+1]; + } else { + dy = d[ns--]; + } + y1 += dy;} + delete[] c; + delete[] d; + + return y1; +} +