File indexing completed on 2024-09-27 07:03:46
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
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
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
0064 if (m_Sigma) sigma = sqrt(sigma*sigma + m_Sigma*m_Sigma);
0065
0066
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 }
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