Math/interface/deltaPhi.h

0001 #ifndef DataFormats_Math_deltaPhi_h
0002 #define DataFormats_Math_deltaPhi_h
0003
0004 /* function to compute deltaPhi
0005  *
0006  * Ported from original code in RecoJets
0007  * by Fedor Ratnikov, FNAL
0008  * stabilize range reduction
0009  */
0010
0011 #include <cmath>
0012 #include "DataFormats/Math/interface/angle_units.h"
0013
0014 namespace reco {
0015
0016   // reduce to [-pi,pi]
0017   template <typename T>
0018   constexpr T reducePhiRange(T x) {
0019     constexpr T o2pi = 1. / (2. * M_PI);
0020     if (std::abs(x) <= T(M_PI))
0021       return x;
0022     T n = std::round(x * o2pi);
0023     return x - n * T(2. * M_PI);
0024   }
0025
0026   constexpr double deltaPhi(double phi1, double phi2) { return reducePhiRange(phi1 - phi2); }
0027
0028   constexpr double deltaPhi(float phi1, double phi2) { return deltaPhi(static_cast<double>(phi1), phi2); }
0029
0030   constexpr double deltaPhi(double phi1, float phi2) { return deltaPhi(phi1, static_cast<double>(phi2)); }
0031
0032   constexpr float deltaPhi(float phi1, float phi2) { return reducePhiRange(phi1 - phi2); }
0033
0034   template <typename T1, typename T2>
0035   constexpr auto deltaPhi(T1 const& t1, T2 const& t2) -> decltype(deltaPhi(t1.phi(), t2.phi())) {
0036     return deltaPhi(t1.phi(), t2.phi());
0037   }
0038
0039   template <typename T>
0040   constexpr T deltaPhi(T phi1, T phi2) {
0041     return reducePhiRange(phi1 - phi2);
0042   }
0043 }  // namespace reco
0044
0045 // lovely!  VI
0046 using reco::deltaPhi;
0047
0048 template <typename T1, typename T2 = T1>
0049 struct DeltaPhi {
0050   constexpr auto operator()(const T1& t1, const T2& t2) -> decltype(reco::deltaPhi(t1, t2)) const {
0051     return reco::deltaPhi(t1, t2);
0052   }
0053 };
0054
0055 namespace angle0to2pi {
0056
0057   using angle_units::operators::operator""_pi;
0058
0059   // make0To2pi constrains an angle to be >= 0 and < 2pi.
0060   // This function is a faster version of reco::reducePhiRange.
0061   // In timing tests, it is almost always faster than reco::reduceRange.
0062   // It also protects against floating-point value drift over repeated calculations.
0063   // This implementation uses multiplication instead of division and avoids
0064   // calling fmod to improve performance.
0065
0066   template <class valType>
0067   inline constexpr valType make0To2pi(valType angle) {
0068     constexpr valType twoPi = 2. * M_PI;
0069     constexpr valType oneOverTwoPi = 1. / (2. * M_PI);
0070     constexpr valType epsilon = 1.e-13;
0071
0072     if ((std::abs(angle) <= epsilon) || (std::abs(twoPi - std::abs(angle)) <= epsilon))
0073       return (0.);
0074     if (std::abs(angle) > twoPi) {
0075       valType nFac = trunc(angle * oneOverTwoPi);
0076       angle -= (nFac * twoPi);
0077       if (std::abs(angle) <= epsilon)
0078         return (0.);
0079     }
0080     if (angle < 0.)
0081       angle += twoPi;
0082     return (angle);
0083   }
0084 }  // namespace angle0to2pi
0085
0086 #endif