#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) |
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_decayChannel_getFirstProduct(), tpia_decayChannel_getNextProduct(), and tpia_kinetics_COMKineticEnergy2LabEnergyAndMomentum().
Referenced by tpia_decayChannel_sampleProductsAtE().
00049 { 00050 00051 tpia_product *pp3 = tpia_decayChannel_getFirstProduct( decayChannel ), *pp4; 00052 double m1 = decayChannel->m1_fullMass_MeV, m2 = decayChannel->m2_fullMass_MeV, m3, m4, mi, mf, Kp, x, beta; 00053 00054 pp4 = tpia_decayChannel_getNextProduct( pp3 ); 00055 m3 = pp3->productID->fullMass_MeV; 00056 m4 = pp4->productID->fullMass_MeV; 00057 mi = m1 + m2; 00058 mf = m3 + m4; 00059 beta = std::sqrt( K * ( K + 2. * m1 ) ) / ( K + mi ); 00060 x = K * m2 / ( mi * mi ); 00061 if( x < 2e-5 ) { /* Kp is the total kinetic energy for m3 and m4 in the COM frame. */ 00062 Kp = mi - mf + K * m2 / mi * ( 1 - 0.5 * x * ( 1 - x ) ); } 00063 else { 00064 Kp = std::sqrt( mi * mi + 2 * K * m2 ) - mf; 00065 } 00066 if( Kp < 0 ) Kp = 0.; /* ???? There needs to be a better test here. */ 00067 outgoingData[0].decayChannel = &(pp3->decayChannel); 00068 outgoingData[1].genre = outgoingData[0].genre; 00069 outgoingData[1].productID = pp4->productID; 00070 outgoingData[1].decayChannel = &(pp4->decayChannel); 00071 return( tpia_kinetics_COMKineticEnergy2LabEnergyAndMomentum( smr, beta, Kp, mu, phi, m3, m4, outgoingData ) ); 00072 }
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_frame_getColumn().
Referenced by tpia_kinetics_2BodyReaction().
00078 { 00079 /* 00080 * beta the velocity/speedOflight of the com frame relative to the lab frame. 00081 * e_kinetic_com Total kinetic energy (K1 + K2) in the COM frame. 00082 * mu std::cos( theta ) in the COM frame. 00083 */ 00084 double x, v_p, p, pp3, pp4, px3, py3, pz3, pz4, pz, p_perp2, E3, E4, gamma, m3cc2 = m3cc * m3cc, m4cc2 = m4cc * m4cc; 00085 00086 p = std::sqrt( e_kinetic_com * ( e_kinetic_com + 2. * m3cc ) * ( e_kinetic_com + 2. * m4cc ) * ( e_kinetic_com + 2. * ( m3cc + m4cc ) ) ) / 00087 ( 2. * ( e_kinetic_com + m3cc + m4cc ) ); 00088 py3 = p * std::sqrt( 1 - mu * mu ); 00089 px3 = py3 * std::cos( phi ); 00090 py3 *= std::sin( phi ); 00091 pz = p * mu; 00092 if( tpia_frame_getColumn( NULL, &(outgoingData[0].frame), 0 ) == tpia_referenceFrame_lab ) { 00093 E3 = std::sqrt( p * p + m3cc2 ); 00094 E4 = std::sqrt( p * p + m4cc2 ); 00095 gamma = std::sqrt( 1. / ( 1. - beta * beta ) ); 00096 pz3 = gamma * ( pz + beta * E3 ); 00097 pz4 = gamma * ( -pz + beta * E4 ); } 00098 else { 00099 pz3 = pz; 00100 pz4 = -pz; 00101 } 00102 outgoingData[1].isVelocity = outgoingData[0].isVelocity; 00103 outgoingData[1].frame = outgoingData[0].frame; 00104 00105 p_perp2 = px3 * px3 + py3 * py3; 00106 00107 outgoingData[0].px_vx = px3; 00108 outgoingData[0].py_vy = py3; 00109 outgoingData[0].pz_vz = pz3; 00110 pp3 = p_perp2 + pz3 * pz3; 00111 x = pp3 / ( 2 * m3cc2 ); 00112 if( x < 1e-5 ) { 00113 outgoingData[0].kineticEnergy = m3cc * x * ( 1 - 0.5 * x * ( 1 - x ) ); } 00114 else { 00115 outgoingData[0].kineticEnergy = std::sqrt( m3cc2 + pp3 ) - m3cc; 00116 } 00117 outgoingData[1].px_vx = -px3; 00118 outgoingData[1].py_vy = -py3; 00119 outgoingData[1].pz_vz = pz4; 00120 pp4 = p_perp2 + pz4 * pz4; 00121 x = pp4 / ( 2 * m4cc2 ); 00122 if( x < 1e-5 ) { 00123 outgoingData[1].kineticEnergy = m4cc * x * ( 1 - 0.5 * x * ( 1 - x ) ); } 00124 else { 00125 outgoingData[1].kineticEnergy = std::sqrt( m4cc2 + pp4 ) - m4cc; 00126 } 00127 00128 if( outgoingData[0].isVelocity ) { 00129 v_p = tpia_speedOfLight_cm_sec / std::sqrt( pp3 + m3cc2 ); 00130 outgoingData[0].px_vx *= v_p; 00131 outgoingData[0].py_vy *= v_p; 00132 outgoingData[0].pz_vz *= v_p; 00133 00134 v_p = tpia_speedOfLight_cm_sec / std::sqrt( pp4 + m4cc2 ); 00135 outgoingData[1].px_vx *= v_p; 00136 outgoingData[1].py_vy *= v_p; 00137 outgoingData[1].pz_vz *= v_p; 00138 } 00139 00140 return( 0 ); 00141 }