TEvtGen/EvtGenModels/EvtBtoXsgammaKagan.cxx

   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 // Module: EvtBtoXsgammaKagan.cc
   9 //
  10 // Description:
  11 //       Routine to perform two-body non-resonant B->Xs,gamma decays.
  12 //       The X_s mass spectrum generated is based on the Kagan-Neubert model.
  13 //       See hep-ph/9805303 for the model details and input parameters.
  14 //
  15 //       The input parameters are 1:fermi_model, 2:mB, 3:mb, 4:mu, 5:lam1,
  16 //       6:delta, 7:z, 8:nIntervalS, 9:nIntervalmH. Choosing fermi_model=1
  17 //       uses an exponential shape function, fermi_model=2 uses a gaussian
  18 //       shape function and fermi_model=3 a roman shape function. The complete mass
  19 //       spectrum for a given set of input parameters is calculated from
  20 //       scratch in bins of nIntervalmH. The s22, s27 and s28 coefficients are calculated
  21 //       in bins of nIntervalS. As the program includes lots of integration, the
  22 //       theoretical hadronic mass spectra is computed for the first time
  23 //       the init method is called. Then, all the other times (eg if we want to decay a B0
  24 //       as well as an anti-B0) the vector mass info stored the first time is used again.
  25 //
  26 // Modification history:
  27 //
  28 //      Jane Tinslay, Francesca Di Lodovico  March 21, 2001       Module created
  29 //------------------------------------------------------------------------
  30 //
  31 #include "EvtGenBase/EvtPatches.hh"
  32
  33 #include <stdlib.h>
  34 #include "EvtGenModels/EvtBtoXsgamma.hh"
  35 #include "EvtGenBase/EvtRandom.hh"
  36 #include "EvtGenModels/EvtBtoXsgammaKagan.hh"
  37 #include <string>
  38 #include "EvtGenBase/EvtConst.hh"
  39 #include "EvtGenBase/EvtParticle.hh"
  40 #include "EvtGenBase/EvtGenKine.hh"
  41 #include "EvtGenBase/EvtPDL.hh"
  42 #include "EvtGenBase/EvtReport.hh"
  43 #include "EvtGenModels/EvtItgSimpsonIntegrator.hh"
  44 #include "EvtGenModels/EvtItgFunction.hh"
  45 #include "EvtGenModels/EvtItgPtrFunction.hh"
  46 #include "EvtGenModels/EvtItgTwoCoeffFcn.hh"
  47 #include "EvtGenModels/EvtItgThreeCoeffFcn.hh"
  48 #include "EvtGenModels/EvtItgFourCoeffFcn.hh"
  49 #include "EvtGenModels/EvtItgAbsIntegrator.hh"
  50 #include "EvtGenModels/EvtBtoXsgammaFermiUtil.hh"
  51
  52 #include <fstream>
  53 using std::endl;
  54 using std::fstream;
  55
  56 bool EvtBtoXsgammaKagan::bbprod = false;
  57 double EvtBtoXsgammaKagan::intervalMH = 0;
  58
  59 EvtBtoXsgammaKagan::~EvtBtoXsgammaKagan(){
  60   delete [] massHad;
  61   delete [] brHad;
  62 }
  63
  64 void EvtBtoXsgammaKagan::init(int nArg, double* args){
  65
  66   if ((nArg) > 12 || (nArg > 1 && nArg <10) || nArg == 11){
  67
  68   report(ERROR,"EvtGen") << "EvtBtoXsgamma generator model "
  69                          << "EvtBtoXsgammaKagan expected "
  70                          << "either 1(default config) or "
  71                          << "10 (default mass range) or "
  72                          << "12 (user range) arguments but found: "
  73                          <<nArg<<endl;
  74   report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
  75     ::abort();
  76   }
  77
  78   if(nArg == 1){
  79     bbprod = true;
  80     getDefaultHadronicMass();
  81   }else{
  82     bbprod = false;
  83     computeHadronicMass(nArg, args);
  84   }
  85
  86   double mHminLimit=0.6373;
  87   double mHmaxLimit=4.5;
  88
  89   if (nArg>10){
  90     _mHmin = args[10];
  91     _mHmax = args[11];
  92     if (_mHmin > _mHmax){
  93       report(ERROR,"EvtGen") << "Minimum hadronic mass exceeds maximum "
  94                              << endl;
  95       report(ERROR,"EvtGen") << "Will terminate execution!" << endl;
  96       ::abort();
  97     }
  98     if (_mHmin < mHminLimit){
  99       report(ERROR,"EvtGen") << "Minimum hadronic mass below K pi threshold"
 100                              << endl;
 101       report(ERROR,"EvtGen") << "Resetting to K pi threshold" << endl;
 102       _mHmin = mHminLimit;
 103     }
 104     if (_mHmax > mHmaxLimit){
 105       report(ERROR,"EvtGen") << "Maximum hadronic mass above 4.5 GeV/c^2"
 106                              << endl;
 107       report(ERROR,"EvtGen") << "Resetting to 4.5 GeV/c^2" << endl;
 108       _mHmax = mHmaxLimit;
 109     }
 110   }else{
 111     _mHmin=mHminLimit; //  usually just above K pi threshold for Xsd/u
 112     _mHmax=mHmaxLimit;
 113   }
 114
 115 }
 116
 117 void EvtBtoXsgammaKagan::getDefaultHadronicMass(){
 118
 119     massHad = new double[81];
 120     brHad = new double[81];
 121
 122     double mass[81] = { 0, 0.0625995, 0.125199, 0.187798, 0.250398, 0.312997, 0.375597, 0.438196, 0.500796, 0.563395, 0.625995, 0.688594, 0.751194, 0.813793, 0.876392, 0.938992, 1.00159, 1.06419, 1.12679, 1.18939, 1.25199, 1.31459, 1.37719, 1.43979, 1.50239, 1.56499, 1.62759, 1.69019, 1.75278, 1.81538, 1.87798, 1.94058, 2.00318, 2.06578, 2.12838, 2.19098, 2.25358, 2.31618, 2.37878, 2.44138, 2.50398, 2.56658, 2.62918, 2.69178, 2.75438, 2.81698, 2.87958, 2.94217, 3.00477, 3.06737, 3.12997, 3.19257, 3.25517, 3.31777, 3.38037, 3.44297, 3.50557, 3.56817, 3.63077, 3.69337, 3.75597, 3.81857, 3.88117, 3.94377, 4.00637, 4.06896, 4.13156, 4.19416, 4.25676, 4.31936, 4.38196, 4.44456, 4.50716, 4.56976, 4.63236, 4.69496, 4.75756, 4.82016, 4.88276, 4.94536, 5.00796};
 123
 124     double br[81] = { 0, 1.03244e-09, 3.0239e-08, 1.99815e-07, 7.29392e-07, 1.93129e-06, 4.17806e-06, 7.86021e-06, 1.33421e-05, 2.09196e-05, 3.07815e-05, 4.29854e-05, 5.74406e-05, 7.3906e-05, 9.2003e-05, 0.000111223, 0.000130977, 0.000150618, 0.000169483, 0.000186934, 0.000202392, 0.000215366, 0.000225491, 0.000232496, 0.000236274, 0.000236835, 0.000234313, 0.000228942, 0.000221042, 0.000210994, 0.000199215, 0.000186137, 0.000172194, 0.000157775, 0.000143255, 0.000128952, 0.000115133, 0.000102012, 8.97451e-05, 7.84384e-05, 6.81519e-05, 5.89048e-05, 5.06851e-05, 4.34515e-05, 3.71506e-05, 3.1702e-05, 2.70124e-05, 2.30588e-05, 1.96951e-05, 1.68596e-05, 1.44909e-05, 1.25102e-05, 1.08596e-05, 9.48476e-06, 8.34013e-06, 7.38477e-06, 6.58627e-06, 5.91541e-06, 5.35022e-06, 4.87047e-06, 4.46249e-06, 4.11032e-06, 3.80543e-06, 3.54051e-06, 3.30967e-06, 3.10848e-06, 2.93254e-06, 2.78369e-06, 2.65823e-06, 2.55747e-06, 2.51068e-06, 2.57179e-06, 2.74684e-06, 3.02719e-06, 3.41182e-06, 3.91387e-06, 4.56248e-06, 5.40862e-06, 6.53915e-06, 8.10867e-06, 1.04167e-05 };
 125
 126   for(int i=0; i<81; i++){
 127     massHad[i] = mass[i];
 128     brHad[i] = br[i];
 129   }
 130   intervalMH=80;
 131 }
 132
 133 void EvtBtoXsgammaKagan::computeHadronicMass(int /*nArg*/, double* args){
 134
 135   //Input parameters
 136   int fermiFunction = (int)args[1];
 137   _mB = args[2];
 138   _mb = args[3];
 139   _mu = args[4];
 140   _lam1 = args[5];
 141   _delta = args[6];
 142   _z = args[7];
 143   _nIntervalS = args[8];
 144   _nIntervalmH = args[9];
 145   std::vector<double> mHVect(int(_nIntervalmH+1.0));
 146   massHad = new double[int(_nIntervalmH+1.0)];
 147   brHad = new double[int(_nIntervalmH+1.0)];
 148   intervalMH=_nIntervalmH;
 149
 150   //Going to have to add a new entry into the data file - takes ages...
 151   report(WARNING,"EvtGen") << "EvtBtoXsgammaKagan: calculating new hadronic mass spectra. This takes a while..." << endl;
 152
 153   //Now need to compute the mHVect vector for
 154   //the current parameters
 155
 156   //A few more parameters
 157   double _mubar = _mu;
 158   _mW = 80.33;
 159   _mt = 175.0;
 160   _alpha = 1./137.036;
 161   _lambdabar = _mB - _mb;
 162   _kappabar = 3.382 - 4.14*(sqrt(_z) - 0.29);
 163   _fz=Fz(_z);
 164   _rer8 = (44./9.) - (8./27.)*pow(EvtConst::pi,2.);
 165   _r7 = (-10./3.) - (8./9.)*pow(EvtConst::pi,2.);
 166   _rer2 = -4.092 + 12.78*(sqrt(_z) -.29);
 167   _gam77 = 32./3.;
 168   _gam27 = 416./81.;
 169   _gam87 = -32./9.;
 170   _lam2 = .12;
 171   _beta0 = 23./3.;
 172   _beta1 = 116./3.;
 173   _alphasmZ = .118;
 174   _mZ = 91.187;
 175   _ms = _mb/50.;
 176
 177   double eGammaMin = 0.5*_mB*(1. - _delta);
 178   double eGammaMax = 0.5*_mB;
 179   double yMin = 2.*eGammaMin/_mB;
 180   double yMax = 2.*eGammaMax/_mB;
 181   double _CKMrat= 0.976;
 182   double Nsl = 1.0;
 183
 184   //Calculate alpha the various scales
 185   _alphasmW = CalcAlphaS(_mW);
 186   _alphasmt = CalcAlphaS(_mt);
 187   _alphasmu = CalcAlphaS(_mu);
 188   _alphasmubar = CalcAlphaS(_mubar);
 189
 190   //Calculate the Wilson Coefficients and Delta
 191   _etamu = _alphasmW/_alphasmu;
 192   _kSLemmu = (12./23.)*((1./_etamu) -1.);
 193   CalcWilsonCoeffs();
 194   CalcDelta();
 195
 196   //Build s22 and s27 vector - saves time because double
 197   //integration is required otherwise
 198   std::vector<double> s22Coeffs(int(_nIntervalS+1.0));
 199   std::vector<double> s27Coeffs(int(_nIntervalS+1.0));
 200   std::vector<double> s28Coeffs(int(_nIntervalS+1.0));
 201
 202   double dy = (yMax - yMin)/_nIntervalS;
 203   double yp = yMin;
 204
 205   std::vector<double> sCoeffs(1);
 206   sCoeffs[0] = _z;
 207
 208   //Define s22 and s27 functions
 209   EvtItgPtrFunction *mys22Func = new EvtItgPtrFunction(&s22Func, 0., yMax+0.1, sCoeffs);
 210   EvtItgPtrFunction *mys27Func = new EvtItgPtrFunction(&s27Func, 0., yMax+0.1, sCoeffs);
 211
 212   //Use a simpson integrator
 213   EvtItgAbsIntegrator *mys22Simp = new EvtItgSimpsonIntegrator(*mys22Func, 1.0e-4, 20);
 214   EvtItgAbsIntegrator *mys27Simp = new EvtItgSimpsonIntegrator(*mys27Func, 1.0e-4, 50);
 215
 216   int i;
 217
 218   for (i=0;i<int(_nIntervalS+1.0);i++) {
 219
 220     s22Coeffs[i] = (16./27.)*mys22Simp->evaluate(1.0e-20,yp);
 221     s27Coeffs[i] = (-8./9.)*_z*mys27Simp->evaluate(1.0e-20,yp);
 222     s28Coeffs[i] = -s27Coeffs[i]/3.;
 223     yp = yp + dy;
 224
 225   }
 226
 227   delete mys22Func;
 228   delete mys27Func;
 229   delete mys22Simp;
 230   delete mys27Simp;
 231
 232   //Define functions and vectors used to calculate mHVect. Each function takes a set
 233   //of vectors which are used as the function coefficients
 234   std::vector<double> FermiCoeffs(6);
 235   std::vector<double> varCoeffs(3);
 236   std::vector<double> DeltaCoeffs(1);
 237   std::vector<double> s88Coeffs(2);
 238   std::vector<double> sInitCoeffs(3);
 239
 240   varCoeffs[0] = _mB;
 241   varCoeffs[1] = _mb;
 242   varCoeffs[2] = 0.;
 243
 244   DeltaCoeffs[0] = _alphasmu;
 245
 246   s88Coeffs[0] = _mb;
 247   s88Coeffs[1] = _ms;
 248
 249   sInitCoeffs[0] = _nIntervalS;
 250   sInitCoeffs[1] = yMin;
 251   sInitCoeffs[2] = yMax;
 252
 253   FermiCoeffs[0]=fermiFunction;
 254   FermiCoeffs[1]=0.0;
 255   FermiCoeffs[2]=0.0;
 256   FermiCoeffs[3]=0.0;
 257   FermiCoeffs[4]=0.0;
 258   FermiCoeffs[5]=0.0;
 259
 260   //Coefficients for gamma function
 261   std::vector<double> gammaCoeffs(6);
 262   gammaCoeffs[0]=76.18009172947146;
 263   gammaCoeffs[1]=-86.50532032941677;
 264   gammaCoeffs[2]=24.01409824083091;
 265   gammaCoeffs[3]=-1.231739572450155;
 266   gammaCoeffs[4]=0.1208650973866179e-2;
 267   gammaCoeffs[5]=-0.5395239384953e-5;
 268
 269   //Calculate quantities for the fermi function to be used
 270   //Distinguish among the different shape functions
 271   if (fermiFunction == 1) {
 272
 273     FermiCoeffs[1]=_lambdabar;
 274     FermiCoeffs[2]=(-3.*pow(_lambdabar,2.)/_lam1) - 1.;
 275     FermiCoeffs[3]=_lam1;
 276     FermiCoeffs[4]=1.0;
 277
 278     EvtItgPtrFunction *myNormFunc = new EvtItgPtrFunction(&EvtBtoXsgammaFermiUtil::FermiExpFunc, -_mb, _mB-_mb, FermiCoeffs);
 279     EvtItgAbsIntegrator *myNormSimp = new EvtItgSimpsonIntegrator(*myNormFunc, 1.0e-4, 40);
 280     FermiCoeffs[4]=myNormSimp->normalisation();
 281     delete myNormFunc; myNormFunc=0;
 282     delete myNormSimp; myNormSimp=0;
 283
 284   } else if (fermiFunction == 2) {
 285
 286     double a = EvtBtoXsgammaFermiUtil::FermiGaussFuncRoot(_lambdabar, _lam1, _mb, gammaCoeffs);
 287     FermiCoeffs[1]=_lambdabar;
 288     FermiCoeffs[2]=a;
 289     FermiCoeffs[3]= EvtBtoXsgammaFermiUtil::Gamma((2.0 + a)/2., gammaCoeffs)/
 290       EvtBtoXsgammaFermiUtil::Gamma((1.0 + a)/2., gammaCoeffs);
 291     FermiCoeffs[4]=1.0;
 292
 293     EvtItgPtrFunction *myNormFunc = new EvtItgPtrFunction(&EvtBtoXsgammaFermiUtil::FermiGaussFunc, -_mb, _mB-_mb, FermiCoeffs);
 294     EvtItgAbsIntegrator *myNormSimp = new EvtItgSimpsonIntegrator(*myNormFunc, 1.0e-4, 40);
 295     FermiCoeffs[4]=myNormSimp->normalisation();
 296     delete myNormFunc; myNormFunc=0;
 297     delete myNormSimp; myNormSimp=0;
 298
 299   }
 300   else if (fermiFunction == 3) {
 301
 302     double rho = EvtBtoXsgammaFermiUtil::FermiRomanFuncRoot(_lambdabar, _lam1);
 303     FermiCoeffs[1]=_mB;
 304     FermiCoeffs[2]=_mb;
 305     FermiCoeffs[3]= rho;
 306     FermiCoeffs[4]=_lambdabar;
 307     FermiCoeffs[5]=1.0;
 308
 309     EvtItgPtrFunction *myNormFunc = new EvtItgPtrFunction(&EvtBtoXsgammaFermiUtil::FermiRomanFunc, -_mb, _mB-_mb, FermiCoeffs);
 310     EvtItgAbsIntegrator *myNormSimp = new EvtItgSimpsonIntegrator(*myNormFunc, 1.0e-4, 40);
 311     FermiCoeffs[5]=myNormSimp->normalisation();
 312     delete myNormFunc; myNormFunc=0;
 313     delete myNormSimp; myNormSimp=0;
 314
 315   }
 316
 317   //Define functions
 318   EvtItgThreeCoeffFcn* myDeltaFermiFunc = new EvtItgThreeCoeffFcn(&DeltaFermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, DeltaCoeffs);
 319   EvtItgThreeCoeffFcn* mys88FermiFunc = new EvtItgThreeCoeffFcn(&s88FermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, s88Coeffs);
 320   EvtItgTwoCoeffFcn* mys77FermiFunc = new EvtItgTwoCoeffFcn(&s77FermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs);
 321   EvtItgTwoCoeffFcn* mys78FermiFunc = new EvtItgTwoCoeffFcn(&s78FermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs);
 322   EvtItgFourCoeffFcn* mys22FermiFunc = new EvtItgFourCoeffFcn(&sFermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, sInitCoeffs, s22Coeffs);
 323   EvtItgFourCoeffFcn* mys27FermiFunc = new EvtItgFourCoeffFcn(&sFermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, sInitCoeffs, s27Coeffs);
 324   EvtItgFourCoeffFcn* mys28FermiFunc = new EvtItgFourCoeffFcn(&sFermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, sInitCoeffs, s28Coeffs);
 325
 326   //Define integrators
 327   EvtItgSimpsonIntegrator* myDeltaFermiSimp =
 328     new EvtItgSimpsonIntegrator(*myDeltaFermiFunc, 1.0e-4, 40);
 329   EvtItgSimpsonIntegrator* mys77FermiSimp =
 330     new EvtItgSimpsonIntegrator(*mys77FermiFunc, 1.0e-4, 40);
 331   EvtItgSimpsonIntegrator* mys88FermiSimp =
 332     new EvtItgSimpsonIntegrator(*mys88FermiFunc, 1.0e-4, 40);
 333   EvtItgSimpsonIntegrator* mys78FermiSimp =
 334     new EvtItgSimpsonIntegrator(*mys78FermiFunc, 1.0e-4, 40);
 335   EvtItgSimpsonIntegrator* mys22FermiSimp =
 336     new EvtItgSimpsonIntegrator(*mys22FermiFunc, 1.0e-4, 40);
 337   EvtItgSimpsonIntegrator* mys27FermiSimp =
 338     new EvtItgSimpsonIntegrator(*mys27FermiFunc, 1.0e-4, 40);
 339   EvtItgSimpsonIntegrator* mys28FermiSimp =
 340     new EvtItgSimpsonIntegrator(*mys28FermiFunc, 1.0e-4, 40);
 341
 342   //Finally calculate mHVect for the range of hadronic masses
 343   double mHmin = sqrt(_mB*_mB - 2.*_mB*eGammaMax);
 344   double mHmax = sqrt(_mB*_mB - 2.*_mB*eGammaMin);
 345   double dmH = (mHmax - mHmin)/_nIntervalmH;
 346
 347   double mH=mHmin;
 348
 349   //Calculating the Branching Fractions
 350   for (i=0;i<int(_nIntervalmH+1.0);i++) {
 351
 352     double ymH = 1. - ((mH*mH)/(_mB*_mB));
 353
 354     //Need to set ymH as one of the input parameters
 355     myDeltaFermiFunc->setCoeff(2, 2, ymH);
 356     mys77FermiFunc->setCoeff(2, 2, ymH);
 357     mys88FermiFunc->setCoeff(2, 2, ymH);
 358     mys78FermiFunc->setCoeff(2, 2, ymH);
 359     mys22FermiFunc->setCoeff(2, 2, ymH);
 360     mys27FermiFunc->setCoeff(2, 2, ymH);
 361     mys28FermiFunc->setCoeff(2, 2, ymH);
 362
 363     //Integrate
 364
 365     double deltaResult = myDeltaFermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
 366     double s77Result = mys77FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
 367     double s88Result = mys88FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
 368     double s78Result = mys78FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
 369     double s22Result = mys22FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
 370     double s27Result = mys27FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
 371     mys28FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
 372
 373     double py = (pow(_CKMrat,2.)*(6./_fz)*(_alpha/EvtConst::pi)*(deltaResult*_cDeltatot  + (_alphasmu/EvtConst::pi)*(s77Result*pow(_c70mu,2.) + s27Result*_c2mu*(_c70mu  - _c80mu/3.) + s78Result*_c70mu*_c80mu + s22Result*_c2mu*_c2mu  + s88Result*_c80mu*_c80mu )  ) );
 374
 375     mHVect[i] = 2.*(mH/(_mB*_mB))*0.105*Nsl*py;
 376
 377     massHad[i] = mH;
 378     brHad[i] =  2.*(mH/(_mB*_mB))*0.105*Nsl*py;
 379
 380     mH = mH+dmH;
 381
 382   }
 383
 384   //Clean up
 385   delete  myDeltaFermiFunc; myDeltaFermiFunc=0;
 386   delete mys88FermiFunc; mys88FermiFunc=0;
 387   delete mys77FermiFunc; mys77FermiFunc=0;
 388   delete mys78FermiFunc; mys78FermiFunc=0;
 389   delete mys22FermiFunc; mys22FermiFunc=0;
 390   delete mys27FermiFunc; mys27FermiFunc=0;
 391   delete mys28FermiFunc; mys28FermiFunc=0;
 392
 393   delete myDeltaFermiSimp; myDeltaFermiSimp=0;
 394   delete mys77FermiSimp; mys77FermiSimp=0;
 395   delete mys88FermiSimp; mys88FermiSimp=0;
 396   delete mys78FermiSimp; mys78FermiSimp=0;
 397   delete mys22FermiSimp; mys22FermiSimp=0;
 398   delete mys27FermiSimp; mys27FermiSimp=0;
 399   delete mys28FermiSimp; mys28FermiSimp=0;
 400
 401 }
 402
 403 double EvtBtoXsgammaKagan::GetMass( int /*Xscode*/ ){
 404
 405 //  Get hadronic mass for the event according to the hadronic mass spectra computed in computeHadronicMass
 406   double mass=0.0;
 407   double min=_mHmin;
 408   if(bbprod)min=1.1;
 409   //  double max=4.5;
 410   double max=_mHmax;
 411   double xbox(0), ybox(0);
 412   double boxheight(0);
 413   double trueHeight(0);
 414   double boxwidth=max-min;
 415   double wgt(0.);
 416
 417   for (int i=0;i<int(intervalMH+1.0);i++) {
 418     if(brHad[i]>boxheight)boxheight=brHad[i];
 419   }
 420   while ((mass > max) || (mass < min)){
 421     xbox = EvtRandom::Flat(boxwidth)+min;
 422     ybox=EvtRandom::Flat(boxheight);
 423     trueHeight=0.0;
 424     // Correction by Peter Richardson
 425     for( int i = 1 ; i < int( intervalMH + 1.0 ) ; ++i ) {
 426       if ( ( massHad[i] >= xbox ) && ( 0.0 == trueHeight ) ) {
 427         wgt=(xbox-massHad[i-1])/(massHad[i]-massHad[i-1]);
 428         trueHeight=brHad[i-1]+wgt*(brHad[i]-brHad[i-1]);
 429       }
 430     }
 431
 432     if (ybox>trueHeight) {
 433       mass=0.0;
 434     } else {
 435       mass=xbox;
 436     }
 437   }
 438
 439   return mass;
 440 }
 441
 442 double EvtBtoXsgammaKagan::CalcAlphaS(double scale) {
 443
 444   double v = 1. -_beta0*(_alphasmZ/(2.*EvtConst::pi))*(log(_mZ/scale));
 445   return (_alphasmZ/v)*(1. - ((_beta1/_beta0)*(_alphasmZ/(4.*EvtConst::pi))*(log(v)/v)));
 446
 447 }
 448
 449 void EvtBtoXsgammaKagan::CalcWilsonCoeffs( ){
 450
 451    double mtatmw=_mt*pow((_alphasmW/_alphasmt),(12./23.))*(1 + (12./23.)*((253./18.) - (116./23.))*((_alphasmW - _alphasmt)/(4.0*EvtConst::pi)) - (4./3.)*(_alphasmt/EvtConst::pi));
 452   double xt=pow(mtatmw,2.)/pow(_mW,2.);
 453
 454
 455
 456   /////LO
 457   _c2mu = .5*pow(_etamu,(-12./23.)) + .5*pow(_etamu,(6./23.));
 458
 459   double c7mWsm = ((3.*pow(xt,3.) - 2.*pow(xt,2.))/(4.*pow((xt - 1.),4.)))*log(xt)
 460     + ((-8.*pow(xt,3.) - 5.*pow(xt,2.) + 7.*xt)/(24.*pow((xt - 1.),3.) )) ;
 461
 462   double c8mWsm =  ((-3.*pow(xt,2.))/(4.*pow((xt - 1.),4.)))*log(xt)
 463     + ((- pow(xt,3.) + 5.*pow(xt,2.) + 2.*xt)/(8.*pow((xt - 1.),3.)));
 464
 465   double c7constmu = (626126./272277.)*pow(_etamu,(14./23.))
 466     - (56281./51730.)*pow(_etamu,(16./23.)) - (3./7.)*pow(_etamu,(6./23.))
 467     - (1./14.)*pow(_etamu,(-12./23.)) - .6494*pow(_etamu,.4086) - .038*pow(_etamu,-.423)
 468     - .0186*pow(_etamu,-.8994) - .0057*pow(_etamu,.1456);
 469
 470   _c70mu = c7mWsm*pow(_etamu,(16./23.)) + (8./3.)*(pow(_etamu,(14./23.))
 471     -pow(_etamu,(16./23.)))*c8mWsm + c7constmu;
 472
 473   double c8constmu =  (313063./363036.)*pow(_etamu,(14./23.))
 474     -.9135*pow(_etamu,.4086) + .0873*pow(_etamu,-.423) - .0571*pow(_etamu,-.8994)
 475     + .0209*pow(_etamu,.1456);
 476
 477   _c80mu = c8mWsm*pow(_etamu,(14./23.)) + c8constmu;
 478
 479  //Compute the dilogarithm (PolyLog(2,x)) with the Simpson integrator
 480  //The dilogarithm is defined as: Li_2(x)=Int_0^x(-log(1.-z)/z)
 481  //however, Mathematica implements it as  Sum[z^k/k^2,{k,1,Infinity}], so, althought the two
 482  //results are similar and both implemented in the program, we prefer to use the
 483  //one closer to the Mathematica implementation as identical to what used by the theorists.
 484
 485  // EvtItgFunction *myDiLogFunc = new EvtItgFunction(&diLogFunc, 0., 1.-1./xt);
 486  //EvtItgAbsIntegrator *myDiLogSimp = new EvtItgSimpsonIntegrator(*myDiLogFunc, 1.0e-4, 50);
 487  //double li2 = myDiLogSimp->evaluate(1.0e-20,1.-1./xt);
 488
 489  double li2=diLogMathematica(1.-1./xt);
 490
 491 double c7mWsm1 = ( (-16. *pow(xt,4.) -122. *pow(xt,3.) + 80. *pow(xt,2.) -8. *xt)/
 492 (9. *pow((xt -1.),4.)) * li2 +
 493 (6. *pow(xt,4.) + 46. *pow(xt,3.) -28. *pow(xt,2.))/(3. *pow((xt-1.),5.)) *pow(log(xt),2.)
 494 + (-102. *pow(xt,5.) -588. *pow(xt,4.) -2262. *pow(xt,3.) + 3244. *pow(xt,2.) -1364. *xt
 495 + 208.)/(81. *pow((xt-1),5.)) *log(xt)
 496 + (1646. *pow(xt,4.) + 12205. *pow(xt,3.) -10740. *pow(xt,2.) + 2509. *xt -436.)/
 497 (486. *pow((xt-1),4.)) );
 498
 499 double c8mWsm1 = ((-4. *pow(xt,4.) + 40. *pow(xt,3.) + 41. *pow(xt,2.) + xt)/
 500 (6. *pow((xt-1.),4.))  * li2
 501 + (-17. *pow(xt,3.) -31. *pow(xt,2.))/(2. *pow((xt-1.),5.) ) *pow(log(xt),2.)
 502 + (-210. *pow(xt,5.) + 1086. *pow(xt,4.) + 4893. *pow(xt,3.) + 2857. *pow(xt,2.)
 503 -1994. *xt + 280.)/(216. *pow((xt-1),5.)) *log(xt)
 504 + (737. *pow(xt,4.) -14102. *pow(xt,3.) -28209. *pow(xt,2.) + 610. *xt -508.)/
 505 (1296. *pow((xt-1),4.)) );
 506
 507 double E1 = (xt *(18. -11. *xt -pow(xt,2.))/(12.*pow( (1. -xt),3.))
 508 + pow(xt,2.)* (15. -16. *xt + 4. *pow(xt,2.))/(6. *pow((1. -xt),4.)) *log(xt)
 509 -2./3. *log(xt) );
 510
 511 double e1 = 4661194./816831.;
 512 double e2 = -8516./2217. ;
 513 double e3 = 0.;
 514 double e4 = 0.;
 515 double e5 = -1.9043;
 516 double e6 = -.1008;
 517 double e7 = .1216;
 518 double e8 = .0183;
 519
 520 double f1 = -17.3023;
 521 double f2 = 8.5027;
 522 double f3 = 4.5508;
 523 double f4 = .7519;
 524 double f5 =  2.004;
 525 double f6 = .7476;
 526 double f7 = -.5385;
 527 double f8 = .0914;
 528
 529 double g1 = 14.8088;
 530 double g2 = -10.809;
 531 double g3 = -.874;
 532 double g4 = .4218;
 533 double g5 = -2.9347;
 534 double g6 = .3971;
 535 double g7 = .1600;
 536 double g8 = .0225;
 537
 538
 539 double c71constmu  = ((e1 *_etamu *E1 + f1 + g1 *_etamu) *pow(_etamu,(14./23.))
 540 + (e2 *_etamu *E1 + f2 + g2 *_etamu) *pow(_etamu,(16./23.))
 541 + (e3 *_etamu *E1 + f3 + g3 *_etamu) *pow(_etamu,(6./23.))
 542 + (e4 *_etamu *E1 + f4 + g4 *_etamu) *pow(_etamu,(-12./23.))
 543 + (e5 *_etamu *E1 + f5 + g5 *_etamu) *pow(_etamu,.4086)
 544 + (e6 *_etamu *E1 + f6 + g6 *_etamu) *pow(_etamu,(-.423))
 545 + (e7 *_etamu *E1 + f7 + g7 *_etamu) *pow(_etamu,(-.8994))
 546 + (e8 *_etamu *E1 + f8 + g8 *_etamu) *pow(_etamu,.1456 ));
 547
 548 double c71pmu = ( ((297664./14283. *pow(_etamu,(16./23.))
 549 -7164416./357075. *pow(_etamu,(14./23.))
 550 + 256868./14283. *pow(_etamu,(37./23.)) - 6698884./357075. *pow(_etamu,(39./23.)))
 551 *(c8mWsm))
 552 + 37208./4761. *(pow(_etamu,(39./23.)) - pow(_etamu,(16./23.))) *(c7mWsm)
 553 + c71constmu );
 554
 555 _c71mu = (_alphasmW/_alphasmu *(pow(_etamu,(16./23.))* c7mWsm1 + 8./3. *(pow(_etamu,(14./23.))
 556 - pow(_etamu,(16./23.)) ) *c8mWsm1 ) + c71pmu);
 557
 558 _c7emmu = ((32./75. *pow(_etamu,(-9./23.)) - 40./69. *pow(_etamu,(-7./23.)) +
 559                88./575. *pow(_etamu,(16./23.))) *c7mWsm + (-32./575. *pow(_etamu,(-9./23.)) +
 560                32./1449. *pow(_etamu,(-7./23.)) + 640./1449.*pow(_etamu,(14./23.)) -
 561                704./1725.*pow(_etamu,(16./23.)) ) *c8mWsm
 562          - 190./8073.*pow(_etamu,(-35./23.))  - 359./3105. *pow(_etamu,(-17./23.)) +
 563          4276./121095. *pow(_etamu,(-12./23.)) + 350531./1009125.*pow(_etamu,(-9./23.))
 564          + 2./4347. *pow(_etamu,(-7./23.)) - 5956./15525. *pow(_etamu,(6./23.)) +
 565          38380./169533. *pow(_etamu,(14./23.))   - 748./8625. *pow(_etamu,(16./23.)));
 566
 567 // Wilson coefficients values as according to Kagan's program
 568 // _c2mu=1.10566;
 569 //_c70mu=-0.314292;
 570 // _c80mu=-0.148954;
 571 // _c71mu=0.480964;
 572 // _c7emmu=0.0323219;
 573
 574 }
 575
 576 void EvtBtoXsgammaKagan::CalcDelta() {
 577
 578   double cDelta77 = (1. + (_alphasmu/(2.*EvtConst::pi)) *(_r7 - (16./3.) + _gam77*log(_mb/_mu)) + ( (pow((1. - _z),4.)/_fz) - 1.)*(6.*_lam2/pow(_mb,2.)) + (_alphasmubar/(2.*EvtConst::pi))*_kappabar )*pow(_c70mu,2.);
 579
 580   double cDelta27 = ((_alphasmu/(2.*EvtConst::pi))*(_rer2 + _gam27*log(_mb/_mu)) - (_lam2/(9.*_z*pow(_mb,2.))))*_c2mu*_c70mu;
 581
 582   double cDelta78 = (_alphasmu/(2.*EvtConst::pi))*(_rer8 + _gam87*log(_mb/_mu))*_c70mu*_c80mu;
 583
 584   _cDeltatot = cDelta77  + cDelta27 + cDelta78 + (_alphasmu/(2.*EvtConst::pi))*_c71mu*_c70mu + (_alpha/_alphasmu)*(2.*_c7emmu*_c70mu - _kSLemmu*pow(_c70mu,2.));
 585
 586 }
 587
 588 double EvtBtoXsgammaKagan::Delta(double y, double alphasMu) {
 589
 590   //Fix for singularity at endpoint
 591   if (y >= 1.0) y = 0.9999999999;
 592
 593   return ( - 4.*(alphasMu/(3.*EvtConst::pi*(1. - y)))*(log(1. - y) + 7./4.)*
 594            exp(-2.*(alphasMu/(3.*EvtConst::pi))*(pow(log(1. - y),2) + (7./2.)*log(1. - y))));
 595
 596 }
 597
 598 double EvtBtoXsgammaKagan::s77(double y) {
 599
 600   //Fix for singularity at endpoint
 601   if (y >= 1.0) y = 0.9999999999;
 602
 603   return ((1./3.)*(7. + y - 2.*pow(y,2) - 2.*(1. + y)*log(1. - y)));
 604 }
 605
 606 double EvtBtoXsgammaKagan::s88(double y, double mb, double ms) {
 607
 608   //Fix for singularity at endpoint
 609   if (y >= 1.0) y = 0.9999999999;
 610
 611   return ((1./27.)*((2.*(2. - 2.*y + pow(y,2))/y)*(log(1. - y) + 2.*log(mb/ms))
 612                     - 2.*pow(y,2) - y - 8.*((1. - y)/y)));
 613 }
 614
 615 double EvtBtoXsgammaKagan::s78(double y) {
 616
 617   //Fix for singularity at endpoint
 618   if (y >= 1.0) y = 0.9999999999;
 619
 620   return ((8./9.)*(((1. - y)/y)*log(1. - y) + 1. + (pow(y,2)/4.)));
 621 }
 622
 623 double EvtBtoXsgammaKagan::ReG(double y) {
 624
 625   if (y < 4.) return -2.*pow(atan(sqrt(y/(4. - y))),2.);
 626   else {
 627     return 2.*(pow(log((sqrt(y) + sqrt(y - 4.))/2.),2.)) - (1./2.)*pow(EvtConst::pi,2.);
 628   }
 629
 630 }
 631
 632 double EvtBtoXsgammaKagan::ImG(double y) {
 633
 634   if (y < 4.) return 0.0;
 635   else {
 636     return (-2.*EvtConst::pi*log((sqrt(y) + sqrt(y - 4.))/2.));
 637   }
 638 }
 639
 640 double EvtBtoXsgammaKagan::s22Func(double y, const std::vector<double> &coeffs) {
 641
 642   //coeffs[0]=z
 643   return (1. - y)*((pow(coeffs[0],2.)/pow(y,2.))*(pow(ReG(y/coeffs[0]),2.) + pow(ImG(y/coeffs[0]),2.)) + (coeffs[0]/y)*ReG(y/coeffs[0]) + (1./4.));
 644
 645 }
 646
 647 double EvtBtoXsgammaKagan::s27Func(double y, const std::vector<double> &coeffs) {
 648
 649   //coeffs[0] = z
 650   return (ReG(y/coeffs[0]) + y/(2.*coeffs[0]));
 651
 652 }
 653
 654 double EvtBtoXsgammaKagan::DeltaFermiFunc(double y, const std::vector<double> &coeffs1,
 655                                         const std::vector<double> &coeffs2, const std::vector<double> &coeffs3) {
 656
 657   //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
 658   //coeffs2[2]=ymH, coeffs3[0]=DeltaCoeff (alphasmu)
 659
 660   return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))*
 661     Delta((coeffs2[0]*coeffs2[2])/(coeffs2[1]+y),coeffs3[0]);
 662
 663 }
 664
 665 double EvtBtoXsgammaKagan::s77FermiFunc(double y, const std::vector<double> &coeffs1,
 666                                       const std::vector<double> &coeffs2) {
 667
 668   //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
 669   //coeffs2[2]=ymH
 670   return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))*
 671     s77((coeffs2[0]*coeffs2[2])/(coeffs2[1]+y));
 672
 673 }
 674
 675 double EvtBtoXsgammaKagan::s88FermiFunc(double y, const std::vector<double> &coeffs1,
 676                                       const std::vector<double> &coeffs2, const std::vector<double> &coeffs3) {
 677
 678   //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
 679   //coeffs2[2]=ymH, coeffs3=s88 coeffs
 680   return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))*
 681    s88((coeffs2[0]*coeffs2[2])/(coeffs2[1]+y),coeffs3[0], coeffs3[1]);
 682
 683 }
 684
 685 double EvtBtoXsgammaKagan::s78FermiFunc(double y, const std::vector<double> &coeffs1,
 686                                       const std::vector<double> &coeffs2) {
 687
 688   //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
 689   //coeffs2[2]=ymH
 690   return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))*
 691     s78((coeffs2[0]*coeffs2[2])/(coeffs2[1]+y));
 692
 693 }
 694
 695 double EvtBtoXsgammaKagan::sFermiFunc(double y, const std::vector<double> &coeffs1,
 696                                       const std::vector<double> &coeffs2, const std::vector<double> &coeffs3,
 697                                       const std::vector<double> &coeffs4) {
 698
 699   //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
 700   //coeffs2[2]=ymH, coeffs3[0]=nIntervals in s22 or s27 array, coeffs3[1]=yMin,
 701   //coeffs3[2]=yMax, coeffs4=s22 or s27 array
 702   return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))*
 703     GetArrayVal(coeffs2[0]*coeffs2[2]/(coeffs2[1]+y), coeffs3[0], coeffs3[1], coeffs3[2], coeffs4);
 704
 705 }
 706
 707 double EvtBtoXsgammaKagan::Fz(double z) {
 708
 709   return (1. -8.*z + 8.*pow(z,3.) - pow(z,4.) - 12.*pow(z,2.)*log(z));
 710 }
 711
 712 double EvtBtoXsgammaKagan::GetArrayVal(double xp, double nInterval, double xMin, double xMax, std::vector<double> array) {
 713
 714   double dx = (xMax - xMin)/nInterval;
 715   int bin1 = int(((xp-xMin)/(xMax - xMin))*nInterval);
 716
 717   double x1 = double(bin1)*dx + xMin;
 718
 719   if (xp == x1) return array[bin1];
 720
 721   int bin2(0);
 722   if (xp > x1) {
 723     bin2 = bin1 + 1;
 724   }
 725   else if (xp < x1) {
 726     bin2 = bin1 - 1;
 727   }
 728
 729   if (bin1 <= 0) {
 730     bin1=0;
 731     bin2 = 1;
 732   }
 733
 734   //If xp is in the last bin, always interpolate between the last two bins
 735   if (bin1 == (int)nInterval){
 736     bin2 = (int)nInterval;
 737     bin1 = (int)nInterval - 1;
 738     x1 = double(bin1)*dx + xMin;
 739   }
 740
 741   double x2 = double(bin2)*dx + xMin;
 742   double y1 = array[bin1];
 743
 744   double y2 = array[bin2];
 745   double m = (y2 - y1)/(x2 - x1);
 746   double c =  y1 - m*x1;
 747   double result = m*xp + c;
 748
 749   return result;
 750
 751 }
 752
 753 double EvtBtoXsgammaKagan::FermiFunc(double y, const std::vector<double> &coeffs) {
 754
 755   //Fermi shape functions :1=exponential, 2=gaussian, 3=roman
 756   if (int(coeffs[0]) == 1) return EvtBtoXsgammaFermiUtil::FermiExpFunc(y, coeffs);
 757   if (int(coeffs[0]) == 2) return EvtBtoXsgammaFermiUtil::FermiGaussFunc(y, coeffs);
 758   if (int(coeffs[0]) == 3) return EvtBtoXsgammaFermiUtil::FermiRomanFunc(y, coeffs);
 759   return 1.;
 760
 761 }
 762
 763 double EvtBtoXsgammaKagan::diLogFunc(double y) {
 764
 765   return -log(fabs(1. - y))/y;
 766
 767 }
 768
 769
 770 double EvtBtoXsgammaKagan::diLogMathematica(double y) {
 771
 772   double li2(0);
 773   for(int i=1; i<1000; i++){ //the value 1000 should actually be Infinite...
 774     li2+=pow(y,i)/(i*i);
 775   }
 776   return li2;
 777 }