Fix for definitions for CINT
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenModels / EvtVubNLO.hh
1 //--------------------------------------------------------------------------
2 //
3 // Environment:
4 //      This software is part of the EvtGen package developed jointly
5 //      for the BaBar and CLEO collaborations.  If you use all or part
6 //      of it, please give an appropriate acknowledgement.
7 //
8 // Copyright Information: See EvtGen/COPYRIGHT
9 //      Copyright (C) 1998      Caltech, UCSB
10 //
11 // Module: EvtGen/EvtVubNLO.hh
12 //
13 // Description:
14 // Class to generate inclusive B to X_u l nu decays according to various
15 // decay models. Implemtented are ACCM, parton-model and a QCD model.
16 //
17 // Modification history:
18 //
19 //   Sven Menke     January 17, 2001         Module created
20 //
21 //------------------------------------------------------------------------
22
23 #ifndef EVTVUBNLO_HH
24 #define EVTVUBNLO_HH
25
26 #include <vector>
27 #include "EvtGenBase/EvtDecayIncoherent.hh"
28
29 class EvtParticle;
30 class RandGeneral;
31
32 class EvtVubNLO:public  EvtDecayIncoherent  {
33
34 public:
35   
36   EvtVubNLO() {}
37   virtual ~EvtVubNLO();
38
39   std::string getName();
40
41   EvtDecayBase* clone();
42
43   void initProbMax();
44
45   void init();
46
47   void decay(EvtParticle *p); 
48
49
50 private:
51
52   // cache
53   double _lbar;
54   double _mupi2;
55
56   double _mb;     // the b-quark pole mass in GeV 
57   double _mB;
58   double _lambdaSF;
59   double _b;      // Parameter for the Fermi Motion 
60   double _kpar;  
61   double _mui; // renormalization scale (preferred value=1.5 GeV)
62   double _SFNorm; // SF normalization 
63   double _dGMax;  // max dGamma*p2 value;
64   int    _nbins;
65   int    _idSF;// which shape function?
66   double * _masses;
67   double * _weights;
68
69   double _gmax;
70   int _ngood,_ntot;
71
72
73   double tripleDiff(double pp, double pl, double pm);
74   double SFNorm(const std::vector<double> &coeffs);
75   static double integrand(double omega, const std::vector<double> &coeffs);
76   double F10(const std::vector<double> &coeffs);
77   static double F1Int(double omega,const std::vector<double> &coeffs);
78   double F20(const std::vector<double> &coeffs);
79   static double F2Int(double omega,const std::vector<double> &coeffs);
80   double F30(const std::vector<double> &coeffs);
81   static double F3Int(double omega,const std::vector<double> &coeffs);
82   static double g1(double y, double z);
83   static double g2(double y, double z);
84   static double g3(double y, double z);
85
86   static double Gamma(double z);// Euler Gamma Function
87   static double dgamma(double t, const std::vector<double> &c){  return pow(t,c[0]-1)*exp(-t);}
88   static double Gamma(double z, double tmax);
89
90   // theory parameters
91   inline double mu_i(){return _mui;} // intermediate scale
92   inline double mu_bar(){return _mui;} 
93   inline double mu_h(){return _mb/sqrt(2.0);} // high scale
94   inline double lambda1(){return -_mupi2;}
95   
96   // expansion coefficients for RGE
97   static double beta0(int nf=4){return 11.-2./3.*nf;}
98   static double beta1(int nf=4){return 34.*3.-38./3.*nf;}
99   static double beta2(int nf=4){return 1428.5-5033./18.*nf+325./54.*nf*nf;}
100   static double gamma0(){return 16./3.;}
101   static double gamma1(int nf=4){return 4./3.*(49.85498-40./9.*nf);}
102   static double gamma2(int nf=4){return 64./3.*(55.07242-8.58691*nf-nf*nf/27.);} /*  zeta3=1.20206 */
103   static double gammap0(){return -20./3.;}
104   static double gammap1(int nf=4){return -32./3.*(6.92653-0.9899*nf);} /* ??  zeta3=1.202 */
105
106
107   // running constants
108
109   static double alphas(double mu) ; 
110   static double C_F(double mu){return (4.0/3.0)*alphas(mu)/4./EvtConst::pi;}
111
112   // Shape Functions
113
114   inline double lambda_SF(){ return _lambdaSF;}
115   double lambda_bar(double omega0);
116   inline double lambda2(){return 0.12;}
117   double mu_pi2(double omega0);
118   inline double lambda(double){ return _mB-_mb;}
119
120   // specail for gaussian SF
121   static double cGaus(double b){return pow(Gamma(1+b/2.)/Gamma((1+b)/2.),2);}
122
123   double M0(double mui,double omega0);
124   static double shapeFunction(double omega, const std::vector<double> &coeffs);
125   static double expShapeFunction(double omega, const std::vector<double> &coeffs);
126   static double gausShapeFunction(double omega, const std::vector<double> &coeffs);
127   // SSF (not yet implemented)
128   double subS(const std::vector<double> &coeffs );
129   double subT(const std::vector<double> &coeffs);
130   double subU(const std::vector<double> &coeffs);
131   double subV(const std::vector<double> &coeffs);
132
133
134   // Sudakov
135
136   inline double S0(double a, double r){return -gamma0()/4/a/pow(beta0(),2)*(1/r-1+log(r));}
137   inline double S1(double /*a*/, double r){return gamma0()/4./pow(beta0(),2)*(
138                                                                   pow(log(r),2)*beta1()/2./beta0()+(gamma1()/gamma0()-beta1()/beta0())*(1.-r+log(r))
139                                                                   );}
140   inline double S2(double a, double r){return gamma0()*a/4./pow(beta0(),2)*(
141                                                                            -0.5*pow((1-r),2)*(
142                                                                                          pow(beta1()/beta0(),2)-beta2()/beta0()-beta1()/beta0()*gamma1()/gamma0()+gamma2()/gamma0()
143                                                                                          )
144                                                                            +(pow(beta1()/beta0(),2)-beta2()/beta0())*(1-r)*log(r)
145                                                                            +(beta1()/beta0()*gamma1()/gamma0()-beta2()/beta0())*(1-r+r*log(r))
146                                                                            );}
147   inline double dSudakovdepsi(double mu1, double mu2){return S2(alphas(mu1)/(4*EvtConst::pi),alphas(mu2)/alphas(mu1));}
148   inline double Sudakov(double mu1, double mu2, double epsi=0){double fp(4*EvtConst::pi);return S0(alphas(mu1)/fp,alphas(mu2)/alphas(mu1))+S1(alphas(mu1)/fp,alphas(mu2)/alphas(mu1))+epsi*dSudakovdepsi(mu1,mu2);}
149
150   // RG 
151   inline double dGdepsi(double mu1, double mu2){return 1./8./EvtConst::pi*(alphas(mu2)-alphas(mu1))*(gamma1()/beta0()-beta1()*gamma0()/pow(beta0(),2));}
152   inline double aGamma(double mu1, double mu2, double epsi=0){return gamma0()/2/beta0()*log(alphas(mu2)/alphas(mu1))+epsi*dGdepsi( mu1, mu2);}
153   inline double dgpdepsi(double mu1, double mu2){return 1./8./EvtConst::pi*(alphas(mu2)-alphas(mu1))*(gammap1()/beta0()-beta1()*gammap0()/pow(beta0(),2));}
154   inline double agammap(double mu1, double mu2, double epsi=0){return gammap0()/2/beta0()*log(alphas(mu2)/alphas(mu1))+epsi*dgpdepsi( mu1, mu2);}
155   inline double U1(double mu1, double mu2, double epsi=0){return exp(2*(Sudakov(mu1,mu2,epsi)-agammap(mu1,mu2,epsi)-aGamma(mu1,mu2,epsi)*log(_mb/mu1)));}
156   inline double U1lo(double mu1, double mu2){return U1(mu1,mu2);}
157   inline double U1nlo(double mu1, double mu2){return U1(mu1,mu2)*(1+2*(dSudakovdepsi(mu1,mu2)-dgpdepsi( mu1, mu2)-log(_mb/mu1)*dGdepsi( mu1, mu2)));}
158   inline double alo(double mu1, double mu2){return -2*aGamma(mu1,mu2);}
159   inline double anlo(double mu1, double mu2){return -2*dGdepsi(mu1,mu2);}
160
161 };
162
163 #endif
164