Changeset View
Changeset View
Standalone View
Standalone View
EvtGenBase/EvtPdfSum.hh
Show All 35 Lines | public: | ||||
virtual ~EvtPdfSum(); | virtual ~EvtPdfSum(); | ||||
EvtPdfSum* clone() const override { return new EvtPdfSum( *this ); } | EvtPdfSum* clone() const override { return new EvtPdfSum( *this ); } | ||||
// Manipulate terms and coefficients | // Manipulate terms and coefficients | ||||
void addTerm( double c, const EvtPdf<T>& pdf ) | void addTerm( double c, const EvtPdf<T>& pdf ) | ||||
{ | { | ||||
assert( c >= 0. ); | assert( c >= 0. ); | ||||
_c.push_back( c ); | m_c.push_back( c ); | ||||
_term.push_back( pdf.clone() ); | m_term.push_back( pdf.clone() ); | ||||
} | } | ||||
void addOwnedTerm( double c, std::unique_ptr<EvtPdf<T>> pdf ) | void addOwnedTerm( double c, std::unique_ptr<EvtPdf<T>> pdf ) | ||||
{ | { | ||||
_c.push_back( c ); | m_c.push_back( c ); | ||||
_term.push_back( pdf.release() ); | m_term.push_back( pdf.release() ); | ||||
} | } | ||||
size_t nTerms() const { return _term.size(); } // number of terms | size_t nTerms() const { return m_term.size(); } // number of terms | ||||
inline double c( int i ) const { return _c[i]; } | inline double c( int i ) const { return m_c[i]; } | ||||
inline EvtPdf<T>* getPdf( int i ) const { return _term[i]; } | inline EvtPdf<T>* getPdf( int i ) const { return m_term[i]; } | ||||
// Integrals | // Integrals | ||||
EvtValError compute_integral() const override; | EvtValError compute_integral() const override; | ||||
EvtValError compute_integral( int N ) const override; | EvtValError compute_integral( int N ) const override; | ||||
T randomPoint() override; | T randomPoint() override; | ||||
protected: | protected: | ||||
double pdf( const T& p ) const override; | double pdf( const T& p ) const override; | ||||
vector<double> _c; // coefficients | vector<double> m_c; // coefficients | ||||
vector<EvtPdf<T>*> _term; // pointers to pdfs | vector<EvtPdf<T>*> m_term; // pointers to pdfs | ||||
}; | }; | ||||
template <class T> | template <class T> | ||||
EvtPdfSum<T>::EvtPdfSum( const EvtPdfSum<T>& other ) : EvtPdf<T>( other ) | EvtPdfSum<T>::EvtPdfSum( const EvtPdfSum<T>& other ) : EvtPdf<T>( other ) | ||||
{ | { | ||||
for ( size_t i = 0; i < other.nTerms(); i++ ) { | for ( size_t i = 0; i < other.nTerms(); i++ ) { | ||||
_c.push_back( other._c[i] ); | m_c.push_back( other.m_c[i] ); | ||||
_term.push_back( other._term[i]->clone() ); | m_term.push_back( other.m_term[i]->clone() ); | ||||
} | } | ||||
} | } | ||||
template <class T> | template <class T> | ||||
EvtPdfSum<T>::~EvtPdfSum() | EvtPdfSum<T>::~EvtPdfSum() | ||||
{ | { | ||||
for ( size_t i = 0; i < _c.size(); i++ ) { | for ( size_t i = 0; i < m_c.size(); i++ ) { | ||||
delete _term[i]; | delete m_term[i]; | ||||
} | } | ||||
} | } | ||||
template <class T> | template <class T> | ||||
double EvtPdfSum<T>::pdf( const T& p ) const | double EvtPdfSum<T>::pdf( const T& p ) const | ||||
{ | { | ||||
double ret = 0.; | double ret = 0.; | ||||
for ( size_t i = 0; i < _c.size(); i++ ) { | for ( size_t i = 0; i < m_c.size(); i++ ) { | ||||
ret += _c[i] * _term[i]->evaluate( p ); | ret += m_c[i] * m_term[i]->evaluate( p ); | ||||
} | } | ||||
return ret; | return ret; | ||||
} | } | ||||
/* | /* | ||||
* Compute the sum integral by summing all term integrals. | * Compute the sum integral by summing all term integrals. | ||||
*/ | */ | ||||
template <class T> | template <class T> | ||||
EvtValError EvtPdfSum<T>::compute_integral() const | EvtValError EvtPdfSum<T>::compute_integral() const | ||||
{ | { | ||||
EvtValError itg( 0.0, 0.0 ); | EvtValError itg( 0.0, 0.0 ); | ||||
for ( size_t i = 0; i < nTerms(); i++ ) { | for ( size_t i = 0; i < nTerms(); i++ ) { | ||||
itg += _c[i] * _term[i]->getItg(); | itg += m_c[i] * m_term[i]->getItg(); | ||||
} | } | ||||
return itg; | return itg; | ||||
} | } | ||||
template <class T> | template <class T> | ||||
EvtValError EvtPdfSum<T>::compute_integral( int N ) const | EvtValError EvtPdfSum<T>::compute_integral( int N ) const | ||||
{ | { | ||||
EvtValError itg( 0.0, 0.0 ); | EvtValError itg( 0.0, 0.0 ); | ||||
for ( size_t i = 0; i < nTerms(); i++ ) | for ( size_t i = 0; i < nTerms(); i++ ) | ||||
itg += _c[i] * _term[i]->getItg( N ); | itg += m_c[i] * m_term[i]->getItg( N ); | ||||
return itg; | return itg; | ||||
} | } | ||||
/* | /* | ||||
* Sample points randomly according to the sum of PDFs. First throw a random number uniformly | * Sample points randomly according to the sum of PDFs. First throw a random number uniformly | ||||
* between zero and the value of the sum integral. Using this random number select one | * between zero and the value of the sum integral. Using this random number select one | ||||
* of the PDFs. The generate a random point according to that PDF. | * of the PDFs. The generate a random point according to that PDF. | ||||
*/ | */ | ||||
template <class T> | template <class T> | ||||
T EvtPdfSum<T>::randomPoint() | T EvtPdfSum<T>::randomPoint() | ||||
{ | { | ||||
if ( !this->_itg.valueKnown() ) | if ( !this->m_itg.valueKnown() ) | ||||
this->_itg = compute_integral(); | this->m_itg = compute_integral(); | ||||
double max = this->_itg.value(); | double max = this->m_itg.value(); | ||||
double rnd = EvtRandom::Flat( 0, max ); | double rnd = EvtRandom::Flat( 0, max ); | ||||
double sum = 0.; | double sum = 0.; | ||||
size_t i; | size_t i; | ||||
for ( i = 0; i < nTerms(); i++ ) { | for ( i = 0; i < nTerms(); i++ ) { | ||||
double itg = _term[i]->getItg().value(); | double itg = m_term[i]->getItg().value(); | ||||
sum += _c[i] * itg; | sum += m_c[i] * itg; | ||||
if ( sum > rnd ) | if ( sum > rnd ) | ||||
break; | break; | ||||
} | } | ||||
return _term[i]->randomPoint(); | return m_term[i]->randomPoint(); | ||||
} | } | ||||
#endif | #endif |