Changeset View
Changeset View
Standalone View
Standalone View
src/EvtGenModels/EvtbTosllAmp.cpp
Show All 28 Lines | |||||
#include "EvtGenBase/EvtParticle.hh" | #include "EvtGenBase/EvtParticle.hh" | ||||
#include "EvtGenBase/EvtPatches.hh" | #include "EvtGenBase/EvtPatches.hh" | ||||
#include "EvtGenBase/EvtReport.hh" | #include "EvtGenBase/EvtReport.hh" | ||||
#include "EvtGenBase/EvtScalarParticle.hh" | #include "EvtGenBase/EvtScalarParticle.hh" | ||||
#include "EvtGenBase/EvtTensor4C.hh" | #include "EvtGenBase/EvtTensor4C.hh" | ||||
#include "EvtGenBase/EvtVector4C.hh" | #include "EvtGenBase/EvtVector4C.hh" | ||||
#include "EvtGenBase/EvtVectorParticle.hh" | #include "EvtGenBase/EvtVectorParticle.hh" | ||||
#include <fstream> | |||||
double getKineWeight( EvtParticle* p ) | |||||
{ | |||||
int nd = p->getNDaug(); | |||||
double prod = 1; | |||||
for ( int i = 0; i < nd - 1; i++ ) { | |||||
EvtVector4R M0; | |||||
for ( int j = i + 1; j < nd; j++ ) | |||||
M0 += p->getDaug( j )->getP4(); | |||||
EvtVector4R M1 = M0 + p->getDaug( i )->getP4(); | |||||
double Mi12 = M1.mass2(), Mi = M0.mass(), | |||||
mi = p->getDaug( i )->getP4().mass(), Mipm = Mi + mi, | |||||
Mimm = Mi - mi; | |||||
// std::cout << Mi12 << std::endl; | |||||
prod *= 0.5 * | |||||
sqrt( ( Mi12 - Mipm * Mipm ) * ( Mi12 - Mimm * Mimm ) / Mi12 ); | |||||
} | |||||
return prod; | |||||
} | |||||
double getKineWeight2( EvtParticle* p ) | |||||
{ | |||||
int nd = p->getNDaug(); | |||||
double prod = 1; | |||||
for ( int i = nd - 1; i > 0; i-- ) { | |||||
EvtVector4R M0; | |||||
for ( int j = i - 1; j >= 0; j-- ) | |||||
M0 += p->getDaug( j )->getP4(); | |||||
EvtVector4R M1 = M0 + p->getDaug( i )->getP4(); | |||||
double Mi12 = M1.mass2(), Mi = M0.mass(), | |||||
mi = p->getDaug( i )->getP4().mass(), Mipm = Mi + mi, | |||||
Mimm = Mi - mi; | |||||
prod *= 0.5 * | |||||
sqrt( ( Mi12 - Mipm * Mipm ) * ( Mi12 - Mimm * Mimm ) / Mi12 ); | |||||
} | |||||
return prod; | |||||
} | |||||
double EvtbTosllAmp::CalcMaxProb( EvtId parent, EvtId meson, EvtId lepton1, | double EvtbTosllAmp::CalcMaxProb( EvtId parent, EvtId meson, EvtId lepton1, | ||||
EvtId lepton2, EvtbTosllFF* FormFactors, | EvtId lepton2, EvtbTosllFF* FormFactors, | ||||
double& poleSize ) | double& poleSize ) | ||||
{ | { | ||||
//This routine takes the arguements parent, meson, and lepton | //This routine takes the arguements parent, meson, and lepton | ||||
//number, and a form factor model, and returns a maximum | //number, and a form factor model, and returns a maximum | ||||
//probability for this semileptonic form factor model. A | //probability for this semileptonic form factor model. A | ||||
//brute force method is used. The 2D cos theta lepton and | //brute force method is used. The 2D cos theta lepton and | ||||
Show All 35 Lines | double EvtbTosllAmp::CalcMaxProb( EvtId parent, EvtId meson, EvtId lepton1, | ||||
lep1 = root_part->getDaug( 1 ); | lep1 = root_part->getDaug( 1 ); | ||||
lep2 = root_part->getDaug( 2 ); | lep2 = root_part->getDaug( 2 ); | ||||
//cludge to avoid generating random numbers! | //cludge to avoid generating random numbers! | ||||
daughter->noLifeTime(); | daughter->noLifeTime(); | ||||
lep1->noLifeTime(); | lep1->noLifeTime(); | ||||
lep2->noLifeTime(); | lep2->noLifeTime(); | ||||
//Initial particle is unpolarized, well it is a scalar so it is | //Initial particle is unpolarized, well it is a scalar so it is trivial | ||||
//trivial | |||||
EvtSpinDensity rho; | EvtSpinDensity rho; | ||||
rho.setDiag( root_part->getSpinStates() ); | rho.setDiag( root_part->getSpinStates() ); | ||||
double mass[3]; | double mass[3]; | ||||
double m = root_part->mass(); | double m = root_part->mass(); | ||||
EvtVector4R p4meson, p4lepton1, p4lepton2, p4w; | EvtVector4R p4meson, p4lepton1, p4lepton2, p4w; | ||||
double q2max; | |||||
double q2, elepton, plepton; | |||||
int i, j; | |||||
double erho, prho, costl; | |||||
double maxfoundprob = 0.0; | double maxfoundprob = 0.0; | ||||
double prob = -10.0; | double prob = -10.0; | ||||
int massiter; | |||||
double maxpole = 0; | double maxpole = 0; | ||||
for ( massiter = 0; massiter < 3; massiter++ ) { | |||||
mass[0] = EvtPDL::getMeanMass( meson ); | |||||
mass[1] = EvtPDL::getMeanMass( lepton1 ); | mass[1] = EvtPDL::getMeanMass( lepton1 ); | ||||
mass[2] = EvtPDL::getMeanMass( lepton2 ); | mass[2] = EvtPDL::getMeanMass( lepton2 ); | ||||
if ( massiter == 1 ) { | std::vector<double> mH; | ||||
mass[0] = EvtPDL::getMinMass( meson ); | mH.push_back( EvtPDL::getMeanMass( meson ) ); | ||||
} | //if the particle is narrow dont bother with changing the mass. | ||||
if ( massiter == 2 ) { | double g = EvtPDL::getWidth( meson ); | ||||
mass[0] = EvtPDL::getMaxMass( meson ); | if ( g > 0 ) { | ||||
if ( ( mass[0] + mass[1] + mass[2] ) > m ) | mH.push_back( EvtPDL::getMinMass( meson ) ); | ||||
mass[0] = m - mass[1] - mass[2] - 0.00001; | mH.push_back( | ||||
std::min( EvtPDL::getMaxMass( meson ), m - mass[1] - mass[2] ) ); | |||||
mH.push_back( mH.front() - g ); | |||||
mH.push_back( mH.front() + g ); | |||||
} | } | ||||
q2max = ( m - mass[0] ) * ( m - mass[0] ); | double q2min = ( mass[1] + mass[2] ) * ( mass[1] + mass[2] ); | ||||
for ( double m0 : mH ) { | |||||
double q2max = ( m - m0 ) * ( m - m0 ); | |||||
//loop over q2 | //loop over q2 | ||||
//cout << "m " << m << "mass[0] " << mass[0] << " q2max "<< q2max << endl; | for ( int i = 0, n = 25; i < n; i++ ) { | ||||
for ( i = 0; i < 25; i++ ) { | |||||
//want to avoid picking up the tail of the photon propagator | // want to avoid picking up the tail of the photon propagator | ||||
q2 = ( ( i + 1.5 ) * q2max ) / 26.0; | double q2 = q2min + ( i + 0.5 ) / n * ( q2max - q2min ); | ||||
// std::cout<<"m0 = "<<m0<<" q2 = "<<q2<<" prob = "<<prob<<std::endl; | |||||
if ( i == 0 ) | if ( i == 0 ) | ||||
q2 = 4 * ( mass[1] * mass[1] ); | q2 = q2min; // pole | ||||
double Erho = ( m * m + m0 * m0 - q2 ) / ( 2.0 * m ), | |||||
Prho = sqrt( ( Erho - m0 ) * ( Erho + m0 ) ); | |||||
erho = ( m * m + mass[0] * mass[0] - q2 ) / ( 2.0 * m ); | p4meson.set( Erho, 0.0, 0.0, -Prho ); | ||||
p4w.set( m - Erho, 0.0, 0.0, Prho ); | |||||
prho = sqrt( erho * erho - mass[0] * mass[0] ); | |||||
p4meson.set( erho, 0.0, 0.0, -1.0 * prho ); | |||||
p4w.set( m - erho, 0.0, 0.0, prho ); | |||||
//This is in the W rest frame | //This is in the W rest frame | ||||
elepton = ( q2 + mass[1] * mass[1] ) / ( 2.0 * sqrt( q2 ) ); | double Elep = sqrt( q2 ) * 0.5, | ||||
plepton = sqrt( elepton * elepton - mass[1] * mass[1] ); | Plep = sqrt( Elep * Elep - mass[1] * mass[1] ); | ||||
double probctl[3]; | |||||
for ( j = 0; j < 3; j++ ) { | const int nj = 3 + 2 + 4 + 8; | ||||
costl = 0.99 * ( j - 1.0 ); | double p[nj], cmin = -1, cmax = 1, dc = ( cmax - cmin ) / ( nj - 1 ); | ||||
for ( int j = 0; j < nj; j++ ) { | |||||
double c = cmin + dc * j; | |||||
//These are in the W rest frame. Need to boost out into the B frame. | |||||
double Py = Plep * sqrt( 1.0 - c * c ), Pz = Plep * c; | |||||
p4lepton1.set( Elep, 0.0, Py, Pz ); | |||||
p4lepton2.set( Elep, 0.0, -Py, -Pz ); | |||||
//These are in the W rest frame. Need to boost out into | p4lepton1 = boostTo( p4lepton1, p4w ); | ||||
//the B frame. | p4lepton2 = boostTo( p4lepton2, p4w ); | ||||
p4lepton1.set( elepton, 0.0, | |||||
plepton * sqrt( 1.0 - costl * costl ), | |||||
plepton * costl ); | |||||
p4lepton2.set( elepton, 0.0, | |||||
-1.0 * plepton * sqrt( 1.0 - costl * costl ), | |||||
-1.0 * plepton * costl ); | |||||
EvtVector4R boost( ( m - erho ), 0.0, 0.0, 1.0 * prho ); | |||||
p4lepton1 = boostTo( p4lepton1, boost ); | |||||
p4lepton2 = boostTo( p4lepton2, boost ); | |||||
//Now initialize the daughters... | //Now initialize the daughters... | ||||
daughter->init( meson, p4meson ); | daughter->init( meson, p4meson ); | ||||
lep1->init( lepton1, p4lepton1 ); | lep1->init( lepton1, p4lepton1 ); | ||||
lep2->init( lepton2, p4lepton2 ); | lep2->init( lepton2, p4lepton2 ); | ||||
// std::cout<<j<<" kw "<<q2<<" "<<getKineWeight(root_part)<<" "<<getKineWeight2(root_part)<<std::endl; | |||||
CalcAmp( root_part, amp, FormFactors ); | CalcAmp( root_part, amp, FormFactors ); | ||||
//Now find the probability at this q2 and cos theta lepton point | //Now find the probability at this q2 and cos(theta) lepton point | ||||
//and compare to maxfoundprob. | //and compare to maxfoundprob. | ||||
//Do a little magic to get the probability!! | //Do a little magic to get the probability!! | ||||
//cout <<"amp:"<<amp.getSpinDensity()<<endl; | p[j] = rho.normalizedProb( amp.getSpinDensity() ); | ||||
// std::cout<<j<<" "<<p[j]<<std::endl; | |||||
prob = rho.normalizedProb( amp.getSpinDensity() ); | |||||
//cout << "prob:"<<q2<<" "<<costl<<" "<<prob<<endl; | |||||
probctl[j] = prob; | |||||
} | } | ||||
//probclt contains prob at ctl=-1,0,1. | for ( int j = 0; j < nj - 2; j++ ) { | ||||
//prob=a+b*ctl+c*ctl^2 | // pm,p0,pp contain probabilities at x=-1,0,1. | ||||
// prob = a + b*x + c*x^2 | |||||
double a = probctl[1]; | double pm = p[j + 0], p0 = p[j + 1], pp = p[j + 2]; | ||||
double b = 0.5 * ( probctl[2] - probctl[0] ); | double a = p0, b = 0.5 * ( pp - pm ), c = 0.5 * ( pp + pm ) - p0; | ||||
double c = 0.5 * ( probctl[2] + probctl[0] ) - probctl[1]; | prob = std::max( p0, std::max( pp, pm ) ); | ||||
prob = probctl[0]; | |||||
if ( probctl[1] > prob ) | |||||
prob = probctl[1]; | |||||
if ( probctl[2] > prob ) | |||||
prob = probctl[2]; | |||||
if ( fabs( c ) > 1e-20 ) { | if ( fabs( c ) > 1e-20 ) { | ||||
double ctlx = -0.5 * b / c; | double x = -0.5 * b / c; | ||||
if ( fabs( ctlx ) < 1.0 ) { | if ( fabs( x ) < 1.0 ) | ||||
double probtmp = a + b * ctlx + c * ctlx * ctlx; | prob = std::max( prob, a + b * x + c * x * x ); | ||||
if ( probtmp > prob ) | |||||
prob = probtmp; | |||||
} | } | ||||
} | } | ||||
//EvtGenReport(EVTGEN_DEBUG,"EvtGen") << "prob,probctl:"<<prob<<" " | |||||
// << probctl[0]<<" " | |||||
// << probctl[1]<<" " | |||||
// << probctl[2]<<endl; | |||||
if ( i == 0 ) { | if ( i == 0 ) { | ||||
maxpole = prob; | maxpole = prob; | ||||
continue; | continue; | ||||
} | } | ||||
if ( prob > maxfoundprob ) { | if ( prob > maxfoundprob ) { | ||||
maxfoundprob = prob; | maxfoundprob = prob; | ||||
} | } | ||||
//cout << "q2,maxfoundprob:"<<q2<<" "<<maxfoundprob<<endl; | |||||
} | |||||
if ( EvtPDL::getWidth( meson ) <= 0.0 ) { | |||||
//if the particle is narrow dont bother with changing the mass. | |||||
massiter = 4; | |||||
} | } | ||||
} | } | ||||
root_part->deleteTree(); | if ( 0 ) { | ||||
std::vector<float> v; | |||||
float x, y; | |||||
std::ifstream INR( "eRe_wide.dat" ); | |||||
while ( INR >> x >> y ) | |||||
v.push_back( x ); | |||||
std::cout << "q2 points read -- " << v.size() << std::endl; | |||||
std::ofstream OUT( "b2sllprob_amp.dat" ); | |||||
double m0 = mH[0]; | |||||
m0 = EvtPDL::getMinMass( meson ); | |||||
double q2max = ( m - m0 ) * ( m - m0 ); | |||||
//loop over q2 | |||||
for ( double q2 : v ) { | |||||
// want to avoid picking up the tail of the photon propagator | |||||
if ( !( q2min <= q2 && q2 < q2max ) ) | |||||
continue; | |||||
double Erho = ( m * m + m0 * m0 - q2 ) / ( 2.0 * m ), | |||||
Prho = sqrt( ( Erho - m0 ) * ( Erho + m0 ) ); | |||||
p4meson.set( Erho, 0.0, 0.0, -Prho ); | |||||
p4w.set( m - Erho, 0.0, 0.0, Prho ); | |||||
poleSize = 0.04 * ( maxpole / maxfoundprob ) * 4 * ( mass[1] * mass[1] ); | //This is in the W rest frame | ||||
double Elep = sqrt( q2 ) * 0.5, | |||||
Plep = sqrt( Elep * Elep - mass[1] * mass[1] ); | |||||
//poleSize=0.002; | const int nj = 3 + 2 + 4 + 8 + 32; | ||||
double cmin = -1, cmax = 1, dc = ( cmax - cmin ) / ( nj - 1 ); | |||||
double maxprob = 0; | |||||
for ( int j = 0; j < nj; j++ ) { | |||||
double c = cmin + dc * j; | |||||
//These are in the W rest frame. Need to boost out into the B frame. | |||||
double Py = Plep * sqrt( 1.0 - c * c ), Pz = Plep * c; | |||||
p4lepton1.set( Elep, 0.0, Py, Pz ); | |||||
p4lepton2.set( Elep, 0.0, -Py, -Pz ); | |||||
//cout <<"maxfoundprob,maxpole,poleSize:"<<maxfoundprob<<" " | p4lepton1 = boostTo( p4lepton1, p4w ); | ||||
// <<maxpole<<" "<<poleSize<<endl; | p4lepton2 = boostTo( p4lepton2, p4w ); | ||||
//Now initialize the daughters... | |||||
daughter->init( meson, p4meson ); | |||||
lep1->init( lepton1, p4lepton1 ); | |||||
lep2->init( lepton2, p4lepton2 ); | |||||
CalcAmp( root_part, amp, FormFactors ); | |||||
double prob = rho.normalizedProb( amp.getSpinDensity() ); | |||||
maxprob = std::max( maxprob, prob ); | |||||
} | |||||
OUT << q2 << " " << maxprob << "\n"; | |||||
} | |||||
} | |||||
root_part->deleteTree(); | |||||
maxfoundprob *= 1.15; | // poleSize = 0.04 * ( maxpole / maxfoundprob ) * q2min; | ||||
poleSize = ( maxpole / maxfoundprob ) * q2min; | |||||
std::cout << "max " << maxfoundprob << " " << maxpole << " " << poleSize | |||||
<< std::endl; | |||||
// maxfoundprob *= 1.5e-3; return maxfoundprob; // with resonances | |||||
// maxfoundprob *= 4e-7; return maxfoundprob; | |||||
// maxfoundprob *= 3.7e-4; return maxfoundprob; // with resonances | |||||
return maxfoundprob; | return maxfoundprob * 3; // without resonances | ||||
return maxfoundprob / 1000 * 3.5; | |||||
} | } | ||||
EvtComplex EvtbTosllAmp::GetC7Eff( double q2, bool nnlo ) | EvtComplex EvtbTosllAmp::GetC7Eff( double q2, bool nnlo ) | ||||
{ | { | ||||
if ( !nnlo ) | if ( !nnlo ) | ||||
return -0.313; | return -0.313; | ||||
double mbeff = 4.8; | double mbeff = 4.8; | ||||
double shat = q2 / mbeff / mbeff; | double shat = q2 / mbeff / mbeff; | ||||
▲ Show 20 Lines • Show All 437 Lines • Show Last 20 Lines |