Geant4  10.01.p01
tpia_kinetics.cc
Go to the documentation of this file.
1 /*
2 # <<BEGIN-copyright>>
3 # Copyright (c) 2010, Lawrence Livermore National Security, LLC.
4 # Produced at the Lawrence Livermore National Laboratory
5 # Written by Bret R. Beck, beck6@llnl.gov.
6 # CODE-461393
7 # All rights reserved.
8 #
9 # This file is part of GIDI. For details, see nuclear.llnl.gov.
10 # Please also read the "Additional BSD Notice" at nuclear.llnl.gov.
11 #
12 # Redistribution and use in source and binary forms, with or without modification,
13 # are permitted provided that the following conditions are met:
14 #
15 # 1) Redistributions of source code must retain the above copyright notice,
16 # this list of conditions and the disclaimer below.
17 # 2) Redistributions in binary form must reproduce the above copyright notice,
18 # this list of conditions and the disclaimer (as noted below) in the
19 # documentation and/or other materials provided with the distribution.
20 # 3) Neither the name of the LLNS/LLNL nor the names of its contributors may be
21 # used to endorse or promote products derived from this software without
22 # specific prior written permission.
23 #
24 # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
25 # EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
26 # OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
27 # SHALL LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR
28 # CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
29 # CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
30 # OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
31 # AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
32 # (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
33 # EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
34 # <<END-copyright>>
35 */
36 #include <iostream>
37 #include <float.h>
38 #include <string.h>
39 #include <cmath>
40 #include <tpia_target.h>
41 
42 #if defined __cplusplus
43 namespace GIDI {
44 using namespace GIDI;
45 #endif
46 
47 /*
48 ************************************************************
49 */
50 int tpia_kinetics_2BodyReaction( statusMessageReporting *smr, tpia_decayChannel *decayChannel, double K, double mu, double phi,
51  tpia_productOutgoingData *outgoingData ) {
52 
53  tpia_product *pp3 = tpia_decayChannel_getFirstProduct( decayChannel ), *pp4;
54  double m1 = decayChannel->m1_fullMass_MeV, m2 = decayChannel->m2_fullMass_MeV, m3, m4, mi, mf, Kp, x, beta;
55 
57  m3 = pp3->productID->fullMass_MeV;
58  m4 = pp4->productID->fullMass_MeV;
59  mi = m1 + m2;
60  mf = m3 + m4;
61  beta = std::sqrt( K * ( K + 2. * m1 ) ) / ( K + mi );
62  x = K * m2 / ( mi * mi );
63  if( x < 2e-5 ) { /* Kp is the total kinetic energy for m3 and m4 in the COM frame. */
64  Kp = mi - mf + K * m2 / mi * ( 1 - 0.5 * x * ( 1 - x ) ); }
65  else {
66  Kp = std::sqrt( mi * mi + 2 * K * m2 ) - mf;
67  }
68  if( Kp < 0 ) Kp = 0.; /* ???? There needs to be a better test here. */
69  outgoingData[0].decayChannel = &(pp3->decayChannel);
70  outgoingData[1].genre = outgoingData[0].genre;
71  outgoingData[1].productID = pp4->productID;
72  outgoingData[1].decayChannel = &(pp4->decayChannel);
73  return( tpia_kinetics_COMKineticEnergy2LabEnergyAndMomentum( smr, beta, Kp, mu, phi, m3, m4, outgoingData ) );
74 }
75 /*
76 ************************************************************
77 */
78 //int tpia_kinetics_COMKineticEnergy2LabEnergyAndMomentum( statusMessageReporting *smr, double beta, double e_kinetic_com, double mu, double phi,
79 int tpia_kinetics_COMKineticEnergy2LabEnergyAndMomentum( statusMessageReporting *, double beta, double e_kinetic_com, double mu, double phi,
80  double m3cc, double m4cc, tpia_productOutgoingData *outgoingData ) {
81 /*
82 * beta the velocity/speedOflight of the com frame relative to the lab frame.
83 * e_kinetic_com Total kinetic energy (K1 + K2) in the COM frame.
84 * mu std::cos( theta ) in the COM frame.
85 */
86  double x, v_p, p, pp3, pp4, px3, py3, pz3, pz4, pz, p_perp2, E3, E4, gamma, m3cc2 = m3cc * m3cc, m4cc2 = m4cc * m4cc;
87 
88  p = std::sqrt( e_kinetic_com * ( e_kinetic_com + 2. * m3cc ) * ( e_kinetic_com + 2. * m4cc ) * ( e_kinetic_com + 2. * ( m3cc + m4cc ) ) ) /
89  ( 2. * ( e_kinetic_com + m3cc + m4cc ) );
90  py3 = p * std::sqrt( 1 - mu * mu );
91  px3 = py3 * std::cos( phi );
92  py3 *= std::sin( phi );
93  pz = p * mu;
94  if( tpia_frame_getColumn( NULL, &(outgoingData[0].frame), 0 ) == tpia_referenceFrame_lab ) {
95  E3 = std::sqrt( p * p + m3cc2 );
96  E4 = std::sqrt( p * p + m4cc2 );
97  gamma = std::sqrt( 1. / ( 1. - beta * beta ) );
98  pz3 = gamma * ( pz + beta * E3 );
99  pz4 = gamma * ( -pz + beta * E4 ); }
100  else {
101  pz3 = pz;
102  pz4 = -pz;
103  }
104  outgoingData[1].isVelocity = outgoingData[0].isVelocity;
105  outgoingData[1].frame = outgoingData[0].frame;
106 
107  p_perp2 = px3 * px3 + py3 * py3;
108 
109  outgoingData[0].px_vx = px3;
110  outgoingData[0].py_vy = py3;
111  outgoingData[0].pz_vz = pz3;
112  pp3 = p_perp2 + pz3 * pz3;
113 //TK140602 Modified for protecting divided by 0 BEGIN
114  if ( m3cc2 != 0 )
115  x = pp3 / ( 2 * m3cc2 );
116  else
117  x = FLT_MIN;
118 //TK140602 Modified for protecting divided by 0 END
119  if( x < 1e-5 ) {
120  outgoingData[0].kineticEnergy = m3cc * x * ( 1 - 0.5 * x * ( 1 - x ) ); }
121  else {
122  outgoingData[0].kineticEnergy = std::sqrt( m3cc2 + pp3 ) - m3cc;
123  }
124  outgoingData[1].px_vx = -px3;
125  outgoingData[1].py_vy = -py3;
126  outgoingData[1].pz_vz = pz4;
127  pp4 = p_perp2 + pz4 * pz4;
128  x = pp4 / ( 2 * m4cc2 );
129  if( x < 1e-5 ) {
130  outgoingData[1].kineticEnergy = m4cc * x * ( 1 - 0.5 * x * ( 1 - x ) ); }
131  else {
132  outgoingData[1].kineticEnergy = std::sqrt( m4cc2 + pp4 ) - m4cc;
133  }
134 
135  if( outgoingData[0].isVelocity ) {
136  v_p = tpia_speedOfLight_cm_sec / std::sqrt( pp3 + m3cc2 );
137  outgoingData[0].px_vx *= v_p;
138  outgoingData[0].py_vy *= v_p;
139  outgoingData[0].pz_vz *= v_p;
140 
141  v_p = tpia_speedOfLight_cm_sec / std::sqrt( pp4 + m4cc2 );
142  outgoingData[1].px_vx *= v_p;
143  outgoingData[1].py_vy *= v_p;
144  outgoingData[1].pz_vz *= v_p;
145  }
146 
147  return( 0 );
148 }
149 
150 #if defined __cplusplus
151 }
152 #endif
static const double m3
Definition: G4SIunits.hh:112
tpia_product * tpia_decayChannel_getNextProduct(tpia_product *product)
int tpia_kinetics_2BodyReaction(statusMessageReporting *smr, tpia_decayChannel *decayChannel, double K, double mu, double phi, tpia_productOutgoingData *outgoingData)
static const double m2
Definition: G4SIunits.hh:111
int tpia_kinetics_COMKineticEnergy2LabEnergyAndMomentum(statusMessageReporting *, double beta, double e_kinetic_com, double mu, double phi, double m3cc, double m4cc, tpia_productOutgoingData *outgoingData)
int tpia_frame_getColumn(statusMessageReporting *smr, tpia_data_frame *frame, int column)
Definition: tpia_frame.cc:214
#define FLT_MIN
Definition: templates.hh:91
tpia_product * tpia_decayChannel_getFirstProduct(tpia_decayChannel *decayChannel)