Page MenuHomeHEPForge
Differential D36 Diff 134 src/EvtGenModels/EvtThreeBodyPhsp.cpp

Changeset View

Standalone View

View Options

src/EvtGenModels/EvtThreeBodyPhsp.cpp

Show All 18 Lines
***********************************************************************/ ***********************************************************************/
#include "EvtGenModels/EvtThreeBodyPhsp.hh" #include "EvtGenModels/EvtThreeBodyPhsp.hh"
#include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtReport.hh" #include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtConst.hh" #include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtRandom.hh" #include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include <iostream> #include <iostream>
#include <algorithm> #include <algorithm>
#include <cmath> #include <cmath>
std::string EvtThreeBodyPhsp::getName() std::string EvtThreeBodyPhsp::getName()
{ {
return "THREEBODYPHSP"; return "THREEBODYPHSP";
▲ Show 20 Lines • Show All 75 Lines • ▼ Show 20 Linesvoid EvtThreeBodyPhsp::decay( EvtParticle* p )
} while ( goodEvent == false && iTrial < nMaxTrials ); } while ( goodEvent == false && iTrial < nMaxTrials );
if ( !goodEvent ) { if ( !goodEvent ) {
EvtGenReport( EVTGEN_WARNING, "EvtThreeBodyPhsp" ) EvtGenReport( EVTGEN_WARNING, "EvtThreeBodyPhsp" )
<< "Failed to generate m12Sq and m23Sq. Taking last m12Sq and midpoint of allowed m23Sq.\n"; << "Failed to generate m12Sq and m23Sq. Taking last m12Sq and midpoint of allowed m23Sq.\n";
m23Sq = 0.5 * ( m23LowLimit + m23HighLimit ); m23Sq = 0.5 * ( m23LowLimit + m23HighLimit );
} }
// At this moment we have valid invariant masses squared // At this moment we have valid invariant masses squared
const double En1 = 0.5 * ( mParent * mParent + mDaug1 * mDaug1 - m23Sq ) / mParent; EvtGenKine::ThreeBodyKine(m12Sq, m23Sq, p);
const double En3 = 0.5 * ( mParent * mParent + mDaug3 * mDaug3 - m12Sq ) / mParent;
const double En2 = mParent - En1 - En3;
const double p1mag = std::sqrt( En1 * En1 - mDaug1 * mDaug1 );
const double p2mag = std::sqrt( En2 * En2 - mDaug2 * mDaug2 );
double cosPhi = 0.5 * ( mDaug1 * mDaug1 + mDaug2 * mDaug2 +
2 * En1 * En2 - m12Sq )/p1mag/p2mag;
if ( cosPhi > 1.0 ) {
EvtGenReport( EVTGEN_WARNING, "EvtThreeBodyPhsp" )
<< " Model: cosine larger than one: " << cosPhi << std::endl;
cosPhi = 1.0;
}
double sinPhi = std::sqrt( 1 - cosPhi * cosPhi );
if ( EvtRandom::Flat( 0., 1. ) > 0.5 ) {
sinPhi *= -1;
}
const double p2x = p2mag * cosPhi;
const double p2y = p2mag * sinPhi;
const double p3x = -p1mag - p2x;
const double p3y = -p2y;
// Construct 4-momenta and rotate them randomly in space
EvtVector4R p1( En1, p1mag, 0., 0. );
EvtVector4R p2( En2, p2x, p2y, 0. );
EvtVector4R p3( En3, p3x, p3y, 0. );
const double euler1 = EvtRandom::Flat( 0., EvtConst::twoPi );
const double euler2 = std::acos( EvtRandom::Flat( -1.0, 1.0 ) );
const double euler3 = EvtRandom::Flat( 0., EvtConst::twoPi );
p1.applyRotateEuler(euler1, euler2, euler3);
p2.applyRotateEuler(euler1, euler2, euler3);
p3.applyRotateEuler(euler1, euler2, euler3);
// set momenta for daughters
daug1->init( getDaug( 0 ), p1 );
daug2->init( getDaug( 1 ), p2 );
daug3->init( getDaug( 2 ), p3 );
return; return;
} }