RKTrackRep.h

/* Copyright 2008-2010, Technische Universitaet Muenchen,
   Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch
 
   This file is part of GENFIT.
 
   GENFIT is free software: you can redistribute it and/or modify
   it under the terms of the GNU Lesser General Public License as published
   by the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.
 
   GENFIT is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU Lesser General Public License for more details.
 
   You should have received a copy of the GNU Lesser General Public License
   along with GENFIT.  If not, see <http://www.gnu.org/licenses/>.
*/
 
#ifndef genfit_RKTrackRep_h
#define genfit_RKTrackRep_h
 
#include "AbsTrackRep.h"
#include "StateOnPlane.h"
#include "RKTools.h"
#include "StepLimits.h"
#include "Material.h"
 
#include <algorithm>
 
namespace genfit {
 
struct RKStep {
  MatStep matStep_; // material properties and stepsize
  M1x7 state7_; // 7D state vector
  StepLimits limits_;
 
  RKStep() {
    std::fill(state7_.begin(), state7_.end(), 0);
  }
};
 
 
struct ExtrapStep {
  M7x7 jac7_; // 5D jacobian of transport
  M7x7 noise7_; // 5D noise matrix
 
  ExtrapStep() {
    std::fill(jac7_.begin(), jac7_.end(), 0);
    std::fill(noise7_.begin(), noise7_.end(), 0);
  }
};
 
 
class RKTrackRep : public AbsTrackRep {
    friend class RKTrackRepTests_momMag_Test;
    friend class RKTrackRepTests_calcForwardJacobianAndNoise_Test;
    friend class RKTrackRepTests_getState7_Test;
    friend class RKTrackRepTests_getState5_Test;
 
 public:
 
  RKTrackRep();
  RKTrackRep(int pdgCode, char propDir = 0);
 
  virtual ~RKTrackRep();
 
  virtual AbsTrackRep* clone() const override {return new RKTrackRep(*this);}
 
  virtual double extrapolateToPlane(StateOnPlane& state,
      const SharedPlanePtr& plane,
      bool stopAtBoundary = false,
      bool calcJacobianNoise = false) const override;
 
  using AbsTrackRep::extrapolateToLine;
 
  virtual double extrapolateToLine(StateOnPlane& state,
      const TVector3& linePoint,
      const TVector3& lineDirection,
      bool stopAtBoundary = false,
      bool calcJacobianNoise = false) const override;
 
  virtual double extrapolateToPoint(StateOnPlane& state,
      const TVector3& point,
      bool stopAtBoundary = false,
      bool calcJacobianNoise = false) const override {
    return extrapToPoint(state, point, nullptr, stopAtBoundary, calcJacobianNoise);
  }
 
  virtual double extrapolateToPoint(StateOnPlane& state,
      const TVector3& point,
      const TMatrixDSym& G, // weight matrix (metric)
      bool stopAtBoundary = false,
      bool calcJacobianNoise = false) const override {
    return extrapToPoint(state, point, &G, stopAtBoundary, calcJacobianNoise);
  }
 
  virtual double extrapolateToCylinder(StateOnPlane& state,
      double radius,
      const TVector3& linePoint = TVector3(0.,0.,0.),
      const TVector3& lineDirection = TVector3(0.,0.,1.),
      bool stopAtBoundary = false,
      bool calcJacobianNoise = false) const override;
 
  
  virtual double extrapolateToCone(StateOnPlane& state,
      double radius,
      const TVector3& linePoint = TVector3(0.,0.,0.),
      const TVector3& lineDirection = TVector3(0.,0.,1.),
      bool stopAtBoundary = false,
      bool calcJacobianNoise = false) const override ;
 
  virtual double extrapolateToSphere(StateOnPlane& state,
      double radius,
      const TVector3& point = TVector3(0.,0.,0.),
      bool stopAtBoundary = false,
      bool calcJacobianNoise = false) const override;
 
  virtual double extrapolateBy(StateOnPlane& state,
      double step,
      bool stopAtBoundary = false,
      bool calcJacobianNoise = false) const override;
 
 
  unsigned int getDim() const override {return 5;}
 
  virtual TVector3 getPos(const StateOnPlane& state) const override;
 
  virtual TVector3 getMom(const StateOnPlane& state) const override;
  virtual void getPosMom(const StateOnPlane& state, TVector3& pos, TVector3& mom) const override;
 
  virtual double getMomMag(const StateOnPlane& state) const override;
  virtual double getMomVar(const MeasuredStateOnPlane& state) const override;
 
  virtual TMatrixDSym get6DCov(const MeasuredStateOnPlane& state) const override;
  virtual void getPosMomCov(const MeasuredStateOnPlane& state, TVector3& pos, TVector3& mom, TMatrixDSym& cov) const override;
  virtual double getCharge(const StateOnPlane& state) const override;
  virtual double getQop(const StateOnPlane& state) const override {return state.getState()(0);}
  double getSpu(const StateOnPlane& state) const;
  double getTime(const StateOnPlane& state) const override;
 
  virtual void getForwardJacobianAndNoise(TMatrixD& jacobian, TMatrixDSym& noise, TVectorD& deltaState) const override;
 
