MathUtils.hh
// -*- C++ -*-
#ifndef RIVET_MathUtils_HH
#define RIVET_MathUtils_HH
#include "Rivet/Math/MathHeader.hh"
#include "Rivet/RivetBoost.hh"
#include <cassert>
namespace Rivet {
inline bool isZero(double val, double tolerance=1E-8) {
return (fabs(val) < tolerance);
}
inline bool isZero(long val, double UNUSED(tolerance)=1E-8) {
return val == 0;
}
inline bool fuzzyEquals(double a, double b, double tolerance=1E-5) {
const double absavg = (fabs(a) + fabs(b))/2.0;
const double absdiff = fabs(a - b);
const bool rtn = (isZero(a) && isZero(b)) || absdiff < tolerance*absavg;
// cout << a << " == " << b << "? " << rtn << endl;
return rtn;
}
inline bool fuzzyEquals(long a, long b, double UNUSED(tolerance)=1E-5) {
return a == b;
}
inline bool fuzzyGtrEquals(double a, double b, double tolerance=1E-5) {
return a > b || fuzzyEquals(a, b, tolerance);
}
inline bool fuzzyGtrEquals(long a, long b, double UNUSED(tolerance)=1E-5) {
return a >= b;
}
inline bool fuzzyLessEquals(double a, double b, double tolerance=1E-5) {
return a < b || fuzzyEquals(a, b, tolerance);
}
inline bool fuzzyLessEquals(long a, long b, double UNUSED(tolerance)=1E-5) {
return a <= b;
}
enum RangeBoundary { OPEN=0, SOFT=0, CLOSED=1, HARD=1 };
template<typename NUM>
inline bool inRange(NUM value, NUM low, NUM high,
RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
if (lowbound == OPEN && highbound == OPEN) {
return (value > low && value < high);
} else if (lowbound == OPEN && highbound == CLOSED) {
return (value > low && fuzzyLessEquals(value, high));
} else if (lowbound == CLOSED && highbound == OPEN) {
return (fuzzyGtrEquals(value, low) && value < high);
} else { // if (lowbound == CLOSED && highbound == CLOSED) {
return (fuzzyGtrEquals(value, low) && fuzzyLessEquals(value, high));
}
}
template<typename NUM>
inline bool inRange(NUM value, pair<NUM, NUM> lowhigh,
RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound);
}
inline bool inRange(int value, int low, int high,
RangeBoundary lowbound=CLOSED, RangeBoundary highbound=CLOSED) {
if (lowbound == OPEN && highbound == OPEN) {
return (value > low && value < high);
} else if (lowbound == OPEN && highbound == CLOSED) {
return (value > low && value <= high);
} else if (lowbound == CLOSED && highbound == OPEN) {
return (value >= low && value < high);
} else { // if (lowbound == CLOSED && highbound == CLOSED) {
return (value >= low && value <= high);
}
}
inline bool inRange(int value, pair<int, int> lowhigh,
RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound);
}
template <typename NUM>
inline NUM sqr(NUM a) {
return a*a;
}
template <typename Num>
inline Num add_quad(Num a, Num b) {
return sqrt(a*a + b*b);
}
template <typename Num>
inline Num add_quad(Num a, Num b, Num c) {
return sqrt(a*a + b*b + c*c);
}
template <typename Num>
inline Num intpow(Num val, unsigned int exp) {
assert(exp >= 0);
if (exp == 0) return (Num) 1;
else if (exp == 1) return val;
return val * intpow(val, exp-1);
}
inline int sign(double val) {
if (isZero(val)) return ZERO;
const int valsign = (val > 0) ? PLUS : MINUS;
return valsign;
}
inline int sign(int val) {
if (val == 0) return ZERO;
return (val > 0) ? PLUS : MINUS;
}
inline int sign(long val) {
if (val == 0) return ZERO;
return (val > 0) ? PLUS : MINUS;
}
inline vector<double> linspace(size_t nbins, double start, double end) {
assert(end >= start);
assert(nbins > 0);
vector<double> rtn;
const double interval = (end-start)/static_cast<double>(nbins);
double edge = start;
while (inRange(edge, start, end, CLOSED, CLOSED)) {
rtn.push_back(edge);
edge += interval;
}
assert(rtn.size() == nbins+1);
return rtn;
}
inline vector<double> logspace(size_t nbins, double start, double end) {
assert(end >= start);
assert(start > 0);
assert(nbins > 0);
const double logstart = std::log(start);
const double logend = std::log(end);
const vector<double> logvals = linspace(nbins, logstart, logend);
assert(logvals.size() == nbins+1);
vector<double> rtn;
foreach (double logval, logvals) {
rtn.push_back(std::exp(logval));
}
assert(rtn.size() == nbins+1);
return rtn;
}
namespace BWHelpers {
inline double CDF(double x, double mu, double gamma) {
// normalize to (0;1) distribution
const double xn = (x - mu)/gamma;
return std::atan(xn)/M_PI + 0.5;
}
inline double antiCDF(double p, double mu, double gamma) {
const double xn = std::tan(M_PI*(p-0.5));
return gamma*xn + mu;
}
}
inline vector<double> BWspace(size_t nbins, double start, double end,
double mu, double gamma) {
assert(end >= start);
assert(nbins > 0);
const double pmin = BWHelpers::CDF(start,mu,gamma);
const double pmax = BWHelpers::CDF(end, mu,gamma);
const vector<double> edges = linspace(nbins, pmin, pmax);
assert(edges.size() == nbins+1);
vector<double> rtn;
foreach (double edge, edges) {
rtn.push_back(BWHelpers::antiCDF(edge,mu,gamma));
}
assert(rtn.size() == nbins+1);
return rtn;
}
template <typename NUM>
inline int index_between(const NUM& val, const vector<NUM>& binedges) {
if (!inRange(val, binedges.front(), binedges.back())) return -1; //< Out of histo range
int index = -1;
for (size_t i = 1; i < binedges.size(); ++i) {
if (val < binedges[i]) {
index = i-1;
break;
}
}
assert(inRange(index, -1, binedges.size()-1));
return index;
}
inline double mean(const vector<int>& sample) {
double mean = 0.0;
for (size_t i=0; i<sample.size(); ++i) {
mean += sample[i];
}
return mean/sample.size();
}
// Calculate the error on the mean, assuming poissonian errors
inline double mean_err(const vector<int>& sample) {
double mean_e = 0.0;
for (size_t i=0; i<sample.size(); ++i) {
mean_e += sqrt(sample[i]);
}
return mean_e/sample.size();
}
inline double covariance(const vector<int>& sample1, const vector<int>& sample2) {
const double mean1 = mean(sample1);
const double mean2 = mean(sample2);
const size_t N = sample1.size();
double cov = 0.0;
for (size_t i = 0; i < N; i++) {
const double cov_i = (sample1[i] - mean1)*(sample2[i] - mean2);
cov += cov_i;
}
if (N > 1) return cov/(N-1);
else return 0.0;
}
inline double covariance_err(const vector<int>& sample1, const vector<int>& sample2) {
const double mean1 = mean(sample1);
const double mean2 = mean(sample2);
const double mean1_e = mean_err(sample1);
const double mean2_e = mean_err(sample2);
const size_t N = sample1.size();
double cov_e = 0.0;
for (size_t i = 0; i < N; i++) {
const double cov_i = (sqrt(sample1[i]) - mean1_e)*(sample2[i] - mean2) +
(sample1[i] - mean1)*(sqrt(sample2[i]) - mean2_e);
cov_e += cov_i;
}
if (N > 1) return cov_e/(N-1);
else return 0.0;
}
inline double correlation(const vector<int>& sample1, const vector<int>& sample2) {
const double cov = covariance(sample1, sample2);
const double var1 = covariance(sample1, sample1);
const double var2 = covariance(sample2, sample2);
const double correlation = cov/sqrt(var1*var2);
const double corr_strength = correlation*sqrt(var2/var1);
return corr_strength;
}
inline double correlation_err(const vector<int>& sample1, const vector<int>& sample2) {
const double cov = covariance(sample1, sample2);
const double var1 = covariance(sample1, sample1);
const double var2 = covariance(sample2, sample2);
const double cov_e = covariance_err(sample1, sample2);
const double var1_e = covariance_err(sample1, sample1);
const double var2_e = covariance_err(sample2, sample2);
// Calculate the correlation
const double correlation = cov/sqrt(var1*var2);
// Calculate the error on the correlation
const double correlation_err = cov_e/sqrt(var1*var2) -
cov/(2*pow(3./2., var1*var2)) * (var1_e * var2 + var1 * var2_e);
// Calculate the error on the correlation strength
const double corr_strength_err = correlation_err*sqrt(var2/var1) +
correlation/(2*sqrt(var2/var1)) * (var2_e/var1 - var2*var1_e/pow(2, var2));
return corr_strength_err;
}
inline double _mapAngleM2PITo2Pi(double angle) {
double rtn = fmod(angle, TWOPI);
if (isZero(rtn)) return 0;
assert(rtn >= -TWOPI && rtn <= TWOPI);
return rtn;
}
inline double mapAngleMPiToPi(double angle) {
double rtn = _mapAngleM2PITo2Pi(angle);
if (isZero(rtn)) return 0;
rtn = (rtn > PI ? rtn-TWOPI :
rtn <= -PI ? rtn+TWOPI : rtn);
assert(rtn > -PI && rtn <= PI);
return rtn;
}
inline double mapAngle0To2Pi(double angle) {
double rtn = _mapAngleM2PITo2Pi(angle);
if (isZero(rtn)) return 0;
if (rtn < 0) rtn += TWOPI;
if (rtn == TWOPI) rtn = 0;
assert(rtn >= 0 && rtn < TWOPI);
return rtn;
}
inline double mapAngle0ToPi(double angle) {
double rtn = fabs(mapAngleMPiToPi(angle));
if (isZero(rtn)) return 0;
assert(rtn > 0 && rtn <= PI);
return rtn;
}
inline double mapAngle(double angle, PhiMapping mapping) {
switch (mapping) {
case MINUSPI_PLUSPI:
return mapAngleMPiToPi(angle);
case ZERO_2PI:
return mapAngle0To2Pi(angle);
case ZERO_PI:
return mapAngle0To2Pi(angle);
default:
throw Rivet::UserError("The specified phi mapping scheme is not implemented");
}
}
inline double deltaPhi(double phi1, double phi2) {
return mapAngle0ToPi(phi1 - phi2);
}
inline double deltaEta(double eta1, double eta2) {
return fabs(eta1 - eta2);
}
inline double deltaR(double rap1, double phi1, double rap2, double phi2) {
const double dphi = deltaPhi(phi1, phi2);
return sqrt( sqr(rap1-rap2) + sqr(dphi) );
}
inline double rapidity(double E, double pz) {
if (isZero(E - pz)) {
throw std::runtime_error("Divergent positive rapidity");
return MAXDOUBLE;
}
if (isZero(E + pz)) {
throw std::runtime_error("Divergent negative rapidity");
return -MAXDOUBLE;
}
return 0.5*log((E+pz)/(E-pz));
}
}
#endif
Rivet
Definition: MC_JetAnalysis.hh:9
Rivet::covariance_err
double covariance_err(const vector< int > &sample1, const vector< int > &sample2)
Calculate the error on the covariance (variance) of two samples, assuming poissonian errors...
Definition: MathUtils.hh:346
Rivet::mean
double mean(const vector< int > &sample)
Calculate the mean of a sample.
Definition: MathUtils.hh:314
Rivet::BWspace
vector< double > BWspace(size_t nbins, double start, double end, double mu, double gamma)
Make a list of nbins + 1 values spaced for equal area Breit-Wigner binning between start and end incl...
Definition: MathUtils.hh:273
Rivet::rapidity
double rapidity(double E, double pz)
Calculate a rapidity value from the supplied energy E and longitudinal momentum pz.
Definition: MathUtils.hh:481
Rivet::inRange
bool inRange(NUM value, NUM low, NUM high, RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN)
Determine if value is in the range low to high, for floating point numbers.
Definition: MathUtils.hh:109
Rivet::mapAngleMPiToPi
double mapAngleMPiToPi(double angle)
Map an angle into the range (-PI, PI].
Definition: MathUtils.hh:413
Rivet::add_quad
Num add_quad(Num a, Num b)
Named number-type addition in quadrature operation.
Definition: MathUtils.hh:167
Rivet::fuzzyLessEquals
bool fuzzyLessEquals(double a, double b, double tolerance=1E-5)
Compare two floating point numbers for <= with a degree of fuzziness.
Definition: MathUtils.hh:76
Rivet::PI
const double PI
A pre-defined value of .
Definition: MathHeader.hh:48
Rivet::correlation_err
double correlation_err(const vector< int > &sample1, const vector< int > &sample2)
Calculate the error of the correlation strength between two samples assuming Poissonian errors...
Definition: MathUtils.hh:374
Rivet::mapAngle
double mapAngle(double angle, PhiMapping mapping)
Map an angle into the enum-specified range.
Definition: MathUtils.hh:441
Rivet::sqr
NUM sqr(NUM a)
Named number-type squaring operation.
Definition: MathUtils.hh:161
Rivet::UserError
Error specialisation for where the problem is between the chair and computer.
Definition: Exceptions.hh:61
Rivet::PhiMapping
PhiMapping
Enum for range of to be mapped into.
Definition: MathHeader.hh:63
Rivet::angle
double angle(const Vector3 &a, const Vector3 &b)
Angle (in radians) between two 3-vectors.
Definition: Vector3.hh:244
Rivet::logspace
vector< double > logspace(size_t nbins, double start, double end)
Make a list of nbins + 1 values exponentially spaced between start and end inclusive.
Definition: MathUtils.hh:235
Rivet::TWOPI
const double TWOPI
A pre-defined value of .
Definition: MathHeader.hh:51
Rivet::isZero
bool isZero(double val, double tolerance=1E-8)
Definition: MathUtils.hh:17
Rivet::RangeBoundary
RangeBoundary
Definition: MathUtils.hh:101
Rivet::intpow
Num intpow(Num val, unsigned int exp)
A more efficient version of pow for raising numbers to integer powers.
Definition: MathUtils.hh:179
Rivet::sign
int sign(double val)
Find the sign of a number.
Definition: MathUtils.hh:187
Rivet::covariance
double covariance(const vector< int > &sample1, const vector< int > &sample2)
Calculate the covariance (variance) between two samples.
Definition: MathUtils.hh:332
Rivet::linspace
vector< double > linspace(size_t nbins, double start, double end)
Make a list of nbins + 1 values equally spaced between start and end inclusive.
Definition: MathUtils.hh:215
Rivet::mapAngle0To2Pi
double mapAngle0To2Pi(double angle)
Map an angle into the range [0, 2PI).
Definition: MathUtils.hh:423
Rivet::mapAngle0ToPi
double mapAngle0ToPi(double angle)
Map an angle into the range [0, PI].
Definition: MathUtils.hh:433
Rivet::index_between
int index_between(const NUM &val, const vector< NUM > &binedges)
Return the bin index of the given value, val, given a vector of bin edges.
Definition: MathUtils.hh:294
Rivet::fuzzyGtrEquals
bool fuzzyGtrEquals(double a, double b, double tolerance=1E-5)
Compare two floating point numbers for >= with a degree of fuzziness.
Definition: MathUtils.hh:57
Rivet::fuzzyEquals
bool fuzzyEquals(double a, double b, double tolerance=1E-5)
Compare two floating point numbers for equality with a degree of fuzziness.
Definition: MathUtils.hh:34
Rivet::correlation
double correlation(const vector< int > &sample1, const vector< int > &sample2)
Calculate the correlation strength between two samples.
Definition: MathUtils.hh:364