Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MCGIDI_kinetics.cc File Reference
#include <string.h>
#include <cmath>
#include "MCGIDI.h"
Include dependency graph for MCGIDI_kinetics.cc:

Go to the source code of this file.

Functions

int MCGIDI_kinetics_2BodyReaction (statusMessageReporting *smr, MCGIDI_angular *angular, double K, double mu, double phi, MCGIDI_sampledProductsData *outgoingData)
 
int MCGIDI_kinetics_COMKineticEnergy2LabEnergyAndMomentum (statusMessageReporting *, double beta, double e_kinetic_com, double mu, double phi, double m3cc, double m4cc, MCGIDI_sampledProductsData *outgoingData)
 
int MCGIDI_kinetics_COM2Lab (statusMessageReporting *smr, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo, double masses[3])
 

Function Documentation

int MCGIDI_kinetics_2BodyReaction ( statusMessageReporting smr,
MCGIDI_angular angular,
double  K,
double  mu,
double  phi,
MCGIDI_sampledProductsData outgoingData 
)

Definition at line 18 of file MCGIDI_kinetics.cc.

19  {
20 
21  double m1 = angular->projectileMass_MeV, m2 = angular->targetMass_MeV, m3 = angular->productMass_MeV, m4 = angular->residualMass_MeV, mi, mf, Kp, x, beta;
22 
23  mi = m1 + m2;
24  mf = m3 + m4;
25  beta = std::sqrt( K * ( K + 2. * m1 ) ) / ( K + mi );
26  x = K * m2 / ( mi * mi );
27  if( x < 2e-5 ) { /* Kp is the total kinetic energy for m3 and m4 in the COM frame. */
28  Kp = mi - mf + K * m2 / mi * ( 1 - 0.5 * x * ( 1 - x ) ); }
29  else {
30  Kp = std::sqrt( mi * mi + 2 * K * m2 ) - mf;
31  }
32  if( Kp < 0 ) Kp = 0.; /* ???? There needs to be a better test here. */
33  return( MCGIDI_kinetics_COMKineticEnergy2LabEnergyAndMomentum( smr, beta, Kp, mu, phi, m3, m4, outgoingData ) );
34 }
double projectileMass_MeV
Definition: MCGIDI.h:323
int MCGIDI_kinetics_COMKineticEnergy2LabEnergyAndMomentum(statusMessageReporting *smr, double beta, double e_kinetic_com, double mu, double phi, double m3cc, double m4cc, MCGIDI_sampledProductsData *outgoingData)
static constexpr double m3
Definition: G4SIunits.hh:131
double residualMass_MeV
Definition: MCGIDI.h:323
double productMass_MeV
Definition: MCGIDI.h:323
double targetMass_MeV
Definition: MCGIDI.h:323
static constexpr double m2
Definition: G4SIunits.hh:130

Here is the call graph for this function:

Here is the caller graph for this function:

int MCGIDI_kinetics_COM2Lab ( statusMessageReporting smr,
MCGIDI_quantitiesLookupModes modes,
MCGIDI_decaySamplingInfo decaySamplingInfo,
double  masses[3] 
)

Definition at line 105 of file MCGIDI_kinetics.cc.

105  {
106 /*
107 * massProjectile = masses[0], massTarget = masses[1], massProduct = masses[2];
108 */
109  double a = masses[0] + masses[1], b, e_in = modes.getProjectileEnergy( ) * masses[0] * masses[2] / ( a * a ), Ep;
110 
111  if( decaySamplingInfo->frame != xDataTOM_frame_centerOfMass ) {
112  smr_setReportError2( smr, smr_unknownID, 1, "bad frame = %d for COM to lab conversion of mu/energy", decaySamplingInfo->frame );
113  return( 1 );
114  }
115  a = std::sqrt( e_in );
116  b = std::sqrt( decaySamplingInfo->Ep );
117  Ep = decaySamplingInfo->Ep + e_in + 2. * decaySamplingInfo->mu * a * b;
118  if( Ep != 0 ) {
119  decaySamplingInfo->mu = ( a + decaySamplingInfo->mu * b ) / std::sqrt( Ep );
120  }
121  decaySamplingInfo->Ep = Ep;
122  decaySamplingInfo->frame = xDataTOM_frame_lab;
123  return( 0 );
124 }
#define smr_setReportError2(smr, libraryID, code, fmt,...)
#define smr_unknownID
enum xDataTOM_frame frame
Definition: MCGIDI.h:256
double getProjectileEnergy(void) const
Definition: MCGIDI.h:97

