2 #ifndef RIVET_MathUtils_HH
3 #define RIVET_MathUtils_HH
5 #include "Rivet/Math/MathHeader.hh"
6 #include "Rivet/RivetBoost.hh"
17 inline bool isZero(
double val,
double tolerance=1E-8) {
18 return (fabs(val) < tolerance);
26 inline bool isZero(
long val,
double UNUSED(tolerance)=1E-8) {
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;
49 inline bool fuzzyEquals(
long a,
long b,
double UNUSED(tolerance)=1E-5) {
108 template<
typename NUM>
109 inline bool inRange(NUM value, NUM low, NUM high,
111 if (lowbound == OPEN && highbound == OPEN) {
112 return (value > low && value < high);
113 }
else if (lowbound == OPEN && highbound == CLOSED) {
115 }
else if (lowbound == CLOSED && highbound == OPEN) {
123 template<
typename NUM>
124 inline bool inRange(NUM value, pair<NUM, NUM> lowhigh,
126 return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound);
134 inline bool inRange(
int value,
int low,
int high,
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);
143 return (value >= low && value <= high);
148 inline bool inRange(
int value, pair<int, int> lowhigh,
150 return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound);
160 template <
typename NUM>
166 template <
typename Num>
168 return sqrt(a*a + b*b);
172 template <
typename Num>
174 return sqrt(a*a + b*b + c*c);
178 template <
typename Num>
179 inline Num
intpow(Num val,
unsigned int exp) {
181 if (exp == 0)
return (Num) 1;
182 else if (exp == 1)
return val;
183 return val *
intpow(val, exp-1);
188 if (
isZero(val))
return ZERO;
189 const int valsign = (val > 0) ? PLUS : MINUS;
195 if (val == 0)
return ZERO;
196 return (val > 0) ? PLUS : MINUS;
201 if (val == 0)
return ZERO;
202 return (val > 0) ? PLUS : MINUS;
215 inline vector<double>
linspace(
size_t nbins,
double start,
double end) {
216 assert(end >= start);
219 const double interval = (end-start)/static_cast<double>(nbins);
221 while (
inRange(edge, start, end, CLOSED, CLOSED)) {
225 assert(rtn.size() == nbins+1);
235 inline vector<double>
logspace(
size_t nbins,
double start,
double end) {
236 assert(end >= start);
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);
244 foreach (
double logval, logvals) {
245 rtn.push_back(std::exp(logval));
247 assert(rtn.size() == nbins+1);
251 namespace BWHelpers {
253 inline double CDF(
double x,
double mu,
double gamma) {
255 const double xn = (x - mu)/gamma;
256 return std::atan(xn)/M_PI + 0.5;
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;
273 inline vector<double>
BWspace(
size_t nbins,
double start,
double end,
274 double mu,
double gamma) {
275 assert(end >= start);
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);
282 foreach (
double edge, edges) {
283 rtn.push_back(BWHelpers::antiCDF(edge,mu,gamma));
285 assert(rtn.size() == nbins+1);
293 template <
typename NUM>
295 if (!
inRange(val, binedges.front(), binedges.back()))
return -1;
297 for (
size_t i = 1; i < binedges.size(); ++i) {
298 if (val < binedges[i]) {
303 assert(
inRange(index, -1, binedges.size()-1));
314 inline double mean(
const vector<int>& sample) {
316 for (
size_t i=0; i<sample.size(); ++i) {
319 return mean/sample.size();
323 inline double mean_err(
const vector<int>& sample) {
325 for (
size_t i=0; i<sample.size(); ++i) {
326 mean_e += sqrt(sample[i]);
328 return mean_e/sample.size();
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();
337 for (
size_t i = 0; i < N; i++) {
338 const double cov_i = (sample1[i] - mean1)*(sample2[i] - mean2);
341 if (N > 1)
return cov/(N-1);
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();
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);
358 if (N > 1)
return cov_e/(N-1);
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);
369 const double corr_strength = correlation*sqrt(var2/var1);
370 return corr_strength;
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);
386 cov/(2*pow(3./2., var1*var2)) * (var1_e * var2 + var1 * var2_e);
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));
393 return corr_strength_err;
405 inline double _mapAngleM2PITo2Pi(
double angle) {
406 double rtn = fmod(angle,
TWOPI);
407 if (
isZero(rtn))
return 0;
414 double rtn = _mapAngleM2PITo2Pi(angle);
415 if (
isZero(rtn))
return 0;
418 assert(rtn > -
PI && rtn <=
PI);
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);
435 if (
isZero(rtn))
return 0;
436 assert(rtn > 0 && rtn <=
PI);
450 throw Rivet::UserError(
"The specified phi mapping scheme is not implemented");
463 inline double deltaPhi(
double phi1,
double phi2) {
469 inline double deltaEta(
double eta1,
double eta2) {
470 return fabs(eta1 - eta2);
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) );
483 throw std::runtime_error(
"Divergent positive rapidity");
487 throw std::runtime_error(
"Divergent negative rapidity");
490 return 0.5*log((E+pz)/(E-pz));
Definition: MC_JetAnalysis.hh:9
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
double mean(const vector< int > &sample)
Calculate the mean of a sample.
Definition: MathUtils.hh:314
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
double rapidity(double E, double pz)
Calculate a rapidity value from the supplied energy E and longitudinal momentum pz.
Definition: MathUtils.hh:481
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
double mapAngleMPiToPi(double angle)
Map an angle into the range (-PI, PI].
Definition: MathUtils.hh:413
Num add_quad(Num a, Num b)
Named number-type addition in quadrature operation.
Definition: MathUtils.hh:167
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
const double PI
A pre-defined value of .
Definition: MathHeader.hh:48
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
double mapAngle(double angle, PhiMapping mapping)
Map an angle into the enum-specified range.
Definition: MathUtils.hh:441
NUM sqr(NUM a)
Named number-type squaring operation.
Definition: MathUtils.hh:161
Error specialisation for where the problem is between the chair and computer.
Definition: Exceptions.hh:61
PhiMapping
Enum for range of to be mapped into.
Definition: MathHeader.hh:63
double angle(const Vector3 &a, const Vector3 &b)
Angle (in radians) between two 3-vectors.
Definition: Vector3.hh:244
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
const double TWOPI
A pre-defined value of .
Definition: MathHeader.hh:51
bool isZero(double val, double tolerance=1E-8)
Definition: MathUtils.hh:17
RangeBoundary
Definition: MathUtils.hh:101
Num intpow(Num val, unsigned int exp)
A more efficient version of pow for raising numbers to integer powers.
Definition: MathUtils.hh:179
int sign(double val)
Find the sign of a number.
Definition: MathUtils.hh:187
double covariance(const vector< int > &sample1, const vector< int > &sample2)
Calculate the covariance (variance) between two samples.
Definition: MathUtils.hh:332
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
double mapAngle0To2Pi(double angle)
Map an angle into the range [0, 2PI).
Definition: MathUtils.hh:423
double mapAngle0ToPi(double angle)
Map an angle into the range [0, PI].
Definition: MathUtils.hh:433
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
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
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
double correlation(const vector< int > &sample1, const vector< int > &sample2)
Calculate the correlation strength between two samples.
Definition: MathUtils.hh:364