Page MenuHomeHEPForge
Differential D58 Diff 259 src/EvtGenBase/EvtBreitWignerPdf.cpp

Changeset View

Standalone View

View Options

src/EvtGenBase/EvtBreitWignerPdf.cpp

Show All 23 Lines
#include "EvtGenBase/EvtPatches.hh" #include "EvtGenBase/EvtPatches.hh"
#include <assert.h> #include <assert.h>
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
EvtBreitWignerPdf::EvtBreitWignerPdf( double min, double max, double m0, EvtBreitWignerPdf::EvtBreitWignerPdf( double min, double max, double m0,
double g0 ) : double g0 ) :
EvtIntegPdf1D( min, max ), _m0( m0 ), _g0( g0 ) EvtIntegPdf1D( min, max ), m_m0( m0 ), m_g0( g0 )
{ {
} }
EvtBreitWignerPdf::EvtBreitWignerPdf( const EvtBreitWignerPdf& other ) : EvtBreitWignerPdf::EvtBreitWignerPdf( const EvtBreitWignerPdf& other ) :
EvtIntegPdf1D( other ), _m0( other._m0 ), _g0( other._g0 ) EvtIntegPdf1D( other ), m_m0( other.m_m0 ), m_g0( other.m_g0 )
{ {
} }
double EvtBreitWignerPdf::pdf( const EvtPoint1D& x ) const double EvtBreitWignerPdf::pdf( const EvtPoint1D& x ) const
{ {
double m = x.value(); double m = x.value();
if ( ( 0 == ( m - _m0 ) ) && ( 0. == _g0 ) ) { if ( ( 0 == ( m - m_m0 ) ) && ( 0. == m_g0 ) ) {
printf( "Delta function Breit-Wigner\n" ); printf( "Delta function Breit-Wigner\n" );
assert( 0 ); assert( 0 );
} }
double ret = _g0 / EvtConst::twoPi / double ret = m_g0 / EvtConst::twoPi /
( ( m - _m0 ) * ( m - _m0 ) + _g0 * _g0 / 4 ); ( ( m - m_m0 ) * ( m - m_m0 ) + m_g0 * m_g0 / 4 );
return ret; return ret;
} }
double EvtBreitWignerPdf::pdfIntegral( double m ) const double EvtBreitWignerPdf::pdfIntegral( double m ) const
{ {
double itg = 0; double itg = 0;
if ( _g0 == 0 ) { if ( m_g0 == 0 ) {
if ( m > _m0 ) if ( m > m_m0 )
itg = 1.; itg = 1.;
else if ( m < _m0 ) else if ( m < m_m0 )
itg = 0.; itg = 0.;
else else
itg = 0.5; itg = 0.5;
} else } else
itg = atan( ( m - _m0 ) / ( _g0 / 2. ) ) / EvtConst::pi + 0.5; itg = atan( ( m - m_m0 ) / ( m_g0 / 2. ) ) / EvtConst::pi + 0.5;
return itg; return itg;
} }
double EvtBreitWignerPdf::pdfIntegralInverse( double x ) const double EvtBreitWignerPdf::pdfIntegralInverse( double x ) const
{ {
if ( x < 0 || x > 1 ) { if ( x < 0 || x > 1 ) {
printf( "Invalid integral value %f\n", x ); printf( "Invalid integral value %f\n", x );
assert( 0 ); assert( 0 );
} }
double m = _m0; double m = m_m0;
if ( _g0 != 0 ) if ( m_g0 != 0 )
m = _m0 + ( _g0 / 2. ) * tan( EvtConst::pi * ( x - 0.5 ) ); m = m_m0 + ( m_g0 / 2. ) * tan( EvtConst::pi * ( x - 0.5 ) );
return m; return m;
} }