Changeset View
Changeset View
Standalone View
Standalone View
src/EvtGenBase/EvtKine.cpp
Show All 22 Lines | |||||
#include "EvtGenBase/EvtConst.hh" | #include "EvtGenBase/EvtConst.hh" | ||||
#include "EvtGenBase/EvtPatches.hh" | #include "EvtGenBase/EvtPatches.hh" | ||||
#include "EvtGenBase/EvtReport.hh" | #include "EvtGenBase/EvtReport.hh" | ||||
#include "EvtGenBase/EvtTensor4C.hh" | #include "EvtGenBase/EvtTensor4C.hh" | ||||
#include "EvtGenBase/EvtVector4C.hh" | #include "EvtGenBase/EvtVector4C.hh" | ||||
#include "EvtGenBase/EvtVector4R.hh" | #include "EvtGenBase/EvtVector4R.hh" | ||||
#include "EvtGenBase/EvtdFunction.hh" | #include "EvtGenBase/EvtdFunction.hh" | ||||
#include <math.h> | #include <cmath> | ||||
double EvtDecayAngle( const EvtVector4R& p, const EvtVector4R& q, | double EvtDecayAngle( const EvtVector4R& p, const EvtVector4R& q, | ||||
const EvtVector4R& d ) | const EvtVector4R& d ) | ||||
{ | { | ||||
double pd = p * d; | double pd = p * d; | ||||
double pq = p * q; | double pq = p * q; | ||||
double qd = q * d; | double qd = q * d; | ||||
double mp2 = p.mass2(); | double mp2 = p.mass2(); | ||||
double mq2 = q.mass2(); | double mq2 = q.mass2(); | ||||
double md2 = d.mass2(); | double md2 = d.mass2(); | ||||
double cost = ( pd * mq2 - pq * qd ) / | double cost = ( pd * mq2 - pq * qd ) / | ||||
sqrt( ( pq * pq - mq2 * mp2 ) * ( qd * qd - mq2 * md2 ) ); | std::sqrt( ( pq * pq - mq2 * mp2 ) * ( qd * qd - mq2 * md2 ) ); | ||||
return cost; | return cost; | ||||
} | } | ||||
double EvtDecayAngleChi( const EvtVector4R& p4_p, const EvtVector4R& p4_d1, | double EvtDecayAngleChi( const EvtVector4R& p4_p, const EvtVector4R& p4_d1, | ||||
const EvtVector4R& p4_d2, const EvtVector4R& p4_h1, | const EvtVector4R& p4_d2, const EvtVector4R& p4_h1, | ||||
const EvtVector4R& p4_h2 ) | const EvtVector4R& p4_h2 ) | ||||
{ | { | ||||
Show All 22 Lines | double EvtDecayAngleChi( const EvtVector4R& p4_p, const EvtVector4R& p4_d1, | ||||
d1_perp = d1_perp / d1_perp.d3mag(); | d1_perp = d1_perp / d1_perp.d3mag(); | ||||
d1_prime = d1_prime / d1_prime.d3mag(); | d1_prime = d1_prime / d1_prime.d3mag(); | ||||
double x, y; | double x, y; | ||||
x = d1_perp.dot( h1_perp ); | x = d1_perp.dot( h1_perp ); | ||||
y = d1_prime.dot( h1_perp ); | y = d1_prime.dot( h1_perp ); | ||||
double chi = atan2( y, x ); | double chi = std::atan2( y, x ); | ||||
if ( chi < 0.0 ) | if ( chi < 0.0 ) | ||||
chi += EvtConst::twoPi; | chi += EvtConst::twoPi; | ||||
return chi; | return chi; | ||||
} | } | ||||
double EvtDecayPlaneNormalAngle( const EvtVector4R& p, const EvtVector4R& q, | double EvtDecayPlaneNormalAngle( const EvtVector4R& p, const EvtVector4R& q, | ||||
const EvtVector4R& d1, const EvtVector4R& d2 ) | const EvtVector4R& d1, const EvtVector4R& d2 ) | ||||
{ | { | ||||
EvtVector4C lc = dual( EvtGenFunctions::directProd( d1, d2 ) ).cont2( q ); | EvtVector4C lc = dual( EvtGenFunctions::directProd( d1, d2 ) ).cont2( q ); | ||||
EvtVector4R l( real( lc.get( 0 ) ), real( lc.get( 1 ) ), | EvtVector4R l( real( lc.get( 0 ) ), real( lc.get( 1 ) ), | ||||
real( lc.get( 2 ) ), real( lc.get( 3 ) ) ); | real( lc.get( 2 ) ), real( lc.get( 3 ) ) ); | ||||
double pq = p * q; | double pq = p * q; | ||||
return q.mass() * ( p * l ) / | return q.mass() * ( p * l ) / | ||||
sqrt( -( pq * pq - p.mass2() * q.mass2() ) * l.mass2() ); | std::sqrt( -( pq * pq - p.mass2() * q.mass2() ) * l.mass2() ); | ||||
} | } | ||||
// Calculate phi using the given 4 vectors (all in the same frame) | // Calculate phi using the given 4 vectors (all in the same frame) | ||||
double EvtDecayAnglePhi( const EvtVector4R& z, const EvtVector4R& p, | double EvtDecayAnglePhi( const EvtVector4R& z, const EvtVector4R& p, | ||||
const EvtVector4R& q, const EvtVector4R& d ) | const EvtVector4R& q, const EvtVector4R& d ) | ||||
{ | { | ||||
double eq = ( p * q ) / p.mass(); | double eq = ( p * q ) / p.mass(); | ||||
double ed = ( p * d ) / p.mass(); | double ed = ( p * d ) / p.mass(); | ||||
double mq = q.mass(); | double mq = q.mass(); | ||||
double q2 = p.mag2r3( q ); | double q2 = p.mag2r3( q ); | ||||
double qd = p.dotr3( q, d ); | double qd = p.dotr3( q, d ); | ||||
double zq = p.dotr3( z, q ); | double zq = p.dotr3( z, q ); | ||||
double zd = p.dotr3( z, d ); | double zd = p.dotr3( z, d ); | ||||
double alpha = ( eq - mq ) / ( q2 * mq ) * qd - ed / mq; | double alpha = ( eq - mq ) / ( q2 * mq ) * qd - ed / mq; | ||||
double y = p.scalartripler3( z, q, d ) + alpha * p.scalartripler3( z, q, q ); | double y = p.scalartripler3( z, q, d ) + alpha * p.scalartripler3( z, q, q ); | ||||
double x = ( zq * ( qd + alpha * q2 ) - q2 * ( zd + alpha * zq ) ) / | double x = ( zq * ( qd + alpha * q2 ) - q2 * ( zd + alpha * zq ) ) / | ||||
sqrt( q2 ); | std::sqrt( q2 ); | ||||
double phi = atan2( y, x ); | double phi = std::atan2( y, x ); | ||||
return phi < 0 ? ( phi + EvtConst::twoPi ) : phi; | return phi < 0 ? ( phi + EvtConst::twoPi ) : phi; | ||||
} | } | ||||
EvtComplex wignerD( int j, int m1, int m2, double phi, double theta, double gamma ) | EvtComplex wignerD( int j, int m1, int m2, double phi, double theta, double gamma ) | ||||
{ | { | ||||
EvtComplex gp( 0.0, -phi * m1 ); | EvtComplex gp( 0.0, -phi * m1 ); | ||||
EvtComplex gm( 0.0, -gamma * m2 ); | EvtComplex gm( 0.0, -gamma * m2 ); | ||||
return exp( gp ) * EvtdFunction::d( j, m1, m2, theta ) * exp( gm ); | return exp( gp ) * EvtdFunction::d( j, m1, m2, theta ) * exp( gm ); | ||||
} | } | ||||
double twoBodyMomentum( const double M, const double m1, const double m2 ) | |||||
jback: This function can be constant. | |||||
Done Inline ActionsThis is not class member, just standalone function, so const does not apply here. kreps: This is not class member, just standalone function, so const does not apply here. | |||||
{ | |||||
const double MSq = M * M; | |||||
const double mSum = m1 + m2; | |||||
const double mDiff = m1 - m2; | |||||
double result = ( MSq - mDiff * mDiff ) * ( MSq - mSum * mSum ); | |||||
if ( result < 0 ) { | |||||
return 0; | |||||
} | |||||
return std::sqrt( result ) / ( 2 * M ); | |||||
} |
This function can be constant.