Changeset View
Changeset View
Standalone View
Standalone View
src/EvtGenModels/EvtVubBLNP.cpp
Show First 20 Lines • Show All 258 Lines • ▼ Show 20 Lines | void EvtVubBLNP::decay( EvtParticle* Bmeson ) | ||||
// 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( EX * EX - sh ); | ptmp = sqrt( EX * EX - sh ); | ||||
double pHB[4] = {EX, ptmp * sttmp * cos( phH ), ptmp * sttmp * sin( phH ), | double pHB[4] = { EX, 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 ); | ||||
bool _storeWhat( true ); | bool _storeWhat( true ); | ||||
if ( _storeWhat ) { | if ( _storeWhat ) { | ||||
// cludge to store the hidden parameter what with the decay; | // cludge to store the hidden parameter what with the decay; | ||||
// the lifetime of the Xu is abused for this purpose. | // the lifetime of the Xu is abused for this purpose. | ||||
// tau = 1 ps corresponds to ctau = 0.3 mm -> in order to | // tau = 1 ps corresponds to ctau = 0.3 mm -> in order to | ||||
// stay well below BaBars sensitivity we take what/(10000 GeV). | // stay well below BaBars sensitivity we take what/(10000 GeV). | ||||
// To extract what back from the StdHepTrk its necessary to get | // To extract what back from the StdHepTrk its necessary to get | ||||
// delta_ctau = Xu->decayVtx()->point().distanceTo(XuDaughter->decayVtx()->point()); | // delta_ctau = Xu->decayVtx()->point().distanceTo(XuDaughter->decayVtx()->point()); | ||||
// | // | ||||
// what = delta_ctau * 100000 * Mass_Xu/Momentum_Xu | // what = delta_ctau * 100000 * Mass_Xu/Momentum_Xu | ||||
// | // | ||||
xuhad->setLifetime( what / 10000. ); | xuhad->setLifetime( what / 10000. ); | ||||
} | } | ||||
// 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] = {mBB - EX, -pHB[1], -pHB[2], -pHB[3]}; | double pWB[4] = { mBB - EX, -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 = mBB * mBB + sh - 2 * mBB * EX; | double mW2 = mBB * mBB + sh - 2 * mBB * EX; | ||||
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 ( j = 0; j < 2; j++ ) | for ( 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 ( j = 0; j < 3; j++ ) | for ( 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 ( j = 0; j < 3; j++ ) | for ( j = 0; j < 3; j++ ) | ||||
Show All 9 Lines | void EvtVubBLNP::decay( EvtParticle* Bmeson ) | ||||
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[4] = {pWB[0] - El, 0, 0, 0}; | double pNB[4] = { pWB[0] - El, 0, 0, 0 }; | ||||
for ( j = 1; j < 4; j++ ) { | for ( 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 20 Lines • Show All 135 Lines • ▼ Show 20 Lines | double EvtVubBLNP::DoneJS( double Pp, double Pm, double /*mui*/ ) | ||||
vars[1] = Pm; | vars[1] = Pm; | ||||
for ( int j = 2; j < 12; j++ ) { | for ( int j = 2; j < 12; j++ ) { | ||||
vars[j] = gvars[j]; | vars[j] = gvars[j]; | ||||
} | } | ||||
double lowerlim = 0.001 * Pp; | double lowerlim = 0.001 * Pp; | ||||
double upperlim = ( 1.0 - 0.001 ) * Pp; | double upperlim = ( 1.0 - 0.001 ) * Pp; | ||||
auto func = EvtItgPtrFunction{&IntJS, lowerlim, upperlim, vars}; | auto func = EvtItgPtrFunction{ &IntJS, lowerlim, upperlim, vars }; | ||||
auto integ = EvtItgSimpsonIntegrator{func, precision, maxLoop}; | auto integ = EvtItgSimpsonIntegrator{ func, precision, maxLoop }; | ||||
return integ.evaluate( lowerlim, upperlim ); | return integ.evaluate( lowerlim, upperlim ); | ||||
} | } | ||||
double EvtVubBLNP::Done1( double Pp, double Pm, double /*mui*/ ) | double EvtVubBLNP::Done1( double Pp, double Pm, double /*mui*/ ) | ||||
{ | { | ||||
std::vector<double> vars( 12 ); | std::vector<double> vars( 12 ); | ||||
vars[0] = Pp; | vars[0] = Pp; | ||||
vars[1] = Pm; | vars[1] = Pm; | ||||
for ( int j = 2; j < 12; j++ ) { | for ( int j = 2; j < 12; j++ ) { | ||||
vars[j] = gvars[j]; | vars[j] = gvars[j]; | ||||
} | } | ||||
double lowerlim = 0.001 * Pp; | double lowerlim = 0.001 * Pp; | ||||
double upperlim = ( 1.0 - 0.001 ) * Pp; | double upperlim = ( 1.0 - 0.001 ) * Pp; | ||||
auto func = EvtItgPtrFunction{&Int1, lowerlim, upperlim, vars}; | auto func = EvtItgPtrFunction{ &Int1, lowerlim, upperlim, vars }; | ||||
auto integ = EvtItgSimpsonIntegrator{func, precision, maxLoop}; | auto integ = EvtItgSimpsonIntegrator{ func, precision, maxLoop }; | ||||
return integ.evaluate( lowerlim, upperlim ); | return integ.evaluate( lowerlim, upperlim ); | ||||
} | } | ||||
double EvtVubBLNP::Done2( double Pp, double Pm, double /*mui*/ ) | double EvtVubBLNP::Done2( double Pp, double Pm, double /*mui*/ ) | ||||
{ | { | ||||
std::vector<double> vars( 12 ); | std::vector<double> vars( 12 ); | ||||
vars[0] = Pp; | vars[0] = Pp; | ||||
vars[1] = Pm; | vars[1] = Pm; | ||||
for ( int j = 2; j < 12; j++ ) { | for ( int j = 2; j < 12; j++ ) { | ||||
vars[j] = gvars[j]; | vars[j] = gvars[j]; | ||||
} | } | ||||
double lowerlim = 0.001 * Pp; | double lowerlim = 0.001 * Pp; | ||||
double upperlim = ( 1.0 - 0.001 ) * Pp; | double upperlim = ( 1.0 - 0.001 ) * Pp; | ||||
auto func = EvtItgPtrFunction{&Int2, lowerlim, upperlim, vars}; | auto func = EvtItgPtrFunction{ &Int2, lowerlim, upperlim, vars }; | ||||
auto integ = EvtItgSimpsonIntegrator{func, precision, maxLoop}; | auto integ = EvtItgSimpsonIntegrator{ func, precision, maxLoop }; | ||||
return integ.evaluate( lowerlim, upperlim ); | return integ.evaluate( lowerlim, upperlim ); | ||||
} | } | ||||
double EvtVubBLNP::Done3( double Pp, double Pm, double /*mui*/ ) | double EvtVubBLNP::Done3( double Pp, double Pm, double /*mui*/ ) | ||||
{ | { | ||||
std::vector<double> vars( 12 ); | std::vector<double> vars( 12 ); | ||||
vars[0] = Pp; | vars[0] = Pp; | ||||
vars[1] = Pm; | vars[1] = Pm; | ||||
for ( int j = 2; j < 12; j++ ) { | for ( int j = 2; j < 12; j++ ) { | ||||
vars[j] = gvars[j]; | vars[j] = gvars[j]; | ||||
} | } | ||||
double lowerlim = 0.001 * Pp; | double lowerlim = 0.001 * Pp; | ||||
double upperlim = ( 1.0 - 0.001 ) * Pp; | double upperlim = ( 1.0 - 0.001 ) * Pp; | ||||
auto func = EvtItgPtrFunction{&Int3, lowerlim, upperlim, vars}; | auto func = EvtItgPtrFunction{ &Int3, lowerlim, upperlim, vars }; | ||||
auto integ = EvtItgSimpsonIntegrator{func, precision, maxLoop}; | auto integ = EvtItgSimpsonIntegrator{ func, precision, maxLoop }; | ||||
return integ.evaluate( lowerlim, upperlim ); | return integ.evaluate( lowerlim, upperlim ); | ||||
} | } | ||||
double EvtVubBLNP::Int1( double what, const std::vector<double>& vars ) | double EvtVubBLNP::Int1( double what, const std::vector<double>& vars ) | ||||
{ | { | ||||
return Shat( what, vars ) * g1( what, vars ); | return Shat( what, vars ) * g1( what, vars ); | ||||
} | } | ||||
▲ Show 20 Lines • Show All 490 Lines • Show Last 20 Lines |