genfit2/trackReps/include/RKTrackRep.h

Go to the documentation of this file.

//****************************************************************************
//*                   This file is part of PandaRoot.                        *
//*                                                                          *
//*            PandaRoot is distributed under the terms of the               *
//*              GNU General Public License (GPL) version 3,                 *
//*                 copied verbatim in the file "LICENSE".                   *
//*                                                                          *
//*  Copyright (C) 2006 - 2024 FAIR GmbH and copyright holders of PandaRoot  *
//*     The copyright holders are listed in the file "COPYRIGHTHOLDERS".     *
//*               The authors are listed in the file "AUTHORS".              *
//****************************************************************************

/* 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"

namespace genfit {

struct RKStep {
  MatStep matStep_; // material properties and stepsize
  M1x7 state7_;     // 7D state vector
  StepLimits limits_;

  RKStep() { memset(state7_, 0x00, 7 * sizeof(double)); }
};

struct ExtrapStep {
  M7x7 jac7_;   // 5D jacobian of transport
  M7x7 noise7_; // 5D noise matrix

  ExtrapStep()
  {
    memset(jac7_, 0, sizeof(M7x7));
    memset(noise7_, 0, sizeof(M7x7));
  }
};

class RKTrackRep : public AbsTrackRep {

 public:
  RKTrackRep();
  RKTrackRep(int pdgCode, char propDir = 0);

  virtual ~RKTrackRep();

  virtual AbsTrackRep *clone() const { return new RKTrackRep(*this); }

  virtual double extrapolateToPlane(StateOnPlane &state, const SharedPlanePtr &plane, bool stopAtBoundary = false, bool calcJacobianNoise = false) const;

  using AbsTrackRep::extrapolateToLine;

  virtual double
  extrapolateToLine(StateOnPlane &state, const TVector3 &linePoint, const TVector3 &lineDirection, bool stopAtBoundary = false, bool calcJacobianNoise = false) const;

  virtual double extrapolateToPoint(StateOnPlane &state, const TVector3 &point, bool stopAtBoundary = false, bool calcJacobianNoise = false) const
  {
    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
  {
    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;

  virtual double
  extrapolateToSphere(StateOnPlane &state, double radius, const TVector3 &point = TVector3(0., 0., 0.), bool stopAtBoundary = false, bool calcJacobianNoise = false) const;

  virtual double extrapolateBy(StateOnPlane &state, double step, bool stopAtBoundary = false, bool calcJacobianNoise = false) const;

  unsigned int getDim() const { return 5; }

  virtual TVector3 getPos(const StateOnPlane &state) const;

  virtual TVector3 getMom(const StateOnPlane &state) const;
  virtual void getPosMom(const StateOnPlane &state, TVector3 &pos, TVector3 &mom) const;

  virtual double getMomMag(const StateOnPlane &state) const;
  virtual double getMomVar(const MeasuredStateOnPlane &state) const;

  virtual TMatrixDSym get6DCov(const MeasuredStateOnPlane &state) const;
  virtual void getPosMomCov(const MeasuredStateOnPlane &state, TVector3 &pos, TVector3 &mom, TMatrixDSym &cov) const;
  virtual double getCharge(const StateOnPlane &state) const;
  virtual double getQop(const StateOnPlane &state) const { return state.getState()(0); }
  double getSpu(const StateOnPlane &state) const;
  double getTime(const StateOnPlane &state) const;

  virtual void getForwardJacobianAndNoise(TMatrixD &jacobian, TMatrixDSym &noise, TVectorD &deltaState) const;

  virtual void getBackwardJacobianAndNoise(TMatrixD &jacobian, TMatrixDSym &noise, TVectorD &deltaState) const;

  std::vector<genfit::MatStep> getSteps() const;

  virtual double getRadiationLenght() const;

  virtual void setPosMom(StateOnPlane &state, const TVector3 &pos, const TVector3 &mom) const;
  virtual void setPosMom(StateOnPlane &state, const TVectorD &state6) const;
  virtual void setPosMomErr(MeasuredStateOnPlane &state, const TVector3 &pos, const TVector3 &mom, const TVector3 &posErr, const TVector3 &momErr) const;
  virtual void setPosMomCov(MeasuredStateOnPlane &state, const TVector3 &pos, const TVector3 &mom, const TMatrixDSym &cov6x6) const;
  virtual void setPosMomCov(MeasuredStateOnPlane &state, const TVectorD &state6, const TMatrixDSym &cov6x6) const;

  virtual void setChargeSign(StateOnPlane &state, double charge) const;
  virtual void setQop(StateOnPlane &state, double qop) const { state.getState()(0) = qop; }

  void setSpu(StateOnPlane &state, double spu) const;
  void setTime(StateOnPlane &state, double time) const;


  double RKPropagate(M1x7 &state7, M7x7 *jacobian, M1x3 &SA, double S, bool varField = true, bool calcOnlyLastRowOfJ = false) const;

  virtual bool isSameType(const AbsTrackRep *other);
  virtual bool isSame(const AbsTrackRep *other);

 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 transformPM7(const MeasuredStateOnPlane &state, M7x7 &out7x7) const;

  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 transformM7P(const M7x7 &in7x7, const M1x7 &state7,
                    MeasuredStateOnPlane &state) const; // plane must already be set!

  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;

  mutable StateOnPlane lastStartState_;         
  mutable StateOnPlane lastEndState_;           
  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:
  ClassDef(RKTrackRep, 1)
};

} /* End of namespace genfit */
#endif // genfit_RKTrackRep_h