GeometryVector/interface/Phi.h

0001 #ifndef GeometryVector_Geom_Phi_h
0002 #define GeometryVector_Geom_Phi_h
0003
0004 #include "DataFormats/GeometryVector/interface/Pi.h"
0005 #include "DataFormats/Math/interface/deltaPhi.h"
0006 #include "DataFormats/Math/interface/angle_units.h"
0007 #include <cmath>
0008
0009 namespace Geom {
0010
0011   /** \class Phi
0012  *  A class for azimuthal angle represantation and algebra.
0013  *  The use of Phi<T> is tranparant due to the implicit conversion to T
0014  *  Constructs like cos(phi) work as with float or double.
0015  *  The difference with respect to built-in types is that
0016  *  Phi is kept in the range (-pi, pi] by default, and this is consistently
0017  *  implemented in aritmetic operations. In other words, Phi implements
0018  *  "modulo(2 pi)" arithmetics.
0019  *  Phi can be instantiated to implement the range 0 to 2pi.
0020  */
0021
0022   using angle_units::operators::operator""_deg;
0023   using angle_units::operators::convertRadToDeg;
0024
0025   struct MinusPiToPi {};  // Dummy struct to indicate -pi to pi range
0026   struct ZeroTo2pi {};    // Dummy struct to indicate 0 to 2pi range
0027
0028   template <typename T1, typename Range>
0029   class NormalizeWrapper {};
0030
0031   template <typename T1>
0032   class NormalizeWrapper<T1, MinusPiToPi> {
0033   public:
0034     static void normalize(T1& value) {  // Reduce range to -pi to pi
0035       if (value > twoPi() || value < -twoPi()) {
0036         value = std::fmod(value, static_cast<T1>(twoPi()));
0037       }
0038       if (value <= -pi())
0039         value += twoPi();
0040       if (value > pi())
0041         value -= twoPi();
0042     }
0043   };
0044
0045   template <typename T1>
0046   class NormalizeWrapper<T1, ZeroTo2pi> {  // Reduce range to 0 to 2pi
0047   public:
0048     static void normalize(T1& value) { value = angle0to2pi::make0To2pi(value); }
0049   };
0050
0051   template <typename T1, typename Range = MinusPiToPi>
0052   class Phi {
0053   public:
0054     /// Default constructor does not initialise - just as double.
0055     Phi() {}
0056
0057     // Constructor from T1.
0058     // Not "explicit" to enable convenient conversions.
0059     // There may be cases of ambiguities because of multiple possible
0060     // conversions, in which case explicit casts must be used.
0061     // The constructor provides range checking and normalization,
0062     // e.g. the value of Phi(2*pi()+1) is 1
0063     Phi(const T1& val) : theValue(val) { normalize(theValue); }
0064
0065     /// conversion operator makes transparent use possible.
0066     operator T1() const { return theValue; }
0067
0068     /// Template argument conversion
0069     template <typename T2, typename Range1>
0070     operator Phi<T2, Range1>() {
0071       return Phi<T2, Range1>(theValue);
0072     }
0073
0074     /// Explicit access to value in case implicit conversion not OK
0075     T1 value() const { return theValue; }
0076
0077     // so that template classes expecting phi() works! (deltaPhi)
0078     T1 phi() const { return theValue; }
0079
0080     /// Standard arithmetics
0081     Phi& operator+=(const T1& a) {
0082       theValue += a;
0083       normalize(theValue);
0084       return *this;
0085     }
0086     Phi& operator+=(const Phi& a) { return operator+=(a.value()); }
0087
0088     Phi& operator-=(const T1& a) {
0089       theValue -= a;
0090       normalize(theValue);
0091       return *this;
0092     }
0093     Phi& operator-=(const Phi& a) { return operator-=(a.value()); }
0094
0095     Phi& operator*=(const T1& a) {
0096       theValue *= a;
0097       normalize(theValue);
0098       return *this;
0099     }
0100
0101     Phi& operator/=(const T1& a) {
0102       theValue /= a;
0103       normalize(theValue);
0104       return *this;
0105     }
0106
0107     T1 degrees() const { return convertRadToDeg(theValue); }
0108
0109     // nearZero() tells whether the angle is close enough to 0 to be considered 0.
0110     // The default tolerance is 1 degree.
0111     inline bool nearZero(float tolerance = 1.0_deg) const { return (std::abs(theValue) - tolerance <= 0.0); }
0112
0113     // nearEqual() tells whether two angles are close enough to be considered equal.
0114     // The default tolerance is 0.001 radian.
0115     inline bool nearEqual(const Phi<T1, Range>& angle, float tolerance = 0.001) const {
0116       return (std::abs(theValue - angle) - tolerance <= 0.0);
0117     }
0118
0119   private:
0120     T1 theValue;
0121
0122     void normalize(T1& value) { NormalizeWrapper<T1, Range>::normalize(value); }
0123   };
0124
0125   /// - operator
0126   template <typename T1, typename Range>
0127   inline Phi<T1, Range> operator-(const Phi<T1, Range>& a) {
0128     return Phi<T1, Range>(-a.value());
0129   }
0130
0131   /// Addition
0132   template <typename T1, typename Range>
0133   inline Phi<T1, Range> operator+(const Phi<T1, Range>& a, const Phi<T1, Range>& b) {
0134     return Phi<T1, Range>(a) += b;
0135   }
0136   /// Addition with scalar, does not change the precision
0137   template <typename T1, typename Range, typename Scalar>
0138   inline Phi<T1, Range> operator+(const Phi<T1, Range>& a, const Scalar& b) {
0139     return Phi<T1, Range>(a) += b;
0140   }
0141   /// Addition with scalar, does not change the precision
0142   template <typename T1, typename Range, typename Scalar>
0143   inline Phi<T1, Range> operator+(const Scalar& a, const Phi<T1, Range>& b) {
0144     return Phi<T1, Range>(b) += a;
0145   }
0146
0147   /// Subtraction
0148   template <typename T1, typename Range>
0149   inline Phi<T1, Range> operator-(const Phi<T1, Range>& a, const Phi<T1, Range>& b) {
0150     return Phi<T1, Range>(a) -= b;
0151   }
0152   /// Subtraction with scalar, does not change the precision
0153   template <typename T1, typename Range, typename Scalar>
0154   inline Phi<T1, Range> operator-(const Phi<T1, Range>& a, const Scalar& b) {
0155     return Phi<T1, Range>(a) -= b;
0156   }
0157   /// Subtraction with scalar, does not change the precision
0158   template <typename T1, typename Range, typename Scalar>
0159   inline Phi<T1, Range> operator-(const Scalar& a, const Phi<T1, Range>& b) {
0160     return Phi<T1, Range>(a - b.value());
0161   }
0162
0163   /// Multiplication with scalar, does not change the precision
0164   template <typename T1, typename Range, typename Scalar>
0165   inline Phi<T1, Range> operator*(const Phi<T1, Range>& a, const Scalar& b) {
0166     return Phi<T1, Range>(a) *= b;
0167   }
0168   /// Multiplication with scalar
0169   template <typename T1, typename Range>
0170   inline Phi<T1, Range> operator*(double a, const Phi<T1, Range>& b) {
0171     return Phi<T1, Range>(b) *= a;
0172   }
0173
0174   /// Division
0175   template <typename T1, typename Range>
0176   inline T1 operator/(const Phi<T1, Range>& a, const Phi<T1, Range>& b) {
0177     return a.value() / b.value();
0178   }
0179   /// Division by scalar
0180   template <typename T1, typename Range>
0181   inline Phi<T1, Range> operator/(const Phi<T1, Range>& a, double b) {
0182     return Phi<T1, Range>(a) /= b;
0183   }
0184
0185   // For convenience
0186   template <typename T>
0187   using Phi0To2pi = Phi<T, ZeroTo2pi>;
0188 }  // namespace Geom
0189
0190 /*
0191 // this a full mess with the above that is a mess in itself
0192 #include "DataFormats/Math/interface/deltaPhi.h"
0193 namespace reco {
0194   template <class T1,class T2>
0195   inline double deltaPhi(const Geom::Phi<T1> phi1, const Geom::Phi<T2> phi2) {
0196     return deltaPhi(static_cast<double>(phi1.value()), static_cast<double>(phi2.value()));
0197   }
0198
0199   template <class T>
0200   inline double deltaPhi(const Geom::Phi<T> phi1, double phi2) {
0201     return deltaPhi(static_cast<double>(phi1.value()), phi2);
0202   }
0203   template <class T>
0204   inline double deltaPhi(const Geom::Phi<T> phi1, float phi2) {
0205     return deltaPhi(static_cast<double>(phi1.value()), static_cast<double>(phi2));
0206   }
0207   template <class T>
0208   inline double deltaPhi(double phi1, const Geom::Phi<T>  phi2) {
0209     return deltaPhi(phi1, static_cast<double>(phi2.value()) );
0210   }
0211   template <class T>
0212   inline double deltaPhi(float phi1, const Geom::Phi<T>  phi2) {
0213     return deltaPhi(static_cast<double>(phi1),static_cast<double>(phi2.value()) );
0214   }
0215 }
0216 */
0217
0218 #endif