Changeset View
Changeset View
Standalone View
Standalone View
src/EvtGenBase/EvtDalitzResPdf.cpp
| Show All 23 Lines | |||||
| #include "EvtGenBase/EvtDalitzCoord.hh" | #include "EvtGenBase/EvtDalitzCoord.hh" | ||||
| #include "EvtGenBase/EvtPatches.hh" | #include "EvtGenBase/EvtPatches.hh" | ||||
| #include "EvtGenBase/EvtRandom.hh" | #include "EvtGenBase/EvtRandom.hh" | ||||
| #include <math.h> | #include <math.h> | ||||
| #include <stdio.h> | #include <stdio.h> | ||||
| using namespace EvtCyclic3; | using namespace EvtCyclic3; | ||||
| EvtDalitzResPdf::EvtDalitzResPdf( const EvtDalitzPlot& dp, double _m0, | EvtDalitzResPdf::EvtDalitzResPdf( const EvtDalitzPlot& dp, double m0, double g0, | ||||
| double _g0, EvtCyclic3::Pair pair ) : | EvtCyclic3::Pair pair ) : | ||||
| EvtPdf<EvtDalitzPoint>(), _dp( dp ), _m0( _m0 ), _g0( _g0 ), _pair( pair ) | EvtPdf<EvtDalitzPoint>(), m_dp( dp ), m_m0( m0 ), m_g0( g0 ), m_pair( pair ) | ||||
| { | { | ||||
| } | } | ||||
| EvtValError EvtDalitzResPdf::compute_integral( int N ) const | EvtValError EvtDalitzResPdf::compute_integral( int N ) const | ||||
| { | { | ||||
| assert( N != 0 ); | assert( N != 0 ); | ||||
| EvtCyclic3::Pair i = _pair; | EvtCyclic3::Pair i = m_pair; | ||||
| EvtCyclic3::Pair j = EvtCyclic3::next( i ); | EvtCyclic3::Pair j = EvtCyclic3::next( i ); | ||||
| // Trapezoidal integral | // Trapezoidal integral | ||||
| double dh = ( _dp.qAbsMax( j ) - _dp.qAbsMin( j ) ) / ( (double)N ); | double dh = ( m_dp.qAbsMax( j ) - m_dp.qAbsMin( j ) ) / ( (double)N ); | ||||
| double sum = 0; | double sum = 0; | ||||
| int ii; | int ii; | ||||
| for ( ii = 1; ii < N; ii++ ) { | for ( ii = 1; ii < N; ii++ ) { | ||||
| double x = _dp.qAbsMin( j ) + ii * dh; | double x = m_dp.qAbsMin( j ) + ii * dh; | ||||
| double min = ( _dp.qMin( i, j, x ) - _m0 * _m0 ) / _m0 / _g0; | double min = ( m_dp.qMin( i, j, x ) - m_m0 * m_m0 ) / m_m0 / m_g0; | ||||
| double max = ( _dp.qMax( i, j, x ) - _m0 * _m0 ) / _m0 / _g0; | double max = ( m_dp.qMax( i, j, x ) - m_m0 * m_m0 ) / m_m0 / m_g0; | ||||
| double itg = 1 / EvtConst::pi * ( atan( max ) - atan( min ) ); | double itg = 1 / EvtConst::pi * ( atan( max ) - atan( min ) ); | ||||
| sum += itg; | sum += itg; | ||||
| } | } | ||||
| EvtValError ret( sum * dh, 0. ); | EvtValError ret( sum * dh, 0. ); | ||||
| return ret; | return ret; | ||||
| } | } | ||||
| EvtDalitzPoint EvtDalitzResPdf::randomPoint() | EvtDalitzPoint EvtDalitzResPdf::randomPoint() | ||||
| { | { | ||||
| // Random point generation must be done in a box encompassing the | // Random point generation must be done in a box encompassing the | ||||
| // Dalitz plot | // Dalitz plot | ||||
| EvtCyclic3::Pair i = _pair; | EvtCyclic3::Pair i = m_pair; | ||||
| EvtCyclic3::Pair j = EvtCyclic3::next( i ); | EvtCyclic3::Pair j = EvtCyclic3::next( i ); | ||||
| double min = 1 / EvtConst::pi * | double min = 1 / EvtConst::pi * | ||||
| atan( ( _dp.qAbsMin( i ) - _m0 * _m0 ) / _m0 / _g0 ); | atan( ( m_dp.qAbsMin( i ) - m_m0 * m_m0 ) / m_m0 / m_g0 ); | ||||
| double max = 1 / EvtConst::pi * | double max = 1 / EvtConst::pi * | ||||
| atan( ( _dp.qAbsMax( i ) - _m0 * _m0 ) / _m0 / _g0 ); | atan( ( m_dp.qAbsMax( i ) - m_m0 * m_m0 ) / m_m0 / m_g0 ); | ||||
| int n = 0; | int n = 0; | ||||
| while ( n++ < 1000 ) { | while ( n++ < 1000 ) { | ||||
| double qj = EvtRandom::Flat( _dp.qAbsMin( j ), _dp.qAbsMax( j ) ); | double qj = EvtRandom::Flat( m_dp.qAbsMin( j ), m_dp.qAbsMax( j ) ); | ||||
| double r = EvtRandom::Flat( min, max ); | double r = EvtRandom::Flat( min, max ); | ||||
| double qi = tan( EvtConst::pi * r ) * _g0 * _m0 + _m0 * _m0; | double qi = tan( EvtConst::pi * r ) * m_g0 * m_m0 + m_m0 * m_m0; | ||||
| EvtDalitzCoord x( i, qi, j, qj ); | EvtDalitzCoord x( i, qi, j, qj ); | ||||
| EvtDalitzPoint ret( _dp, x ); | EvtDalitzPoint ret( m_dp, x ); | ||||
| if ( ret.isValid() ) | if ( ret.isValid() ) | ||||
| return ret; | return ret; | ||||
| } | } | ||||
| // All generated points turned out to be outside of the Dalitz plot | // All generated points turned out to be outside of the Dalitz plot | ||||
| // (in the outer box) | // (in the outer box) | ||||
| printf( "No point generated for dalitz plot after 1000 tries\n" ); | printf( "No point generated for dalitz plot after 1000 tries\n" ); | ||||
| return EvtDalitzPoint( 0., 0., 0., 0., 0., 0. ); | return EvtDalitzPoint( 0., 0., 0., 0., 0., 0. ); | ||||
| } | } | ||||
| double EvtDalitzResPdf::pdf( const EvtDalitzPoint& x ) const | double EvtDalitzResPdf::pdf( const EvtDalitzPoint& x ) const | ||||
| { | { | ||||
| EvtCyclic3::Pair i = _pair; | EvtCyclic3::Pair i = m_pair; | ||||
| double dq = x.q( i ) - _m0 * _m0; | double dq = x.q( i ) - m_m0 * m_m0; | ||||
| return 1 / EvtConst::pi * _g0 * _m0 / ( dq * dq + _g0 * _g0 * _m0 * _m0 ); | return 1 / EvtConst::pi * m_g0 * m_m0 / | ||||
| ( dq * dq + m_g0 * m_g0 * m_m0 * m_m0 ); | |||||
| } | } | ||||
| double EvtDalitzResPdf::pdfMaxValue() const | double EvtDalitzResPdf::pdfMaxValue() const | ||||
| { | { | ||||
| return 1 / ( EvtConst::pi * _g0 * _m0 ); | return 1 / ( EvtConst::pi * m_g0 * m_m0 ); | ||||
| } | } | ||||