LlrUpmu.cxx
Go to the documentation of this file.
1 #include "Eval/LlrUpmu.h"
2 #include <boost/math/special_functions.hpp>
3 
4 /******
5 The above header is here because: the `genreflex` tool (the one that's being run during the `codegen` step) is built on the older C++03 library. our version of boost uses newer C++11 stuff (including that `std::__atomic_thread_fence` object). But `genreflex` only needs to look at the `.h` for a class it's going to produce a dictionary for: it doesn't care about the implementation (.cxx). So by pushing the `boost` stuff into the `.cxx`, `genreflex` never sees it, and there's no problem.
6 *******/
7 
8 namespace upmuana
9 {
10 
12  :
13  outliers_(), P_up_(1e-30), P_dn_(1e-30),
14  chi2_(1e6), slope_(-1e6), llr_(0.)
15  {}
16 
18  {}
19 
20  inline double LlrUpmu::getErr(double PE) {
21  //return 231594.0/(2267.16+pow(PE,2.01011)) + 9.55689;
22  return 165143/(1882.9+pow(PE,2.11447)) + 10.4321;
23  }
24 
25  inline double LlrUpmu::getDist(double x1, double y1, double z1,
26  double x2, double y2, double z2) {
27  return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
28  }
29 
30  void LlrUpmu::setLLR(std::set< std::pair<rb::CellHit,rb::RecoHit>,
31  bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
32  std::pair<rb::CellHit,rb::RecoHit>)>
33  sortedHits)
34  {
35  std::vector<rb::RecoHit> outliers;
36  double P_up;
37  double P_dn;
38  double chi2;
39  double slope;
40 
41  if (sortedHits.size() < 2) {
42  P_up = 1e-30;
43  P_dn = 1e-30;
44  chi2 = 1e6;
45  slope = -1e6;
46  return;
47  }
48 
49  std::vector<double> measuredTimes;
50  std::vector<double> expectedTimes;
51  std::vector<double> wgts;
52 
53  // first hit
54  measuredTimes.push_back(0.0);
55  expectedTimes.push_back(0.0);
56  double err = getErr(sortedHits.begin()->first.PE());
57  wgts.push_back(1.0/err/err);
58 
59  double startX = sortedHits.begin()->second.X(),
60  startY = sortedHits.begin()->second.Y(),
61  startZ = sortedHits.begin()->second.Z(),
62  startT = sortedHits.begin()->second.T();
63 
64  std::vector<rb::RecoHit> in;
65 
66  for (auto i_hit = ++sortedHits.begin();
67  i_hit != sortedHits.end();
68  ++i_hit) {
69  in.push_back(i_hit->second);
70  measuredTimes.push_back(i_hit->second.T() - startT);
71  double dist = getDist(i_hit->second.X(), i_hit->second.Y(),
72  i_hit->second.Z(), startX, startY, startZ);
73  expectedTimes.push_back(dist/29.97);
74  err = getErr(i_hit->first.PE());
75  wgts.push_back(1.0/err/err);
76  }
77 
78  LLR(expectedTimes, measuredTimes, wgts, slope, chi2, P_up, P_dn,
79  in, outliers);
80 
82  P_up_ = P_up;
83  P_dn_ = P_dn;
84  chi2_ = chi2;
85  slope_ = slope;
86  llr_ = log(P_up/P_dn);
87 
88  }
89 
90  void LlrUpmu::LLR(std::vector<double>& eT,
91  std::vector<double>& mT,
92  std::vector<double>& wgts, double& slope,
93  double& chi2, double& P_up, double& P_dn,
94  std::vector<rb::RecoHit>& in,
95  std::vector<rb::RecoHit>& outliers)
96  {
97  // eT - x, mT - y
98  size_t n = mT.size();
99 
100  // y = a + bx
101  double a, b;
102  if (eT.size() >= 2)
103  chi2 = util::LinFit(eT, mT, wgts, b, a);
104  else {
105  chi2 = 1e6;
106  slope = -1e6;
107  P_up = 1e-30;
108  P_dn = 1e-30;
109  return;
110  }
111 
112  size_t totAllowOutlier = n/10;
113  size_t nOutliers = 0;
114  std::vector<double> x_filt, y_filt, w_filt;
115  for (size_t i=0; i < n; i++) {
116  double y_fit = a + b*eT.at(i);
117  if ( fabs(mT.at(i) - y_fit) < 5*(1.0/sqrt(wgts.at(i)))
118  || nOutliers>totAllowOutlier)
119  {
120  x_filt.push_back(eT.at(i));
121  y_filt.push_back(mT.at(i));
122  w_filt.push_back(wgts.at(i));
123  }
124  else {
125  ++nOutliers;
126  outliers.push_back(in[i]);
127  }
128  }
129 
130 
131  if (x_filt.size() >= 2)
132  chi2 = util::LinFit(x_filt, y_filt, w_filt, b, a);
133  else {
134  chi2 = 1e6;
135  slope = -1e6;
136  P_up = 1e-30;
137  P_dn = 1e-30;
138  }
139  n = x_filt.size();
140  if (n < 5) {
141  slope = 0;
142  chi2 = 999;
143  P_up = 1e-30;
144  P_dn = 1;
145  return;
146  }
147 
148  slope = b;
149  chi2 /= (double)(n-2);
150 
151  double one_over_ee = 0.0,
152  x_over_ee = 0.0,
153  y_over_ee = 0.0;
154 
155  for (size_t i=0; i<n; ++i) {
156  one_over_ee += w_filt.at(i);
157  x_over_ee += x_filt.at(i)*w_filt.at(i);
158  y_over_ee += y_filt.at(i)*w_filt.at(i);
159  }
160 
161  double up_inter = (y_over_ee-x_over_ee)/one_over_ee;
162  double dn_inter = (y_over_ee+x_over_ee)/one_over_ee;
163 
164  double up_chi2 = 0.0, dn_chi2 = 0.0;
165  for (size_t i=0; i<n; ++i) {
166  double e = 1.0/sqrt(w_filt.at(i));
167  double up_expected = up_inter + x_filt.at(i);
168  double dn_expected = dn_inter - x_filt.at(i);
169  up_chi2 += pow((y_filt.at(i)-up_expected) / e, 2.0);
170  dn_chi2 += pow((y_filt.at(i)-dn_expected) / e, 2.0);
171  }
172 
173  double prob_up = boost::math::gamma_q((double)(n-2)/2.0,up_chi2/2.0),
174  prob_dn = boost::math::gamma_q((double)(n-2)/2.0,dn_chi2/2.0);
175 
176  if (prob_up < 1e-30) prob_up = 1e-30;
177  if (prob_dn < 1e-30) prob_dn = 1e-30;
178 
179  P_up = prob_up;
180  P_dn = prob_dn;
181 
182  }
183 
184 }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
void setLLR(std::set< std::pair< rb::CellHit, rb::RecoHit >, bool(*)(std::pair< rb::CellHit, rb::RecoHit >, std::pair< rb::CellHit, rb::RecoHit >)>)
Definition: LlrUpmu.cxx:30
double getDist(double x1, double y1, double z1, double x2, double y2, double z2)
Definition: LlrUpmu.cxx:25
double P_dn_
Definition: LlrUpmu.h:42
double chi2() const
Definition: LlrUpmu.h:33
double dist
Definition: runWimpSim.h:113
double P_up() const
Definition: LlrUpmu.h:31
double slope() const
Definition: LlrUpmu.h:34
double llr_
Definition: LlrUpmu.h:45
const double a
double P_up_
Definition: LlrUpmu.h:41
std::vector< rb::RecoHit > outliers_
Definition: LlrUpmu.h:40
double getErr(double PE)
Definition: LlrUpmu.cxx:20
std::vector< rb::RecoHit > outliers() const
Definition: LlrUpmu.h:36
ifstream in
Definition: comparison.C:7
const hit & b
Definition: hits.cxx:21
double LinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &m, double &c)
Find the best-fit line to a collection of points in 2-D by minimizing the squared vertical distance f...
Definition: MathUtil.cxx:36
double P_dn() const
Definition: LlrUpmu.h:32
double chi2_
Definition: LlrUpmu.h:43
double slope_
Definition: LlrUpmu.h:44
void LLR(std::vector< double > &eT, std::vector< double > &mT, std::vector< double > &mTErr, double &slope, double &chi2, double &P_up, double &P_dn, std::vector< rb::RecoHit > &in, std::vector< rb::RecoHit > &outliers)
Definition: LlrUpmu.cxx:90
Float_t e
Definition: plot.C:35
fvar< T > gamma_q(const fvar< T > &x1, const fvar< T > &x2)
Definition: gamma_q.hpp:12