phihelper.h

Go to the documentation of this file.

/*
  Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
*/
#ifndef CXXUTILS_PHIHELPER_H
#define CXXUTILS_PHIHELPER_H
 
#include <climits>
#include <numbers>
#include <concepts>
#include <cmath>
 
namespace CxxUtils {

  template <std::floating_point T>
  inline constexpr T wrapToPi(T phi)
  {
    constexpr T TWOPI = 2*std::numbers::pi_v<T>;
    constexpr T INV2PI = std::numbers::inv_pi_v<T>/2;
    T x = phi * INV2PI;
 
    // Round x to the nearest integer.
    // https://stackoverflow.com/questions/17035464/a-fast-method-to-round-a-double-to-a-32-bit-int-explained
    static_assert (std::numeric_limits<T>::digits/2 <= sizeof (long int) * CHAR_BIT);
    // Break up the exponent into two pieces to avoid overflowing a long int
    // for an ARM long double with 113 mantissa bits.
    constexpr T TOINT = 0x1.8p0 *
      (1ull<<(std::numeric_limits<T>::digits/2)) *
      (1ull<<(std::numeric_limits<T>::digits-1-std::numeric_limits<T>::digits/2));
    T ix = (x + TOINT) - TOINT;
 
    // Above gives banker's rounding; that is, halves round to even integers.
    // However, to get the cases of phi=Npi right, we want halves to round
    // towards zero.  Fix up the rounding in that case.  Is there a
    // better way of doing this?
    T diff = ix - x;
    if (std::abs(diff) == 0.5) {
      if (ix > 0) --ix;
      else if (ix < 0) ++ix;
    }
 
    return phi - TWOPI*ix;
  }

  template <std::floating_point T>
  inline constexpr T deltaPhi(T phiA, T phiB)
  {
    return wrapToPi(phiA - phiB);
  }

  template <std::floating_point T>
  inline constexpr T phiMean(T phiA, T phiB)
  {
    const T diff = wrapToPi(phiA - phiB);
    return wrapToPi(phiB + 0.5 * diff);
  }

  template <std::floating_point T>
  inline constexpr T phiBisect(T phiA, T phiB)
  {
    T phi = 0.5 * (phiA + phiB);
    if (phiA > phiB) phi += std::numbers::pi;
    return wrapToPi(phi);
  }
 
} // namespace CxxUtils
 
#endif