HelixSwimmer.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#pragma once
 
#include <Math/Vector3D.h>
#include <Math/Point3D.h>
#include <Math/Transform3D.h>
#include <cmath>
#include <algorithm>
 
namespace Belle2 {
  namespace TOP {
 
    class HelixSwimmer {
 
    public:
 
      HelixSwimmer()
      {}
 
      void set(const ROOT::Math::XYZPoint& position, const ROOT::Math::XYZVector& momentum, double charge, double Bz);
 
      void setTransformation(const ROOT::Math::Transform3D& transform) {m_transformInv = transform.Inverse();}
 
      void moveReferencePosition(double length);
 
      ROOT::Math::XYZPoint getPosition(double length) const;
 
      ROOT::Math::XYZVector getDirection(double length) const;
 
      double getDistanceToPlane(const ROOT::Math::XYZPoint& point, const ROOT::Math::XYZVector& normal) const;
 
    private:
 
      double shortestDistance(double cosAlpha, double phi) const;
 
      double withinSingleTurn(double t) const;
 
      double m_R = 0; 
      double m_xc = 0; 
      double m_yc = 0; 
      double m_omega = 0; 
      double m_kz = 0; 
      double m_phi0 = 0; 
      double m_z0 = 0; 
      double m_T0 = 0; 
      double& x0 = m_xc; 
      double& y0 = m_yc; 
      double& z0 = m_z0; 
      double& kx = m_R;  
      double& ky = m_phi0; 
      double& kz = m_kz; 
      ROOT::Math::Transform3D m_transformInv; 
    };
 
 
    inline double HelixSwimmer::shortestDistance(double cosAlpha, double phi) const
    {
      double alpha = acos(cosAlpha);
      double t1 = withinSingleTurn((phi - m_phi0 - alpha) / m_omega);
      double t2 = withinSingleTurn((phi - m_phi0 + alpha) / m_omega);
      double ta = std::min(t1, t2);
      double tb = std::max(t1 - m_T0, t2 - m_T0);
      double t = (std::abs(ta) < std::abs(tb)) ? ta : tb;
      return t;
    }
 
    inline double HelixSwimmer::withinSingleTurn(double t) const
    {
      t = fmod(t, m_T0);
      if (t < 0) t += m_T0;
      return t;
    }
 
  } // namespace TOP
} // namespace Belle2