Rivet  1.8.3
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 
215  inline vector<double> linspace(size_t nbins, double start, double end) {
216  assert(end >= start);
217  assert(nbins > 0);
218  vector<double> rtn;
219  const double interval = (end-start)/static_cast<double>(nbins);
220  double edge = start;
221  while (inRange(edge, start, end, CLOSED, CLOSED)) {
222  rtn.push_back(edge);
223  edge += interval;
224  }
225  assert(rtn.size() == nbins+1);
226  return rtn;
227  }
228 
229 
235  inline vector<double> logspace(size_t nbins, double start, double end) {
236  assert(end >= start);
237  assert(start > 0);
238  assert(nbins > 0);
239  const double logstart = std::log(start);
240  const double logend = std::log(end);
241  const vector<double> logvals = linspace(nbins, logstart, logend);
242  assert(logvals.size() == nbins+1);
243  vector<double> rtn;
244  foreach (double logval, logvals) {
245  rtn.push_back(std::exp(logval));
246  }
247  assert(rtn.size() == nbins+1);
248  return rtn;
249  }
250 
251  namespace BWHelpers {
253  inline double CDF(double x, double mu, double gamma) {
254  // normalize to (0;1) distribution
255  const double xn = (x - mu)/gamma;
256  return std::atan(xn)/M_PI + 0.5;
257  }
258 
260  inline double antiCDF(double p, double mu, double gamma) {
261  const double xn = std::tan(M_PI*(p-0.5));
262  return gamma*xn + mu;
263  }
264  }
265 
273  inline vector<double> BWspace(size_t nbins, double start, double end,
274  double mu, double gamma) {
275  assert(end >= start);
276  assert(nbins > 0);
277  const double pmin = BWHelpers::CDF(start,mu,gamma);
278  const double pmax = BWHelpers::CDF(end, mu,gamma);
279  const vector<double> edges = linspace(nbins, pmin, pmax);
280  assert(edges.size() == nbins+1);
281  vector<double> rtn;
282  foreach (double edge, edges) {
283  rtn.push_back(BWHelpers::antiCDF(edge,mu,gamma));
284  }
285  assert(rtn.size() == nbins+1);
286  return rtn;
287  }
288 
289 
293  template <typename NUM>
294  inline int index_between(const NUM& val, const vector<NUM>& binedges) {
295  if (!inRange(val, binedges.front(), binedges.back())) return -1; //< Out of histo range
296  int index = -1;
297  for (size_t i = 1; i < binedges.size(); ++i) {
298  if (val < binedges[i]) {
299  index = i-1;
300  break;
301  }
302  }
303  assert(inRange(index, -1, binedges.size()-1));
304  return index;
305  }
306 
308 
309 
311 
312 
314  inline double mean(const vector<int>& sample) {
315  double mean = 0.0;
316  for (size_t i=0; i<sample.size(); ++i) {
317  mean += sample[i];
318  }
319  return mean/sample.size();
320  }
321 
322  // Calculate the error on the mean, assuming poissonian errors
323  inline double mean_err(const vector<int>& sample) {
324  double mean_e = 0.0;
325  for (size_t i=0; i<sample.size(); ++i) {
326  mean_e += sqrt(sample[i]);
327  }
328  return mean_e/sample.size();
329  }
330 
332  inline double covariance(const vector<int>& sample1, const vector<int>& sample2) {
333  const double mean1 = mean(sample1);
334  const double mean2 = mean(sample2);
335  const size_t N = sample1.size();
336  double cov = 0.0;
337  for (size_t i = 0; i < N; i++) {
338  const double cov_i = (sample1[i] - mean1)*(sample2[i] - mean2);
339  cov += cov_i;
340  }
341  if (N > 1) return cov/(N-1);
342  else return 0.0;
343  }
344 
346  inline double covariance_err(const vector<int>& sample1, const vector<int>& sample2) {
347  const double mean1 = mean(sample1);
348  const double mean2 = mean(sample2);
349  const double mean1_e = mean_err(sample1);
350  const double mean2_e = mean_err(sample2);
351  const size_t N = sample1.size();
352  double cov_e = 0.0;
353  for (size_t i = 0; i < N; i++) {
354  const double cov_i = (sqrt(sample1[i]) - mean1_e)*(sample2[i] - mean2) +
355  (sample1[i] - mean1)*(sqrt(sample2[i]) - mean2_e);
356  cov_e += cov_i;
357  }
358  if (N > 1) return cov_e/(N-1);
359  else return 0.0;
360  }
361 
362 
364  inline double correlation(const vector<int>& sample1, const vector<int>& sample2) {
365  const double cov = covariance(sample1, sample2);
366  const double var1 = covariance(sample1, sample1);
367  const double var2 = covariance(sample2, sample2);
368  const double correlation = cov/sqrt(var1*var2);
369  const double corr_strength = correlation*sqrt(var2/var1);
370  return corr_strength;
371  }
372 
374  inline double correlation_err(const vector<int>& sample1, const vector<int>& sample2) {
375  const double cov = covariance(sample1, sample2);
376  const double var1 = covariance(sample1, sample1);
377  const double var2 = covariance(sample2, sample2);
378  const double cov_e = covariance_err(sample1, sample2);
379  const double var1_e = covariance_err(sample1, sample1);
380  const double var2_e = covariance_err(sample2, sample2);
381 
382  // Calculate the correlation
383  const double correlation = cov/sqrt(var1*var2);
384  // Calculate the error on the correlation
385  const double correlation_err = cov_e/sqrt(var1*var2) -
386  cov/(2*pow(3./2., var1*var2)) * (var1_e * var2 + var1 * var2_e);
387 
388 
389  // Calculate the error on the correlation strength
390  const double corr_strength_err = correlation_err*sqrt(var2/var1) +
391  correlation/(2*sqrt(var2/var1)) * (var2_e/var1 - var2*var1_e/pow(2, var2));
392 
393  return corr_strength_err;
394  }
396 
397 
399 
400 
405  inline double _mapAngleM2PITo2Pi(double angle) {
406  double rtn = fmod(angle, TWOPI);
407  if (isZero(rtn)) return 0;
408  assert(rtn >= -TWOPI && rtn <= TWOPI);
409  return rtn;
410  }
411 
413  inline double mapAngleMPiToPi(double angle) {
414  double rtn = _mapAngleM2PITo2Pi(angle);
415  if (isZero(rtn)) return 0;
416  rtn = (rtn > PI ? rtn-TWOPI :
417  rtn <= -PI ? rtn+TWOPI : rtn);
418  assert(rtn > -PI && rtn <= PI);
419  return rtn;
420  }
421 
423  inline double mapAngle0To2Pi(double angle) {
424  double rtn = _mapAngleM2PITo2Pi(angle);
425  if (isZero(rtn)) return 0;
426  if (rtn < 0) rtn += TWOPI;
427  if (rtn == TWOPI) rtn = 0;
428  assert(rtn >= 0 && rtn < TWOPI);
429  return rtn;
430  }
431 
433  inline double mapAngle0ToPi(double angle) {
434  double rtn = fabs(mapAngleMPiToPi(angle));
435  if (isZero(rtn)) return 0;
436  assert(rtn > 0 && rtn <= PI);
437  return rtn;
438  }
439 
441  inline double mapAngle(double angle, PhiMapping mapping) {
442  switch (mapping) {
443  case MINUSPI_PLUSPI:
444  return mapAngleMPiToPi(angle);
445  case ZERO_2PI:
446  return mapAngle0To2Pi(angle);
447  case ZERO_PI:
448  return mapAngle0To2Pi(angle);
449  default:
450  throw Rivet::UserError("The specified phi mapping scheme is not implemented");
451  }
452  }
453 
455 
456 
458 
459 
463  inline double deltaPhi(double phi1, double phi2) {
464  return mapAngle0ToPi(phi1 - phi2);
465  }
466 
469  inline double deltaEta(double eta1, double eta2) {
470  return fabs(eta1 - eta2);
471  }
472 
475  inline double deltaR(double rap1, double phi1, double rap2, double phi2) {
476  const double dphi = deltaPhi(phi1, phi2);
477  return sqrt( sqr(rap1-rap2) + sqr(dphi) );
478  }
479 
481  inline double rapidity(double E, double pz) {
482  if (isZero(E - pz)) {
483  throw std::runtime_error("Divergent positive rapidity");
484  return MAXDOUBLE;
485  }
486  if (isZero(E + pz)) {
487  throw std::runtime_error("Divergent negative rapidity");
488  return -MAXDOUBLE;
489  }
490  return 0.5*log((E+pz)/(E-pz));
491  }
492 
494 
495 
496 }
497 
498 
499 #endif