include/EvtGenBase/EvtPdfSum.hh

0001
0002 /***********************************************************************
0003 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors       *
0004 *                                                                      *
0005 * This file is part of EvtGen.                                         *
0006 *                                                                      *
0007 * EvtGen is free software: you can redistribute it and/or modify       *
0008 * it under the terms of the GNU General Public License as published by *
0009 * the Free Software Foundation, either version 3 of the License, or    *
0010 * (at your option) any later version.                                  *
0011 *                                                                      *
0012 * EvtGen is distributed in the hope that it will be useful,            *
0013 * but WITHOUT ANY WARRANTY; without even the implied warranty of       *
0014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        *
0015 * GNU General Public License for more details.                         *
0016 *                                                                      *
0017 * You should have received a copy of the GNU General Public License    *
0018 * along with EvtGen.  If not, see <https://www.gnu.org/licenses/>.     *
0019 ***********************************************************************/
0020
0021 #ifndef EVT_PDF_SUM_HH
0022 #define EVT_PDF_SUM_HH
0023
0024 #include <stdio.h>
0025 #include <vector>
0026 using std::vector;
0027 #include "EvtGenBase/EvtPdf.hh"
0028
0029 // Sum of PDF functions.
0030
0031 template <class T>
0032 class EvtPdfSum : public EvtPdf<T> {
0033   public:
0034     EvtPdfSum() {}
0035     EvtPdfSum( const EvtPdfSum<T>& other );
0036     virtual ~EvtPdfSum();
0037     EvtPdfSum* clone() const override { return new EvtPdfSum( *this ); }
0038
0039     // Manipulate terms and coefficients
0040
0041     void addTerm( double c, const EvtPdf<T>& pdf )
0042     {
0043         assert( c >= 0. );
0044         _c.push_back( c );
0045         _term.push_back( pdf.clone() );
0046     }
0047
0048     void addOwnedTerm( double c, std::unique_ptr<EvtPdf<T>> pdf )
0049     {
0050         _c.push_back( c );
0051         _term.push_back( pdf.release() );
0052     }
0053
0054     size_t nTerms() const { return _term.size(); }    // number of terms
0055
0056     inline double c( int i ) const { return _c[i]; }
0057     inline EvtPdf<T>* getPdf( int i ) const { return _term[i]; }
0058
0059     // Integrals
0060
0061     EvtValError compute_integral() const override;
0062     EvtValError compute_integral( int N ) const override;
0063     T randomPoint() override;
0064
0065   protected:
0066     double pdf( const T& p ) const override;
0067
0068     vector<double> _c;           // coefficients
0069     vector<EvtPdf<T>*> _term;    // pointers to pdfs
0070 };
0071
0072 template <class T>
0073 EvtPdfSum<T>::EvtPdfSum( const EvtPdfSum<T>& other ) : EvtPdf<T>( other )
0074 {
0075     for ( size_t i = 0; i < other.nTerms(); i++ ) {
0076         _c.push_back( other._c[i] );
0077         _term.push_back( other._term[i]->clone() );
0078     }
0079 }
0080
0081 template <class T>
0082 EvtPdfSum<T>::~EvtPdfSum()
0083 {
0084     for ( size_t i = 0; i < _c.size(); i++ ) {
0085         delete _term[i];
0086     }
0087 }
0088
0089 template <class T>
0090 double EvtPdfSum<T>::pdf( const T& p ) const
0091 {
0092     double ret = 0.;
0093     for ( size_t i = 0; i < _c.size(); i++ ) {
0094         ret += _c[i] * _term[i]->evaluate( p );
0095     }
0096     return ret;
0097 }
0098
0099 /*
0100  * Compute the sum integral by summing all term integrals.
0101  */
0102
0103 template <class T>
0104 EvtValError EvtPdfSum<T>::compute_integral() const
0105 {
0106     EvtValError itg( 0.0, 0.0 );
0107     for ( size_t i = 0; i < nTerms(); i++ ) {
0108         itg += _c[i] * _term[i]->getItg();
0109     }
0110     return itg;
0111 }
0112
0113 template <class T>
0114 EvtValError EvtPdfSum<T>::compute_integral( int N ) const
0115 {
0116     EvtValError itg( 0.0, 0.0 );
0117     for ( size_t i = 0; i < nTerms(); i++ )
0118         itg += _c[i] * _term[i]->getItg( N );
0119     return itg;
0120 }
0121
0122 /*
0123  * Sample points randomly according to the sum of PDFs. First throw a random number uniformly
0124  * between zero and the value of the sum integral. Using this random number select one
0125  * of the PDFs. The generate a random point according to that PDF.
0126  */
0127
0128 template <class T>
0129 T EvtPdfSum<T>::randomPoint()
0130 {
0131     if ( !this->_itg.valueKnown() )
0132         this->_itg = compute_integral();
0133
0134     double max = this->_itg.value();
0135     double rnd = EvtRandom::Flat( 0, max );
0136
0137     double sum = 0.;
0138     size_t i;
0139     for ( i = 0; i < nTerms(); i++ ) {
0140         double itg = _term[i]->getItg().value();
0141         sum += _c[i] * itg;
0142         if ( sum > rnd )
0143             break;
0144     }
0145
0146     return _term[i]->randomPoint();
0147 }
0148
0149 #endif