Differential D36 Diff 156 src/EvtGenModels/EvtThreeBodyPhsp.cpp

Changeset View

Standalone View

src/EvtGenModels/EvtThreeBodyPhsp.cpp

Show All 14 Lines
* GNU General Public License for more details. *		* GNU General Public License for more details. *
* *		* *
* You should have received a copy of the GNU General Public License *		* You should have received a copy of the GNU General Public License *
* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *		* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
***********************************************************************/		***********************************************************************/

#include "EvtGenModels/EvtThreeBodyPhsp.hh"		#include "EvtGenModels/EvtThreeBodyPhsp.hh"

#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtConst.hh"		#include "EvtGenBase/EvtConst.hh"
		#include "EvtGenBase/EvtGenKine.hh"
		#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtRandom.hh"		#include "EvtGenBase/EvtRandom.hh"
		#include "EvtGenBase/EvtReport.hh"

#include <iostream>
#include <algorithm>		#include <algorithm>
#include <cmath>		#include <cmath>
		#include <iostream>

std::string EvtThreeBodyPhsp::getName()		std::string EvtThreeBodyPhsp::getName()
{		{
return "THREEBODYPHSP";		return "THREEBODYPHSP";
}		}

EvtDecayBase* EvtThreeBodyPhsp::clone()		EvtDecayBase* EvtThreeBodyPhsp::clone()
{		{
Show All 23 Lines	void EvtThreeBodyPhsp::initProbMax()
noProbMax();		noProbMax();
}		}

void EvtThreeBodyPhsp::decay( EvtParticle* p )		void EvtThreeBodyPhsp::decay( EvtParticle* p )
{		{
p->makeDaughters( getNDaug(), getDaugs() );		p->makeDaughters( getNDaug(), getDaugs() );
p->generateMassTree();		p->generateMassTree();
const double mParent = p->mass();		const double mParent = p->mass();
EvtParticle* daug1 = p->getDaug(0);		EvtParticle* daug1 = p->getDaug( 0 );
EvtParticle* daug2 = p->getDaug(1);		EvtParticle* daug2 = p->getDaug( 1 );
EvtParticle* daug3 = p->getDaug(2);		EvtParticle* daug3 = p->getDaug( 2 );
const double mDaug1 = daug1->mass();		const double mDaug1 = daug1->mass();
const double mDaug2 = daug2->mass();		const double mDaug2 = daug2->mass();
const double mDaug3 = daug3->mass();		const double mDaug3 = daug3->mass();

const double m12borderMin = mDaug1 + mDaug2;		const double m12borderMin = mDaug1 + mDaug2;
const double m12borderMax = mParent - mDaug3;		const double m12borderMax = mParent - mDaug3;
const double m12Min = std::max( m_m12SqMin, m12borderMin * m12borderMin );		const double m12Min = std::max( m_m12SqMin, m12borderMin * m12borderMin );
const double m12Max = std::min( m_m12SqMax, m12borderMax * m12borderMax );		const double m12Max = std::min( m_m12SqMax, m12borderMax * m12borderMax );

const double m23borderMin = mDaug2 + mDaug3;		const double m23borderMin = mDaug2 + mDaug3;
Show All 27 Lines	void 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;
}		}