CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Quantile.h
Go to the documentation of this file.
1 #ifndef QUANTILE_H
2 #define QUANTILE_H
3 
4 #include "TH1.h"
5 #include <vector>
6 #include <iostream>
7 
8 struct Quantile {
9  typedef std::pair<double,double> pair;
10  typedef std::vector<pair> array;
11 
12  pair operator()(const double frac) const {return fromHead(frac);}
13  pair operator[](const double frac) const {return fromTail(frac);}
14 
15  Quantile(const TH1* h) :
16  N( 1 + h->GetNbinsX()),
17  Total(h->Integral(0,N))
18  { for(int i=0;i<N; i++) {
19  const double H = h->GetBinContent(i) + (head.size()?head.back().second:0);
20  const double T = h->GetBinContent(N-i) + (tail.size()?tail.back().second:0);
21  if(H) head.push_back( pair( h->GetBinWidth(i) + h->GetBinLowEdge(i) , H));
22  if(T) tail.push_back( pair( h->GetBinLowEdge(N-i),T));
23  }
24  }
25 
26  pair fromHead(const double frac) const {return calculateQ(frac,true);}
27  pair fromTail(const double frac) const {return calculateQ(frac,false);}
28 
29 private:
30 
31  pair calculateQ(const double frac, const bool fromHead) const {
32  const double f = frac<0.5 ? frac : 1-frac ;
33  array::const_iterator
34  begin( ( (frac<0.5) == fromHead ) ? head.begin() : tail.begin()),
35  end( ( (frac<0.5) == fromHead ) ? head.end() : tail.end()),
36  bin(begin);
37 
38  while( bin->second < f*Total ) bin++;
39 //dk if( bin==begin ) return pair(sqrt(-1),0);
40  if( bin==begin ) return pair(-1,0);
41 
42  array::const_iterator
43  binNext( next(bin,end)),
44  binPrev( prev(bin,begin)),
45  binPPrev( prev(binPrev,begin));
46 
47  const double
48  DX( binNext->first - binPPrev->first ),
49  DY( (binNext->second - binPPrev->second)/Total ),
50 
51  dX( bin->first - binPrev->first ),
52  dY( (bin->second - binPrev->second)/Total ),
53 
54  avgX( ( bin->first + binPrev->first) /2 ),
55  avgY( ( bin->second + binPrev->second) /(2*Total) ),
56 
57  x_q( avgX + dX/dY * ( f - avgY ) ),
58  xerr_q( std::max(fabs(DX/DY),fabs(dX/dY)) * sqrt(f*(1-f)/Total) );
59 
60  return pair(x_q,xerr_q);
61  }
62 
63 
64  template<class T> T prev(T bin, T begin) const {
65  T binPrev = bin;
66  while( binPrev > begin &&
67  (binPrev-1)->second == (bin-1)->second )
68  binPrev--;
69  return binPrev;
70  }
71 
72  template<class T> T next(T bin, T end) const {
73  T binNext = bin;
74  while( binNext<end-1 &&
75  (++binNext)->second == bin->second)
76  ;
77  return binNext;
78  }
79 
80  const int N;
81  const double Total;
83 };
84 #endif
int i
Definition: DBlmapReader.cc:9
T next(T bin, T end) const
Definition: Quantile.h:72
pair operator[](const double frac) const
Definition: Quantile.h:13
pair operator()(const double frac) const
Definition: Quantile.h:12
std::pair< double, double > pair
Definition: Quantile.h:9
U second(std::pair< T, U > const &p)
const T & max(const T &a, const T &b)
Quantile(const TH1 *h)
Definition: Quantile.h:15
T sqrt(T t)
Definition: SSEVec.h:48
array head
Definition: Quantile.h:82
double f[11][100]
#define end
Definition: vmac.h:37
const double Total
Definition: Quantile.h:81
pair fromHead(const double frac) const
Definition: Quantile.h:26
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
#define begin
Definition: vmac.h:30
array tail
Definition: Quantile.h:82
const int N
Definition: Quantile.h:80
pair calculateQ(const double frac, const bool fromHead) const
Definition: Quantile.h:31
T prev(T bin, T begin) const
Definition: Quantile.h:64
long double T
pair fromTail(const double frac) const
Definition: Quantile.h:27
std::vector< pair > array
Definition: Quantile.h:10