Rivet
1.8.0
|
00001 // -*- C++ -*- 00002 #ifndef RIVET_MathUtils_HH 00003 #define RIVET_MathUtils_HH 00004 00005 #include "Rivet/Math/MathHeader.hh" 00006 #include "Rivet/RivetBoost.hh" 00007 #include <cassert> 00008 00009 namespace Rivet { 00010 00011 00013 00014 00017 inline bool isZero(double val, double tolerance=1E-8) { 00018 return (fabs(val) < tolerance); 00019 } 00020 00026 inline bool isZero(long val, double UNUSED(tolerance)=1E-8) { 00027 return val == 0; 00028 } 00029 00030 00034 inline bool fuzzyEquals(double a, double b, double tolerance=1E-5) { 00035 const double absavg = (fabs(a) + fabs(b))/2.0; 00036 const double absdiff = fabs(a - b); 00037 const bool rtn = (isZero(a) && isZero(b)) || absdiff < tolerance*absavg; 00038 // cout << a << " == " << b << "? " << rtn << endl; 00039 return rtn; 00040 } 00041 00049 inline bool fuzzyEquals(long a, long b, double UNUSED(tolerance)=1E-5) { 00050 return a == b; 00051 } 00052 00053 00057 inline bool fuzzyGtrEquals(double a, double b, double tolerance=1E-5) { 00058 return a > b || fuzzyEquals(a, b, tolerance); 00059 } 00060 00068 inline bool fuzzyGtrEquals(long a, long b, double UNUSED(tolerance)=1E-5) { 00069 return a >= b; 00070 } 00071 00072 00076 inline bool fuzzyLessEquals(double a, double b, double tolerance=1E-5) { 00077 return a < b || fuzzyEquals(a, b, tolerance); 00078 } 00079 00087 inline bool fuzzyLessEquals(long a, long b, double UNUSED(tolerance)=1E-5) { 00088 return a <= b; 00089 } 00090 00092 00093 00095 00096 00101 enum RangeBoundary { OPEN=0, SOFT=0, CLOSED=1, HARD=1 }; 00102 00103 00108 template<typename NUM> 00109 inline bool inRange(NUM value, NUM low, NUM high, 00110 RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) { 00111 if (lowbound == OPEN && highbound == OPEN) { 00112 return (value > low && value < high); 00113 } else if (lowbound == OPEN && highbound == CLOSED) { 00114 return (value > low && fuzzyLessEquals(value, high)); 00115 } else if (lowbound == CLOSED && highbound == OPEN) { 00116 return (fuzzyGtrEquals(value, low) && value < high); 00117 } else { // if (lowbound == CLOSED && highbound == CLOSED) { 00118 return (fuzzyGtrEquals(value, low) && fuzzyLessEquals(value, high)); 00119 } 00120 } 00121 00123 template<typename NUM> 00124 inline bool inRange(NUM value, pair<NUM, NUM> lowhigh, 00125 RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) { 00126 return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound); 00127 } 00128 00129 00134 inline bool inRange(int value, int low, int high, 00135 RangeBoundary lowbound=CLOSED, RangeBoundary highbound=CLOSED) { 00136 if (lowbound == OPEN && highbound == OPEN) { 00137 return (value > low && value < high); 00138 } else if (lowbound == OPEN && highbound == CLOSED) { 00139 return (value > low && value <= high); 00140 } else if (lowbound == CLOSED && highbound == OPEN) { 00141 return (value >= low && value < high); 00142 } else { // if (lowbound == CLOSED && highbound == CLOSED) { 00143 return (value >= low && value <= high); 00144 } 00145 } 00146 00148 inline bool inRange(int value, pair<int, int> lowhigh, 00149 RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) { 00150 return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound); 00151 } 00152 00154 00155 00157 00158 00160 template <typename NUM> 00161 inline NUM sqr(NUM a) { 00162 return a*a; 00163 } 00164 00166 template <typename Num> 00167 inline Num add_quad(Num a, Num b) { 00168 return sqrt(a*a + b*b); 00169 } 00170 00172 template <typename Num> 00173 inline Num add_quad(Num a, Num b, Num c) { 00174 return sqrt(a*a + b*b + c*c); 00175 } 00176 00178 template <typename Num> 00179 inline Num intpow(Num val, unsigned int exp) { 00180 assert(exp >= 0); 00181 if (exp == 0) return (Num) 1; 00182 else if (exp == 1) return val; 00183 return val * intpow(val, exp-1); 00184 } 00185 00187 inline int sign(double val) { 00188 if (isZero(val)) return ZERO; 00189 const int valsign = (val > 0) ? PLUS : MINUS; 00190 return valsign; 00191 } 00192 00194 inline int sign(int val) { 00195 if (val == 0) return ZERO; 00196 return (val > 0) ? PLUS : MINUS; 00197 } 00198 00200 inline int sign(long val) { 00201 if (val == 0) return ZERO; 00202 return (val > 0) ? PLUS : MINUS; 00203 } 00204 00206 00207 00209 00210 00212 inline vector<double> linspace(double start, double end, size_t nbins) { 00213 assert(end >= start); 00214 assert(nbins > 0); 00215 vector<double> rtn; 00216 const double interval = (end-start)/static_cast<double>(nbins); 00217 double edge = start; 00218 while (inRange(edge, start, end, CLOSED, CLOSED)) { 00219 rtn.push_back(edge); 00220 edge += interval; 00221 } 00222 assert(rtn.size() == nbins+1); 00223 return rtn; 00224 } 00225 00226 00228 inline vector<double> logspace(double start, double end, size_t nbins) { 00229 assert(end >= start); 00230 assert(start > 0); 00231 assert(nbins > 0); 00232 const double logstart = std::log(start); 00233 const double logend = std::log(end); 00234 const vector<double> logvals = linspace(logstart, logend, nbins); 00235 assert(logvals.size() == nbins+1); 00236 vector<double> rtn; 00237 foreach (double logval, logvals) { 00238 rtn.push_back(std::exp(logval)); 00239 } 00240 assert(rtn.size() == nbins+1); 00241 return rtn; 00242 } 00243 00244 00248 template <typename NUM> 00249 inline int index_between(const NUM& val, const vector<NUM>& binedges) { 00250 if (!inRange(val, binedges.front(), binedges.back())) return -1; //< Out of histo range 00251 int index = -1; 00252 for (size_t i = 1; i < binedges.size(); ++i) { 00253 if (val < binedges[i]) { 00254 index = i-1; 00255 break; 00256 } 00257 } 00258 assert(inRange(index, -1, binedges.size()-1)); 00259 return index; 00260 } 00261 00263 00264 00266 00267 00269 inline double mean(const vector<int>& sample) { 00270 double mean = 0.0; 00271 for (size_t i=0; i<sample.size(); ++i) { 00272 mean += sample[i]; 00273 } 00274 return mean/sample.size(); 00275 } 00276 00277 // Calculate the error on the mean, assuming poissonian errors 00278 inline double mean_err(const vector<int>& sample) { 00279 double mean_e = 0.0; 00280 for (size_t i=0; i<sample.size(); ++i) { 00281 mean_e += sqrt(sample[i]); 00282 } 00283 return mean_e/sample.size(); 00284 } 00285 00287 inline double covariance(const vector<int>& sample1, const vector<int>& sample2) { 00288 const double mean1 = mean(sample1); 00289 const double mean2 = mean(sample2); 00290 const size_t N = sample1.size(); 00291 double cov = 0.0; 00292 for (size_t i = 0; i < N; i++) { 00293 const double cov_i = (sample1[i] - mean1)*(sample2[i] - mean2); 00294 cov += cov_i; 00295 } 00296 if (N > 1) return cov/(N-1); 00297 else return 0.0; 00298 } 00299 00301 inline double covariance_err(const vector<int>& sample1, const vector<int>& sample2) { 00302 const double mean1 = mean(sample1); 00303 const double mean2 = mean(sample2); 00304 const double mean1_e = mean_err(sample1); 00305 const double mean2_e = mean_err(sample2); 00306 const size_t N = sample1.size(); 00307 double cov_e = 0.0; 00308 for (size_t i = 0; i < N; i++) { 00309 const double cov_i = (sqrt(sample1[i]) - mean1_e)*(sample2[i] - mean2) + 00310 (sample1[i] - mean1)*(sqrt(sample2[i]) - mean2_e); 00311 cov_e += cov_i; 00312 } 00313 if (N > 1) return cov_e/(N-1); 00314 else return 0.0; 00315 } 00316 00317 00319 inline double correlation(const vector<int>& sample1, const vector<int>& sample2) { 00320 const double cov = covariance(sample1, sample2); 00321 const double var1 = covariance(sample1, sample1); 00322 const double var2 = covariance(sample2, sample2); 00323 const double correlation = cov/sqrt(var1*var2); 00324 const double corr_strength = correlation*sqrt(var2/var1); 00325 return corr_strength; 00326 } 00327 00329 inline double correlation_err(const vector<int>& sample1, const vector<int>& sample2) { 00330 const double cov = covariance(sample1, sample2); 00331 const double var1 = covariance(sample1, sample1); 00332 const double var2 = covariance(sample2, sample2); 00333 const double cov_e = covariance_err(sample1, sample2); 00334 const double var1_e = covariance_err(sample1, sample1); 00335 const double var2_e = covariance_err(sample2, sample2); 00336 00337 // Calculate the correlation 00338 const double correlation = cov/sqrt(var1*var2); 00339 // Calculate the error on the correlation 00340 const double correlation_err = cov_e/sqrt(var1*var2) - 00341 cov/(2*pow(3./2., var1*var2)) * (var1_e * var2 + var1 * var2_e); 00342 00343 00344 // Calculate the error on the correlation strength 00345 const double corr_strength_err = correlation_err*sqrt(var2/var1) + 00346 correlation/(2*sqrt(var2/var1)) * (var2_e/var1 - var2*var1_e/pow(2, var2)); 00347 00348 return corr_strength_err; 00349 } 00351 00352 00354 00355 00360 inline double _mapAngleM2PITo2Pi(double angle) { 00361 double rtn = fmod(angle, TWOPI); 00362 if (isZero(rtn)) return 0; 00363 assert(rtn >= -TWOPI && rtn <= TWOPI); 00364 return rtn; 00365 } 00366 00368 inline double mapAngleMPiToPi(double angle) { 00369 double rtn = _mapAngleM2PITo2Pi(angle); 00370 if (isZero(rtn)) return 0; 00371 rtn = (rtn > PI ? rtn-TWOPI : 00372 rtn <= -PI ? rtn+TWOPI : rtn); 00373 assert(rtn > -PI && rtn <= PI); 00374 return rtn; 00375 } 00376 00378 inline double mapAngle0To2Pi(double angle) { 00379 double rtn = _mapAngleM2PITo2Pi(angle); 00380 if (isZero(rtn)) return 0; 00381 if (rtn < 0) rtn += TWOPI; 00382 if (rtn == TWOPI) rtn = 0; 00383 assert(rtn >= 0 && rtn < TWOPI); 00384 return rtn; 00385 } 00386 00388 inline double mapAngle0ToPi(double angle) { 00389 double rtn = fabs(mapAngleMPiToPi(angle)); 00390 if (isZero(rtn)) return 0; 00391 assert(rtn > 0 && rtn <= PI); 00392 return rtn; 00393 } 00394 00396 00397 00399 00400 00404 inline double deltaPhi(double phi1, double phi2) { 00405 return mapAngle0ToPi(phi1 - phi2); 00406 } 00407 00410 inline double deltaEta(double eta1, double eta2) { 00411 return fabs(eta1 - eta2); 00412 } 00413 00416 inline double deltaR(double rap1, double phi1, double rap2, double phi2) { 00417 const double dphi = deltaPhi(phi1, phi2); 00418 return sqrt( sqr(rap1-rap2) + sqr(dphi) ); 00419 } 00420 00422 inline double rapidity(double E, double pz) { 00423 if (isZero(E - pz)) { 00424 throw std::runtime_error("Divergent positive rapidity"); 00425 return MAXDOUBLE; 00426 } 00427 if (isZero(E + pz)) { 00428 throw std::runtime_error("Divergent negative rapidity"); 00429 return -MAXDOUBLE; 00430 } 00431 return 0.5*log((E+pz)/(E-pz)); 00432 } 00433 00435 00436 00437 } 00438 00439 00440 #endif