interface/SusyDQM/Quantile.h

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     //dk    if( bin==begin ) return pair(sqrt(-1),0);
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