|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
|
#ifndef QUANTILE_H
#define QUANTILE_H
#include "TH1.h"
#include <vector>
#include <iostream>
struct Quantile {
typedef std::pair<double, double> pair;
typedef std::vector<pair> array;
pair operator()(const double frac) const { return fromHead(frac); }
pair operator[](const double frac) const { return fromTail(frac); }
Quantile(const TH1* h) : N(1 + h->GetNbinsX()), Total(h->Integral(0, N)) {
for (int i = 0; i < N; i++) {
const double H = h->GetBinContent(i) + (!head.empty() ? head.back().second : 0);
const double T = h->GetBinContent(N - i) + (!tail.empty() ? tail.back().second : 0);
if (H)
head.push_back(pair(h->GetBinWidth(i) + h->GetBinLowEdge(i), H));
if (T)
tail.push_back(pair(h->GetBinLowEdge(N - i), T));
}
}
pair fromHead(const double frac) const { return calculateQ(frac, true); }
pair fromTail(const double frac) const { return calculateQ(frac, false); }
private:
pair calculateQ(const double frac, const bool fromHead) const {
const double f = frac < 0.5 ? frac : 1 - frac;
array::const_iterator begin(((frac < 0.5) == fromHead) ? head.begin() : tail.begin()),
end(((frac < 0.5) == fromHead) ? head.end() : tail.end()), bin(begin);
while (bin->second < f * Total)
bin++;
//dk if( bin==begin ) return pair(sqrt(-1),0);
if (bin == begin)
return pair(-1, 0);
array::const_iterator binNext(next(bin, end)), binPrev(prev(bin, begin)), binPPrev(prev(binPrev, begin));
const double DX(binNext->first - binPPrev->first), DY((binNext->second - binPPrev->second) / Total),
dX(bin->first - binPrev->first), dY((bin->second - binPrev->second) / Total),
avgX((bin->first + binPrev->first) / 2), avgY((bin->second + binPrev->second) / (2 * Total)),
x_q(avgX + dX / dY * (f - avgY)), xerr_q(std::max(fabs(DX / DY), fabs(dX / dY)) * sqrt(f * (1 - f) / Total));
return pair(x_q, xerr_q);
}
template <class T>
T prev(T bin, T begin) const {
T binPrev = bin;
while (binPrev > begin && (binPrev - 1)->second == (bin - 1)->second)
binPrev--;
return binPrev;
}
template <class T>
T next(T bin, T end) const {
T binNext = bin;
while (binNext < end - 1 && (++binNext)->second == bin->second)
;
return binNext;
}
const int N;
const double Total;
array head, tail;
};
#endif
|