Here is the call graph for this function:

Here is the caller graph for this function:

int MCGIDI_kinetics_COMKineticEnergy2LabEnergyAndMomentum ( statusMessageReporting ,
double  beta,
double  e_kinetic_com,
double  mu,
double  phi,
double  m3cc,
double  m4cc,
MCGIDI_sampledProductsData outgoingData 
)

Definition at line 38 of file MCGIDI_kinetics.cc.

39  {
40 /*
41 * beta the velocity/speedOflight of the com frame relative to the lab frame.
42 * e_kinetic_com Total kinetic energy (K1 + K2) in the COM frame.
43 * mu cos( theta ) in the COM frame.
44 */
45  double x, v_p, p, pp3, pp4, px3, py3, pz3, pz4, pz, p_perp2, E3, E4, gamma, m3cc2 = m3cc * m3cc, m4cc2 = m4cc * m4cc;
46 
47  p = std::sqrt( e_kinetic_com * ( e_kinetic_com + 2. * m3cc ) * ( e_kinetic_com + 2. * m4cc ) * ( e_kinetic_com + 2. * ( m3cc + m4cc ) ) ) /
48  ( 2. * ( e_kinetic_com + m3cc + m4cc ) );
49  py3 = p * std::sqrt( 1 - mu * mu );
50  px3 = py3 * std::cos( phi );
51  py3 *= std::sin( phi );
52  pz = p * mu;
53  if( 1 ) { /* ????????? Assuming the answer is wanted in the lab frame for now. */
54  E3 = std::sqrt( p * p + m3cc2 );
55  E4 = std::sqrt( p * p + m4cc2 );
56  gamma = std::sqrt( 1. / ( 1. - beta * beta ) );
57  pz3 = gamma * ( pz + beta * E3 );
58  pz4 = gamma * ( -pz + beta * E4 ); }
59  else { /* COM frame. */
60  pz3 = pz;
61  pz4 = -pz;
62  }
63  outgoingData[1].isVelocity = outgoingData[0].isVelocity;
64 
65  p_perp2 = px3 * px3 + py3 * py3;
66 
67  outgoingData[0].px_vx = px3;
68  outgoingData[0].py_vy = py3;
69  outgoingData[0].pz_vz = pz3;
70  pp3 = p_perp2 + pz3 * pz3;
71  x = ( m3cc > 0 ) ? pp3 / ( 2 * m3cc2 ) : 1.;
72  if( x < 1e-5 ) {
73  outgoingData[0].kineticEnergy = m3cc * x * ( 1 - 0.5 * x * ( 1 - x ) ); }
74  else {
75  outgoingData[0].kineticEnergy = std::sqrt( m3cc2 + pp3 ) - m3cc;
76  }
77  outgoingData[1].px_vx = -px3;
78  outgoingData[1].py_vy = -py3;
79  outgoingData[1].pz_vz = pz4;
80  pp4 = p_perp2 + pz4 * pz4;
81  x = ( m4cc > 0 ) ? pp4 / ( 2 * m4cc2 ) : 1.;
82  if( x < 1e-5 ) {
83  outgoingData[1].kineticEnergy = m4cc * x * ( 1 - 0.5 * x * ( 1 - x ) ); }
84  else {
85  outgoingData[1].kineticEnergy = std::sqrt( m4cc2 + pp4 ) - m4cc;
86  }
87 
88  if( outgoingData[0].isVelocity ) {
89  v_p = MCGIDI_speedOfLight_cm_sec / std::sqrt( pp3 + m3cc2 );
90  outgoingData[0].px_vx *= v_p;
91  outgoingData[0].py_vy *= v_p;
92  outgoingData[0].pz_vz *= v_p;
93 
94  v_p = MCGIDI_speedOfLight_cm_sec / std::sqrt( pp4 + m4cc2 );
95  outgoingData[1].px_vx *= v_p;
96  outgoingData[1].py_vy *= v_p;
97  outgoingData[1].pz_vz *= v_p;
98  }
99 
100  return( 0 );
101 }
#define MCGIDI_speedOfLight_cm_sec
Definition: MCGIDI.h:183
const char * p
Definition: xmltok.h:285

Here is the caller graph for this function: