WireLine.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 <tracking/trackFindingCDC/geometry/Vector3D.h>
#include <tracking/trackFindingCDC/geometry/Vector2D.h>
 
#include <cmath>
 
namespace Belle2 {
  namespace TrackFindingCDC {
    class WireLine  {
 
    public:
      WireLine()
      {}

      WireLine(const Vector3D& forward, const Vector3D& backward, double sagCoeff);

      WireLine movedBy(const Vector3D& offset) const;

      WireLine movedBy(const Vector2D& offset) const;

      Vector3D nominalPos3DAtZ(const double z) const
      { return Vector3D(nominalPos2DAtZ(z), z); }

      Vector2D nominalPos2DAtZ(const double z) const
      { return refPos2D() + nominalMovePerZ() * z; }

      Vector3D sagPos3DAtZ(const double z) const
      { return Vector3D(sagPos2DAtZ(z), z); }

      Vector2D sagPos2DAtZ(const double z) const
      {
        Vector2D pos2D = nominalPos2DAtZ(z);
        double xzFactorSquared = 1 + nominalMovePerZ().x() * nominalMovePerZ().x();
        double sagY = - sagCoeff() * xzFactorSquared * (forwardZ() - z) * (z - backwardZ());
        pos2D.setY(pos2D.y() + sagY);
        return pos2D;
      }

      const Vector2D& nominalMovePerZ() const
      { return m_nominalMovePerZ; }

      Vector2D sagMovePerZ(const double z) const
      {
        Vector2D movePerZ = nominalMovePerZ();
        double xzFactorSquared = 1 + nominalMovePerZ().x() * nominalMovePerZ().x();
        double sagDYDZ = - sagCoeff() * xzFactorSquared * (forwardZ() + backwardZ() - 2 * z);
        movePerZ.setY(movePerZ.y() + sagDYDZ);
        return movePerZ;
      }

      double nominalDistance(const Vector3D& pos3D) const
      { return (pos3D - refPos3D()).orthogonalComp(Vector3D(nominalMovePerZ(), 1)); }

      double sagDistance(const Vector3D& pos3D) const
      {
        Vector3D wirePos3D = sagPos3DAtZ(pos3D.z());
        Vector3D movePerZ(sagMovePerZ(pos3D.z()), 1);
        return (pos3D - wirePos3D).orthogonalComp(movePerZ);
      }

      Vector3D nominalClosest3D(const Vector3D& point) const
      { return  refPos3D() - (point - refPos3D()).parallelVector(Vector3D(nominalMovePerZ(), 1)); }

      Vector3D sagClosest3D(const Vector3D& point) const
      {
        Vector3D wirePos3D = sagPos3DAtZ(point.z());
        Vector3D movePerZ(sagMovePerZ(point.z()), 1);
        return  wirePos3D - (point - wirePos3D).parallelVector(movePerZ);
      }

      Vector3D forward3D() const
      { return nominalPos3DAtZ(forwardZ()); }

      Vector2D forward2D() const
      { return nominalPos2DAtZ(forwardZ()); }

      Vector3D backward3D() const
      { return nominalPos3DAtZ(backwardZ()); }

      Vector2D backward2D() const
      { return nominalPos2DAtZ(backwardZ()); }

      Vector3D wireVector() const
      { return Vector3D(nominalMovePerZ() * deltaZ(), deltaZ()); }

      double forwardZ() const
      { return m_forwardZ; }

      double backwardZ() const
      { return m_backwardZ; }

      double deltaZ() const
      { return forwardZ() - backwardZ(); }

      double outOfZBoundsFactor(double z) const
      { return std::fmax(backwardZ() - z, z - forwardZ()) / deltaZ(); }

      double forwardPhi() const
      { return forward2D().phi(); }

      double backwardPhi() const
      { return backward2D().phi(); }

      double forwardCylindricalR() const
      { return forward2D().cylindricalR(); }

      double backwardCylindricalR() const
      { return backward2D().cylindricalR(); }

      double forwardPhiToRef() const
      { return forward2D().angleWith(refPos2D()); }

      double backwardPhiToRef() const
      { return backward2D().angleWith(refPos2D()); }

      /* backwardToForwardAngle means how far the backward position has to be rotated in the xy projection
         in the mathematical positive sense that it seems to be coaligned with the forward position. */
      double backwardToForwardAngle() const
      { return backward2D().angleWith(forward2D()) ; }

      double tanLambda() const
      { return 1 / nominalMovePerZ().norm(); }

      double lambda() const
      { return std::atan(tanLambda()); }

      double tanTheta() const
      { return std::atan(nominalMovePerZ().norm()); }

      double theta() const
      { return std::atan(nominalMovePerZ().norm()); }

      double nominalPerigeeZ() const
      { return -refPos2D().dot(nominalMovePerZ()) / nominalMovePerZ().normSquared(); }

      Vector3D nominalPerigee3D() const
      { return nominalPos3DAtZ(nominalPerigeeZ()); }

      Vector2D nominalPerigee2D() const
      { return refPos2D().orthogonalVector(nominalMovePerZ()); }

      double refX() const
      { return m_refPos3D.x(); }

      double refY() const
      { return m_refPos3D.y(); }

      double refZ() const
      { return m_refPos3D.z(); }

      double refCylindricalRSquared() const
      { return m_refPos3D.cylindricalRSquared(); }

      const Vector2D& refPos2D() const
      { return m_refPos3D.xy(); }

      const Vector3D& refPos3D() const
      { return m_refPos3D; }

      double sagCoeff() const
      { return m_sagCoeff; }
 
    private:
      Vector3D m_refPos3D;

      Vector2D m_nominalMovePerZ;

      double m_forwardZ = 0.0;

      double m_backwardZ = 0.0;

      double m_sagCoeff = 0.0;
 
    };
  }
}