include/IRT2/SinglePDF.h

0001 #pragma once
0002
0003 #include <vector>
0004 #include <math.h>
0005
0006 #include <TObject.h>
0007
0008 namespace IRT2 {
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 Reset( void ) {
0091     for(auto member: m_Members)
0092       delete member;
0093
0094     m_Members.clear();
0095   };
0096
0097   void AddMember(UniformPDF *pdf) { m_Members.push_back(pdf); };
0098   unsigned GetMemberCount( void ) const { return m_Members.size(); };
0099   unsigned GetWithinRangeCount(double x, double dx = 0.0) const {
0100     unsigned ret = 0;
0101
0102     for(auto member: m_Members)
0103       if (member->WithinRange(x, dx))
0104     ret++;
0105
0106     return ret;
0107   };
0108   double GetValue(double x) const {
0109     double ret = 0.0;
0110
0111     for(auto member: m_Members)
0112       ret += member->GetValue(x);
0113
0114     return ret;
0115   };
0116   double GetAverage( void ) const {
0117     double avg = 0.0, wtsum = 0.0;
0118
0119     for(auto member: m_Members) {
0120       avg   += member->GetWeight() * member->GetAverage();
0121       wtsum += member->GetWeight();
0122     } //for member
0123
0124     return (wtsum ? avg / wtsum : 0.0);
0125   };
0126   double GetAverage(double min, double max) const {
0127     double avg = 0.0, wtsum = 0.0;
0128
0129     for(auto member: m_Members) {
0130       double value = member->GetAverage();
0131       if (value < min || value > max) continue;
0132
0133       avg   += member->GetWeight() * value;//member->GetAverage();
0134       wtsum += member->GetWeight();
0135     } //for member
0136
0137     return (wtsum ? avg / wtsum : 0.0);
0138   };
0139   inline double GetRangeIntegral(double x0, double x1) const {
0140     double ret = 0.0;
0141
0142     for(auto member: m_Members)
0143       ret += member->GetRangeIntegral(x0, x1);
0144
0145     return ret;
0146   };
0147   inline double GetGaussianIntegral(double mu, double sigma) const {
0148     double ret = 0.0;
0149
0150     for(auto member: m_Members)
0151       ret += member->GetGaussianIntegral(mu, sigma);
0152
0153     return ret;
0154   };
0155
0156   UniformPDF *GetMember(unsigned id) const {
0157     return (id < m_Members.size() ? m_Members[id] : 0);
0158   };
0159   const std::vector<UniformPDF*> &GetMembers( void ) const { return m_Members; };
0160
0161  private:
0162   std::vector<UniformPDF*> m_Members;
0163
0164 #ifndef DISABLE_ROOT_IO
0165   ClassDef(VectorPDF, 1);
0166 #endif
0167 };
0168
0169 } // namespace IRT2