Page MenuHomeHEPForge
Differential D87 Diff 389 src/EvtGenModels/EvtPi0Dalitz.cpp

Changeset View

Standalone View

View Options

src/EvtGenModels/EvtPi0Dalitz.cpp

Show First 20 Lines • Show All 46 Lines • ▼ Show 20 Lines
void EvtPi0Dalitz::initProbMax() void EvtPi0Dalitz::initProbMax()
{ {
// Search for maximum probability. In order to avoid setting up all // Search for maximum probability. In order to avoid setting up all
// particles with spinors and four-momenta, we use result after all // particles with spinors and four-momenta, we use result after all
// contractions, which is: // contractions, which is:
// 1/((m_R^2-q^2)^2+m_R^2 Gamma_R^2) 1/(2q^2) (M^2-q^2)^2 beta_l // 1/((m_R^2-q^2)^2+m_R^2 Gamma_R^2) 1/(2q^2) (M^2-q^2)^2 beta_l
// (1+cos(theta)^2) where we set cos(theta)=1 // (1+cos(theta)^2) where we set cos(theta)=1
auto daughter1 = getDaug( 0 ); double q2Min = EvtPDL::getMass( getDaug( 0 ) ) +
auto daughter2 = getDaug( 1 ); EvtPDL::getMass( getDaug( 1 ) );
double q2Min = EvtPDL::getMass( daughter1 ) + EvtPDL::getMass( daughter2 );
q2Min *= q2Min; q2Min *= q2Min;
double q2Max = EvtPDL::getMass( getParentId() ); double q2Max = EvtPDL::getMass( getParentId() );
q2Max *= q2Max; q2Max *= q2Max;
const int steps = 20000; const int steps = 100;
const double step = ( q2Max - q2Min ) / steps; const double step = expm1( log( q2Max / q2Min ) / steps );
double maxProb = 0; double maxProb = 0, q2 = q2Min;
for ( int ii = 0; ii < steps; ++ii ) { do {
double q2 = q2Min + ii * step; double prob = ( q2Max - q2 ) * ( q2Max - q2 ) * ( q2 - q2Min );
const double mSqDiff = m_m0Sq - q2; double m2d = m_m0Sq - q2;
const double q2Sq = q2 * q2; prob /= q2 * ( m2d * m2d + m_m0SqG0Sq ) * ( q2 + m_poleSize );
double prob = ( q2Max - q2 ) * ( q2Max - q2 ) * ( q2 - q2Min ) / ( q2Sq ); if ( prob > maxProb )
prob *= ( 1.0 / ( mSqDiff * mSqDiff + m_m0SqG0Sq ) );
// When generating events, we do not start from phase-space, but
// add some pole to it, weight of which is taken into account
// elsewhere
prob /= 1.0 + m_poleSize / ( q2Sq );
if ( prob > maxProb ) {
maxProb = prob; maxProb = prob;
} } while ( ( q2 += q2 * step ) < q2Max );
} setProbMax( maxProb * 1.01 );
setProbMax( maxProb * 1.05 ); EvtGenReport( EVTGEN_INFO, "EvtGenProbMax" )
<< " " << sqrt( q2 ) << " " << maxProb << std::endl;
} }
void EvtPi0Dalitz::init() void EvtPi0Dalitz::init()
{ {
// check that there are 0 arguments // check that there are 0 arguments
checkNArg( 0 ); checkNArg( 0 );
checkNDaug( 3 ); checkNDaug( 3 );
checkSpinParent( EvtSpinType::SCALAR ); checkSpinParent( EvtSpinType::SCALAR );
checkSpinDaughter( 0, EvtSpinType::DIRAC ); checkSpinDaughter( 0, EvtSpinType::DIRAC );
checkSpinDaughter( 1, EvtSpinType::DIRAC ); checkSpinDaughter( 1, EvtSpinType::DIRAC );
checkSpinDaughter( 2, EvtSpinType::PHOTON ); checkSpinDaughter( 2, EvtSpinType::PHOTON );
// Rescale pole size to improve efficiency. Not sure about exact // Rescale pole size to improve efficiency. Not sure about exact
// factor, but this seem to be best simple rescaling for // factor, but this seem to be best simple rescaling for
// eta-->e+e-gamma. // eta-->e+e-gamma.
const double parentMass = EvtPDL::getMass( getParentId() ); const double parentMass = EvtPDL::getMass( getParentId() );
m_poleSize *= parentMass * parentMass / ( 0.135 * 0.135 ); m_poleSize *= parentMass * parentMass / ( 0.135 * 0.135 );
} }
void EvtPi0Dalitz::decay( EvtParticle* p ) void EvtPi0Dalitz::decay( EvtParticle* p )
{ {
EvtParticle *ep, *em, *gamma; using EvtGenFunctions::directProd;
setWeight( p->initializePhaseSpace( getNDaug(), getDaugs(), false, double ew = p->initializePhaseSpace( getNDaug(), getDaugs(), false,
m_poleSize, 0, 1 ) ); m_poleSize, 0, 1 );
ep = p->getDaug( 0 ); setWeight( ew );
em = p->getDaug( 1 );
gamma = p->getDaug( 2 ); const EvtVector4R& pp = p->getDaug( 0 )->getP4();
const EvtVector4R& pm = p->getDaug( 1 )->getP4();
// the next four lines generates events with a weight such that const EvtVector4R& pg = p->getDaug( 2 )->getP4();
// the efficiency for selecting them is good. The parameter below of EvtVector4R q = pp + pm;
// 0.1 is the size of the peak at low q^2 (in arbitrary units). double q2 = q.mass2(), t = pg * q;
// The value of 0.1 is appropriate for muons.
// when you use this remember to remove the cut on q^2!
//ep em invariant mass^2
double m2 = ( ep->getP4() + em->getP4() ).mass2();
EvtVector4R q = ep->getP4() + em->getP4();
//Just use the prob summed over spins...
EvtTensor4C w, v; //Just use the prob summed over spins...
//EvtTensor4C v = (2.0 * t) * directProd( q, pg ) - t2 * EvtTensor4C::g() - q2 * directProd( pg, pg );
//EvtTensor4C w = 4.0 * ( directProd( pp, pm ) + directProd( pm, pp ) - EvtTensor4C::g() * ( pp*pm - pp.mass2() ) );
//double prob = real( cont( v, w ) ) / ( q2 * q2 );
EvtTensor4C v = directProd( q, pg )
.addScaled( -0.5 * t, EvtTensor4C::g() )
.addScaled( -0.5 * q2 / t, directProd( pg, pg ) );
EvtTensor4C w = directProd( pp, pm ).addDirProd( pm, pp ).addScaled(
pp.mass2() - pp * pm, EvtTensor4C::g() );
double prob = 8 * t * real( cont( v, w ) ) / ( q2 * q2 );
double m2d = m_m0Sq - q2;
prob /= m2d * m2d + m_m0SqG0Sq;
v = 2.0 * ( gamma->getP4() * q ) * // EvtGenReport(EVTGEN_INFO,"EvtGenDDD") <<" "<<q2<<" "<<prob<<" "<<ew<<std::endl;
EvtGenFunctions::directProd( q, gamma->getP4() ) -
( gamma->getP4() * q ) * ( gamma->getP4() * q ) * EvtTensor4C::g() -
m2 * EvtGenFunctions::directProd( gamma->getP4(), gamma->getP4() );
w = 4.0 * ( EvtGenFunctions::directProd( ep->getP4(), em->getP4() ) +
EvtGenFunctions::directProd( em->getP4(), ep->getP4() ) -
EvtTensor4C::g() *
( ep->getP4() * em->getP4() - ep->getP4().mass2() ) );
double prob = ( real( cont( v, w ) ) ) / ( m2 * m2 );
const double m2Diff = m_m0Sq - m2;
prob *= ( 1.0 / ( m2Diff * m2Diff + m_m0SqG0Sq ) );
// EvtGenReport(EVTGEN_INFO,"EvtGen") << "prob is "<<prob<<endl;
setProb( prob ); setProb( prob );
return; return;
} }