File indexing completed on 2024-04-06 12:09:25
0001 #ifndef QUANTILE_H
0002 #define QUANTILE_H
0003
0004 #include "TH1.h"
0005 #include <vector>
0006 #include <iostream>
0007
0008 struct Quantile {
0009 typedef std::pair<double, double> pair;
0010 typedef std::vector<pair> array;
0011
0012 pair operator()(const double frac) const { return fromHead(frac); }
0013 pair operator[](const double frac) const { return fromTail(frac); }
0014
0015 Quantile(const TH1* h) : N(1 + h->GetNbinsX()), Total(h->Integral(0, N)) {
0016 for (int i = 0; i < N; i++) {
0017 const double H = h->GetBinContent(i) + (!head.empty() ? head.back().second : 0);
0018 const double T = h->GetBinContent(N - i) + (!tail.empty() ? tail.back().second : 0);
0019 if (H)
0020 head.push_back(pair(h->GetBinWidth(i) + h->GetBinLowEdge(i), H));
0021 if (T)
0022 tail.push_back(pair(h->GetBinLowEdge(N - i), T));
0023 }
0024 }
0025
0026 pair fromHead(const double frac) const { return calculateQ(frac, true); }
0027 pair fromTail(const double frac) const { return calculateQ(frac, false); }
0028
0029 private:
0030 pair calculateQ(const double frac, const bool fromHead) const {
0031 const double f = frac < 0.5 ? frac : 1 - frac;
0032 array::const_iterator begin(((frac < 0.5) == fromHead) ? head.begin() : tail.begin()),
0033 end(((frac < 0.5) == fromHead) ? head.end() : tail.end()), bin(begin);
0034
0035 while (bin->second < f * Total)
0036 bin++;
0037
0038 if (bin == begin)
0039 return pair(-1, 0);
0040
0041 array::const_iterator binNext(next(bin, end)), binPrev(prev(bin, begin)), binPPrev(prev(binPrev, begin));
0042
0043 const double DX(binNext->first - binPPrev->first), DY((binNext->second - binPPrev->second) / Total),
0044
0045 dX(bin->first - binPrev->first), dY((bin->second - binPrev->second) / Total),
0046
0047 avgX((bin->first + binPrev->first) / 2), avgY((bin->second + binPrev->second) / (2 * Total)),
0048
0049 x_q(avgX + dX / dY * (f - avgY)), xerr_q(std::max(fabs(DX / DY), fabs(dX / dY)) * sqrt(f * (1 - f) / Total));
0050
0051 return pair(x_q, xerr_q);
0052 }
0053
0054 template <class T>
0055 T prev(T bin, T begin) const {
0056 T binPrev = bin;
0057 while (binPrev > begin && (binPrev - 1)->second == (bin - 1)->second)
0058 binPrev--;
0059 return binPrev;
0060 }
0061
0062 template <class T>
0063 T next(T bin, T end) const {
0064 T binNext = bin;
0065 while (binNext < end - 1 && (++binNext)->second == bin->second)
0066 ;
0067 return binNext;
0068 }
0069
0070 const int N;
0071 const double Total;
0072 array head, tail;
0073 };
0074 #endif