| Geant4 Cross Reference |
1 /* 1 /*
2 # <<BEGIN-copyright>> 2 # <<BEGIN-copyright>>
3 # <<END-copyright>> 3 # <<END-copyright>>
4 */ 4 */
5 #include <string.h> 5 #include <string.h>
6 #include <cmath> 6 #include <cmath>
7 7
8 #include "MCGIDI.h" 8 #include "MCGIDI.h"
9 9
10 #if defined __cplusplus 10 #if defined __cplusplus
11 namespace GIDI { 11 namespace GIDI {
12 using namespace GIDI; 12 using namespace GIDI;
13 #endif 13 #endif
14 14
15 /* 15 /*
16 ********************************************** 16 ************************************************************
17 */ 17 */
18 int MCGIDI_kinetics_2BodyReaction( statusMessa 18 int MCGIDI_kinetics_2BodyReaction( statusMessageReporting *smr, MCGIDI_angular *angular, double K, double mu, double phi,
19 MCGIDI_sampledProductsData *outgoingDa 19 MCGIDI_sampledProductsData *outgoingData ) {
20 20
21 double m1 = angular->projectileMass_MeV, m 21 double m1 = angular->projectileMass_MeV, m2 = angular->targetMass_MeV, m3 = angular->productMass_MeV, m4 = angular->residualMass_MeV, mi, mf, Kp, x, beta;
22 22
23 mi = m1 + m2; 23 mi = m1 + m2;
24 mf = m3 + m4; 24 mf = m3 + m4;
25 beta = std::sqrt( K * ( K + 2. * m1 ) ) / 25 beta = std::sqrt( K * ( K + 2. * m1 ) ) / ( K + mi );
26 x = K * m2 / ( mi * mi ); 26 x = K * m2 / ( mi * mi );
27 if( x < 2e-5 ) { 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 28 Kp = mi - mf + K * m2 / mi * ( 1 - 0.5 * x * ( 1 - x ) ); }
29 else { 29 else {
30 Kp = std::sqrt( mi * mi + 2 * K * m2 ) 30 Kp = std::sqrt( mi * mi + 2 * K * m2 ) - mf;
31 } 31 }
32 if( Kp < 0 ) Kp = 0.; /* ???? Th 32 if( Kp < 0 ) Kp = 0.; /* ???? There needs to be a better test here. */
33 return( MCGIDI_kinetics_COMKineticEnergy2L 33 return( MCGIDI_kinetics_COMKineticEnergy2LabEnergyAndMomentum( smr, beta, Kp, mu, phi, m3, m4, outgoingData ) );
34 } 34 }
35 /* 35 /*
36 ********************************************** 36 ************************************************************
37 */ 37 */
38 int MCGIDI_kinetics_COMKineticEnergy2LabEnergy 38 int MCGIDI_kinetics_COMKineticEnergy2LabEnergyAndMomentum( statusMessageReporting * /*smr*/, double beta, double e_kinetic_com, double mu, double phi,
39 double m3cc, double m4cc, MCGIDI_sampl 39 double m3cc, double m4cc, MCGIDI_sampledProductsData *outgoingData ) {
40 /* 40 /*
41 * beta the velocity/speedOflight 41 * beta the velocity/speedOflight of the com frame relative to the lab frame.
42 * e_kinetic_com Total kinetic energy (K1 + 42 * e_kinetic_com Total kinetic energy (K1 + K2) in the COM frame.
43 * mu cos( theta ) in the COM fr 43 * mu cos( theta ) in the COM frame.
44 */ 44 */
45 double x, v_p, p, pp3, pp4, px3, py3, pz3, 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 46
47 p = std::sqrt( e_kinetic_com * ( e_kinetic 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 + m4 48 ( 2. * ( e_kinetic_com + m3cc + m4cc ) );
49 py3 = p * std::sqrt( 1 - mu * mu ); 49 py3 = p * std::sqrt( 1 - mu * mu );
50 px3 = py3 * std::cos( phi ); 50 px3 = py3 * std::cos( phi );
51 py3 *= std::sin( phi ); 51 py3 *= std::sin( phi );
52 pz = p * mu; 52 pz = p * mu;
53 if( 1 ) { /* ??? 53 if( 1 ) { /* ????????? Assuming the answer is wanted in the lab frame for now. */
54 E3 = std::sqrt( p * p + m3cc2 ); 54 E3 = std::sqrt( p * p + m3cc2 );
55 E4 = std::sqrt( p * p + m4cc2 ); 55 E4 = std::sqrt( p * p + m4cc2 );
56 gamma = std::sqrt( 1. / ( 1. - beta * 56 gamma = std::sqrt( 1. / ( 1. - beta * beta ) );
57 pz3 = gamma * ( pz + beta * E3 ); 57 pz3 = gamma * ( pz + beta * E3 );
58 pz4 = gamma * ( -pz + beta * E4 ); } 58 pz4 = gamma * ( -pz + beta * E4 ); }
59 else { /* COM 59 else { /* COM frame. */
60 pz3 = pz; 60 pz3 = pz;
61 pz4 = -pz; 61 pz4 = -pz;
62 } 62 }
63 outgoingData[1].isVelocity = outgoingData[ 63 outgoingData[1].isVelocity = outgoingData[0].isVelocity;
64 64
65 p_perp2 = px3 * px3 + py3 * py3; 65 p_perp2 = px3 * px3 + py3 * py3;
66 66
67 outgoingData[0].px_vx = px3; 67 outgoingData[0].px_vx = px3;
68 outgoingData[0].py_vy = py3; 68 outgoingData[0].py_vy = py3;
69 outgoingData[0].pz_vz = pz3; 69 outgoingData[0].pz_vz = pz3;
70 pp3 = p_perp2 + pz3 * pz3; 70 pp3 = p_perp2 + pz3 * pz3;
71 x = ( m3cc > 0 ) ? pp3 / ( 2 * m3cc2 ) : 1 71 x = ( m3cc > 0 ) ? pp3 / ( 2 * m3cc2 ) : 1.;
72 if( x < 1e-5 ) { 72 if( x < 1e-5 ) {
73 outgoingData[0].kineticEnergy = m3cc * 73 outgoingData[0].kineticEnergy = m3cc * x * ( 1 - 0.5 * x * ( 1 - x ) ); }
74 else { 74 else {
75 outgoingData[0].kineticEnergy = std::s 75 outgoingData[0].kineticEnergy = std::sqrt( m3cc2 + pp3 ) - m3cc;
76 } 76 }
77 outgoingData[1].px_vx = -px3; 77 outgoingData[1].px_vx = -px3;
78 outgoingData[1].py_vy = -py3; 78 outgoingData[1].py_vy = -py3;
79 outgoingData[1].pz_vz = pz4; 79 outgoingData[1].pz_vz = pz4;
80 pp4 = p_perp2 + pz4 * pz4; 80 pp4 = p_perp2 + pz4 * pz4;
81 x = ( m4cc > 0 ) ? pp4 / ( 2 * m4cc2 ) : 1 81 x = ( m4cc > 0 ) ? pp4 / ( 2 * m4cc2 ) : 1.;
82 if( x < 1e-5 ) { 82 if( x < 1e-5 ) {
83 outgoingData[1].kineticEnergy = m4cc * 83 outgoingData[1].kineticEnergy = m4cc * x * ( 1 - 0.5 * x * ( 1 - x ) ); }
84 else { 84 else {
85 outgoingData[1].kineticEnergy = std::s 85 outgoingData[1].kineticEnergy = std::sqrt( m4cc2 + pp4 ) - m4cc;
86 } 86 }
87 87
88 if( outgoingData[0].isVelocity ) { 88 if( outgoingData[0].isVelocity ) {
89 v_p = MCGIDI_speedOfLight_cm_sec / std 89 v_p = MCGIDI_speedOfLight_cm_sec / std::sqrt( pp3 + m3cc2 );
90 outgoingData[0].px_vx *= v_p; 90 outgoingData[0].px_vx *= v_p;
91 outgoingData[0].py_vy *= v_p; 91 outgoingData[0].py_vy *= v_p;
92 outgoingData[0].pz_vz *= v_p; 92 outgoingData[0].pz_vz *= v_p;
93 93
94 v_p = MCGIDI_speedOfLight_cm_sec / std 94 v_p = MCGIDI_speedOfLight_cm_sec / std::sqrt( pp4 + m4cc2 );
95 outgoingData[1].px_vx *= v_p; 95 outgoingData[1].px_vx *= v_p;
96 outgoingData[1].py_vy *= v_p; 96 outgoingData[1].py_vy *= v_p;
97 outgoingData[1].pz_vz *= v_p; 97 outgoingData[1].pz_vz *= v_p;
98 } 98 }
99 99
100 return( 0 ); 100 return( 0 );
101 } 101 }
102 /* 102 /*
103 ********************************************** 103 ************************************************************
104 */ 104 */
105 int MCGIDI_kinetics_COM2Lab( statusMessageRepo 105 int MCGIDI_kinetics_COM2Lab( statusMessageReporting *smr, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo, double masses[3] ) {
106 /* 106 /*
107 * massProjectile = masses[0], massTarget 107 * massProjectile = masses[0], massTarget = masses[1], massProduct = masses[2];
108 */ 108 */
109 double a = masses[0] + masses[1], b, e_in 109 double a = masses[0] + masses[1], b, e_in = modes.getProjectileEnergy( ) * masses[0] * masses[2] / ( a * a ), Ep;
110 110
111 if( decaySamplingInfo->frame != xDataTOM_f 111 if( decaySamplingInfo->frame != xDataTOM_frame_centerOfMass ) {
112 smr_setReportError2( smr, smr_unknownI 112 smr_setReportError2( smr, smr_unknownID, 1, "bad frame = %d for COM to lab conversion of mu/energy", decaySamplingInfo->frame );
113 return( 1 ); 113 return( 1 );
114 } 114 }
115 a = std::sqrt( e_in ); 115 a = std::sqrt( e_in );
116 b = std::sqrt( decaySamplingInfo->Ep ); 116 b = std::sqrt( decaySamplingInfo->Ep );
117 Ep = decaySamplingInfo->Ep + e_in + 2. * d 117 Ep = decaySamplingInfo->Ep + e_in + 2. * decaySamplingInfo->mu * a * b;
118 if( Ep != 0 ) { 118 if( Ep != 0 ) {
119 decaySamplingInfo->mu = ( a + decaySam 119 decaySamplingInfo->mu = ( a + decaySamplingInfo->mu * b ) / std::sqrt( Ep );
120 } 120 }
121 decaySamplingInfo->Ep = Ep; 121 decaySamplingInfo->Ep = Ep;
122 decaySamplingInfo->frame = xDataTOM_frame_ 122 decaySamplingInfo->frame = xDataTOM_frame_lab;
123 return( 0 ); 123 return( 0 );
124 } 124 }
125 125
126 #if defined __cplusplus 126 #if defined __cplusplus
127 } 127 }
128 #endif 128 #endif
129 129