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

Changeset View

Standalone View

View Options

src/EvtGenModels/EvtVub.cpp

Show First 20 Lines • Show All 268 Lines • ▼ Show 20 Linesvoid EvtVub::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 );
if ( _storeQplus ) { if ( _storeQplus ) {
// cludge to store the hidden parameter q+ with the decay; // cludge to store the hidden parameter q+ 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 q+/(10000 GeV) which // stay well below BaBars sensitivity we take q+/(10000 GeV) which
// goes up to 0.0005 in the most extreme cases as ctau in mm. // goes up to 0.0005 in the most extreme cases as ctau in mm.
// To extract q+ back from the StdHepTrk its necessary to get // To extract q+ back from the StdHepTrk its necessary to get
// delta_ctau = Xu->anyDaughter->getVertexTime()-Xu->getVertexTime() // delta_ctau = Xu->anyDaughter->getVertexTime()-Xu->getVertexTime()
// where these pseudo calls refere to the StdHep time stored at // where these pseudo calls refere to the StdHep time stored at
// the production vertex in the lab for each particle. The boost // the production vertex in the lab for each particle. The boost
// has to be reversed and the result is: // has to be reversed and the result is:
// //
// q+ = delta_ctau * 10000 GeV/mm * Mass_Xu/Energy_Xu // q+ = delta_ctau * 10000 GeV/mm * Mass_Xu/Energy_Xu
// //
xuhad->setLifetime( qplus / 10000. ); xuhad->setLifetime( qplus / 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] = {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;
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++ )
pLW[j + 1] = sttmp * cos( phL ) * ptmp * xW[j] + pLW[j + 1] = sttmp * cos( phL ) * ptmp * xW[j] +
sttmp * sin( phL ) * ptmp * yW[j] + ctL * ptmp * zW[j]; sttmp * sin( phL ) * ptmp * yW[j] + ctL * ptmp * zW[j];
double apLW = ptmp; double apLW = ptmp;
// boost them back in the B Meson restframe // boost them back in the B Meson restframe
double appLB = beta * gamma * pLW[0] + gamma * ctL * apLW; double appLB = beta * gamma * pLW[0] + gamma * ctL * apLW;
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 All 38 Lines