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;
212 inline vector<double>
linspace(
double start,
double end,
size_t nbins) {
213 assert(end >= start);
216 const double interval = (end-start)/static_cast<double>(nbins);
218 while (
inRange(edge, start, end, CLOSED, CLOSED)) {
222 assert(rtn.size() == nbins+1);
228 inline vector<double>
logspace(
double start,
double end,
size_t nbins) {
229 assert(end >= start);
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);
237 foreach (
double logval, logvals) {
238 rtn.push_back(std::exp(logval));
240 assert(rtn.size() == nbins+1);
248 template <
typename NUM>
250 if (!
inRange(val, binedges.front(), binedges.back()))
return -1;
252 for (
size_t i = 1; i < binedges.size(); ++i) {
253 if (val < binedges[i]) {
258 assert(
inRange(index, -1, binedges.size()-1));
269 inline double mean(
const vector<int>& sample) {
271 for (
size_t i=0; i<sample.size(); ++i) {
274 return mean/sample.size();
278 inline double mean_err(
const vector<int>& sample) {
280 for (
size_t i=0; i<sample.size(); ++i) {
281 mean_e += sqrt(sample[i]);
283 return mean_e/sample.size();
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();
292 for (
size_t i = 0; i < N; i++) {
293 const double cov_i = (sample1[i] - mean1)*(sample2[i] - mean2);
296 if (N > 1)
return cov/(N-1);
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();
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);
313 if (N > 1)
return cov_e/(N-1);
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);
324 const double corr_strength = correlation*sqrt(var2/var1);
325 return corr_strength;
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);
341 cov/(2*pow(3./2., var1*var2)) * (var1_e * var2 + var1 * var2_e);
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));
348 return corr_strength_err;
360 inline double _mapAngleM2PITo2Pi(
double angle) {
361 double rtn = fmod(angle,
TWOPI);
362 if (
isZero(rtn))
return 0;
369 double rtn = _mapAngleM2PITo2Pi(angle);
370 if (
isZero(rtn))
return 0;
373 assert(rtn > -
PI && rtn <=
PI);
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);
390 if (
isZero(rtn))
return 0;
391 assert(rtn > 0 && rtn <=
PI);
404 inline double deltaPhi(
double phi1,
double phi2) {
410 inline double deltaEta(
double eta1,
double eta2) {
411 return fabs(eta1 - eta2);
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) );
424 throw std::runtime_error(
"Divergent positive rapidity");
428 throw std::runtime_error(
"Divergent negative rapidity");
431 return 0.5*log((E+pz)/(E-pz));