Page MenuHomeHEPForge
Differential D90 Diff 388 src/EvtGenModels/EvtVubNLO.cpp

Changeset View

Standalone View

View Options

src/EvtGenModels/EvtVubNLO.cpp

Show First 20 Lines • Show All 255 Lines • ▼ Show 20 Linesvoid EvtVubNLO::decay( EvtParticle* p )
// now compute the four vectors in the B Meson restframe // now compute the four vectors in the B Meson restframe
double ptmp, sttmp; double ptmp, sttmp;
// calculate the hadron 4 vector in the B Meson restframe // calculate the hadron 4 vector in the B Meson restframe
sttmp = sqrt( 1 - ctH * ctH ); sttmp = sqrt( 1 - ctH * ctH );
ptmp = sqrt( Eh * Eh - sh ); ptmp = sqrt( Eh * Eh - sh );
double pHB[4] = {Eh, ptmp * sttmp * cos( phH ), ptmp * sttmp * sin( phH ), double pHB[4] = { Eh, ptmp * sttmp * cos( phH ), ptmp * sttmp * sin( phH ),
ptmp * ctH}; ptmp * ctH };
p4.set( pHB[0], pHB[1], pHB[2], pHB[3] ); p4.set( pHB[0], pHB[1], pHB[2], pHB[3] );
xuhad->init( getDaug( 0 ), p4 ); xuhad->init( getDaug( 0 ), p4 );
// calculate the W 4 vector in the B Meson restrframe // calculate the W 4 vector in the B Meson restrframe
double apWB = ptmp; double apWB = ptmp;
double pWB[4] = {_mB - Eh, -pHB[1], -pHB[2], -pHB[3]}; double pWB[4] = { _mB - Eh, -pHB[1], -pHB[2], -pHB[3] };
// first go in the W restframe and calculate the lepton and // first go in the W restframe and calculate the lepton and
// the neutrino in the W frame // the neutrino in the W frame
double mW2 = _mB * _mB + sh - 2 * _mB * Eh; double mW2 = _mB * _mB + sh - 2 * _mB * Eh;
// if(mW2<0.1){ // if(mW2<0.1){
// cout <<" low Q2! "<<pp<<" "<<epp<<" "<<x<<" "<<y<<endl; // cout <<" low Q2! "<<pp<<" "<<epp<<" "<<x<<" "<<y<<endl;
//} //}
double beta = ptmp / pWB[0]; double beta = ptmp / pWB[0];
double gamma = pWB[0] / sqrt( mW2 ); double gamma = pWB[0] / sqrt( mW2 );
double pLW[4]; double pLW[4];
ptmp = ( mW2 - ml * ml ) / 2 / sqrt( mW2 ); ptmp = ( mW2 - ml * ml ) / 2 / sqrt( mW2 );
pLW[0] = sqrt( ml * ml + ptmp * ptmp ); pLW[0] = sqrt( ml * ml + ptmp * ptmp );
double ctL = ( El - gamma * pLW[0] ) / beta / gamma / ptmp; double ctL = ( El - gamma * pLW[0] ) / beta / gamma / ptmp;
if ( ctL < -1 ) if ( ctL < -1 )
ctL = -1; ctL = -1;
if ( ctL > 1 ) if ( ctL > 1 )
ctL = 1; ctL = 1;
sttmp = sqrt( 1 - ctL * ctL ); sttmp = sqrt( 1 - ctL * ctL );
// eX' = eZ x eW // eX' = eZ x eW
double xW[3] = {-pWB[2], pWB[1], 0}; double xW[3] = { -pWB[2], pWB[1], 0 };
// eZ' = eW // eZ' = eW
double zW[3] = {pWB[1] / apWB, pWB[2] / apWB, pWB[3] / apWB}; double zW[3] = { pWB[1] / apWB, pWB[2] / apWB, pWB[3] / apWB };
double lx = sqrt( xW[0] * xW[0] + xW[1] * xW[1] ); double lx = sqrt( xW[0] * xW[0] + xW[1] * xW[1] );
for ( int j = 0; j < 2; j++ ) for ( int j = 0; j < 2; j++ )
xW[j] /= lx; xW[j] /= lx;
// eY' = eZ' x eX' // eY' = eZ' x eX'
double yW[3] = {-pWB[1] * pWB[3], -pWB[2] * pWB[3], double yW[3] = { -pWB[1] * pWB[3], -pWB[2] * pWB[3],
pWB[1] * pWB[1] + pWB[2] * pWB[2]}; pWB[1] * pWB[1] + pWB[2] * pWB[2] };
double ly = sqrt( yW[0] * yW[0] + yW[1] * yW[1] + yW[2] * yW[2] ); double ly = sqrt( yW[0] * yW[0] + yW[1] * yW[1] + yW[2] * yW[2] );
for ( int j = 0; j < 3; j++ ) for ( int j = 0; j < 3; j++ )
yW[j] /= ly; yW[j] /= ly;
// p_lep = |p_lep| * ( sin(Theta) * cos(Phi) * eX' // p_lep = |p_lep| * ( sin(Theta) * cos(Phi) * eX'
// + sin(Theta) * sin(Phi) * eY' // + sin(Theta) * sin(Phi) * eY'
// + cos(Theta) * eZ') // + cos(Theta) * eZ')
for ( int j = 0; j < 3; j++ ) for ( int j = 0; j < 3; j++ )
Show All 9 Linesvoid EvtVubNLO::decay( EvtParticle* p )
ptmp = sqrt( El * El - ml * ml ); ptmp = sqrt( El * El - ml * ml );
double ctLL = appLB / ptmp; double ctLL = appLB / ptmp;
if ( ctLL > 1 ) if ( ctLL > 1 )
ctLL = 1; ctLL = 1;
if ( ctLL < -1 ) if ( ctLL < -1 )
ctLL = -1; ctLL = -1;
double pLB[4] = {El, 0, 0, 0}; double pLB[4] = { El, 0, 0, 0 };
double pNB[8] = {pWB[0] - El, 0, 0, 0}; double pNB[8] = { pWB[0] - El, 0, 0, 0 };
for ( int j = 1; j < 4; j++ ) { for ( int j = 1; j < 4; j++ ) {
pLB[j] = pLW[j] + ( ctLL * ptmp - ctL * apLW ) / apWB * pWB[j]; pLB[j] = pLW[j] + ( ctLL * ptmp - ctL * apLW ) / apWB * pWB[j];
pNB[j] = pWB[j] - pLB[j]; pNB[j] = pWB[j] - pLB[j];
} }
p4.set( pLB[0], pLB[1], pLB[2], pLB[3] ); p4.set( pLB[0], pLB[1], pLB[2], pLB[3] );
lepton->init( getDaug( 1 ), p4 ); lepton->init( getDaug( 1 ), p4 );
Show All 22 Linesdouble EvtVubNLO::tripleDiff( double pp, double pl, double pm )
double c1 = ( _mB + pl - pp - pm ) * ( pm - pl ); double c1 = ( _mB + pl - pp - pm ) * ( pm - pl );
double c2 = 2 * ( pl - pp ) * ( pm - pl ); double c2 = 2 * ( pl - pp ) * ( pm - pl );
double c3 = ( _mB - pm ) * ( pm - pp ); double c3 = ( _mB - pm ) * ( pm - pp );
double aF1 = F10( sCoeffs ); double aF1 = F10( sCoeffs );
double aF2 = F20( sCoeffs ); double aF2 = F20( sCoeffs );
double aF3 = F30( sCoeffs ); double aF3 = F30( sCoeffs );
double td0 = c1 * aF1 + c2 * aF2 + c3 * aF3; double td0 = c1 * aF1 + c2 * aF2 + c3 * aF3;
auto func = EvtItgPtrFunction{&integrand, 0., _mB, sCoeffs}; auto func = EvtItgPtrFunction{ &integrand, 0., _mB, sCoeffs };
auto jetSF = EvtItgSimpsonIntegrator{func, 0.01, 25}; auto jetSF = EvtItgSimpsonIntegrator{ func, 0.01, 25 };
double smallfrac = double smallfrac =
0.000001; // stop a bit before the end to avoid problems with numerical integration 0.000001; // stop a bit before the end to avoid problems with numerical integration
double tdInt = jetSF.evaluate( 0, pp * ( 1 - smallfrac ) ); double tdInt = jetSF.evaluate( 0, pp * ( 1 - smallfrac ) );
double SU = U1lo( mu_h(), mu_i() ) * double SU = U1lo( mu_h(), mu_i() ) *
pow( ( pm - pp ) / ( _mB - pp ), alo( mu_h(), mu_i() ) ); pow( ( pm - pp ) / ( _mB - pp ), alo( mu_h(), mu_i() ) );
double TD = ( _mB - pp ) * SU * ( td0 + tdInt ); double TD = ( _mB - pp ) * SU * ( td0 + tdInt );
▲ Show 20 Lines • Show All 316 Lines • ▼ Show 20 Linesdouble EvtVubNLO::expShapeFunction( double omega,
return pow( wL, b - 1 ) * exp( -b * wL ); return pow( wL, b - 1 ) * exp( -b * wL );
} }
double EvtVubNLO::Gamma( double z ) double EvtVubNLO::Gamma( double z )
{ {
std::array<double, 6> gammaCoeffs{ std::array<double, 6> gammaCoeffs{
76.18009172947146, -86.50532032941677, 24.01409824083091, 76.18009172947146, -86.50532032941677, 24.01409824083091,
-1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5}; -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5 };
//Lifted from Numerical Recipies in C //Lifted from Numerical Recipies in C
double y = z; double y = z;
double x = z; double x = z;
double tmp = x + 5.5; double tmp = x + 5.5;
tmp = tmp - ( x + 0.5 ) * log( tmp ); tmp = tmp - ( x + 0.5 ) * log( tmp );
double ser = 1.000000000190015; double ser = 1.000000000190015;
for ( const auto& gammaCoeff : gammaCoeffs ) { for ( const auto& gammaCoeff : gammaCoeffs ) {
y += 1.0; y += 1.0;
ser += gammaCoeff / y; ser += gammaCoeff / y;
} }
return exp( -tmp + log( 2.5066282746310005 * ser / x ) ); return exp( -tmp + log( 2.5066282746310005 * ser / x ) );
} }
double EvtVubNLO::Gamma( double z, double tmin ) double EvtVubNLO::Gamma( double z, double tmin )
{ {
std::vector<double> c( 1 ); std::vector<double> c( 1 );
c[0] = z; c[0] = z;
auto func = EvtItgPtrFunction{&dgamma, tmin, 100., c}; auto func = EvtItgPtrFunction{ &dgamma, tmin, 100., c };
auto jetSF = EvtItgSimpsonIntegrator{func, 0.001}; auto jetSF = EvtItgSimpsonIntegrator{ func, 0.001 };
return jetSF.evaluate( tmin, 100. ); return jetSF.evaluate( tmin, 100. );
} }