Page MenuHomeHEPForge
Files F19245099

EvtBtoXsgammaKagan.cpp
No OneTemporary
Actions

Size
34 KB
Referenced Files
None
Subscribers
None

EvtBtoXsgammaKagan.cpp
View Options

/***********************************************************************
* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
* *
* This file is part of EvtGen. *
* *
* EvtGen is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* EvtGen is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
***********************************************************************/
#include "EvtGenModels/EvtBtoXsgammaKagan.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenModels/EvtBtoXsgamma.hh"
#include "EvtGenModels/EvtBtoXsgammaFermiUtil.hh"
#include "EvtGenModels/EvtItgAbsIntegrator.hh"
#include "EvtGenModels/EvtItgFourCoeffFcn.hh"
#include "EvtGenModels/EvtItgFunction.hh"
#include "EvtGenModels/EvtItgPtrFunction.hh"
#include "EvtGenModels/EvtItgSimpsonIntegrator.hh"
#include "EvtGenModels/EvtItgThreeCoeffFcn.hh"
#include "EvtGenModels/EvtItgTwoCoeffFcn.hh"
#include <fstream>
#include <stdlib.h>
#include <string>
using std::endl;
using std::fstream;
bool EvtBtoXsgammaKagan::bbprod = false;
double EvtBtoXsgammaKagan::intervalMH = 0;
void EvtBtoXsgammaKagan::init( int nArg, double* args )
{
if ( ( nArg ) > 12 || ( nArg > 1 && nArg < 10 ) || nArg == 11 ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "EvtBtoXsgamma generator model "
<< "EvtBtoXsgammaKagan expected "
<< "either 1(default config) or "
<< "10 (default mass range) or "
<< "12 (user range) arguments but found: " << nArg << endl;
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Will terminate execution!" << endl;
::abort();
}
if ( nArg == 1 ) {
bbprod = true;
getDefaultHadronicMass();
} else {
bbprod = false;
computeHadronicMass( nArg, args );
}
double mHminLimit = 0.6373;
double mHmaxLimit = 4.5;
if ( nArg > 10 ) {
_mHmin = args[10];
_mHmax = args[11];
if ( _mHmin > _mHmax ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Minimum hadronic mass exceeds maximum " << endl;
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Will terminate execution!" << endl;
::abort();
}
if ( _mHmin < mHminLimit ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Minimum hadronic mass below K pi threshold" << endl;
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Resetting to K pi threshold" << endl;
_mHmin = mHminLimit;
}
if ( _mHmax > mHmaxLimit ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Maximum hadronic mass above 4.5 GeV/c^2" << endl;
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Resetting to 4.5 GeV/c^2" << endl;
_mHmax = mHmaxLimit;
}
} else {
_mHmin = mHminLimit; // usually just above K pi threshold for Xsd/u
_mHmax = mHmaxLimit;
}
}
void EvtBtoXsgammaKagan::getDefaultHadronicMass()
{
massHad = {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};
brHad = {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};
massHad.resize( 81 );
brHad.resize( 81 );
intervalMH = 80;
}
void EvtBtoXsgammaKagan::computeHadronicMass( int /*nArg*/, double* args )
{
//Input parameters
int fermiFunction = (int)args[1];
_mB = args[2];
_mb = args[3];
_mu = args[4];
_lam1 = args[5];
_delta = args[6];
_z = args[7];
_nIntervalS = args[8];
_nIntervalmH = args[9];
std::vector<double> mHVect( int( _nIntervalmH + 1.0 ) );
massHad.clear();
massHad.resize( int( _nIntervalmH + 1.0 ) );
brHad.clear();
brHad.resize( int( _nIntervalmH + 1.0 ) );
intervalMH = _nIntervalmH;
//Going to have to add a new entry into the data file - takes ages...
EvtGenReport( EVTGEN_WARNING, "EvtGen" )
<< "EvtBtoXsgammaKagan: calculating new hadronic mass spectra. This takes a while..."
<< endl;
//Now need to compute the mHVect vector for
//the current parameters
//A few more parameters
double _mubar = _mu;
_mW = 80.33;
_mt = 175.0;
_alpha = 1. / 137.036;
_lambdabar = _mB - _mb;
_kappabar = 3.382 - 4.14 * ( sqrt( _z ) - 0.29 );
_fz = Fz( _z );
_rer8 = ( 44. / 9. ) - ( 8. / 27. ) * pow( EvtConst::pi, 2. );
_r7 = ( -10. / 3. ) - ( 8. / 9. ) * pow( EvtConst::pi, 2. );
_rer2 = -4.092 + 12.78 * ( sqrt( _z ) - .29 );
_gam77 = 32. / 3.;
_gam27 = 416. / 81.;
_gam87 = -32. / 9.;
_lam2 = .12;
_beta0 = 23. / 3.;
_beta1 = 116. / 3.;
_alphasmZ = .118;
_mZ = 91.187;
_ms = _mb / 50.;
double eGammaMin = 0.5 * _mB * ( 1. - _delta );
double eGammaMax = 0.5 * _mB;
double yMin = 2. * eGammaMin / _mB;
double yMax = 2. * eGammaMax / _mB;
double _CKMrat = 0.976;
double Nsl = 1.0;
//Calculate alpha the various scales
_alphasmW = CalcAlphaS( _mW );
_alphasmt = CalcAlphaS( _mt );
_alphasmu = CalcAlphaS( _mu );
_alphasmubar = CalcAlphaS( _mubar );
//Calculate the Wilson Coefficients and Delta
_etamu = _alphasmW / _alphasmu;
_kSLemmu = ( 12. / 23. ) * ( ( 1. / _etamu ) - 1. );
CalcWilsonCoeffs();
CalcDelta();
//Build s22 and s27 vector - saves time because double
//integration is required otherwise
std::vector<double> s22Coeffs( int( _nIntervalS + 1.0 ) );
std::vector<double> s27Coeffs( int( _nIntervalS + 1.0 ) );
std::vector<double> s28Coeffs( int( _nIntervalS + 1.0 ) );
double dy = ( yMax - yMin ) / _nIntervalS;
double yp = yMin;
std::vector<double> sCoeffs( 1 );
sCoeffs[0] = _z;
//Define s22 and s27 functions
auto mys22Func = EvtItgPtrFunction{&s22Func, 0., yMax + 0.1, sCoeffs};
auto mys27Func = EvtItgPtrFunction{&s27Func, 0., yMax + 0.1, sCoeffs};
//Use a simpson integrator
auto mys22Simp = EvtItgSimpsonIntegrator{mys22Func, 1.0e-4, 20};
auto mys27Simp = EvtItgSimpsonIntegrator{mys27Func, 1.0e-4, 50};
int i;
for ( i = 0; i < int( _nIntervalS + 1.0 ); i++ ) {
s22Coeffs[i] = ( 16. / 27. ) * mys22Simp.evaluate( 1.0e-20, yp );
s27Coeffs[i] = ( -8. / 9. ) * _z * mys27Simp.evaluate( 1.0e-20, yp );
s28Coeffs[i] = -s27Coeffs[i] / 3.;
yp = yp + dy;
}
//Define functions and vectors used to calculate mHVect. Each function takes a set
//of vectors which are used as the function coefficients
std::vector<double> FermiCoeffs( 6 );
std::vector<double> varCoeffs( 3 );
std::vector<double> DeltaCoeffs( 1 );
std::vector<double> s88Coeffs( 2 );
std::vector<double> sInitCoeffs( 3 );
varCoeffs[0] = _mB;
varCoeffs[1] = _mb;
varCoeffs[2] = 0.;
DeltaCoeffs[0] = _alphasmu;
s88Coeffs[0] = _mb;
s88Coeffs[1] = _ms;
sInitCoeffs[0] = _nIntervalS;
sInitCoeffs[1] = yMin;
sInitCoeffs[2] = yMax;
FermiCoeffs[0] = fermiFunction;
FermiCoeffs[1] = 0.0;
FermiCoeffs[2] = 0.0;
FermiCoeffs[3] = 0.0;
FermiCoeffs[4] = 0.0;
FermiCoeffs[5] = 0.0;
//Coefficients for gamma function
std::vector<double> gammaCoeffs( 6 );
gammaCoeffs[0] = 76.18009172947146;
gammaCoeffs[1] = -86.50532032941677;
gammaCoeffs[2] = 24.01409824083091;
gammaCoeffs[3] = -1.231739572450155;
gammaCoeffs[4] = 0.1208650973866179e-2;
gammaCoeffs[5] = -0.5395239384953e-5;
//Calculate quantities for the fermi function to be used
//Distinguish among the different shape functions
if ( fermiFunction == 1 ) {
FermiCoeffs[1] = _lambdabar;
FermiCoeffs[2] = ( -3. * pow( _lambdabar, 2. ) / _lam1 ) - 1.;
FermiCoeffs[3] = _lam1;
FermiCoeffs[4] = 1.0;
auto myNormFunc = std::make_unique<EvtItgPtrFunction>(
&EvtBtoXsgammaFermiUtil::FermiExpFunc, -_mb, _mB - _mb, FermiCoeffs );
auto myNormSimp =
std::make_unique<EvtItgSimpsonIntegrator>( *myNormFunc, 1.0e-4, 40 );
FermiCoeffs[4] = myNormSimp->normalisation();
} else if ( fermiFunction == 2 ) {
double a = EvtBtoXsgammaFermiUtil::FermiGaussFuncRoot( _lambdabar, _lam1,
_mb, gammaCoeffs );
FermiCoeffs[1] = _lambdabar;
FermiCoeffs[2] = a;
FermiCoeffs[3] =
EvtBtoXsgammaFermiUtil::Gamma( ( 2.0 + a ) / 2., gammaCoeffs ) /
EvtBtoXsgammaFermiUtil::Gamma( ( 1.0 + a ) / 2., gammaCoeffs );
FermiCoeffs[4] = 1.0;
auto myNormFunc = std::make_unique<EvtItgPtrFunction>(
&EvtBtoXsgammaFermiUtil::FermiGaussFunc, -_mb, _mB - _mb,
FermiCoeffs );
auto myNormSimp =
std::make_unique<EvtItgSimpsonIntegrator>( *myNormFunc, 1.0e-4, 40 );
FermiCoeffs[4] = myNormSimp->normalisation();
} else if ( fermiFunction == 3 ) {
double rho = EvtBtoXsgammaFermiUtil::FermiRomanFuncRoot( _lambdabar,
_lam1 );
FermiCoeffs[1] = _mB;
FermiCoeffs[2] = _mb;
FermiCoeffs[3] = rho;
FermiCoeffs[4] = _lambdabar;
FermiCoeffs[5] = 1.0;
auto myNormFunc = std::make_unique<EvtItgPtrFunction>(
&EvtBtoXsgammaFermiUtil::FermiRomanFunc, -_mb, _mB - _mb,
FermiCoeffs );
auto myNormSimp =
std::make_unique<EvtItgSimpsonIntegrator>( *myNormFunc, 1.0e-4, 40 );
FermiCoeffs[5] = myNormSimp->normalisation();
}
//Define functions
auto myDeltaFermiFunc = EvtItgThreeCoeffFcn{&DeltaFermiFunc, -_mb,
_mB - _mb, FermiCoeffs,
varCoeffs, DeltaCoeffs};
auto mys88FermiFunc = EvtItgThreeCoeffFcn{&s88FermiFunc, -_mb,
_mB - _mb, FermiCoeffs,
varCoeffs, s88Coeffs};
auto mys77FermiFunc = EvtItgTwoCoeffFcn{&s77FermiFunc, -_mb, _mB - _mb,
FermiCoeffs, varCoeffs};
auto mys78FermiFunc = EvtItgTwoCoeffFcn{&s78FermiFunc, -_mb, _mB - _mb,
FermiCoeffs, varCoeffs};
auto mys22FermiFunc = EvtItgFourCoeffFcn{&sFermiFunc, -_mb,
_mB - _mb, FermiCoeffs,
varCoeffs, sInitCoeffs,
s22Coeffs};
auto mys27FermiFunc = EvtItgFourCoeffFcn{&sFermiFunc, -_mb,
_mB - _mb, FermiCoeffs,
varCoeffs, sInitCoeffs,
s27Coeffs};
auto mys28FermiFunc = EvtItgFourCoeffFcn{&sFermiFunc, -_mb,
_mB - _mb, FermiCoeffs,
varCoeffs, sInitCoeffs,
s28Coeffs};
//Define integrators
auto myDeltaFermiSimp = EvtItgSimpsonIntegrator{myDeltaFermiFunc, 1.0e-4, 40};
auto mys77FermiSimp = EvtItgSimpsonIntegrator{mys77FermiFunc, 1.0e-4, 40};
auto mys88FermiSimp = EvtItgSimpsonIntegrator{mys88FermiFunc, 1.0e-4, 40};
auto mys78FermiSimp = EvtItgSimpsonIntegrator{mys78FermiFunc, 1.0e-4, 40};
auto mys22FermiSimp = EvtItgSimpsonIntegrator{mys22FermiFunc, 1.0e-4, 40};
auto mys27FermiSimp = EvtItgSimpsonIntegrator{mys27FermiFunc, 1.0e-4, 40};
auto mys28FermiSimp = EvtItgSimpsonIntegrator{mys28FermiFunc, 1.0e-4, 40};
//Finally calculate mHVect for the range of hadronic masses
double mHmin = sqrt( _mB * _mB - 2. * _mB * eGammaMax );
double mHmax = sqrt( _mB * _mB - 2. * _mB * eGammaMin );
double dmH = ( mHmax - mHmin ) / _nIntervalmH;
double mH = mHmin;
//Calculating the Branching Fractions
for ( i = 0; i < int( _nIntervalmH + 1.0 ); i++ ) {
double ymH = 1. - ( ( mH * mH ) / ( _mB * _mB ) );
//Need to set ymH as one of the input parameters
myDeltaFermiFunc.setCoeff( 2, 2, ymH );
mys77FermiFunc.setCoeff( 2, 2, ymH );
mys88FermiFunc.setCoeff( 2, 2, ymH );
mys78FermiFunc.setCoeff( 2, 2, ymH );
mys22FermiFunc.setCoeff( 2, 2, ymH );
mys27FermiFunc.setCoeff( 2, 2, ymH );
mys28FermiFunc.setCoeff( 2, 2, ymH );
//Integrate
double deltaResult = myDeltaFermiSimp.evaluate( ( _mB * ymH - _mb ),
_mB - _mb );
double s77Result = mys77FermiSimp.evaluate( ( _mB * ymH - _mb ),
_mB - _mb );
double s88Result = mys88FermiSimp.evaluate( ( _mB * ymH - _mb ),
_mB - _mb );
double s78Result = mys78FermiSimp.evaluate( ( _mB * ymH - _mb ),
_mB - _mb );
double s22Result = mys22FermiSimp.evaluate( ( _mB * ymH - _mb ),
_mB - _mb );
double s27Result = mys27FermiSimp.evaluate( ( _mB * ymH - _mb ),
_mB - _mb );
mys28FermiSimp.evaluate( ( _mB * ymH - _mb ), _mB - _mb );
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 ) ) );
mHVect[i] = 2. * ( mH / ( _mB * _mB ) ) * 0.105 * Nsl * py;
massHad[i] = mH;
brHad[i] = 2. * ( mH / ( _mB * _mB ) ) * 0.105 * Nsl * py;
mH = mH + dmH;
}
}
double EvtBtoXsgammaKagan::GetMass( int /*Xscode*/ )
{
// Get hadronic mass for the event according to the hadronic mass spectra computed in computeHadronicMass
double mass = 0.0;
double min = _mHmin;
if ( bbprod )
min = 1.1;
// double max=4.5;
double max = _mHmax;
double xbox( 0 ), ybox( 0 );
double boxheight( 0 );
double trueHeight( 0 );
double boxwidth = max - min;
double wgt( 0. );
for ( int i = 0; i < int( intervalMH + 1.0 ); i++ ) {
if ( brHad[i] > boxheight )
boxheight = brHad[i];
}
while ( ( mass > max ) || ( mass < min ) ) {
xbox = EvtRandom::Flat( boxwidth ) + min;
ybox = EvtRandom::Flat( boxheight );
trueHeight = 0.0;
// Correction by Peter Richardson
for ( int i = 1; i < int( intervalMH + 1.0 ); ++i ) {
if ( ( massHad[i] >= xbox ) && ( 0.0 == trueHeight ) ) {
wgt = ( xbox - massHad[i - 1] ) / ( massHad[i] - massHad[i - 1] );
trueHeight = brHad[i - 1] + wgt * ( brHad[i] - brHad[i - 1] );
}
}
if ( ybox > trueHeight ) {
mass = 0.0;
} else {
mass = xbox;
}
}
return mass;
}
double EvtBtoXsgammaKagan::CalcAlphaS( double scale )
{
double v = 1. - _beta0 * ( _alphasmZ / ( 2. * EvtConst::pi ) ) *
( log( _mZ / scale ) );
return ( _alphasmZ / v ) *
( 1. - ( ( _beta1 / _beta0 ) *
( _alphasmZ / ( 4. * EvtConst::pi ) ) * ( log( v ) / v ) ) );
}
void EvtBtoXsgammaKagan::CalcWilsonCoeffs()
{
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 ) );
double xt = pow( mtatmw, 2. ) / pow( _mW, 2. );
/////LO
_c2mu = .5 * pow( _etamu, ( -12. / 23. ) ) + .5 * pow( _etamu, ( 6. / 23. ) );
double c7mWsm = ( ( 3. * pow( xt, 3. ) - 2. * pow( xt, 2. ) ) /
( 4. * pow( ( xt - 1. ), 4. ) ) ) *
log( xt ) +
( ( -8. * pow( xt, 3. ) - 5. * pow( xt, 2. ) + 7. * xt ) /
( 24. * pow( ( xt - 1. ), 3. ) ) );
double c8mWsm = ( ( -3. * pow( xt, 2. ) ) / ( 4. * pow( ( xt - 1. ), 4. ) ) ) *
log( xt ) +
( ( -pow( xt, 3. ) + 5. * pow( xt, 2. ) + 2. * xt ) /
( 8. * pow( ( xt - 1. ), 3. ) ) );
double c7constmu = ( 626126. / 272277. ) * pow( _etamu, ( 14. / 23. ) ) -
( 56281. / 51730. ) * pow( _etamu, ( 16. / 23. ) ) -
( 3. / 7. ) * pow( _etamu, ( 6. / 23. ) ) -
( 1. / 14. ) * pow( _etamu, ( -12. / 23. ) ) -
.6494 * pow( _etamu, .4086 ) -
.038 * pow( _etamu, -.423 ) -
.0186 * pow( _etamu, -.8994 ) -
.0057 * pow( _etamu, .1456 );
_c70mu = c7mWsm * pow( _etamu, ( 16. / 23. ) ) +
( 8. / 3. ) *
( pow( _etamu, ( 14. / 23. ) ) - pow( _etamu, ( 16. / 23. ) ) ) *
c8mWsm +
c7constmu;
double c8constmu = ( 313063. / 363036. ) * pow( _etamu, ( 14. / 23. ) ) -
.9135 * pow( _etamu, .4086 ) +
.0873 * pow( _etamu, -.423 ) -
.0571 * pow( _etamu, -.8994 ) +
.0209 * pow( _etamu, .1456 );
_c80mu = c8mWsm * pow( _etamu, ( 14. / 23. ) ) + c8constmu;
//Compute the dilogarithm (PolyLog(2,x)) with the Simpson integrator
//The dilogarithm is defined as: Li_2(x)=Int_0^x(-log(1.-z)/z)
//however, Mathematica implements it as Sum[z^k/k^2,{k,1,Infinity}], so, althought the two
//results are similar and both implemented in the program, we prefer to use the
//one closer to the Mathematica implementation as identical to what used by the theorists.
// EvtItgFunction *myDiLogFunc = new EvtItgFunction(&diLogFunc, 0., 1.-1./xt);
//EvtItgAbsIntegrator *myDiLogSimp = new EvtItgSimpsonIntegrator(*myDiLogFunc, 1.0e-4, 50);
//double li2 = myDiLogSimp->evaluate(1.0e-20,1.-1./xt);
double li2 = diLogMathematica( 1. - 1. / xt );
double c7mWsm1 =
( ( -16. * pow( xt, 4. ) - 122. * pow( xt, 3. ) + 80. * pow( xt, 2. ) -
8. * xt ) /
( 9. * pow( ( xt - 1. ), 4. ) ) * li2 +
( 6. * pow( xt, 4. ) + 46. * pow( xt, 3. ) - 28. * pow( xt, 2. ) ) /
( 3. * pow( ( xt - 1. ), 5. ) ) * pow( log( xt ), 2. ) +
( -102. * pow( xt, 5. ) - 588. * pow( xt, 4. ) -
2262. * pow( xt, 3. ) + 3244. * pow( xt, 2. ) - 1364. * xt + 208. ) /
( 81. * pow( ( xt - 1 ), 5. ) ) * log( xt ) +
( 1646. * pow( xt, 4. ) + 12205. * pow( xt, 3. ) -
10740. * pow( xt, 2. ) + 2509. * xt - 436. ) /
( 486. * pow( ( xt - 1 ), 4. ) ) );
double c8mWsm1 =
( ( -4. * pow( xt, 4. ) + 40. * pow( xt, 3. ) + 41. * pow( xt, 2. ) + xt ) /
( 6. * pow( ( xt - 1. ), 4. ) ) * li2 +
( -17. * pow( xt, 3. ) - 31. * pow( xt, 2. ) ) /
( 2. * pow( ( xt - 1. ), 5. ) ) * pow( log( xt ), 2. ) +
( -210. * pow( xt, 5. ) + 1086. * pow( xt, 4. ) +
4893. * pow( xt, 3. ) + 2857. * pow( xt, 2. ) - 1994. * xt + 280. ) /
( 216. * pow( ( xt - 1 ), 5. ) ) * log( xt ) +
( 737. * pow( xt, 4. ) - 14102. * pow( xt, 3. ) -
28209. * pow( xt, 2. ) + 610. * xt - 508. ) /
( 1296. * pow( ( xt - 1 ), 4. ) ) );
double E1 = ( xt * ( 18. - 11. * xt - pow( xt, 2. ) ) /
( 12. * pow( ( 1. - xt ), 3. ) ) +
pow( xt, 2. ) * ( 15. - 16. * xt + 4. * pow( xt, 2. ) ) /
( 6. * pow( ( 1. - xt ), 4. ) ) * log( xt ) -
2. / 3. * log( xt ) );
double e1 = 4661194. / 816831.;
double e2 = -8516. / 2217.;
double e3 = 0.;
double e4 = 0.;
double e5 = -1.9043;
double e6 = -.1008;
double e7 = .1216;
double e8 = .0183;
double f1 = -17.3023;
double f2 = 8.5027;
double f3 = 4.5508;
double f4 = .7519;
double f5 = 2.004;
double f6 = .7476;
double f7 = -.5385;
double f8 = .0914;
double g1 = 14.8088;
double g2 = -10.809;
double g3 = -.874;
double g4 = .4218;
double g5 = -2.9347;
double g6 = .3971;
double g7 = .1600;
double g8 = .0225;
double c71constmu =
( ( e1 * _etamu * E1 + f1 + g1 * _etamu ) * pow( _etamu, ( 14. / 23. ) ) +
( e2 * _etamu * E1 + f2 + g2 * _etamu ) * pow( _etamu, ( 16. / 23. ) ) +
( e3 * _etamu * E1 + f3 + g3 * _etamu ) * pow( _etamu, ( 6. / 23. ) ) +
( e4 * _etamu * E1 + f4 + g4 * _etamu ) * pow( _etamu, ( -12. / 23. ) ) +
( e5 * _etamu * E1 + f5 + g5 * _etamu ) * pow( _etamu, .4086 ) +
( e6 * _etamu * E1 + f6 + g6 * _etamu ) * pow( _etamu, ( -.423 ) ) +
( e7 * _etamu * E1 + f7 + g7 * _etamu ) * pow( _etamu, ( -.8994 ) ) +
( e8 * _etamu * E1 + f8 + g8 * _etamu ) * pow( _etamu, .1456 ) );
double c71pmu =
( ( ( 297664. / 14283. * pow( _etamu, ( 16. / 23. ) ) -
7164416. / 357075. * pow( _etamu, ( 14. / 23. ) ) +
256868. / 14283. * pow( _etamu, ( 37. / 23. ) ) -
6698884. / 357075. * pow( _etamu, ( 39. / 23. ) ) ) *
( c8mWsm ) ) +
37208. / 4761. *
( pow( _etamu, ( 39. / 23. ) ) - pow( _etamu, ( 16. / 23. ) ) ) *
( c7mWsm ) +
c71constmu );
_c71mu = ( _alphasmW / _alphasmu *
( pow( _etamu, ( 16. / 23. ) ) * c7mWsm1 +
8. / 3. *
( pow( _etamu, ( 14. / 23. ) ) -
pow( _etamu, ( 16. / 23. ) ) ) *
c8mWsm1 ) +
c71pmu );
_c7emmu = ( ( 32. / 75. * pow( _etamu, ( -9. / 23. ) ) -
40. / 69. * pow( _etamu, ( -7. / 23. ) ) +
88. / 575. * pow( _etamu, ( 16. / 23. ) ) ) *
c7mWsm +
( -32. / 575. * pow( _etamu, ( -9. / 23. ) ) +
32. / 1449. * pow( _etamu, ( -7. / 23. ) ) +
640. / 1449. * pow( _etamu, ( 14. / 23. ) ) -
704. / 1725. * pow( _etamu, ( 16. / 23. ) ) ) *
c8mWsm -
190. / 8073. * pow( _etamu, ( -35. / 23. ) ) -
359. / 3105. * pow( _etamu, ( -17. / 23. ) ) +
4276. / 121095. * pow( _etamu, ( -12. / 23. ) ) +
350531. / 1009125. * pow( _etamu, ( -9. / 23. ) ) +
2. / 4347. * pow( _etamu, ( -7. / 23. ) ) -
5956. / 15525. * pow( _etamu, ( 6. / 23. ) ) +
38380. / 169533. * pow( _etamu, ( 14. / 23. ) ) -
748. / 8625. * pow( _etamu, ( 16. / 23. ) ) );
// Wilson coefficients values as according to Kagan's program
// _c2mu=1.10566;
//_c70mu=-0.314292;
// _c80mu=-0.148954;
// _c71mu=0.480964;
// _c7emmu=0.0323219;
}
void EvtBtoXsgammaKagan::CalcDelta()
{
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. );
double cDelta27 = ( ( _alphasmu / ( 2. * EvtConst::pi ) ) *
( _rer2 + _gam27 * log( _mb / _mu ) ) -
( _lam2 / ( 9. * _z * pow( _mb, 2. ) ) ) ) *
_c2mu * _c70mu;
double cDelta78 = ( _alphasmu / ( 2. * EvtConst::pi ) ) *
( _rer8 + _gam87 * log( _mb / _mu ) ) * _c70mu * _c80mu;
_cDeltatot = cDelta77 + cDelta27 + cDelta78 +
( _alphasmu / ( 2. * EvtConst::pi ) ) * _c71mu * _c70mu +
( _alpha / _alphasmu ) *
( 2. * _c7emmu * _c70mu - _kSLemmu * pow( _c70mu, 2. ) );
}
double EvtBtoXsgammaKagan::Delta( double y, double alphasMu )
{
//Fix for singularity at endpoint
if ( y >= 1.0 )
y = 0.9999999999;
return ( -4. * ( alphasMu / ( 3. * EvtConst::pi * ( 1. - y ) ) ) *
( log( 1. - y ) + 7. / 4. ) *
exp( -2. * ( alphasMu / ( 3. * EvtConst::pi ) ) *
( pow( log( 1. - y ), 2 ) + ( 7. / 2. ) * log( 1. - y ) ) ) );
}
double EvtBtoXsgammaKagan::s77( double y )
{
//Fix for singularity at endpoint
if ( y >= 1.0 )
y = 0.9999999999;
return ( ( 1. / 3. ) *
( 7. + y - 2. * pow( y, 2 ) - 2. * ( 1. + y ) * log( 1. - y ) ) );
}
double EvtBtoXsgammaKagan::s88( double y, double mb, double ms )
{
//Fix for singularity at endpoint
if ( y >= 1.0 )
y = 0.9999999999;
return ( ( 1. / 27. ) * ( ( 2. * ( 2. - 2. * y + pow( y, 2 ) ) / y ) *
( log( 1. - y ) + 2. * log( mb / ms ) ) -
2. * pow( y, 2 ) - y - 8. * ( ( 1. - y ) / y ) ) );
}
double EvtBtoXsgammaKagan::s78( double y )
{
//Fix for singularity at endpoint
if ( y >= 1.0 )
y = 0.9999999999;
return ( ( 8. / 9. ) * ( ( ( 1. - y ) / y ) * log( 1. - y ) + 1. +
( pow( y, 2 ) / 4. ) ) );
}
double EvtBtoXsgammaKagan::ReG( double y )
{
if ( y < 4. )
return -2. * pow( atan( sqrt( y / ( 4. - y ) ) ), 2. );
else {
return 2. * ( pow( log( ( sqrt( y ) + sqrt( y - 4. ) ) / 2. ), 2. ) ) -
( 1. / 2. ) * pow( EvtConst::pi, 2. );
}
}
double EvtBtoXsgammaKagan::ImG( double y )
{
if ( y < 4. )
return 0.0;
else {
return ( -2. * EvtConst::pi * log( ( sqrt( y ) + sqrt( y - 4. ) ) / 2. ) );
}
}
double EvtBtoXsgammaKagan::s22Func( double y, const std::vector<double>& coeffs )
{
//coeffs[0]=z
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. ) );
}
double EvtBtoXsgammaKagan::s27Func( double y, const std::vector<double>& coeffs )
{
//coeffs[0] = z
return ( ReG( y / coeffs[0] ) + y / ( 2. * coeffs[0] ) );
}
double EvtBtoXsgammaKagan::DeltaFermiFunc( double y,
const std::vector<double>& coeffs1,
const std::vector<double>& coeffs2,
const std::vector<double>& coeffs3 )
{
//coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
//coeffs2[2]=ymH, coeffs3[0]=DeltaCoeff (alphasmu)
return FermiFunc( y, coeffs1 ) * ( coeffs2[0] / ( coeffs2[1] + y ) ) *
Delta( ( coeffs2[0] * coeffs2[2] ) / ( coeffs2[1] + y ), coeffs3[0] );
}
double EvtBtoXsgammaKagan::s77FermiFunc( double y,
const std::vector<double>& coeffs1,
const std::vector<double>& coeffs2 )
{
//coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
//coeffs2[2]=ymH
return FermiFunc( y, coeffs1 ) * ( coeffs2[0] / ( coeffs2[1] + y ) ) *
s77( ( coeffs2[0] * coeffs2[2] ) / ( coeffs2[1] + y ) );
}
double EvtBtoXsgammaKagan::s88FermiFunc( double y,
const std::vector<double>& coeffs1,
const std::vector<double>& coeffs2,
const std::vector<double>& coeffs3 )
{
//coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
//coeffs2[2]=ymH, coeffs3=s88 coeffs
return FermiFunc( y, coeffs1 ) * ( coeffs2[0] / ( coeffs2[1] + y ) ) *
s88( ( coeffs2[0] * coeffs2[2] ) / ( coeffs2[1] + y ), coeffs3[0],
coeffs3[1] );
}
double EvtBtoXsgammaKagan::s78FermiFunc( double y,
const std::vector<double>& coeffs1,
const std::vector<double>& coeffs2 )
{
//coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
//coeffs2[2]=ymH
return FermiFunc( y, coeffs1 ) * ( coeffs2[0] / ( coeffs2[1] + y ) ) *
s78( ( coeffs2[0] * coeffs2[2] ) / ( coeffs2[1] + y ) );
}
double EvtBtoXsgammaKagan::sFermiFunc( double y,
const std::vector<double>& coeffs1,
const std::vector<double>& coeffs2,
const std::vector<double>& coeffs3,
const std::vector<double>& coeffs4 )
{
//coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
//coeffs2[2]=ymH, coeffs3[0]=nIntervals in s22 or s27 array, coeffs3[1]=yMin,
//coeffs3[2]=yMax, coeffs4=s22 or s27 array
return FermiFunc( y, coeffs1 ) * ( coeffs2[0] / ( coeffs2[1] + y ) ) *
GetArrayVal( coeffs2[0] * coeffs2[2] / ( coeffs2[1] + y ),
coeffs3[0], coeffs3[1], coeffs3[2], coeffs4 );
}
double EvtBtoXsgammaKagan::Fz( double z )
{
return ( 1. - 8. * z + 8. * pow( z, 3. ) - pow( z, 4. ) -
12. * pow( z, 2. ) * log( z ) );
}
double EvtBtoXsgammaKagan::GetArrayVal( double xp, double nInterval, double xMin,
double xMax, std::vector<double> array )
{
double dx = ( xMax - xMin ) / nInterval;
int bin1 = int( ( ( xp - xMin ) / ( xMax - xMin ) ) * nInterval );
double x1 = double( bin1 ) * dx + xMin;
if ( xp == x1 )
return array[bin1];
int bin2( 0 );
if ( xp > x1 ) {
bin2 = bin1 + 1;
} else if ( xp < x1 ) {
bin2 = bin1 - 1;
}
if ( bin1 <= 0 ) {
bin1 = 0;
bin2 = 1;
}
//If xp is in the last bin, always interpolate between the last two bins
if ( bin1 == (int)nInterval ) {
bin2 = (int)nInterval;
bin1 = (int)nInterval - 1;
x1 = double( bin1 ) * dx + xMin;
}
double x2 = double( bin2 ) * dx + xMin;
double y1 = array[bin1];
double y2 = array[bin2];
double m = ( y2 - y1 ) / ( x2 - x1 );
double c = y1 - m * x1;
double result = m * xp + c;
return result;
}
double EvtBtoXsgammaKagan::FermiFunc( double y, const std::vector<double>& coeffs )
{
//Fermi shape functions :1=exponential, 2=gaussian, 3=roman
if ( int( coeffs[0] ) == 1 )
return EvtBtoXsgammaFermiUtil::FermiExpFunc( y, coeffs );
if ( int( coeffs[0] ) == 2 )
return EvtBtoXsgammaFermiUtil::FermiGaussFunc( y, coeffs );
if ( int( coeffs[0] ) == 3 )
return EvtBtoXsgammaFermiUtil::FermiRomanFunc( y, coeffs );
return 1.;
}
double EvtBtoXsgammaKagan::diLogFunc( double y )
{
return -log( fabs( 1. - y ) ) / y;
}
double EvtBtoXsgammaKagan::diLogMathematica( double y )
{
double li2( 0 );
for ( int i = 1; i < 1000;
i++ ) { //the value 1000 should actually be Infinite...
li2 += pow( y, i ) / ( i * i );
}
return li2;
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 4:48 AM (5 h, 52 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6548199
Default Alt Text
EvtBtoXsgammaKagan.cpp (34 KB)

Event Timeline

HEPForge.org · .