CDCTrajectory2D.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 <cdc/topology/ISuperLayer.h>
 
#include <tracking/trackingUtilities/geometry/UncertainPerigeeCircle.h>
#include <tracking/trackingUtilities/geometry/PerigeeCircle.h>
#include <tracking/trackingUtilities/geometry/PerigeeParameters.h>
 
#include <tracking/trackingUtilities/numerics/EForwardBackward.h>
#include <tracking/trackingUtilities/numerics/ERightLeft.h>
#include <tracking/trackingUtilities/numerics/ESign.h>
 
#include <Math/Vector3D.h>
#include <Math/Vector2D.h>
 
#include <TMath.h>
#include <cmath>
#include <cstddef>
#include <iosfwd>
 
namespace Belle2 {
  namespace CDC {
    class WireLine;
  }
  namespace TrackingUtilities {

    class CDCTrajectory2D  {
 
    public:
      CDCTrajectory2D();

      explicit CDCTrajectory2D(const UncertainPerigeeCircle& perigeeCircle);

      CDCTrajectory2D(const ROOT::Math::XYVector& localOrigin,
                      const UncertainPerigeeCircle& localPerigeeCircle,
                      double flightTime = NAN);

      CDCTrajectory2D(const ROOT::Math::XYVector& pos2D,
                      double time,
                      const ROOT::Math::XYVector& mom2D,
                      double charge,
                      double bZ);

      CDCTrajectory2D(const ROOT::Math::XYVector& pos2D, double time, const ROOT::Math::XYVector& mom2D, double charge);
 
    public:
      bool isFitted() const;

      void clear();
 
    public:
      void reverse();

      CDCTrajectory2D reversed() const;
 
    public:
      std::array<double, 2> reconstructBothZ(const CDC::WireLine& wireLine, double distance = 0.0, double z = 0) const;

      double reconstructZ(const CDC::WireLine& wireLine, double distance = 0.0, double z = 0) const;

      std::array<ROOT::Math::XYZVector, 2> reconstructBoth3D(const CDC::WireLine& wireLine, double distance = 0.0, double z = 0) const;

      ROOT::Math::XYZVector reconstruct3D(const CDC::WireLine& wireLine, double distance = 0.0, double z = 0) const;

      ROOT::Math::XYVector getClosest(const ROOT::Math::XYVector& point) const;
 
    private:
      CDC::ISuperLayer getISuperLayerAfter(CDC::ISuperLayer iSuperLayer, bool movingOutward) const;

      CDC::ISuperLayer getISuperLayerAfterStart(bool movingOutward) const;

      CDC::ISuperLayer getISuperLayerAfterStart(EForwardBackward forwardBackwardInfo) const;

      CDC::ISuperLayer getAxialISuperLayerAfterStart(EForwardBackward forwardBackwardInfo) const;
 
    public:
      CDC::ISuperLayer getNextISuperLayer() const;

      CDC::ISuperLayer getPreviousISuperLayer() const;

      CDC::ISuperLayer getNextAxialISuperLayer() const;

      CDC::ISuperLayer getPreviousAxialISuperLayer() const;

      CDC::ISuperLayer getMaximalISuperLayer() const;

      CDC::ISuperLayer getStartISuperLayer() const;

      CDC::ISuperLayer getMinimalISuperLayer() const;
 
    public:
      template <class AHits>
      EForwardBackward isForwardOrBackwardTo(const AHits& hits) const
      {
        return static_cast<EForwardBackward>(sign(getTotalArcLength2D(hits)));
      }

      template <class AFromHits, class AToHits>
      double getArcLength2DGap(const AFromHits& fromHits, const AToHits& toHits) const
      {
        const ROOT::Math::XYVector& fromRecoPos2D = fromHits.back().getRecoPos2D();
        const ROOT::Math::XYVector& toRecoPos2D = toHits.front().getRecoPos2D();
        return calcArcLength2DBetween(fromRecoPos2D, toRecoPos2D);
      }

      template <class AFromHits, class AToHits>
      double getArcLength2DFrontOffset(const AFromHits& fromHits, const AToHits& toHits) const
      {
        const ROOT::Math::XYVector& fromRecoPos2D = fromHits.front().getRecoPos2D();
        const ROOT::Math::XYVector& toRecoPos2D = toHits.front().getRecoPos2D();
        return calcArcLength2DBetween(fromRecoPos2D, toRecoPos2D);
      }

      template <class AFromHits, class AToHits>
      double getArcLength2DBackOffset(const AFromHits& fromHits, const AToHits& toHits) const
      {
        const ROOT::Math::XYVector& fromRecoPos2D = fromHits.back().getRecoPos2D();
        const ROOT::Math::XYVector& toRecoPos2D = toHits.back().getRecoPos2D();
        return calcArcLength2DBetween(fromRecoPos2D, toRecoPos2D);
      }

      template <class AHits>
      double getTotalArcLength2D(const AHits& hits) const
      {
        ROOT::Math::XYVector frontRecoPos2D = hits.front().getRecoPos2D();
        ROOT::Math::XYVector backRecoPos2D = hits.back().getRecoPos2D();
        return calcArcLength2DBetween(frontRecoPos2D, backRecoPos2D);
      }

