Geant4 Cross Reference |
1 /* 1 2 # <<BEGIN-copyright>> 3 # <<END-copyright>> 4 */ 5 #include <string.h> 6 #include <cmath> 7 8 #include "MCGIDI.h" 9 10 #if defined __cplusplus 11 namespace GIDI { 12 using namespace GIDI; 13 #endif 14 15 /* 16 ********************************************** 17 */ 18 int MCGIDI_kinetics_2BodyReaction( statusMessa 19 MCGIDI_sampledProductsData *outgoingDa 20 21 double m1 = angular->projectileMass_MeV, m 22 23 mi = m1 + m2; 24 mf = m3 + m4; 25 beta = std::sqrt( K * ( K + 2. * m1 ) ) / 26 x = K * m2 / ( mi * mi ); 27 if( x < 2e-5 ) { 28 Kp = mi - mf + K * m2 / mi * ( 1 - 0.5 29 else { 30 Kp = std::sqrt( mi * mi + 2 * K * m2 ) 31 } 32 if( Kp < 0 ) Kp = 0.; /* ???? Th 33 return( MCGIDI_kinetics_COMKineticEnergy2L 34 } 35 /* 36 ********************************************** 37 */ 38 int MCGIDI_kinetics_COMKineticEnergy2LabEnergy 39 double m3cc, double m4cc, MCGIDI_sampl 40 /* 41 * beta the velocity/speedOflight 42 * e_kinetic_com Total kinetic energy (K1 + 43 * mu cos( theta ) in the COM fr 44 */ 45 double x, v_p, p, pp3, pp4, px3, py3, pz3, 46 47 p = std::sqrt( e_kinetic_com * ( e_kinetic 48 ( 2. * ( e_kinetic_com + m3cc + m4 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 ) { /* ??? 54 E3 = std::sqrt( p * p + m3cc2 ); 55 E4 = std::sqrt( p * p + m4cc2 ); 56 gamma = std::sqrt( 1. / ( 1. - beta * 57 pz3 = gamma * ( pz + beta * E3 ); 58 pz4 = gamma * ( -pz + beta * E4 ); } 59 else { /* COM 60 pz3 = pz; 61 pz4 = -pz; 62 } 63 outgoingData[1].isVelocity = outgoingData[ 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 * 74 else { 75 outgoingData[0].kineticEnergy = std::s 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 * 84 else { 85 outgoingData[1].kineticEnergy = std::s 86 } 87 88 if( outgoingData[0].isVelocity ) { 89 v_p = MCGIDI_speedOfLight_cm_sec / std 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 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 } 102 /* 103 ********************************************** 104 */ 105 int MCGIDI_kinetics_COM2Lab( statusMessageRepo 106 /* 107 * massProjectile = masses[0], massTarget 108 */ 109 double a = masses[0] + masses[1], b, e_in 110 111 if( decaySamplingInfo->frame != xDataTOM_f 112 smr_setReportError2( smr, smr_unknownI 113 return( 1 ); 114 } 115 a = std::sqrt( e_in ); 116 b = std::sqrt( decaySamplingInfo->Ep ); 117 Ep = decaySamplingInfo->Ep + e_in + 2. * d 118 if( Ep != 0 ) { 119 decaySamplingInfo->mu = ( a + decaySam 120 } 121 decaySamplingInfo->Ep = Ep; 122 decaySamplingInfo->frame = xDataTOM_frame_ 123 return( 0 ); 124 } 125 126 #if defined __cplusplus 127 } 128 #endif 129