Changeset View
Changeset View
Standalone View
Standalone View
EvtGenBase/EvtPdf.hh
Show First 20 Lines • Show All 66 Lines • ▼ Show 20 Lines | |||||
class EvtPdfPred; | class EvtPdfPred; | ||||
template <class T> | template <class T> | ||||
class EvtPdfGen; | class EvtPdfGen; | ||||
template <class T> | template <class T> | ||||
class EvtPdf { | class EvtPdf { | ||||
public: | public: | ||||
EvtPdf() {} | EvtPdf() {} | ||||
EvtPdf( const EvtPdf& other ) : _itg( other._itg ) {} | EvtPdf( const EvtPdf& other ) : m_itg( other.m_itg ) {} | ||||
virtual ~EvtPdf() {} | virtual ~EvtPdf() {} | ||||
virtual EvtPdf<T>* clone() const = 0; | virtual EvtPdf<T>* clone() const = 0; | ||||
double evaluate( const T& p ) const | double evaluate( const T& p ) const | ||||
{ | { | ||||
if ( p.isValid() ) | if ( p.isValid() ) | ||||
return pdf( p ); | return pdf( p ); | ||||
else | else | ||||
return 0.; | return 0.; | ||||
} | } | ||||
// Find PDF maximum. Points are sampled according to pc | // Find PDF maximum. Points are sampled according to pc | ||||
EvtPdfMax<T> findMax( const EvtPdf<T>& pc, int N ); | EvtPdfMax<T> findMax( const EvtPdf<T>& pc, int N ); | ||||
// Find generation efficiency. | // Find generation efficiency. | ||||
EvtValError findGenEff( const EvtPdf<T>& pc, int N, int nFindMax ); | EvtValError findGenEff( const EvtPdf<T>& pc, int N, int nFindMax ); | ||||
// Analytic integration. Calls cascade down until an overridden | // Analytic integration. Calls cascade down until an overridden | ||||
// method is called. | // method is called. | ||||
void setItg( EvtValError itg ) { _itg = itg; } | void setItg( EvtValError itg ) { m_itg = itg; } | ||||
EvtValError getItg() const | EvtValError getItg() const | ||||
{ | { | ||||
if ( !_itg.valueKnown() ) | if ( !m_itg.valueKnown() ) | ||||
_itg = compute_integral(); | m_itg = compute_integral(); | ||||
return _itg; | return m_itg; | ||||
} | } | ||||
EvtValError getItg( int N ) const | EvtValError getItg( int N ) const | ||||
{ | { | ||||
if ( !_itg.valueKnown() ) | if ( !m_itg.valueKnown() ) | ||||
_itg = compute_integral( N ); | m_itg = compute_integral( N ); | ||||
return _itg; | return m_itg; | ||||
} | } | ||||
virtual EvtValError compute_integral() const | virtual EvtValError compute_integral() const | ||||
{ | { | ||||
printf( "Analytic integration of PDF is not defined\n" ); | printf( "Analytic integration of PDF is not defined\n" ); | ||||
assert( 0 ); | assert( 0 ); | ||||
return EvtValError{}; | return EvtValError{}; | ||||
} | } | ||||
Show All 12 Lines | public: | ||||
EvtPredGen<EvtPdfGen<T>, EvtPdfPred<T>> accRejGen( const EvtPdf<T>& pc, | EvtPredGen<EvtPdfGen<T>, EvtPdfPred<T>> accRejGen( const EvtPdf<T>& pc, | ||||
int nMax, | int nMax, | ||||
double factor = 1. ); | double factor = 1. ); | ||||
virtual T randomPoint(); | virtual T randomPoint(); | ||||
protected: | protected: | ||||
virtual double pdf( const T& ) const = 0; | virtual double pdf( const T& ) const = 0; | ||||
mutable EvtValError _itg; | mutable EvtValError m_itg; | ||||
}; | }; | ||||
template <class T> | template <class T> | ||||
class EvtPdfGen { | class EvtPdfGen { | ||||
public: | public: | ||||
typedef T result_type; | typedef T result_type; | ||||
EvtPdfGen() : _pdf( 0 ) {} | EvtPdfGen() : m_pdf( 0 ) {} | ||||
EvtPdfGen( const EvtPdfGen<T>& other ) : | EvtPdfGen( const EvtPdfGen<T>& other ) : | ||||
_pdf( other._pdf ? other._pdf->clone() : 0 ) | m_pdf( other.m_pdf ? other.m_pdf->clone() : 0 ) | ||||
{ | { | ||||
} | } | ||||
EvtPdfGen( const EvtPdf<T>& pdf ) : _pdf( pdf.clone() ) {} | EvtPdfGen( const EvtPdf<T>& pdf ) : m_pdf( pdf.clone() ) {} | ||||
~EvtPdfGen() { delete _pdf; } | ~EvtPdfGen() { delete m_pdf; } | ||||
result_type operator()() { return _pdf->randomPoint(); } | result_type operator()() { return m_pdf->randomPoint(); } | ||||
private: | private: | ||||
EvtPdf<T>* _pdf; | EvtPdf<T>* m_pdf; | ||||
}; | }; | ||||
template <class T> | template <class T> | ||||
class EvtPdfPred { | class EvtPdfPred { | ||||
public: | public: | ||||
typedef T argument_type; | typedef T argument_type; | ||||
typedef bool result_type; | typedef bool result_type; | ||||
EvtPdfPred() {} | EvtPdfPred() {} | ||||
EvtPdfPred( const EvtPdf<T>& thePdf ) : itsPdf( thePdf.clone() ) {} | EvtPdfPred( const EvtPdf<T>& thePdf ) : m_pdf( thePdf.clone() ) {} | ||||
EvtPdfPred( const EvtPdfPred& other ) : | EvtPdfPred( const EvtPdfPred& other ) : | ||||
COPY_PTR( itsPdf ), COPY_MEM( itsPdfMax ) | COPY_PTR( m_pdf ), COPY_MEM( m_pdfMax ) | ||||
{ | { | ||||
} | } | ||||
~EvtPdfPred() { delete itsPdf; } | ~EvtPdfPred() { delete m_pdf; } | ||||
result_type operator()( argument_type p ) | result_type operator()( argument_type p ) | ||||
{ | { | ||||
assert( itsPdf ); | assert( m_pdf ); | ||||
assert( itsPdfMax.valueKnown() ); | assert( m_pdfMax.valueKnown() ); | ||||
double random = EvtRandom::Flat( 0., itsPdfMax.value() ); | double random = EvtRandom::Flat( 0., m_pdfMax.value() ); | ||||
return ( random <= itsPdf->evaluate( p ) ); | return ( random <= m_pdf->evaluate( p ) ); | ||||
} | } | ||||
EvtPdfMax<T> getMax() const { return itsPdfMax; } | EvtPdfMax<T> getMax() const { return m_pdfMax; } | ||||
void setMax( const EvtPdfMax<T>& max ) { itsPdfMax = max; } | void setMax( const EvtPdfMax<T>& max ) { m_pdfMax = max; } | ||||
template <class InputIterator> | template <class InputIterator> | ||||
void compute_max( InputIterator it, InputIterator end, double factor = 1. ) | void compute_max( InputIterator it, InputIterator end, double factor = 1. ) | ||||
{ | { | ||||
T p = *it++; | T p = *it++; | ||||
itsPdfMax = EvtPdfMax<T>( p, itsPdf->evaluate( p ) * factor ); | m_pdfMax = EvtPdfMax<T>( p, m_pdf->evaluate( p ) * factor ); | ||||
while ( !( it == end ) ) { | while ( !( it == end ) ) { | ||||
T p = *it++; | T pp = *it++; | ||||
double val = itsPdf->evaluate( p ) * factor; | double val = m_pdf->evaluate( pp ) * factor; | ||||
if ( val > itsPdfMax.value() ) | if ( val > m_pdfMax.value() ) | ||||
itsPdfMax = EvtPdfMax<T>( p, val ); | m_pdfMax = EvtPdfMax<T>( pp, val ); | ||||
} | } | ||||
} | } | ||||
private: | private: | ||||
EvtPdf<T>* itsPdf; | EvtPdf<T>* m_pdf; | ||||
EvtPdfMax<T> itsPdfMax; | EvtPdfMax<T> m_pdfMax; | ||||
}; | }; | ||||
template <class T> | template <class T> | ||||
class EvtPdfUnary { | class EvtPdfUnary { | ||||
public: | public: | ||||
typedef double result_type; | typedef double result_type; | ||||
typedef T argument_type; | typedef T argument_type; | ||||
EvtPdfUnary() {} | EvtPdfUnary() {} | ||||
EvtPdfUnary( const EvtPdf<T>& thePdf ) : itsPdf( thePdf.clone() ) {} | EvtPdfUnary( const EvtPdf<T>& thePdf ) : m_pdf( thePdf.clone() ) {} | ||||
EvtPdfUnary( const EvtPdfUnary& other ) : COPY_PTR( itsPdf ) {} | EvtPdfUnary( const EvtPdfUnary& other ) : COPY_PTR( m_pdf ) {} | ||||
~EvtPdfUnary() { delete itsPdf; } | ~EvtPdfUnary() { delete m_pdf; } | ||||
result_type operator()( argument_type p ) | result_type operator()( argument_type p ) | ||||
{ | { | ||||
assert( itsPdf ); | assert( m_pdf ); | ||||
double ret = itsPdf->evaluate( p ); | double ret = m_pdf->evaluate( p ); | ||||
return ret; | return ret; | ||||
} | } | ||||
private: | private: | ||||
EvtPdf<T>* itsPdf; | EvtPdf<T>* m_pdf; | ||||
}; | }; | ||||
template <class T> | template <class T> | ||||
class EvtPdfDiv : public EvtPdf<T> { | class EvtPdfDiv : public EvtPdf<T> { | ||||
public: | public: | ||||
EvtPdfDiv() : itsNum( 0 ), itsDen( 0 ) {} | EvtPdfDiv() : m_num( 0 ), m_den( 0 ) {} | ||||
EvtPdfDiv( const EvtPdf<T>& theNum, const EvtPdf<T>& theDen ) : | EvtPdfDiv( const EvtPdf<T>& theNum, const EvtPdf<T>& theDen ) : | ||||
EvtPdf<T>(), itsNum( theNum.clone() ), itsDen( theDen.clone() ) | EvtPdf<T>(), m_num( theNum.clone() ), m_den( theDen.clone() ) | ||||
{ | { | ||||
} | } | ||||
EvtPdfDiv( const EvtPdfDiv<T>& other ) : | EvtPdfDiv( const EvtPdfDiv<T>& other ) : | ||||
EvtPdf<T>( other ), COPY_PTR( itsNum ), COPY_PTR( itsDen ) | EvtPdf<T>( other ), COPY_PTR( m_num ), COPY_PTR( m_den ) | ||||
{ | { | ||||
} | } | ||||
virtual ~EvtPdfDiv() | virtual ~EvtPdfDiv() | ||||
{ | { | ||||
delete itsNum; | delete m_num; | ||||
delete itsDen; | delete m_den; | ||||
} | } | ||||
EvtPdf<T>* clone() const override { return new EvtPdfDiv( *this ); } | EvtPdf<T>* clone() const override { return new EvtPdfDiv( *this ); } | ||||
double pdf( const T& p ) const override | double pdf( const T& p ) const override | ||||
{ | { | ||||
double num = itsNum->evaluate( p ); | double num = m_num->evaluate( p ); | ||||
double den = itsDen->evaluate( p ); | double den = m_den->evaluate( p ); | ||||
assert( den != 0 ); | assert( den != 0 ); | ||||
return num / den; | return num / den; | ||||
} | } | ||||
private: | private: | ||||
EvtPdf<T>* itsNum; // numerator | EvtPdf<T>* m_num; // numerator | ||||
EvtPdf<T>* itsDen; // denominator | EvtPdf<T>* m_den; // denominator | ||||
}; | }; | ||||
template <class T> | template <class T> | ||||
EvtPdfMax<T> EvtPdf<T>::findMax( const EvtPdf<T>& pc, int N ) | EvtPdfMax<T> EvtPdf<T>::findMax( const EvtPdf<T>& pc, int N ) | ||||
{ | { | ||||
EvtPdfPred<T> pred( *this ); | EvtPdfPred<T> pred( *this ); | ||||
EvtPdfGen<T> gen( pc ); | EvtPdfGen<T> gen( pc ); | ||||
pred.compute_max( iter( gen, N ), iter( gen ) ); | pred.compute_max( iter( gen, N ), iter( gen ) ); | ||||
▲ Show 20 Lines • Show All 42 Lines • ▼ Show 20 Lines | if ( N > 0 ) { | ||||
// Due to numerical precision dev2 may sometimes be negative | // Due to numerical precision dev2 may sometimes be negative | ||||
if ( dev2 < 0. ) | if ( dev2 < 0. ) | ||||
dev2 = 0.; | dev2 = 0.; | ||||
double error = sqrt( dev2 / ( (double)N ) ); | double error = sqrt( dev2 / ( (double)N ) ); | ||||
x = EvtValError( av, error ); | x = EvtValError( av, error ); | ||||
} else | } else | ||||
x = EvtValError( av ); | x = EvtValError( av ); | ||||
} | } | ||||
_itg = x * pc.getItg(); | m_itg = x * pc.getItg(); | ||||
return _itg; | return m_itg; | ||||
} | } | ||||
template <class T> | template <class T> | ||||
T EvtPdf<T>::randomPoint() | T EvtPdf<T>::randomPoint() | ||||
{ | { | ||||
printf( "Function defined for analytic PDFs only\n" ); | printf( "Function defined for analytic PDFs only\n" ); | ||||
assert( 0 ); | assert( 0 ); | ||||
T temp; | T temp; | ||||
Show All 16 Lines |