trackReps/include/RKTrackRep.h

0001 /* Copyright 2008-2010, Technische Universitaet Muenchen,
0002    Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch
0003
0004    This file is part of GENFIT.
0005
0006    GENFIT is free software: you can redistribute it and/or modify
0007    it under the terms of the GNU Lesser General Public License as published
0008    by the Free Software Foundation, either version 3 of the License, or
0009    (at your option) any later version.
0010
0011    GENFIT is distributed in the hope that it will be useful,
0012    but WITHOUT ANY WARRANTY; without even the implied warranty of
0013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0014    GNU Lesser General Public License for more details.
0015
0016    You should have received a copy of the GNU Lesser General Public License
0017    along with GENFIT.  If not, see <http://www.gnu.org/licenses/>.
0018 */
0019
0020 /** @addtogroup RKTrackRep
0021  * @{
0022  */
0023
0024 #ifndef genfit_RKTrackRep_h
0025 #define genfit_RKTrackRep_h
0026
0027 #include "AbsTrackRep.h"
0028 #include "StateOnPlane.h"
0029 #include "RKTools.h"
0030 #include "StepLimits.h"
0031 #include "Material.h"
0032
0033 #include <algorithm>
0034
0035 namespace genfit {
0036
0037 /**
0038  * @brief Helper for RKTrackRep
0039  */
0040 struct RKStep {
0041   MatStep matStep_; // material properties and stepsize
0042   M1x7 state7_; // 7D state vector
0043   StepLimits limits_;
0044
0045   RKStep() {
0046     std::fill(state7_.begin(), state7_.end(), 0);
0047   }
0048 };
0049
0050
0051 /**
0052  * @brief Helper for RKTrackRep
0053  */
0054 struct ExtrapStep {
0055   M7x7 jac7_; // 5D jacobian of transport
0056   M7x7 noise7_; // 5D noise matrix
0057
0058   ExtrapStep() {
0059     std::fill(jac7_.begin(), jac7_.end(), 0);
0060     std::fill(noise7_.begin(), noise7_.end(), 0);
0061   }
0062 };
0063
0064
0065 /**
0066  * @brief AbsTrackRep with 5D track parameterization in plane coordinates: (q/p, u', v', u, v)
0067  *
0068  * q/p is charge over momentum.
0069  * u' and v' are direction tangents.
0070  * u and v are positions on a DetPlane.
0071  */
0072 class RKTrackRep : public AbsTrackRep {
0073     friend class RKTrackRepTests_momMag_Test;
0074     friend class RKTrackRepTests_calcForwardJacobianAndNoise_Test;
0075     friend class RKTrackRepTests_getState7_Test;
0076     friend class RKTrackRepTests_getState5_Test;
0077
0078  public:
0079
0080   RKTrackRep();
0081   RKTrackRep(int pdgCode, char propDir = 0);
0082
0083   virtual ~RKTrackRep();
0084
0085   virtual AbsTrackRep* clone() const override {return new RKTrackRep(*this);}
0086
0087   virtual double extrapolateToPlane(StateOnPlane& state,
0088       const SharedPlanePtr& plane,
0089       bool stopAtBoundary = false,
0090       bool calcJacobianNoise = false) const override;
0091
0092   using AbsTrackRep::extrapolateToLine;
0093
0094   virtual double extrapolateToLine(StateOnPlane& state,
0095       const TVector3& linePoint,
0096       const TVector3& lineDirection,
0097       bool stopAtBoundary = false,
0098       bool calcJacobianNoise = false) const override;
0099
0100   virtual double extrapolateToPoint(StateOnPlane& state,
0101       const TVector3& point,
0102       bool stopAtBoundary = false,
0103       bool calcJacobianNoise = false) const override {
0104     return extrapToPoint(state, point, nullptr, stopAtBoundary, calcJacobianNoise);
0105   }
0106
0107   virtual double extrapolateToPoint(StateOnPlane& state,
0108       const TVector3& point,
0109       const TMatrixDSym& G, // weight matrix (metric)
0110       bool stopAtBoundary = false,
0111       bool calcJacobianNoise = false) const override {
0112     return extrapToPoint(state, point, &G, stopAtBoundary, calcJacobianNoise);
0113   }
0114
0115   virtual double extrapolateToCylinder(StateOnPlane& state,
0116       double radius,
0117       const TVector3& linePoint = TVector3(0.,0.,0.),
0118       const TVector3& lineDirection = TVector3(0.,0.,1.),
0119       bool stopAtBoundary = false,
0120       bool calcJacobianNoise = false) const override;
0121
0122
0123   virtual double extrapolateToCone(StateOnPlane& state,
0124       double radius,
0125       const TVector3& linePoint = TVector3(0.,0.,0.),
0126       const TVector3& lineDirection = TVector3(0.,0.,1.),
0127       bool stopAtBoundary = false,
0128       bool calcJacobianNoise = false) const override ;
0129
0130   virtual double extrapolateToSphere(StateOnPlane& state,
0131       double radius,
0132       const TVector3& point = TVector3(0.,0.,0.),
0133       bool stopAtBoundary = false,
0134       bool calcJacobianNoise = false) const override;
0135
0136   virtual double extrapolateBy(StateOnPlane& state,
0137       double step,
0138       bool stopAtBoundary = false,
0139       bool calcJacobianNoise = false) const override;
0140
0141
0142   unsigned int getDim() const override {return 5;}
0143
0144   virtual TVector3 getPos(const StateOnPlane& state) const override;
0145
0146   virtual TVector3 getMom(const StateOnPlane& state) const override;
0147   virtual void getPosMom(const StateOnPlane& state, TVector3& pos, TVector3& mom) const override;
0148
0149   virtual double getMomMag(const StateOnPlane& state) const override;
0150   virtual double getMomVar(const MeasuredStateOnPlane& state) const override;
0151
0152   virtual TMatrixDSym get6DCov(const MeasuredStateOnPlane& state) const override;
0153   virtual void getPosMomCov(const MeasuredStateOnPlane& state, TVector3& pos, TVector3& mom, TMatrixDSym& cov) const override;
0154   virtual double getCharge(const StateOnPlane& state) const override;
0155   virtual double getQop(const StateOnPlane& state) const override {return state.getState()(0);}
0156   double getSpu(const StateOnPlane& state) const;
0157   double getTime(const StateOnPlane& state) const override;
0158
0159   virtual void getForwardJacobianAndNoise(TMatrixD& jacobian, TMatrixDSym& noise, TVectorD& deltaState) const override;
0160
0161   virtual void getBackwardJacobianAndNoise(TMatrixD& jacobian, TMatrixDSym& noise, TVectorD& deltaState) const override;
0162
0163   std::vector<genfit::MatStep> getSteps() const override;
0164
0165   virtual double getRadiationLenght() const override;
0166
0167   virtual void setPosMom(StateOnPlane& state, const TVector3& pos, const TVector3& mom) const override;
0168   virtual void setPosMom(StateOnPlane& state, const TVectorD& state6) const override;
0169   virtual void setPosMomErr(MeasuredStateOnPlane& state, const TVector3& pos, const TVector3& mom, const TVector3& posErr, const TVector3& momErr) const override;
0170   virtual void setPosMomCov(MeasuredStateOnPlane& state, const TVector3& pos, const TVector3& mom, const TMatrixDSym& cov6x6) const override;
0171   virtual void setPosMomCov(MeasuredStateOnPlane& state, const TVectorD& state6, const TMatrixDSym& cov6x6) const override;
0172
0173   virtual void setChargeSign(StateOnPlane& state, double charge) const override;
0174   virtual void setQop(StateOnPlane& state, double qop) const override {state.getState()(0) = qop;}
0175
0176   void setSpu(StateOnPlane& state, double spu) const;
0177   void setTime(StateOnPlane& state, double time) const override;
0178
0179   //! The actual Runge Kutta propagation
0180   /** propagate state7 with step S. Fills SA (Start directions derivatives dA/S).
0181    *  This is a single Runge-Kutta step.
0182    *  If jacobian is nullptr, only the state is propagated,
0183    *  otherwise also the 7x7 jacobian is calculated.
0184    *  If varField is false, the magnetic field will only be evaluated at the starting position.
0185    *  The return value is an estimation on how good the extrapolation is, and it is usually fine if it is > 1.
0186    *  It gives a suggestion how you must scale S so that the quality will be sufficient.
0187    */
0188   virtual double RKPropagate(M1x7& state7,
0189                              M7x7* jacobian,
0190                              M1x3& SA,
0191                              double S,
0192                              bool varField = true,
0193                              bool calcOnlyLastRowOfJ = false) const;
0194
0195   virtual bool isSameType(const AbsTrackRep* other) override;
0196   virtual bool isSame(const AbsTrackRep* other) override;
0197
0198  private:
0199
0200   void initArrays() const;
0201
0202   virtual double extrapToPoint(StateOnPlane& state,
0203       const TVector3& point,
0204       const TMatrixDSym* G = nullptr, // weight matrix (metric)
0205       bool stopAtBoundary = false,
0206       bool calcJacobianNoise = false) const;
0207
0208   void getState7(const StateOnPlane& state, M1x7& state7) const;
0209   void getState5(StateOnPlane& state, const M1x7& state7) const; // state7 must already lie on plane of state!
0210
0211   void calcJ_pM_5x7(M5x7& J_pM, const TVector3& U, const TVector3& V, const M1x3& pTilde, double spu) const;
0212
0213   void transformPM6(const MeasuredStateOnPlane& state,
0214                     M6x6& out6x6) const;
0215
0216   void calcJ_Mp_7x5(M7x5& J_Mp, const TVector3& U, const TVector3& V, const TVector3& W, const M1x3& A) const;
0217
0218   void calcForwardJacobianAndNoise(const M1x7& startState7, const DetPlane& startPlane,
0219                    const M1x7& destState7, const DetPlane& destPlane) const;
0220
0221   void transformM6P(const M6x6& in6x6,
0222                     const M1x7& state7,
0223                     MeasuredStateOnPlane& state) const; // plane and charge must already be set!
0224
0225   //! Propagates the particle through the magnetic field.
0226   /** If the propagation is successful and the plane is reached, the function returns true.
0227     * Propagated state and the jacobian of the extrapolation are written to state7 and jacobianT.
0228     * The jacobian is only calculated if jacobianT != nullptr.
0229     * In the main loop of the Runge Kutta algorithm, the estimateStep() is called
0230     * and may reduce the estimated stepsize so that a maximum momentum loss will not be exceeded,
0231     * and stop at material boundaries.
0232     * If this is the case, RKutta() will only propagate the reduced distance and then return. This is to ensure that
0233     * material effects, which are calculated after the propagation, are taken into account properly.
0234     */
0235   bool RKutta(const M1x4& SU,
0236               const DetPlane& plane,
0237               double charge,
0238               double mass,
0239               M1x7& state7,
0240               M7x7* jacobianT,
0241               M1x7* J_MMT_unprojected_lastRow,
0242               double& coveredDistance, // signed
0243               double& flightTime,
0244               bool& checkJacProj,
0245               M7x7& noiseProjection,
0246               StepLimits& limits,
0247               bool onlyOneStep = false,
0248               bool calcOnlyLastRowOfJ = false) const;
0249
0250   double estimateStep(const M1x7& state7,
0251                       const M1x4& SU,
0252                       const DetPlane& plane,
0253                       const double& charge,
0254                       double& relMomLoss,
0255                       StepLimits& limits) const;
0256
0257   TVector3 pocaOnLine(const TVector3& linePoint,
0258                      const TVector3& lineDirection,
0259                      const TVector3& point) const;
0260
0261   //! Handles propagation and material effects
0262   /** #extrapolateToPlane(), #extrapolateToPoint() and #extrapolateToLine() etc. call this function.
0263     * #Extrap() needs a plane as an argument, hence #extrapolateToPoint() and #extrapolateToLine() create virtual detector planes.
0264     * In this function, #RKutta() is called and the resulting points and point paths are filtered
0265     * so that the direction doesn't change and tiny steps are filtered out.
0266     * After the propagation the material effects are called via the MaterialEffects singleton.
0267     * #Extrap() will loop until the plane is reached, unless the propagation fails or the maximum number of
0268     * iterations is exceeded.
0269     */
0270   double Extrap(const DetPlane& startPlane, // plane where Extrap starts
0271                 const DetPlane& destPlane, // plane where Extrap has to extrapolate to
0272                 double charge,
0273                 double mass,
0274                 bool& isAtBoundary,
0275                 M1x7& state7,
0276                 double& flightTime,
0277                 bool fillExtrapSteps,
0278                 TMatrixDSym* cov = nullptr,
0279                 bool onlyOneStep = false,
0280                 bool stopAtBoundary = false,
0281                 double maxStep = 1.E99) const;
0282
0283   void checkCache(const StateOnPlane& state, const SharedPlanePtr* plane) const;
0284
0285   double momMag(const M1x7& state7) const;
0286
0287  protected:
0288
0289   mutable StateOnPlane lastStartState_; //! state where the last extrapolation has started
0290   mutable StateOnPlane lastEndState_; //! state where the last extrapolation has ended
0291
0292  private:
0293
0294   mutable std::vector<RKStep> RKSteps_; //! RungeKutta steps made in the last extrapolation
0295   mutable int RKStepsFXStart_; //!
0296   mutable int RKStepsFXStop_; //!
0297   mutable std::vector<ExtrapStep> ExtrapSteps_; //! steps made in Extrap during last extrapolation
0298
0299   mutable TMatrixD fJacobian_; //!
0300   mutable TMatrixDSym fNoise_; //!
0301
0302   mutable bool useCache_; //! use cached RKSteps_ for extrapolation
0303   mutable unsigned int cachePos_; //!
0304
0305   // auxiliary variables and arrays
0306   // needed in Extrap()
0307   mutable StepLimits limits_; //!
0308   mutable M7x7 noiseArray_; //! noise matrix of the last extrapolation
0309   mutable M7x7 noiseProjection_; //!
0310   mutable M7x7 J_MMT_; //!
0311
0312  public:
0313
0314   ClassDefOverride(RKTrackRep, 1)
0315
0316 };
0317
0318 } /* End of namespace genfit */
0319 /** @} */
0320
0321 #endif // genfit_RKTrackRep_h