tpia_kinetics.cc File Reference

#include <string.h>
#include <cmath>
#include <tpia_target.h>

Go to the source code of this file.

Functions
int	tpia_kinetics_2BodyReaction (statusMessageReporting *smr, tpia_decayChannel *decayChannel, double K, double mu, double phi, tpia_productOutgoingData *outgoingData)
int	tpia_kinetics_COMKineticEnergy2LabEnergyAndMomentum (statusMessageReporting *, double beta, double e_kinetic_com, double mu, double phi, double m3cc, double m4cc, tpia_productOutgoingData *outgoingData)

Function Documentation

int tpia_kinetics_2BodyReaction	(	statusMessageReporting *	smr,
		tpia_decayChannel *	decayChannel,
		double	K,
		double	mu,
		double	phi,
		tpia_productOutgoingData *	outgoingData
	)

Definition at line 48 of file tpia_kinetics.cc.

References tpia_productOutgoingData_s::decayChannel, tpia_product_s::decayChannel, tpia_particle_s::fullMass_MeV, tpia_productOutgoingData_s::genre, tpia_decayChannel_s::m1_fullMass_MeV, python.hepunit::m2, tpia_decayChannel_s::m2_fullMass_MeV, python.hepunit::m3, tpia_productOutgoingData_s::productID, tpia_product_s::productID, tpia_decayChannel_getFirstProduct(), tpia_decayChannel_getNextProduct(), tpia_kinetics_COMKineticEnergy2LabEnergyAndMomentum(), and test::x.

Referenced by tpia_decayChannel_sampleProductsAtE().

                                                 {

    tpia_product *pp3 = tpia_decayChannel_getFirstProduct( decayChannel ), *pp4;
    double m1 = decayChannel->m1_fullMass_MeV, m2 = decayChannel->m2_fullMass_MeV, m3, m4, mi, mf, Kp, x, beta;

    pp4 = tpia_decayChannel_getNextProduct( pp3 );
    m3 = pp3->productID->fullMass_MeV;
    m4 = pp4->productID->fullMass_MeV;
    mi = m1 + m2;
    mf = m3 + m4;
    beta = std::sqrt( K * ( K + 2. * m1 ) ) / ( K + mi );
    x = K * m2 / ( mi * mi );
    if( x < 2e-5 ) {                                        /* Kp is the total kinetic energy for m3 and m4 in the COM frame. */
        Kp = mi - mf + K * m2 / mi * ( 1 - 0.5 * x * ( 1 - x ) ); }
    else {
        Kp = std::sqrt( mi * mi + 2 * K * m2 ) - mf;
    }
    if( Kp < 0 ) Kp = 0.;           /* ???? There needs to be a better test here. */
    outgoingData[0].decayChannel = &(pp3->decayChannel);
    outgoingData[1].genre = outgoingData[0].genre;
    outgoingData[1].productID = pp4->productID;
    outgoingData[1].decayChannel = &(pp4->decayChannel);
    return( tpia_kinetics_COMKineticEnergy2LabEnergyAndMomentum( smr, beta, Kp, mu, phi, m3, m4, outgoingData ) );
}

int tpia_kinetics_COMKineticEnergy2LabEnergyAndMomentum	(	statusMessageReporting *	,
		double	beta,
		double	e_kinetic_com,
		double	mu,
		double	phi,
		double	m3cc,
		double	m4cc,
		tpia_productOutgoingData *	outgoingData
	)

Definition at line 77 of file tpia_kinetics.cc.

References tpia_productOutgoingData_s::frame, tpia_productOutgoingData_s::isVelocity, tpia_productOutgoingData_s::kineticEnergy, tpia_productOutgoingData_s::px_vx, tpia_productOutgoingData_s::py_vy, tpia_productOutgoingData_s::pz_vz, tpia_frame_getColumn(), tpia_referenceFrame_lab, tpia_speedOfLight_cm_sec, and test::x.

Referenced by tpia_kinetics_2BodyReaction().

                                                                           {
/*
*   beta            the velocity/speedOflight of the com frame relative to the lab frame.
*   e_kinetic_com   Total kinetic energy (K1 + K2) in the COM frame.
*   mu              std::cos( theta ) in the COM frame.
*/
    double x, v_p, p, pp3, pp4, px3, py3, pz3, pz4, pz, p_perp2, E3, E4, gamma, m3cc2 = m3cc * m3cc, m4cc2 = m4cc * m4cc;

    p = std::sqrt( e_kinetic_com * ( e_kinetic_com + 2. * m3cc ) * ( e_kinetic_com + 2. * m4cc )  * ( e_kinetic_com + 2. * ( m3cc + m4cc ) ) ) /
            ( 2. * ( e_kinetic_com + m3cc + m4cc ) );
    py3 = p * std::sqrt( 1 - mu * mu );
    px3 = py3 * std::cos( phi );
    py3 *= std::sin( phi );
    pz = p * mu;
    if( tpia_frame_getColumn( NULL, &(outgoingData[0].frame), 0 ) == tpia_referenceFrame_lab ) {
        E3 = std::sqrt( p * p + m3cc2 );
        E4 = std::sqrt( p * p + m4cc2 );
        gamma = std::sqrt( 1. / ( 1. - beta * beta ) );
        pz3 = gamma * (  pz + beta * E3 );
        pz4 = gamma * ( -pz + beta * E4 ); }
    else {
        pz3 = pz;
        pz4 = -pz;
    }
    outgoingData[1].isVelocity = outgoingData[0].isVelocity;
    outgoingData[1].frame = outgoingData[0].frame;

    p_perp2 = px3 * px3 + py3 * py3;

    outgoingData[0].px_vx = px3;
    outgoingData[0].py_vy = py3;
    outgoingData[0].pz_vz = pz3;
    pp3 = p_perp2 + pz3 * pz3;
    x = pp3 / ( 2 * m3cc2 );
    if( x < 1e-5 ) {
        outgoingData[0].kineticEnergy = m3cc * x  * ( 1 - 0.5 * x * ( 1 - x ) ); }
    else {
        outgoingData[0].kineticEnergy = std::sqrt( m3cc2 + pp3 ) - m3cc;
    }
    outgoingData[1].px_vx = -px3;
    outgoingData[1].py_vy = -py3;
    outgoingData[1].pz_vz = pz4;
    pp4 = p_perp2 + pz4 * pz4;
    x = pp4 / ( 2 * m4cc2 );
    if( x < 1e-5 ) {
        outgoingData[1].kineticEnergy = m4cc * x  * ( 1 - 0.5 * x * ( 1 - x ) ); }
    else {
        outgoingData[1].kineticEnergy = std::sqrt( m4cc2 + pp4 ) - m4cc;
    }

    if( outgoingData[0].isVelocity ) {
        v_p = tpia_speedOfLight_cm_sec / std::sqrt( pp3 + m3cc2 );
        outgoingData[0].px_vx *= v_p;
        outgoingData[0].py_vy *= v_p;
        outgoingData[0].pz_vz *= v_p;

        v_p = tpia_speedOfLight_cm_sec / std::sqrt( pp4 + m4cc2 );
        outgoingData[1].px_vx *= v_p;
        outgoingData[1].py_vy *= v_p;
        outgoingData[1].pz_vz *= v_p;
    }

    return( 0 );
}