TEvtGen/EvtGenModels/EvtBtoXsgammaFermiUtil.cxx

   1 //--------------------------------------------------------------------------
   2 //
   3 //
   4 // Copyright Information: See EvtGen/COPYRIGHT
   5 //
   6 // Environment:
   7 //      This software is part of the EvtGen package developed jointly
   8 //      for the BaBar and CLEO collaborations.  If you use all or part
   9 //      of it, please give an appropriate acknowledgement.
  10 //
  11 // Module: EvtBtoXsgammaFermiUtil.cc
  12 //
  13 // Description:
  14 //      Class to hold various fermi functions and their helper functions. The
  15 //      fermi functions are used in EvtBtoXsgammaKagan.
  16 //
  17 // Modification history:
  18 //
  19 //      Jane Tinslay       March 21, 2001       Module created
  20 //------------------------------------------------------------------------
  21 #include "EvtGenBase/EvtPatches.hh"
  22
  23 //-----------------------
  24 // This Class's Header --
  25 //-----------------------
  26 #include "EvtGenModels/EvtBtoXsgammaFermiUtil.hh"
  27 #include "EvtGenModels/EvtItgTwoCoeffFcn.hh"
  28 #include "EvtGenModels/EvtBtoXsgammaRootFinder.hh"
  29 #include "EvtGenModels/EvtItgFunction.hh"
  30 #include "EvtGenBase/EvtConst.hh"
  31 #include "EvtGenBase/EvtReport.hh"
  32
  33 //---------------
  34 // C++ Headers --
  35 //---------------
  36 #include <iostream>
  37 #include <math.h>
  38 using std::endl;
  39
  40 double EvtBtoXsgammaFermiUtil::FermiExpFunc(double y, const std::vector<double> &coeffs) {
  41
  42   //coeffs: 1 = lambdabar, 2 = a, 3 = lam1, 4 = norm
  43   // report(INFO,"EvtGen")<<coeffs[4]<<endl;
  44   return (pow(1. - (y/coeffs[1]),coeffs[2])*exp((-3.*pow(coeffs[1],2.)/coeffs[3])*y/coeffs[1]))/coeffs[4];
  45
  46 }
  47
  48 double EvtBtoXsgammaFermiUtil::FermiGaussFunc(double y, const std::vector<double> &coeffs) {
  49
  50   //coeffs: 1 = lambdabar, 2 = a, 3 = c, 4 = norm
  51   return (pow(1. - (y/coeffs[1]),coeffs[2])*exp(-pow(coeffs[3],2.)*pow(1. - (y/coeffs[1]),2.)))/coeffs[4];
  52
  53 }
  54
  55 double EvtBtoXsgammaFermiUtil::FermiGaussFuncRoot(double lambdabar, double lam1, double mb, std::vector<double> &gammaCoeffs) {
  56
  57   std::vector<double> coeffs1(3);
  58   std::vector<double> coeffs2(3);
  59
  60   coeffs1[0]=0.2;
  61   coeffs1[1]=lambdabar;
  62   coeffs1[2]=0.0;
  63
  64   coeffs2[0]=0.2;
  65   coeffs2[1]=lambdabar;
  66   coeffs2[2]=-lam1/3.;
  67
  68   EvtItgTwoCoeffFcn *lhFunc = new EvtItgTwoCoeffFcn(&FermiGaussRootFcnA, -mb, lambdabar, coeffs1, gammaCoeffs);
  69   EvtItgTwoCoeffFcn *rhFunc = new EvtItgTwoCoeffFcn(&FermiGaussRootFcnB, -mb, lambdabar, coeffs2, gammaCoeffs);
  70
  71   EvtBtoXsgammaRootFinder *rootFinder = new EvtBtoXsgammaRootFinder();
  72
  73   double root = rootFinder->GetGaussIntegFcnRoot(lhFunc, rhFunc, 1.0e-4, 1.0e-4, 40, 40, -mb, lambdabar, 0.2, 0.4, 1.0e-6);
  74
  75   delete rootFinder; rootFinder=0;
  76   delete lhFunc; lhFunc=0;
  77   delete rhFunc; rhFunc=0;
  78   return root;
  79
  80 }
  81
  82 double EvtBtoXsgammaFermiUtil::FermiGaussRootFcnA(double y, const std::vector<double> &coeffs1, const std::vector<double> &coeffs2) {
  83
  84
  85   //coeffs1: 0=ap, 1=lambdabar, coeffs2=gamma function coeffs
  86   double cp = Gamma((2.0 + coeffs1[0])/2., coeffs2)/Gamma((1.0 + coeffs1[0])/2., coeffs2);
  87
  88   return (y*y)*pow((1. - (y/coeffs1[1])),coeffs1[0])*exp(-pow(cp,2)*pow((1.-(y/coeffs1[1])),2.));
  89
  90 }
  91
  92 double EvtBtoXsgammaFermiUtil::FermiGaussRootFcnB(double y, const std::vector<double> &coeffs1, const std::vector<double> &coeffs2) {
  93
  94   //coeffs1: 0=ap, 1=lambdabar, coeffs2=gamma function coeffs
  95   double cp = Gamma((2.0 + coeffs1[0])/2., coeffs2)/Gamma((1.0 + coeffs1[0])/2., coeffs2);
  96   return pow((1. - (y/coeffs1[1])),coeffs1[0])*exp(-pow(cp,2)*pow((1.-(y/coeffs1[1])),2.));
  97
  98 }
  99
 100 double EvtBtoXsgammaFermiUtil::Gamma(double z, const std::vector<double> &coeffs) {
 101
 102   //Lifted from Numerical Recipies in C
 103   double x, y, tmp, ser;
 104
 105   int j;
 106   y = z;
 107   x = z;
 108
 109   tmp = x + 5.5;
 110   tmp = tmp - (x+0.5)*log(tmp);
 111   ser=1.000000000190015;
 112
 113   for (j=0;j<6;j++) {
 114     y = y +1.0;
 115     ser = ser + coeffs[j]/y;
 116   }
 117
 118   return exp(-tmp+log(2.5066282746310005*ser/x));
 119
 120 }
 121
 122 double EvtBtoXsgammaFermiUtil::BesselK1(double x) {
 123
 124   //Lifted from Numerical Recipies in C : Returns the modified Bessel
 125   //function K_1(x) for positive real x
 126   if (x<0.0) report(INFO,"EvtGen") <<"x is negative !"<<endl;
 127
 128   double y, ans;
 129
 130   if (x <= 2.0) {
 131     y=x*x/4.0;
 132     ans = (log(x/2.0)*BesselI1(x))+(1.0/x)*(1.0+y*(0.15443144+y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1+y*(-0.110404e-2+y*(-0.4686e-4)))))));
 133   }
 134   else {
 135     y=2.0/x;
 136     ans=(exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619+y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2+y*(0.325614e-2+y*(-0.68245e-3)))))));
 137   }
 138   return ans;
 139
 140 }
 141
 142 double EvtBtoXsgammaFermiUtil::BesselI1(double x) {
 143    //Lifted from Numerical Recipies in C : Returns the modified Bessel
 144   //function I_1(x) for any real x
 145
 146   double ax, ans;
 147   double y;
 148
 149   ax=fabs(x);
 150   if ( ax  < 3.75) {
 151     y=x/3.75;
 152     y*=y;
 153     ans=ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934+y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
 154   }
 155   else {
 156     y=3.75/ax;
 157     ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1 -y*0.420059e-2));
 158     ans=0.398914228+y*(-0.3988024e-1+y*(-0.362018e-2+y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
 159     ans*=(exp(ax)/sqrt(ax));
 160   }
 161     return x < 0.0 ? -ans:ans;
 162 }
 163
 164 double EvtBtoXsgammaFermiUtil::FermiRomanFuncRoot(double lambdabar, double lam1) {
 165
 166   EvtItgFunction *lhFunc = new EvtItgFunction(&FermiRomanRootFcnA, -1.e-6, 1.e6);
 167
 168   EvtBtoXsgammaRootFinder *rootFinder = new EvtBtoXsgammaRootFinder();
 169   double rhSide = 1.0 - (lam1/(3.0*lambdabar*lambdabar));
 170
 171   double rho = rootFinder->GetRootSingleFunc(lhFunc, rhSide, 0.1, 0.4, 1.0e-6);
 172   //rho=0.250353;
 173   report(INFO,"EvtGen")<<"rho/2 "<<rho/2.<<" bessel "<<BesselK1(rho/2.)<<endl;
 174   double pF = lambdabar*sqrt(EvtConst::pi)/(rho*exp(rho/2.)*BesselK1(rho/2.));
 175   report(INFO,"EvtGen")<<"rho "<<rho<<" pf "<<pF<<endl;
 176
 177   delete lhFunc; lhFunc=0;
 178   delete rootFinder; rootFinder=0;
 179   return rho;
 180
 181 }
 182
 183 double EvtBtoXsgammaFermiUtil::FermiRomanRootFcnA(double y) {
 184
 185    return EvtConst::pi*(2. + y)*pow(y,-2.)*exp(-y)*pow(BesselK1(y/2.),-2.);
 186
 187 }
 188 double EvtBtoXsgammaFermiUtil::FermiRomanFunc(double y, const std::vector<double> &coeffs) {
 189   if (y == (coeffs[1]-coeffs[2])) y=0.99999999*(coeffs[1]-coeffs[2]);
 190
 191   //coeffs: 1 = mB, 2=mb, 3=rho, 4=lambdabar, 5=norm
 192   double pF = coeffs[4]*sqrt(EvtConst::pi)/(coeffs[3]*exp(coeffs[3]/2.)*BesselK1(coeffs[3]/2.));
 193   //  report(INFO,"EvtGen")<<" pf "<<y<<" "<<pF<<" "<<coeffs[1]<<" "<<coeffs[2]<<" "<<coeffs[3]<<" "<<coeffs[4]<<" "<<coeffs[5]<<endl;
 194   //double pF=0.382533;
 195
 196   //report(INFO,"EvtGen")<<(coeffs[1]-coeffs[2])*(1./(sqrt(EvtConst::pi)*pF))<<endl;
 197   //report(INFO,"EvtGen")<<(1.-y/(coeffs[1]-coeffs[2]))<<endl;
 198   //report(INFO,"EvtGen")<<(coeffs[1]-coeffs[2])<<endl;
 199   //report(INFO,"EvtGen")<<(coeffs[1]-coeffs[2])*(1.-y/(coeffs[1]-coeffs[2]))<<endl;
 200
 201   //report(INFO,"EvtGen")<<" "<<pF*coeffs[3]/((coeffs[1]-coeffs[2])*(1.-y/(coeffs[1]-coeffs[2])))<<endl;
 202   // report(INFO,"EvtGen")<<" "<<((coeffs[1]-coeffs[2])/pF)*(1. -y/(coeffs[1]-coeffs[2]))<<endl;
 203
 204   //report(INFO,"EvtGen")<<"result "<<(coeffs[1]-coeffs[2])*(1./(sqrt(EvtConst::pi)*pF))*exp(-(1./4.)*pow(pF*(coeffs[3]/((coeffs[1]-coeffs[2])*(1.-y/(coeffs[1]-coeffs[2])))) - ((coeffs[1]-coeffs[2])/pF)*(1. -y/(coeffs[1]-coeffs[2])),2.))/coeffs[5];
 205
 206   //report(INFO,"EvtGen")<<"leaving"<<endl;
 207   return (coeffs[1]-coeffs[2])*(1./(sqrt(EvtConst::pi)*pF))*exp(-(1./4.)*pow(pF*(coeffs[3]/((coeffs[1]-coeffs[2])*(1.-y/(coeffs[1]-coeffs[2])))) - ((coeffs[1]-coeffs[2])/pF)*(1. -y/(coeffs[1]-coeffs[2])),2.))/coeffs[5];
 208
 209
 210 }