irt/include/SinglePDF.h

0001
0002 #include <vector>
0003 #include <math.h>
0004
0005 #include <TObject.h>
0006
0007 #ifndef _SINGLE_PDF_
0008 #define _SINGLE_PDF_
0009
0010 class SinglePDF: public TObject {
0011  public:
0012   SinglePDF() {};
0013   ~SinglePDF() {};
0014
0015   virtual double GetValue   (double x)                  const = 0;
0016   virtual double GetAverage(void)                       const = 0;
0017   virtual bool   WithinRange(double x, double dx = 0.0) const = 0;
0018
0019 #ifndef DISABLE_ROOT_IO
0020   ClassDef(SinglePDF, 1);
0021 #endif
0022 };
0023
0024 class UniformPDF: public SinglePDF {
0025  public:
0026  UniformPDF(double x0 = 0.0, double x1 = 0.0, double weight = 1.0, double sigma = 0.0):
0027   m_X0(x0), m_X1(x1), m_Weight(weight), m_Sigma(sigma) {
0028     if (m_X0 > m_X1) std::swap(m_X0, m_X1);
0029
0030     // FIXME: sanity check;
0031     m_Norm = 1.0/(m_X1 - m_X0);
0032   };
0033   ~UniformPDF() {};
0034
0035   inline bool WithinRange(double x, double dx = 0.0) const {
0036     if (!dx) return (x >= m_X0 && x <= m_X1);
0037
0038     double xl = x - dx, xr = x + dx;
0039     if (xl > m_X1 || xr < m_X0) return false;
0040
0041     return true;
0042   };
0043   inline double GetValue(double x) const {
0044     return WithinRange(x) ? m_Norm*m_Weight : 0.0;
0045   };
0046   double GetAverage(void) const {
0047     return (m_X0 + m_X1)/2;
0048   };
0049   double GetWeight(void) const { return m_Weight; };
0050
0051   inline double GetRangeIntegral(double x0, double x1) const {
0052     if (x0 > x1) std::swap(x0, x1);
0053
0054     if (x0 > m_X1) return 0.0;
0055     if (x1 < m_X0) return 0.0;
0056
0057     //printf("%f %f\n", m_Norm, m_Weight);
0058     double l = (x0 > m_X0 ? x0 : m_X0), r = (x1 < m_X1 ? x1 : m_X1);
0059
0060     return m_Norm*m_Weight*(r - l);
0061   };
0062   inline double GetGaussianIntegral(double mu, double sigma) const {
0063     // FIXME: this is indeed not clean, but kind of makes sense;
0064     if (m_Sigma) sigma = sqrt(sigma*sigma + m_Sigma*m_Sigma);
0065
0066     // FIXME: need to think about normalization;
0067     double q0 = (1 - erf((m_X0 - mu)/(sqrt(2.)*sigma)))/2;
0068     double q1 = (1 - erf((m_X1 - mu)/(sqrt(2.)*sigma)))/2;
0069
0070     return m_Norm*m_Weight*(q0 - q1);
0071   };
0072
0073  private:
0074   double m_X0, m_X1, m_Weight, m_Norm, m_Sigma;
0075
0076 #ifndef DISABLE_ROOT_IO
0077   ClassDef(UniformPDF, 1);
0078 #endif
0079 };
0080
0081 class VectorPDF: public TObject {
0082  public:
0083   VectorPDF() {};
0084   ~VectorPDF() {
0085     for(auto pdf : m_Members) {
0086       delete pdf;
0087     }
0088   };
0089
0090   void AddMember(UniformPDF *pdf) { m_Members.push_back(pdf); };
0091   unsigned GetWithinRangeCount(double x, double dx = 0.0) const {
0092     unsigned ret = 0;
0093
0094     for(auto member: m_Members)
0095       if (member->WithinRange(x, dx))
0096     ret++;
0097
0098     return ret;
0099   };
0100   double GetValue(double x) const {
0101     double ret = 0.0;
0102
0103     for(auto member: m_Members)
0104       ret += member->GetValue(x);
0105
0106     return ret;
0107   };
0108   double GetAverage( void ) const {
0109     double avg = 0.0, wtsum = 0.0;
0110
0111     for(auto member: m_Members) {
0112       avg   += member->GetWeight() * member->GetAverage();
0113       wtsum += member->GetWeight();
0114     } //for member
0115
0116     return (wtsum ? avg / wtsum : 0.0);
0117   };
0118   inline double GetRangeIntegral(double x0, double x1) const {
0119     double ret = 0.0;
0120
0121     for(auto member: m_Members)
0122       ret += member->GetRangeIntegral(x0, x1);
0123
0124     return ret;
0125   };
0126   inline double GetGaussianIntegral(double mu, double sigma) const {
0127     double ret = 0.0;
0128
0129     for(auto member: m_Members)
0130       ret += member->GetGaussianIntegral(mu, sigma);
0131
0132     return ret;
0133   };
0134
0135  private:
0136   std::vector<UniformPDF*> m_Members;
0137
0138 #ifndef DISABLE_ROOT_IO
0139   ClassDef(VectorPDF, 1);
0140 #endif
0141 };
0142
0143
0144 #endif