Helix.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/datastore/RelationsObject.h>
 
#include <Math/Vector3D.h>
#include <TMatrixDfwd.h>
 
#include <cmath>
 
namespace Belle2 {
 
  namespace HelixParameterIndex {
    enum HelixParameterIndex {
      iD0 = 0,
      iPhi0 = 1,
      iOmega = 2,
      iZ0 = 3,
      iTanLambda = 4,
      iArcLength2D = 5,
    };
  }

  class Helix : public RelationsObject {
  public:

    Helix();

    Helix(const ROOT::Math::XYZVector& position,
          const ROOT::Math::XYZVector& momentum,
          const short int charge,
          const double bZ);

    Helix(const double& d0,
          const double& phi0,
          const double& omega,
          const double& z0,
          const double& tanLambda);

    friend std::ostream& operator<<(std::ostream& output, const Helix& helix);
 
 
    //---------------------------------------------------------------------------------------------------------------------------
    //---------------------------------------------------------------------------------------------------------------------------
  public:
    double getPerigeeX() const;

    double getPerigeeY() const;

    double getPerigeeZ() const;

    ROOT::Math::XYZVector getPerigee() const;

    double getMomentumX(const double bZ) const;

    double getMomentumY(const double bZ) const;

    double getMomentumZ(const double bZ) const;

    ROOT::Math::XYZVector getMomentum(const double bZ) const;

    ROOT::Math::XYZVector getDirection() const;

    double getTransverseMomentum(const double bZ) const;

    double getKappa(const double bZ) const;

    static double getAlpha(const double bZ);

    short getChargeSign() const;
 
 
    //---------------------------------------------------------------------------------------------------------------------------
    //---------------------------------------------------------------------------------------------------------------------------
  public:
    double getArcLength2DAtCylindricalR(const double& cylindricalR) const;

    double getArcLength2DAtXY(const double& x, const double& y) const;

    double getArcLength2DAtNormalPlane(const double& x, const double& y,
                                       const double& nX, const double& nY) const;

    ROOT::Math::XYZVector getPositionAtArcLength2D(const double& arcLength2D) const;

    ROOT::Math::XYZVector getTangentialAtArcLength2D(const double& arcLength2D) const;

    ROOT::Math::XYZVector getUnitTangentialAtArcLength2D(const double& arcLength2D) const;

    ROOT::Math::XYZVector getMomentumAtArcLength2D(const double& arcLength2D, const double& bz) const;

    double passiveMoveBy(const ROOT::Math::XYZVector& by)
    { return passiveMoveBy(by.X(), by.Y(), by.Z()); }

    double passiveMoveBy(const double& byX,
                         const double& byY,
                         const double& byZ);
 

    TMatrixD calcPassiveMoveByJacobian(const ROOT::Math::XYZVector& by, const double expandBelowChi = M_PI / 8) const;
 
 
  protected:
  public:
    void calcPassiveMoveByJacobian(const double& byX,
                                   const double& byY,
                                   TMatrixD& jacobian,
                                   const double expandBelowChi = M_PI / 8) const;
 
  public:
    void reverse();

    static double reversePhi(const double& phi);

    double calcArcLength2DFromSecantLength(const double& secantLength2D) const;

    double calcSecantLengthToArcLength2DFactor(const double& secantLength2D) const;
 
  protected:
    static double calcASinXDividedByX(const double& x);

    static double calcATanXDividedByX(const double& x);

    static double calcDerivativeOfATanXDividedByX(const double& x);

    void calcArcLength2DAndDrAtXY(const double& x, const double& y, double& arcLength2D, double& dr) const;

    double calcArcLength2DAtDeltaCylindricalRAndDr(const double& deltaCylindricalR, const double& dr) const;
 
 
    //---------------------------------------------------------------------------------------------------------------------------
    //---------------------------------------------------------------------------------------------------------------------------
  public:
    double getD0() const { return m_d0; }

    double getPhi0() const { return m_phi0; }

    double getCosPhi0() const { return std::cos(getPhi0()); }

    double getSinPhi0() const { return std::sin(getPhi0()); }

    double getOmega() const { return m_omega; }

    double getZ0() const { return m_z0; }

    double getTanLambda() const { return m_tanLambda; }

    double getCotTheta() const { return m_tanLambda; }

  private:
    //---------------------------------------------------------------------------------------------------------------------------
    //--- Functions for internal conversions between cartesian and perigee helix parameters
    //--- This can be placed in a separate header which handles all the conversion stuff
    //---------------------------------------------------------------------------------------------------------------------------

    void setCartesian(const ROOT::Math::XYZVector& position,
                      const ROOT::Math::XYZVector& momentum,
                      const short int charge,
                      const double bZ);

    Double32_t m_d0;

    Double32_t m_phi0;

    Double32_t m_omega;

    Double32_t m_z0;

    Double32_t m_tanLambda;
 
    ClassDef(Helix, 2); 
  };

}