  virtual void getBackwardJacobianAndNoise(TMatrixD& jacobian, TMatrixDSym& noise, TVectorD& deltaState) const override;
 
  std::vector<genfit::MatStep> getSteps() const override;
 
  virtual double getRadiationLenght() const override;
 
  virtual void setPosMom(StateOnPlane& state, const TVector3& pos, const TVector3& mom) const override;
  virtual void setPosMom(StateOnPlane& state, const TVectorD& state6) const override;
  virtual void setPosMomErr(MeasuredStateOnPlane& state, const TVector3& pos, const TVector3& mom, const TVector3& posErr, const TVector3& momErr) const override;
  virtual void setPosMomCov(MeasuredStateOnPlane& state, const TVector3& pos, const TVector3& mom, const TMatrixDSym& cov6x6) const override;
  virtual void setPosMomCov(MeasuredStateOnPlane& state, const TVectorD& state6, const TMatrixDSym& cov6x6) const override;
 
  virtual void setChargeSign(StateOnPlane& state, double charge) const override;
  virtual void setQop(StateOnPlane& state, double qop) const override {state.getState()(0) = qop;}
 
  void setSpu(StateOnPlane& state, double spu) const;
  void setTime(StateOnPlane& state, double time) const override;
 
 
  virtual double RKPropagate(M1x7& state7,
                             M7x7* jacobian,
                             M1x3& SA,
                             double S,
                             bool varField = true,
                             bool calcOnlyLastRowOfJ = false) const;
 
  virtual bool isSameType(const AbsTrackRep* other) override;
  virtual bool isSame(const AbsTrackRep* other) override;
 
 private:
 
  void initArrays() const;
 
  virtual double extrapToPoint(StateOnPlane& state,
      const TVector3& point,
      const TMatrixDSym* G = nullptr, // weight matrix (metric)
      bool stopAtBoundary = false,
      bool calcJacobianNoise = false) const;
 
  void getState7(const StateOnPlane& state, M1x7& state7) const;
  void getState5(StateOnPlane& state, const M1x7& state7) const; // state7 must already lie on plane of state!
 
  void calcJ_pM_5x7(M5x7& J_pM, const TVector3& U, const TVector3& V, const M1x3& pTilde, double spu) const;
 
  void transformPM6(const MeasuredStateOnPlane& state,
                    M6x6& out6x6) const;
 
  void calcJ_Mp_7x5(M7x5& J_Mp, const TVector3& U, const TVector3& V, const TVector3& W, const M1x3& A) const;
 
  void calcForwardJacobianAndNoise(const M1x7& startState7, const DetPlane& startPlane,
           const M1x7& destState7, const DetPlane& destPlane) const;
 
  void transformM6P(const M6x6& in6x6,
                    const M1x7& state7,
                    MeasuredStateOnPlane& state) const; // plane and charge must already be set!
 
 
  bool RKutta(const M1x4& SU,
              const DetPlane& plane,
              double charge,
              double mass,
              M1x7& state7,
              M7x7* jacobianT,
              M1x7* J_MMT_unprojected_lastRow,
              double& coveredDistance, // signed
              double& flightTime,
              bool& checkJacProj,
              M7x7& noiseProjection,
              StepLimits& limits,
              bool onlyOneStep = false,
              bool calcOnlyLastRowOfJ = false) const;
 
  double estimateStep(const M1x7& state7,
                      const M1x4& SU,
                      const DetPlane& plane,
                      const double& charge,
                      double& relMomLoss,
                      StepLimits& limits) const;
 
  TVector3 pocaOnLine(const TVector3& linePoint,
                     const TVector3& lineDirection,
                     const TVector3& point) const;
 
 
  double Extrap(const DetPlane& startPlane, // plane where Extrap starts
                const DetPlane& destPlane, // plane where Extrap has to extrapolate to
                double charge,
                double mass,
                bool& isAtBoundary,
                M1x7& state7,
                double& flightTime,
                bool fillExtrapSteps,
                TMatrixDSym* cov = nullptr,
                bool onlyOneStep = false,
                bool stopAtBoundary = false,
                double maxStep = 1.E99) const;
 
  void checkCache(const StateOnPlane& state, const SharedPlanePtr* plane) const;
 
  double momMag(const M1x7& state7) const;
 
 protected:
 
  mutable StateOnPlane lastStartState_; 
  mutable StateOnPlane lastEndState_; 
 
 private:
 
  mutable std::vector<RKStep> RKSteps_; 
  mutable int RKStepsFXStart_; 
  mutable int RKStepsFXStop_; 
  mutable std::vector<ExtrapStep> ExtrapSteps_; 
 
  mutable TMatrixD fJacobian_; 
  mutable TMatrixDSym fNoise_; 
 
  mutable bool useCache_; 
  mutable unsigned int cachePos_; 
 
  // auxiliary variables and arrays
  // needed in Extrap()
  mutable StepLimits limits_; 
  mutable M7x7 noiseArray_; 
  mutable M7x7 noiseProjection_; 
  mutable M7x7 J_MMT_; 
 
 public:
 
  ClassDefOverride(RKTrackRep, 1)
 
};
 
} /* End of namespace genfit */
#endif // genfit_RKTrackRep_h