Rivet  1.8.0
MathUtils.hh
1 // -*- C++ -*-
2 #ifndef RIVET_MathUtils_HH
3 #define RIVET_MathUtils_HH
4 
5 #include "Rivet/Math/MathHeader.hh"
6 #include "Rivet/RivetBoost.hh"
7 #include <cassert>
8 
9 namespace Rivet {
10 
11 
13 
14 
17  inline bool isZero(double val, double tolerance=1E-8) {
18  return (fabs(val) < tolerance);
19  }
20 
26  inline bool isZero(long val, double UNUSED(tolerance)=1E-8) {
27  return val == 0;
28  }
29 
30 
34  inline bool fuzzyEquals(double a, double b, double tolerance=1E-5) {
35  const double absavg = (fabs(a) + fabs(b))/2.0;
36  const double absdiff = fabs(a - b);
37  const bool rtn = (isZero(a) && isZero(b)) || absdiff < tolerance*absavg;
38  // cout << a << " == " << b << "? " << rtn << endl;
39  return rtn;
40  }
41 
49  inline bool fuzzyEquals(long a, long b, double UNUSED(tolerance)=1E-5) {
50  return a == b;
51  }
52 
53 
57  inline bool fuzzyGtrEquals(double a, double b, double tolerance=1E-5) {
58  return a > b || fuzzyEquals(a, b, tolerance);
59  }
60 
68  inline bool fuzzyGtrEquals(long a, long b, double UNUSED(tolerance)=1E-5) {
69  return a >= b;
70  }
71 
72 
76  inline bool fuzzyLessEquals(double a, double b, double tolerance=1E-5) {
77  return a < b || fuzzyEquals(a, b, tolerance);
78  }
79 
87  inline bool fuzzyLessEquals(long a, long b, double UNUSED(tolerance)=1E-5) {
88  return a <= b;
89  }
90 
92 
93 
95 
96 
101  enum RangeBoundary { OPEN=0, SOFT=0, CLOSED=1, HARD=1 };
102 
103 
108  template<typename NUM>
109  inline bool inRange(NUM value, NUM low, NUM high,
110  RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
111  if (lowbound == OPEN && highbound == OPEN) {
112  return (value > low && value < high);
113  } else if (lowbound == OPEN && highbound == CLOSED) {
114  return (value > low && fuzzyLessEquals(value, high));
115  } else if (lowbound == CLOSED && highbound == OPEN) {
116  return (fuzzyGtrEquals(value, low) && value < high);
117  } else { // if (lowbound == CLOSED && highbound == CLOSED) {
118  return (fuzzyGtrEquals(value, low) && fuzzyLessEquals(value, high));
119  }
120  }
121 
123  template<typename NUM>
124  inline bool inRange(NUM value, pair<NUM, NUM> lowhigh,
125  RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
126  return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound);
127  }
128 
129 
134  inline bool inRange(int value, int low, int high,
135  RangeBoundary lowbound=CLOSED, RangeBoundary highbound=CLOSED) {
136  if (lowbound == OPEN && highbound == OPEN) {
137  return (value > low && value < high);
138  } else if (lowbound == OPEN && highbound == CLOSED) {
139  return (value > low && value <= high);
140  } else if (lowbound == CLOSED && highbound == OPEN) {
141  return (value >= low && value < high);
142  } else { // if (lowbound == CLOSED && highbound == CLOSED) {
143  return (value >= low && value <= high);
144  }
145  }
146 
148  inline bool inRange(int value, pair<int, int> lowhigh,
149  RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
150  return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound);
151  }
152 
154 
155 
157 
158 
160  template <typename NUM>
161  inline NUM sqr(NUM a) {
162  return a*a;
163  }
164 
166  template <typename Num>
167  inline Num add_quad(Num a, Num b) {
168  return sqrt(a*a + b*b);
169  }
170 
172  template <typename Num>
173  inline Num add_quad(Num a, Num b, Num c) {
174  return sqrt(a*a + b*b + c*c);
175  }
176 
178  template <typename Num>
179  inline Num intpow(Num val, unsigned int exp) {
180  assert(exp >= 0);
181  if (exp == 0) return (Num) 1;
182  else if (exp == 1) return val;
183  return val * intpow(val, exp-1);
184  }
185 
187  inline int sign(double val) {
188  if (isZero(val)) return ZERO;
189  const int valsign = (val > 0) ? PLUS : MINUS;
190  return valsign;
191  }
192 
194  inline int sign(int val) {
195  if (val == 0) return ZERO;
196  return (val > 0) ? PLUS : MINUS;
197  }
198 
200  inline int sign(long val) {
201  if (val == 0) return ZERO;
202  return (val > 0) ? PLUS : MINUS;
203  }
204 
206 
207 
209 
210 
212  inline vector<double> linspace(double start, double end, size_t nbins) {
213  assert(end >= start);
214  assert(nbins > 0);
215  vector<double> rtn;
216  const double interval = (end-start)/static_cast<double>(nbins);
217  double edge = start;
218  while (inRange(edge, start, end, CLOSED, CLOSED)) {
219  rtn.push_back(edge);
220  edge += interval;
221  }
222  assert(rtn.size() == nbins+1);
223  return rtn;
224  }
225 
226 
228  inline vector<double> logspace(double start, double end, size_t nbins) {
229  assert(end >= start);
230  assert(start > 0);
231  assert(nbins > 0);
232  const double logstart = std::log(start);
233  const double logend = std::log(end);
234  const vector<double> logvals = linspace(logstart, logend, nbins);
235  assert(logvals.size() == nbins+1);
236  vector<double> rtn;
237  foreach (double logval, logvals) {
238  rtn.push_back(std::exp(logval));
239  }
240  assert(rtn.size() == nbins+1);
241  return rtn;
242  }
243 
244 
248  template <typename NUM>
249  inline int index_between(const NUM& val, const vector<NUM>& binedges) {
250  if (!inRange(val, binedges.front(), binedges.back())) return -1; //< Out of histo range
251  int index = -1;
252  for (size_t i = 1; i < binedges.size(); ++i) {
253  if (val < binedges[i]) {
254  index = i-1;
255  break;
256  }
257  }
258  assert(inRange(index, -1, binedges.size()-1));
259  return index;
260  }
261 
263 
264 
266 
267 
269  inline double mean(const vector<int>& sample) {
270  double mean = 0.0;
271  for (size_t i=0; i<sample.size(); ++i) {
272  mean += sample[i];
273  }
274  return mean/sample.size();
275  }
276 
277  // Calculate the error on the mean, assuming poissonian errors
278  inline double mean_err(const vector<int>& sample) {
279  double mean_e = 0.0;
280  for (size_t i=0; i<sample.size(); ++i) {
281  mean_e += sqrt(sample[i]);
282  }
283  return mean_e/sample.size();
284  }
285 
287  inline double covariance(const vector<int>& sample1, const vector<int>& sample2) {
288  const double mean1 = mean(sample1);
289  const double mean2 = mean(sample2);
290  const size_t N = sample1.size();
291  double cov = 0.0;
292  for (size_t i = 0; i < N; i++) {
293  const double cov_i = (sample1[i] - mean1)*(sample2[i] - mean2);
294  cov += cov_i;
295  }
296  if (N > 1) return cov/(N-1);
297  else return 0.0;
298  }
299 
301  inline double covariance_err(const vector<int>& sample1, const vector<int>& sample2) {
302  const double mean1 = mean(sample1);
303  const double mean2 = mean(sample2);
304  const double mean1_e = mean_err(sample1);
305  const double mean2_e = mean_err(sample2);
306  const size_t N = sample1.size();
307  double cov_e = 0.0;
308  for (size_t i = 0; i < N; i++) {
309  const double cov_i = (sqrt(sample1[i]) - mean1_e)*(sample2[i] - mean2) +
310  (sample1[i] - mean1)*(sqrt(sample2[i]) - mean2_e);
311  cov_e += cov_i;
312  }
313  if (N > 1) return cov_e/(N-1);
314  else return 0.0;
315  }
316 
317 
319  inline double correlation(const vector<int>& sample1, const vector<int>& sample2) {
320  const double cov = covariance(sample1, sample2);
321  const double var1 = covariance(sample1, sample1);
322  const double var2 = covariance(sample2, sample2);
323  const double correlation = cov/sqrt(var1*var2);
324  const double corr_strength = correlation*sqrt(var2/var1);
325  return corr_strength;
326  }
327 
329  inline double correlation_err(const vector<int>& sample1, const vector<int>& sample2) {
330  const double cov = covariance(sample1, sample2);
331  const double var1 = covariance(sample1, sample1);
332  const double var2 = covariance(sample2, sample2);
333  const double cov_e = covariance_err(sample1, sample2);
334  const double var1_e = covariance_err(sample1, sample1);
335  const double var2_e = covariance_err(sample2, sample2);
336 
337  // Calculate the correlation
338  const double correlation = cov/sqrt(var1*var2);
339  // Calculate the error on the correlation
340  const double correlation_err = cov_e/sqrt(var1*var2) -
341  cov/(2*pow(3./2., var1*var2)) * (var1_e * var2 + var1 * var2_e);
342 
343 
344  // Calculate the error on the correlation strength
345  const double corr_strength_err = correlation_err*sqrt(var2/var1) +
346  correlation/(2*sqrt(var2/var1)) * (var2_e/var1 - var2*var1_e/pow(2, var2));
347 
348  return corr_strength_err;
349  }
351 
352 
354 
355 
360  inline double _mapAngleM2PITo2Pi(double angle) {
361  double rtn = fmod(angle, TWOPI);
362  if (isZero(rtn)) return 0;
363  assert(rtn >= -TWOPI && rtn <= TWOPI);
364  return rtn;
365  }
366 
368  inline double mapAngleMPiToPi(double angle) {
369  double rtn = _mapAngleM2PITo2Pi(angle);
370  if (isZero(rtn)) return 0;
371  rtn = (rtn > PI ? rtn-TWOPI :
372  rtn <= -PI ? rtn+TWOPI : rtn);
373  assert(rtn > -PI && rtn <= PI);
374  return rtn;
375  }
376 
378  inline double mapAngle0To2Pi(double angle) {
379  double rtn = _mapAngleM2PITo2Pi(angle);
380  if (isZero(rtn)) return 0;
381  if (rtn < 0) rtn += TWOPI;
382  if (rtn == TWOPI) rtn = 0;
383  assert(rtn >= 0 && rtn < TWOPI);
384  return rtn;
385  }
386 
388  inline double mapAngle0ToPi(double angle) {
389  double rtn = fabs(mapAngleMPiToPi(angle));
390  if (isZero(rtn)) return 0;
391  assert(rtn > 0 && rtn <= PI);
392  return rtn;
393  }
394 
396 
397 
399 
400 
404  inline double deltaPhi(double phi1, double phi2) {
405  return mapAngle0ToPi(phi1 - phi2);
406  }
407 
410  inline double deltaEta(double eta1, double eta2) {
411  return fabs(eta1 - eta2);
412  }
413 
416  inline double deltaR(double rap1, double phi1, double rap2, double phi2) {
417  const double dphi = deltaPhi(phi1, phi2);
418  return sqrt( sqr(rap1-rap2) + sqr(dphi) );
419  }
420 
422  inline double rapidity(double E, double pz) {
423  if (isZero(E - pz)) {
424  throw std::runtime_error("Divergent positive rapidity");
425  return MAXDOUBLE;
426  }
427  if (isZero(E + pz)) {
428  throw std::runtime_error("Divergent negative rapidity");
429  return -MAXDOUBLE;
430  }
431  return 0.5*log((E+pz)/(E-pz));
432  }
433 
435 
436 
437 }
438 
439 
440 #endif