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 <framework/geometry/VectorUtil.h>
#include <Math/Vector2D.h>
#include <Math/Vector3D.h>
#include <Math/VectorUtil.h>
 
#include <cmath>
 
namespace Belle2 {
  namespace CDC {
    class WireLine  {
 
    public:
      WireLine()
      {}

      WireLine(const ROOT::Math::XYZVector& forward, const ROOT::Math::XYZVector& backward, double sagCoeff);

      WireLine movedBy(const ROOT::Math::XYZVector& offset) const;

      WireLine movedBy(const ROOT::Math::XYVector& offset) const;

      ROOT::Math::XYZVector nominalPos3DAtZ(const double z) const
      {
        const ROOT::Math::XYVector& tmp = nominalPos2DAtZ(z);
        return ROOT::Math::XYZVector(tmp.X(), tmp.Y(), z);
      }

      ROOT::Math::XYVector nominalPos2DAtZ(const double z) const
      { return refPos2D() + nominalMovePerZ() * z; }

      ROOT::Math::XYZVector sagPos3DAtZ(const double z) const
      {
        const ROOT::Math::XYVector& tmp = sagPos2DAtZ(z);
        return ROOT::Math::XYZVector(tmp.X(), tmp.Y(), z);
      }

      ROOT::Math::XYVector sagPos2DAtZ(const double z) const
      {
        ROOT::Math::XYVector 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 ROOT::Math::XYVector& nominalMovePerZ() const
      { return m_nominalMovePerZ; }

      ROOT::Math::XYVector sagMovePerZ(const double z) const
      {
        ROOT::Math::XYVector 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 ROOT::Math::XYZVector& pos3D) const
      {
        return ROOT::Math::VectorUtil::Perp((pos3D - refPos3D()), ROOT::Math::XYZVector(nominalMovePerZ().X(), nominalMovePerZ().Y(), 1));
      }

      double sagDistance(const ROOT::Math::XYZVector& pos3D) const
      {
        const ROOT::Math::XYZVector& wirePos3D = sagPos3DAtZ(pos3D.z());
        const ROOT::Math::XYVector& tmp = sagMovePerZ(pos3D.z());
        const ROOT::Math::XYZVector movePerZ(tmp.X(), tmp.Y(), 1);
        return ROOT::Math::VectorUtil::Perp((pos3D - wirePos3D), movePerZ);
      }

      ROOT::Math::XYZVector nominalClosest3D(const ROOT::Math::XYZVector& point) const
      {
        return refPos3D() - ROOT::Math::VectorUtil::ProjVector((point - refPos3D()), ROOT::Math::XYZVector(nominalMovePerZ().X(),
                                                               nominalMovePerZ().Y(), 1));
      }

      ROOT::Math::XYZVector sagClosest3D(const ROOT::Math::XYZVector& point) const
      {
        const ROOT::Math::XYZVector& wirePos3D = sagPos3DAtZ(point.z());
        const ROOT::Math::XYVector& tmp = sagMovePerZ(point.z());
        const ROOT::Math::XYZVector movePerZ(tmp.X(), tmp.Y(), 1);
        return  wirePos3D - ROOT::Math::VectorUtil::ProjVector((point - wirePos3D), movePerZ);
      }

      ROOT::Math::XYZVector forward3D() const
      { return nominalPos3DAtZ(forwardZ()); }

      ROOT::Math::XYVector forward2D() const
      { return nominalPos2DAtZ(forwardZ()); }

      ROOT::Math::XYZVector backward3D() const
      { return nominalPos3DAtZ(backwardZ()); }

      ROOT::Math::XYVector backward2D() const
      { return nominalPos2DAtZ(backwardZ()); }

      ROOT::Math::XYZVector wireVector() const
      {
        const ROOT::Math::XYVector& tmp = nominalMovePerZ() * deltaZ();
        return ROOT::Math::XYZVector(tmp.X(), tmp.Y(), 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().R(); }

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

      double forwardPhiToRef() const
      { return VectorUtil::Angle(forward2D(), refPos2D()); }

      double backwardPhiToRef() const
      { return VectorUtil::Angle(backward2D(), 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 VectorUtil::Angle(backward2D(), forward2D()); }

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

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

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

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

      double nominalPerigeeZ() const
      { return -refPos2D().Dot(nominalMovePerZ()) / nominalMovePerZ().Mag2(); }

      ROOT::Math::XYZVector nominalPerigee3D() const
      { return nominalPos3DAtZ(nominalPerigeeZ()); }

      ROOT::Math::XYVector nominalPerigee2D() const
      { return VectorUtil::orthogonalVector(refPos2D(), 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.Perp2(); }

      const ROOT::Math::XYVector& refPos2D() const
      { return m_refPos2D; }

      const ROOT::Math::XYZVector& refPos3D() const
      { return m_refPos3D; }

      double sagCoeff() const
      { return m_sagCoeff; }
 
    private:
      ROOT::Math::XYZVector m_refPos3D;

      ROOT::Math::XYVector m_refPos2D;

      ROOT::Math::XYVector m_nominalMovePerZ;

      double m_forwardZ = 0.0;

      double m_backwardZ = 0.0;

      double m_sagCoeff = 0.0;
 
    };
  }
}