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 | ||||