Differential D47 Diff 191 src/EvtGenBase/EvtKine.cpp

Changeset 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 )
		jbackUnsubmitted Done Inline Actions This function can be constant. jback: This function can be constant.
		krepsAuthorUnsubmitted Done Inline Actions This 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 );
		}