tpia_kinetics.cc

Go to the documentation of this file.
00001 /*
00002 # <<BEGIN-copyright>>
00003 # Copyright (c) 2010, Lawrence Livermore National Security, LLC. 
00004 # Produced at the Lawrence Livermore National Laboratory 
00005 # Written by Bret R. Beck, beck6@llnl.gov. 
00006 # CODE-461393
00007 # All rights reserved. 
00008 #  
00009 # This file is part of GIDI. For details, see nuclear.llnl.gov. 
00010 # Please also read the "Additional BSD Notice" at nuclear.llnl.gov. 
00011 # 
00012 # Redistribution and use in source and binary forms, with or without modification, 
00013 # are permitted provided that the following conditions are met: 
00014 #
00015 #      1) Redistributions of source code must retain the above copyright notice, 
00016 #         this list of conditions and the disclaimer below.
00017 #      2) Redistributions in binary form must reproduce the above copyright notice, 
00018 #         this list of conditions and the disclaimer (as noted below) in the 
00019 #          documentation and/or other materials provided with the distribution.
00020 #      3) Neither the name of the LLNS/LLNL nor the names of its contributors may be 
00021 #         used to endorse or promote products derived from this software without 
00022 #         specific prior written permission. 
00023 #
00024 # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY 
00025 # EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 
00026 # OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT 
00027 # SHALL LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR 
00028 # CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
00029 # CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 
00030 # OR SERVICES;  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED 
00031 # AND ON  ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
00032 # (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, 
00033 # EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
00034 # <<END-copyright>>
00035 */
00036 #include <string.h>
00037 #include <cmath>
00038 #include <tpia_target.h>
00039 
00040 #if defined __cplusplus
00041 namespace GIDI {
00042 using namespace GIDI;
00043 #endif 
00044 
00045 /*
00046 ************************************************************
00047 */
00048 int tpia_kinetics_2BodyReaction( statusMessageReporting *smr, tpia_decayChannel *decayChannel, double K, double mu, double phi, 
00049         tpia_productOutgoingData *outgoingData ) {
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 }
00073 /*
00074 ************************************************************
00075 */
00076 //int tpia_kinetics_COMKineticEnergy2LabEnergyAndMomentum( statusMessageReporting *smr, double beta, double e_kinetic_com, double mu, double phi, 
00077 int tpia_kinetics_COMKineticEnergy2LabEnergyAndMomentum( statusMessageReporting *, double beta, double e_kinetic_com, double mu, double phi, 
00078         double m3cc, double m4cc, tpia_productOutgoingData *outgoingData ) {
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 }
00142 
00143 #if defined __cplusplus
00144 }
00145 #endif

Generated on Mon May 27 17:50:35 2013 for Geant4 by  doxygen 1.4.7