MathUtil.cxx
Go to the documentation of this file.
1 #include "Utilities/func/MathUtil.h"
2 
3 #include <cassert>
4 #include <algorithm>
5 
6 namespace util
7 {
8  //......................................................................
9  void LinFitUnweighted(const std::vector<double>& x,
10  const std::vector<double>& y,
11  double& m, double& c)
12  {
13  // Before going ahead, make sure we have sensible arrays
14  assert(x.size() == y.size());
15  assert(x.size() >= 2);
16 
17  // Accumulate the sums for the fit
18  double Sx = 0;
19  double Sy = 0;
20  double Sxy = 0;
21  double Sy2 = 0;
22  double Sx2 = 0;
23  const unsigned int I = x.size();
24  for(unsigned int i = 0; i < I; ++i) {
25  Sx += x[i];
26  Sy += y[i];
27  Sx2 += x[i]*x[i];
28  Sxy += x[i]*y[i];
29  Sy2 += y[i]*y[i];
30  }
31  const double d = I*Sx2 - Sx*Sx;
32  m = (I*Sxy - Sx*Sy)/d;
33  c = (Sy*Sx2 - Sx*Sxy)/d;
34  }
35 
36  //......................................................................
37  double LinFit(const std::vector<double>& x,
38  const std::vector<double>& y,
39  const std::vector<double>& w,
40  double& m, double& c)
41  {
42  // Before going ahead, make sure we have sensible arrays
43  assert(x.size() == y.size());
44  assert(x.size() == w.size());
45  assert(x.size() >= 2);
46 
47  // Accumulate the sums for the fit
48  double Sw = 0;
49  double Swx = 0;
50  double Swy = 0;
51  double Swxy = 0;
52  double Swy2 = 0;
53  double Swx2 = 0;
54  for(unsigned int i = 0; i < w.size(); ++i) {
55  Sw += w[i];
56  Swx += w[i]*x[i];
57  Swy += w[i]*y[i];
58  Swx2 += w[i]*x[i]*x[i];
59  Swxy += w[i]*x[i]*y[i];
60  Swy2 += w[i]*y[i]*y[i];
61  }
62  const double d = Sw*Swx2 - Swx*Swx;
63  m = (Sw*Swxy - Swx*Swy)/d;
64  c = (Swy*Swx2 - Swx*Swxy)/d;
65 
66  const double chi2 =
67  Swy2 - 2.0*m*Swxy - 2.0*c*Swy + 2.0*m*c*Swx +
68  c*c*Sw + m*m*Swx2;
69 
70  return chi2;
71  }
72 
73  //......................................................................
74 
75  static bool sort_pair_by_first(const std::pair<double,double>& a,
76  const std::pair<double,double>& b)
77  {
78  return (a.first<b.first);
79  }
80 
81  //......................................................................
82 
83  static double find_median(const std::vector< std::pair<double,double> >& a)
84  {
85  unsigned int i;
86  for (i=0; i<a.size(); ++i) {
87  if (a[i].second>0.5) break;
88  }
89  if (i>1) {
90  //
91  // Linearly interpolate along weights to find 0.5
92  //
93  double x1 = a[i-1].second;
94  double x2 = a[i]. second;
95  double y1 = a[i-1].first;
96  double y2 = a[i]. first;
97  return y1 + (0.5-x1)*(y2-y1)/(x2-x1);
98  }
99  return a[i].first;
100  }
101 
102  //......................................................................
103 
104  void LinFitTS(const std::vector<double>& x,
105  const std::vector<double>& y,
106  const std::vector<double>& wt,
107  double& m,
108  double& b,
109  int wmode)
110  {
111  unsigned int i, j;
112 
113  //
114  // Check that we have at least two points and that the x and y
115  // vectors sizes match
116  //
117  unsigned int n = x.size();
118  assert(n>1);
119  assert(n==y.size());
120 
121  //
122  // Reserve a vector of slope,weight pairs. Reserve the space we need
123  // up front to improve efficiency
124  //
125  std::vector<std::pair<double,double> > mw;
126  mw.reserve(n*(n-1)-1);
127 
128  //
129  // Find the slopes of all pairs of points x,y. Depending on calling
130  // options, weight the pairs by their separation distance.
131  //
132  double dx, dy, w;
133  double wnorm = 0.0;
134  for (i=0; (i+1)<n; ++i) {
135  for (j=i+1; j<n; ++j) {
136  dx = x[j]-x[i];
137  dy = y[j]-y[i];
138  if (dx!=0.0) {
139  m = dy/dx;
140  switch (wmode) {
141  case 1: w = sqrt(wt[i]*wt[j]); break;
142  case 2: w = sqrt(wt[i]*wt[j]*(dx*dx+dy*dy)); break;
143  default: w = 1.0; break;
144  }
145  if (w>0.0) {
146  mw.push_back( std::pair<double,double>(m,w) );
147  wnorm += w;
148  }
149  }
150  }
151  }
152  std::sort(mw.begin(), mw.end(), sort_pair_by_first);
153 
154  //
155  // Normalize weights to accumulate along the interval between 0 and 1
156  //
157  double sumw = 0.0;
158  wnorm = 1.0/wnorm;
159  for (i=0; i<mw.size(); ++i) {
160  mw[i].second = wnorm*mw[i].second + sumw;
161  sumw = mw[i].second;
162  }
163 
164  m = find_median(mw);
165 
166  //
167  // Build a vector or intercept, weight pairs
168  //
169  std::vector<std::pair<double,double> > bw;
170  bw.reserve(n);
171 
172  //
173  // Now compute the intercepts
174  //
175  wnorm = 0.0;
176  w = 1.0;
177  for (i=0; i<n; ++i) {
178  b = y[i] - m*x[i];
179  bw.push_back( std::pair<double,double>(b,w) );
180  wnorm += w;
181  }
182  std::sort(bw.begin(), bw.end(), sort_pair_by_first);
183 
184  //
185  // Normalize the weights to accumulate along a range between 0 and
186  // 1
187  //
188  sumw = 0.0;
189  wnorm = 1.0/wnorm;
190  for (i=0; i<bw.size(); ++i) {
191  bw[i].second = wnorm*bw[i].second + sumw;
192  sumw = bw[i].second;
193  }
194 
195  b = find_median(bw);
196  }
197 
198 } // namespace
static bool sort_pair_by_first(const std::pair< double, double > &a, const std::pair< double, double > &b)
Definition: MathUtil.cxx:75
Filter events based on their run/event numbers.
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
double chi2()
double dy[NP][NC]
double dx[NP][NC]
const double a
Float_t d
Definition: plot.C:236
const double j
Definition: BetheBloch.cxx:29
static double find_median(const std::vector< std::pair< double, double > > &a)
Definition: MathUtil.cxx:83
void LinFitUnweighted(const std::vector< double > &x, const std::vector< double > &y, double &m, double &c)
Simplified version of LinFit.
Definition: MathUtil.cxx:8
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
assert(nhit_max >=nhit_nbins)
Float_t w
Definition: plot.C:20
void LinFitTS(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &wt, double &m, double &b, int wmode)
Best fit line to points using Theil-Sens median method.
Definition: MathUtil.cxx:104