SinEqLine.h

/**************************************************************************
 * BASF2 (Belle Analysis Framework 2)                                     *
 * Copyright(C) 2014 - Belle II Collaboration                             *
 *                                                                        *
 * Author: The Belle II Collaboration                                     *
 * Contributors: Oliver Frost                                             *
 *                                                                        *
 * This software is provided "as is" without any warranty.                *
 **************************************************************************/
#pragma once
 
#include <tracking/trackFindingCDC/numerics/Modulo.h>
 
#include <tracking/trackFindingCDC/geometry/Line2D.h>
#include <tracking/trackFindingCDC/geometry/Vector2D.h>
 
#include <tracking/trackFindingCDC/numerics/EIncDec.h>
 
#include <cmath>
 
namespace Belle2 {
  namespace TrackFindingCDC {
 
    class SinEqLine  {
 
    public:
      SinEqLine() :
        m_slope(0.0),
        m_intercept(0.0)
      {}
 
 
      explicit SinEqLine(const Line2D& line2D) :
        m_slope(line2D.slope()),
        m_intercept(line2D.intercept())
      {}
 
      SinEqLine(const double slope, const double intercept) :
        m_slope(slope),
        m_intercept(intercept)
      {}
 
 
      double map(const double x) const
      { return sin(x) - getSlope() * x - getIntercept(); }
 
      double gradient(const double x) const
      { return cos(x) - getSlope(); }
 
      int getIHalfPeriod(const double x) const
      { return floor(x / M_PI); }
 
      /* Computes the positive solution that has the smallest value of x.
      The additional parameter serves as a criterion to abbort the search if the solutions is further away than the specified half period.
      */
      double computeSmallestPositiveRoot(int maxIHalfPeriod = 5) const;
 
 
      double computeRootLargerThanExtemumInHalfPeriod(int iHalfPeriod) const;
 
      double computeRootForLargeSlope() const;
 
      double computeRootInInterval(double lowerX, double upperX) const;
 
    private:
      double newtonX(const Vector2D& pos) const;
 
      static double secantX(const Vector2D& lower, const Vector2D& upper);
 
      static double middleX(const Vector2D& lower, const Vector2D& upper);
 
      static bool updateBounds(Vector2D& lower, Vector2D& upper, const Vector2D& next);
 
      static bool isBetween(const Vector2D& lower, const Vector2D& next, const Vector2D& upper)
      { return lower.x() < next.x() and next.x() < upper.x(); }
 
      static bool isConverged(const Vector2D& lower, const Vector2D& upper)
      {
        return fabs(lower.y()) < 10e-7 or fabs(upper.y()) < 10e-7;
      }
 
      static double getConvergedBound(const Vector2D& lower, const Vector2D& upper)
      {
        if (not std::isfinite(lower.y()) or not std::isfinite(upper.y())) {
          return NAN;
        }
 
        if (fabs(lower.y()) <= fabs(upper.y())) {
          return lower.x();
        }
 
        if (fabs(lower.y()) > fabs(upper.y())) {
          return upper.x();
        }
 
        return NAN;
      }
 
    public:
      static bool changesSign(const Vector2D& lower, const Vector2D& upper)
      { return (lower.y() > 0 and upper.y() < 0) or (lower.y() < 0 and upper.y() > 0); }
 
      static EIncDec getEIncDec(const Vector2D& lower, const Vector2D& upper)
      {
        if (lower.y() < upper.y()) {
          return EIncDec::c_Increasing;
        } else if (lower.y() > upper.y()) {
          return EIncDec::c_Decreasing;
        } else if (lower.y() == upper.y()) {
          return EIncDec::c_Constant;
        } else {
          return EIncDec::c_Invalid;
        }
      }
 
    public:
      double computeExtremumXInHalfPeriod(int iHalfPeriod) const;
 
      static int getIPeriodFromIHalfPeriod(int iHalfPeriod)
      { return isEven(iHalfPeriod) ? iHalfPeriod / 2 : (iHalfPeriod - 1) / 2; }
 
    public:
      bool hasLargeSlope() const
      { return fabs(getSlope()) >= 1; }
 
      double getSlope() const
      { return m_slope; }
 
      double getIntercept() const
      { return m_intercept; }
 
    private:
      double m_slope;
 
      double m_intercept;
 
 
    };
 
  }
}