File indexing completed on 2023-03-17 10:50:35
0001 #ifndef DataFormats_Math_deltaPhi_h
0002 #define DataFormats_Math_deltaPhi_h
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include <cmath>
0012 #include "DataFormats/Math/interface/angle_units.h"
0013
0014 namespace reco {
0015
0016
0017 template <typename T>
0018 constexpr T reduceRange(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 reduceRange(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 reduceRange(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 reduceRange(phi1 - phi2);
0042 }
0043 }
0044
0045
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
0060
0061
0062
0063
0064
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 }
0085
0086 #endif