      double calcArcLength2D(const ROOT::Math::XYVector& point) const
      {
        return getLocalCircle()->arcLengthBetween(ROOT::Math::XYVector(0.0, 0.0), point - getLocalOrigin());
      }

      double calcArcLength2DBetween(const ROOT::Math::XYVector& fromPoint, const ROOT::Math::XYVector& toPoint) const
      {
        return getLocalCircle()->arcLengthBetween(fromPoint - getLocalOrigin(),
                                                  toPoint - getLocalOrigin());
      }

      double getArcLength2DPeriod() const
      {
        return getLocalCircle()->arcLengthPeriod();
      }
 
    public:
      void setPosMom2D(const ROOT::Math::XYVector& pos2D, const ROOT::Math::XYVector& mom2D, double charge);

      ESign getChargeSign() const;

      double getAbsMom2D(double bZ) const;

      double getAbsMom2D() const;

      ROOT::Math::XYVector getMom2DAtSupport(const double bZ) const
      {
        return getFlightDirection2DAtSupport() * getAbsMom2D(bZ);
      }

      ROOT::Math::XYVector getMom2DAtSupport() const
      {
        return getFlightDirection2DAtSupport() * getAbsMom2D();
      }

      ROOT::Math::XYVector getFlightDirection2D(const ROOT::Math::XYVector& point) const
      {
        return getLocalCircle()->tangential(point - getLocalOrigin());
      }

      ROOT::Math::XYVector getFlightDirection2DAtSupport() const
      {
        return getLocalCircle()->tangential();
      }

      bool isMovingOutward() const
      {
        return getFlightDirection2DAtSupport().Dot(getSupport()) > 0;
      }

      ROOT::Math::XYVector getPos2DAtArcLength2D(double arcLength2D)
      {
        return getLocalOrigin() + getLocalCircle()->atArcLength(arcLength2D);
      }

      ROOT::Math::XYVector getSupport() const
      {
        return getLocalCircle()->perigee() + getLocalOrigin();
      }

      ROOT::Math::XYVector getGlobalPerigee() const
      {
        return getLocalCircle()->closest(-m_localOrigin) + m_localOrigin;
      }

      ROOT::Math::XYVector getGlobalCenter() const
      {
        return getLocalCircle()->center() + m_localOrigin;
      }

      ROOT::Math::XYVector getOuterExit(double factor = 1) const;

      ROOT::Math::XYVector getInnerExit() const;

      ROOT::Math::XYVector getExit() const;

      bool isCurler(double factor = 1) const;

      bool isOriginer(double factor = 1) const;

      double getMaximalCylindricalR() const
      {
        return std::fabs(getGlobalImpact() + 2 * getLocalCircle()->radius());
      }

      double getMinimalCylindricalR() const
      {
        return std::fabs(getGlobalImpact());
      }

      double getGlobalImpact() const
      {
        return getLocalCircle()->distance(-m_localOrigin);
      }

      double getDist2D(const ROOT::Math::XYVector& point) const
      {
        return getLocalCircle()->distance(point - getLocalOrigin());
      }

      ERightLeft isRightOrLeft(const ROOT::Math::XYVector& point) const
      {
        return getLocalCircle()->isRightOrLeft(point - getLocalOrigin());
      }
 
    public:
      double getCurvature() const
      {
        return getLocalCircle()->curvature();
      }

      double getLocalCovariance(EPerigeeParameter iRow, EPerigeeParameter iCol) const
      {
        return getLocalCircle().covariance(iRow, iCol);
      }

      double getLocalVariance(EPerigeeParameter i) const
      {
        return getLocalCircle().variance(i);
      }

      PerigeeCircle getGlobalCircle() const
      {
        PerigeeCircle result = getLocalCircle();
        result.passiveMoveBy(-getLocalOrigin());
        return result;
      }

      void setGlobalCircle(const UncertainPerigeeCircle& perigeeCircle)
      {
        m_localPerigeeCircle = perigeeCircle;
        m_localPerigeeCircle.passiveMoveBy(getLocalOrigin());
      }

      const UncertainPerigeeCircle& getLocalCircle() const
      {
        return m_localPerigeeCircle;
      }

      double getPValue() const
      {
        return TMath::Prob(getChi2(), getNDF());
      }

      double getChi2() const
      {
        return getLocalCircle().chi2();
      }

      void setChi2(const double chi2)
      {
        return m_localPerigeeCircle.setChi2(chi2);
      }

      size_t getNDF() const
      {
        return getLocalCircle().ndf();
      }

      void setNDF(std::size_t ndf)
      {
        return m_localPerigeeCircle.setNDF(ndf);
      }

      void setLocalCircle(const UncertainPerigeeCircle& localPerigeeCircle)
      {
        m_localPerigeeCircle = localPerigeeCircle;
      }

      const ROOT::Math::XYVector& getLocalOrigin() const
      {
        return m_localOrigin;
      }

      double setLocalOrigin(const ROOT::Math::XYVector& localOrigin);

      double getFlightTime() const
      {
        return m_flightTime;
      }

      void setFlightTime(double flightTime)
      {
        m_flightTime = flightTime;
      }
 
    private:
      ROOT::Math::XYVector m_localOrigin;

      UncertainPerigeeCircle m_localPerigeeCircle;

      double m_flightTime = NAN;
    };

    std::ostream& operator<<(std::ostream& output, const CDCTrajectory2D& trajectory2D);
  }
}