Loading...
Searching...
No Matches
genfit::RKTrackRep Class Reference
AbsTrackRep with 5D track parameterization in plane coordinates: (q/p, u', v', u, v) More...
#include <RKTrackRep.h>
Inheritance diagram for genfit::RKTrackRep:
Collaboration diagram for genfit::RKTrackRep:
Public Member Functions | |
| RKTrackRep () | |
| RKTrackRep (int pdgCode, char propDir=0) | |
| virtual | ~RKTrackRep () |
| virtual AbsTrackRep * | clone () const |
| Clone the trackRep. | |
| virtual double | extrapolateToPlane (StateOnPlane &state, const SharedPlanePtr &plane, bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| Extrapolates the state to plane, and returns the extrapolation length and, via reference, the extrapolated state. | |
| virtual double | extrapolateToLine (StateOnPlane &state, const TVector3 &linePoint, const TVector3 &lineDirection, bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| Extrapolates the state to the POCA to a line, and returns the extrapolation length and, via reference, the extrapolated state. | |
| virtual double | extrapolateToPoint (StateOnPlane &state, const TVector3 &point, bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| Extrapolates the state to the POCA to a point, and returns the extrapolation length and, via reference, the extrapolated state. | |
| virtual double | extrapolateToPoint (StateOnPlane &state, const TVector3 &point, const TMatrixDSym &G, bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| Extrapolates the state to the POCA to a point in the metric of G, and returns the extrapolation length and, via reference, the extrapolated state. | |
| 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 |
| Extrapolates the state to the cylinder surface, and returns the extrapolation length and, via reference, the extrapolated state. | |
| virtual double | extrapolateToSphere (StateOnPlane &state, double radius, const TVector3 &point=TVector3(0., 0., 0.), bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| Extrapolates the state to the sphere surface, and returns the extrapolation length and, via reference, the extrapolated state. | |
| virtual double | extrapolateBy (StateOnPlane &state, double step, bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| Extrapolates the state by step (cm) and returns the extrapolation length and, via reference, the extrapolated state. | |
| unsigned int | getDim () const |
| Get the dimension of the state vector used by the track representation. | |
| virtual TVector3 | getPos (const StateOnPlane &state) const |
| Get the cartesian position of a state. | |
| virtual TVector3 | getMom (const StateOnPlane &state) const |
| Get the cartesian momentum vector of a state. | |
| virtual void | getPosMom (const StateOnPlane &state, TVector3 &pos, TVector3 &mom) const |
| Get cartesian position and momentum vector of a state. | |
| virtual double | getMomMag (const StateOnPlane &state) const |
| get the magnitude of the momentum in GeV. | |
| virtual double | getMomVar (const MeasuredStateOnPlane &state) const |
| get the variance of the absolute value of the momentum . | |
| virtual TMatrixDSym | get6DCov (const MeasuredStateOnPlane &state) const |
| Get the 6D covariance. | |
| virtual void | getPosMomCov (const MeasuredStateOnPlane &state, TVector3 &pos, TVector3 &mom, TMatrixDSym &cov) const |
| Translates MeasuredStateOnPlane into 3D position, momentum and 6x6 covariance. | |
| virtual double | getCharge (const StateOnPlane &state) const |
| Get the (fitted) charge of a state. This is not always equal the pdg charge (e.g. if the charge sign was flipped during the fit). | |
| virtual double | getQop (const StateOnPlane &state) const |
| Get charge over momentum. | |
| double | getSpu (const StateOnPlane &state) const |
| virtual void | getForwardJacobianAndNoise (TMatrixD &jacobian, TMatrixDSym &noise, TVectorD &deltaState) const |
| Get the jacobian and noise matrix of the last extrapolation. | |
| virtual void | getBackwardJacobianAndNoise (TMatrixD &jacobian, TMatrixDSym &noise, TVectorD &deltaState) const |
| Get the jacobian and noise matrix of the last extrapolation if it would have been done in opposite direction. | |
| std::vector< genfit::MatStep > | getSteps () const |
| Get stepsizes and material properties of crossed materials of the last extrapolation. | |
| virtual double | getRadiationLenght () const |
| Get the accumulated X/X0 (path / radiation length) of the material crossed in the last extrapolation. | |
| virtual double | getTOF () const |
| Get the time of flight [ns] of the last extrapolation. | |
| virtual void | setPosMom (StateOnPlane &state, const TVector3 &pos, const TVector3 &mom) const |
| Set position and momentum of state. | |
| virtual void | setPosMom (StateOnPlane &state, const TVectorD &state6) const |
| Set position and momentum of state. | |
| virtual void | setPosMomErr (MeasuredStateOnPlane &state, const TVector3 &pos, const TVector3 &mom, const TVector3 &posErr, const TVector3 &momErr) const |
| Set position and momentum and error of state. | |
| virtual void | setPosMomCov (MeasuredStateOnPlane &state, const TVector3 &pos, const TVector3 &mom, const TMatrixDSym &cov6x6) const |
| Set position, momentum and covariance of state. | |
| virtual void | setPosMomCov (MeasuredStateOnPlane &state, const TVectorD &state6, const TMatrixDSym &cov6x6) const |
| Set position, momentum and covariance of state. | |
| virtual void | setChargeSign (StateOnPlane &state, double charge) const |
| Set the sign of the charge according to charge. | |
| virtual void | setQop (StateOnPlane &state, double qop) const |
| Set charge/momentum. | |
| void | setSpu (StateOnPlane &state, double spu) const |
| double | RKPropagate (M1x7 &state7, M7x7 *jacobian, M1x3 &SA, double S, bool varField=true, bool calcOnlyLastRowOfJ=false) const |
| The actual Runge Kutta propagation. | |
| virtual bool | isSameType (const AbsTrackRep *other) |
| check if other is of same type (e.g. RKTrackRep). | |
| virtual bool | isSame (const AbsTrackRep *other) |
| check if other is of same type (e.g. RKTrackRep) and has same pdg code. | |
| RKTrackRep () | |
| RKTrackRep (int pdgCode, char propDir=0) | |
| virtual | ~RKTrackRep () |
| virtual AbsTrackRep * | clone () const |
| Clone the trackRep. | |
| SharedPlanePtr | makePlane (const TVector3 &o, const TVector3 &u, const TVector3 &v) |
| virtual double | extrapolateToPlane (StateOnPlane &state, const SharedPlanePtr &plane, bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| Extrapolates the state to plane, and returns the extrapolation length and, via reference, the extrapolated state. | |
| virtual double | extrapolateToPlane (StateOnPlane &state, const TVector3 &point, const TVector3 &dir, bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| virtual double | extrapolateToLine (StateOnPlane &state, const TVector3 &linePoint, const TVector3 &lineDirection, bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| Extrapolates the state to the POCA to a line, and returns the extrapolation length and, via reference, the extrapolated state. | |
| virtual double | extrapolateToPoint (StateOnPlane &state, const TVector3 &point, bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| Extrapolates the state to the POCA to a point, and returns the extrapolation length and, via reference, the extrapolated state. | |
| virtual double | extrapolateToPoint (StateOnPlane &state, const TVector3 &point, const TMatrixDSym &G, bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| Extrapolates the state to the POCA to a point in the metric of G, and returns the extrapolation length and, via reference, the extrapolated state. | |
| 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 |
| Extrapolates the state to the cylinder surface, and returns the extrapolation length and, via reference, the extrapolated state. | |
| virtual double | extrapolateToSphere (StateOnPlane &state, double radius, const TVector3 &point=TVector3(0., 0., 0.), bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| Extrapolates the state to the sphere surface, and returns the extrapolation length and, via reference, the extrapolated state. | |
| virtual double | extrapolateBy (StateOnPlane &state, double step, bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| Extrapolates the state by step (cm) and returns the extrapolation length and, via reference, the extrapolated state. | |
| unsigned int | getDim () const |
| Get the dimension of the state vector used by the track representation. | |
| virtual TVector3 | getPos (const StateOnPlane &state) const |
| Get the cartesian position of a state. | |
| virtual TVector3 | getMom (const StateOnPlane &state) const |
| Get the cartesian momentum vector of a state. | |
| virtual void | getPosMom (const StateOnPlane &state, TVector3 &pos, TVector3 &mom) const |
| Get cartesian position and momentum vector of a state. | |
| virtual double | getMomMag (const StateOnPlane &state) const |
| get the magnitude of the momentum in GeV. | |
| virtual double | getMomVar (const MeasuredStateOnPlane &state) const |
| get the variance of the absolute value of the momentum . | |
| virtual TMatrixDSym | get6DCov (const MeasuredStateOnPlane &state) const |
| Get the 6D covariance. | |
| virtual void | getPosMomCov (const MeasuredStateOnPlane &state, TVector3 &pos, TVector3 &mom, TMatrixDSym &cov) const |
| Translates MeasuredStateOnPlane into 3D position, momentum and 6x6 covariance. | |
| virtual double | getCharge (const StateOnPlane &state) const |
| Get the (fitted) charge of a state. This is not always equal the pdg charge (e.g. if the charge sign was flipped during the fit). | |
| virtual double | getQop (const StateOnPlane &state) const |
| Get charge over momentum. | |
| double | getSpu (const StateOnPlane &state) const |
| virtual void | getForwardJacobianAndNoise (TMatrixD &jacobian, TMatrixDSym &noise, TVectorD &deltaState) const |
| Get the jacobian and noise matrix of the last extrapolation. | |
| virtual void | getBackwardJacobianAndNoise (TMatrixD &jacobian, TMatrixDSym &noise, TVectorD &deltaState) const |
| Get the jacobian and noise matrix of the last extrapolation if it would have been done in opposite direction. | |
| std::vector< genfit::MatStep > | getSteps () const |
| Get stepsizes and material properties of crossed materials of the last extrapolation. | |
| virtual double | getRadiationLenght () const |
| Get the accumulated X/X0 (path / radiation length) of the material crossed in the last extrapolation. | |
| virtual double | getTOF () const |
| Get the time of flight [ns] of the last extrapolation. | |
| virtual void | setPosMom (StateOnPlane &state, const TVector3 &pos, const TVector3 &mom) const |
| Set position and momentum of state. | |
| virtual void | setPosMom (StateOnPlane &state, const TVectorD &state6) const |
| Set position and momentum of state. | |
| virtual void | setPosMomErr (MeasuredStateOnPlane &state, const TVector3 &pos, const TVector3 &mom, const TVector3 &posErr, const TVector3 &momErr) const |
| Set position and momentum and error of state. | |
| virtual void | setPosMomCov (MeasuredStateOnPlane &state, const TVector3 &pos, const TVector3 &mom, const TMatrixDSym &cov6x6) const |
| Set position, momentum and covariance of state. | |
| virtual void | setPosMomCov (MeasuredStateOnPlane &state, const TVectorD &state6, const TMatrixDSym &cov6x6) const |
| Set position, momentum and covariance of state. | |
| virtual void | setChargeSign (StateOnPlane &state, double charge) const |
| Set the sign of the charge according to charge. | |
| virtual void | setQop (StateOnPlane &state, double qop) const |
| Set charge/momentum. | |
| void | setSpu (StateOnPlane &state, double spu) const |
| double | RKPropagate (M1x7 &state7, M7x7 *jacobian, M1x3 &SA, double S, bool varField=true, bool calcOnlyLastRowOfJ=false) const |
| The actual Runge Kutta propagation. | |
| virtual bool | isSameType (const AbsTrackRep *other) |
| check if other is of same type (e.g. RKTrackRep). | |
| virtual bool | isSame (const AbsTrackRep *other) |
| check if other is of same type (e.g. RKTrackRep) and has same pdg code. | |
| virtual double | extrapolateToLine (StateOnPlane &state, const TVector3 &linePoint, const TVector3 &lineDirection, bool stopAtBoundary=false, bool calcJacobianNoise=false) const=0 |
| Extrapolates the state to the POCA to a line, and returns the extrapolation length and, via reference, the extrapolated state. | |
| virtual double | extrapolateToLine (StateOnPlane &state, const TVector3 &point1, const TVector3 &point2, TVector3 &poca, TVector3 &dirInPoca, TVector3 &poca_onwire, bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| Resembles the interface of GFAbsTrackRep in old versions of genfit. | |
| virtual double | extrapolateToLine (StateOnPlane &state, const TVector3 &linePoint, const TVector3 &lineDirection, bool stopAtBoundary=false, bool calcJacobianNoise=false) const=0 |
| Extrapolates the state to the POCA to a line, and returns the extrapolation length and, via reference, the extrapolated state. | |
| virtual double | extrapolateToLine (StateOnPlane &state, const TVector3 &point1, const TVector3 &point2, TVector3 &poca, TVector3 &dirInPoca, TVector3 &poca_onwire, bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| Resembles the interface of GFAbsTrackRep in old versions of genfit. | |
 Public Member Functions inherited from genfit::AbsTrackRep | |
| AbsTrackRep () | |
| AbsTrackRep (int pdgCode, char propDir=0) | |
| virtual | ~AbsTrackRep () |
| double | extrapolateToMeasurement (StateOnPlane &state, const AbsMeasurement *measurement, bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| extrapolate to an AbsMeasurement | |
| TVector3 | getDir (const StateOnPlane &state) const |
| Get the direction vector of a state. | |
| void | getPosDir (const StateOnPlane &state, TVector3 &pos, TVector3 &dir) const |
| Get cartesian position and direction vector of a state. | |
| virtual TVectorD | get6DState (const StateOnPlane &state) const |
| Get the 6D state vector (x, y, z, p_x, p_y, p_z). | |
| virtual void | get6DStateCov (const MeasuredStateOnPlane &state, TVectorD &stateVec, TMatrixDSym &cov) const |
| Translates MeasuredStateOnPlane into 6D state vector (x, y, z, p_x, p_y, p_z) and 6x6 covariance. | |
| int | getPDG () const |
| Get the pdg code. | |
| double | getPDGCharge () const |
| Get the charge of the particle of the pdg code. | |
| double | getMass (const StateOnPlane &state) const |
| Get tha particle mass in GeV/c^2. | |
| char | getPropDir () const |
| Get propagation direction. (-1, 0, 1) -> (backward, auto, forward). | |
| void | calcJacobianNumerically (const genfit::StateOnPlane &origState, const genfit::SharedPlanePtr destPlane, TMatrixD &jacobian) const |
| Calculate Jacobian of transportation numerically. Slow but accurate. Can be used to validate (semi)analytic calculations. | |
| bool | switchPDGSign () |
| try to multiply pdg code with -1. (Switch from particle to anti-particle and vice versa). | |
| void | setPropDir (int dir) |
| Set propagation direction. (-1, 0, 1) -> (backward, auto, forward). | |
| void | switchPropDir () |
| Switch propagation direction. Has no effect if propDir_ is set to 0. | |
| virtual void | setDebugLvl (unsigned int lvl=1) |
| virtual void | Print (const Option_t *="") const |
Private Member Functions | |
| void | initArrays () const |
| virtual double | extrapToPoint (StateOnPlane &state, const TVector3 &point, const TMatrixDSym *G=NULL, bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| void | getState7 (const StateOnPlane &state, M1x7 &state7) const |
| void | getState5 (StateOnPlane &state, const M1x7 &state7) const |
| 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 |
| 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 |
| bool | RKutta (const M1x4 &SU, const DetPlane &plane, double charge, M1x7 &state7, M7x7 *jacobianT, double &coveredDistance, bool &checkJacProj, M7x7 &noiseProjection, StepLimits &limits, bool onlyOneStep=false, bool calcOnlyLastRowOfJ=false) const |
| Propagates the particle through the magnetic field. | |
| 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, const DetPlane &destPlane, double charge, bool &isAtBoundary, M1x7 &state7, bool fillExtrapSteps, TMatrixDSym *cov=nullptr, bool onlyOneStep=false, bool stopAtBoundary=false, double maxStep=1.E99) const |
| Handles propagation and material effects. | |
| void | checkCache (const StateOnPlane &state, const SharedPlanePtr *plane) const |
| double | momMag (const M1x7 &state7) const |
| void | initArrays () const |
| virtual double | extrapToPoint (StateOnPlane &state, const TVector3 &point, const TMatrixDSym *G=NULL, bool stopAtBoundary=false, bool calcJacobianNoise=false) const |
| void | getState7 (const StateOnPlane &state, M1x7 &state7) const |
| void | getState5 (StateOnPlane &state, const M1x7 &state7) const |
| 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 |
| 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 |
| bool | RKutta (const M1x4 &SU, const DetPlane &plane, double charge, M1x7 &state7, M7x7 *jacobianT, double &coveredDistance, bool &checkJacProj, M7x7 &noiseProjection, StepLimits &limits, bool onlyOneStep=false, bool calcOnlyLastRowOfJ=false) const |
| Propagates the particle through the magnetic field. | |
| 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, const DetPlane &destPlane, double charge, bool &isAtBoundary, M1x7 &state7, bool fillExtrapSteps, TMatrixDSym *cov=nullptr, bool onlyOneStep=false, bool stopAtBoundary=false, double maxStep=1.E99) const |
| Handles propagation and material effects. | |
| void | checkCache (const StateOnPlane &state, const SharedPlanePtr *plane) const |
| double | momMag (const M1x7 &state7) const |
Private Attributes | |
| StateOnPlane | lastStartState_ |
| StateOnPlane | lastEndState_ |
| state where the last extrapolation has started | |
| std::vector< RKStep > | RKSteps_ |
| state where the last extrapolation has ended | |
| int | RKStepsFXStart_ |
| RungeKutta steps made in the last extrapolation. | |
| int | RKStepsFXStop_ |
| std::vector< ExtrapStep > | ExtrapSteps_ |
| TMatrixD | fJacobian_ |
| steps made in Extrap during last extrapolation | |
| TMatrixDSym | fNoise_ |
| bool | useCache_ |
| unsigned int | cachePos_ |
| use cached RKSteps_ for extrapolation | |
| StepLimits | limits_ |
| M7x7 | noiseArray_ |
| M7x7 | noiseProjection_ |
| noise matrix of the last extrapolation | |
| M7x7 | J_MMT_ |
Additional Inherited Members | |
| Protected Member Functions inherited from genfit::AbsTrackRep | |
| AbsTrackRep (const AbsTrackRep &) | |
| protect from calling copy c'tor from outside the class. Use clone() if you want a copy! | |
| AbsTrackRep & | operator= (const AbsTrackRep &) |
| protect from calling assignment operator from outside the class. Use clone() instead! | |
| Protected Attributes inherited from genfit::AbsTrackRep | |
| int | pdgCode_ |
| Particle code. | |
| char | propDir_ |
| propagation direction (-1, 0, 1) -> (backward, auto, forward) | |
| unsigned int | debugLvl_ |
Detailed Description
AbsTrackRep with 5D track parameterization in plane coordinates: (q/p, u', v', u, v)
q/p is charge over momentum. u' and v' are direction tangents. u and v are positions on a DetPlane.
Definition at line 70 of file RKTrackRep.h.
Constructor & Destructor Documentation
◆ RKTrackRep() [1/4]
| genfit::RKTrackRep::RKTrackRep | ( | ) |
Definition at line 45 of file RKTrackRep.cc.
45 :
46 AbsTrackRep(),
49 RKStepsFXStart_(0),
50 RKStepsFXStop_(0),
51 fJacobian_(5,5),
52 fNoise_(5),
54 cachePos_(0)
55{
56 initArrays();
57}
◆ RKTrackRep() [2/4]
| genfit::RKTrackRep::RKTrackRep | ( | int | pdgCode, |
| char | propDir = 0 |
||
| ) |
Definition at line 60 of file RKTrackRep.cc.
60 :
61 AbsTrackRep(pdgCode, propDir),
64 RKStepsFXStart_(0),
65 RKStepsFXStop_(0),
66 fJacobian_(5,5),
67 fNoise_(5),
69 cachePos_(0)
70{
71 initArrays();
72}
◆ ~RKTrackRep() [1/2]
|
virtual |
Definition at line 75 of file RKTrackRep.cc.
75 {
76 ;
77}
◆ RKTrackRep() [3/4]
| genfit::RKTrackRep::RKTrackRep | ( | ) |
◆ RKTrackRep() [4/4]
| genfit::RKTrackRep::RKTrackRep | ( | int | pdgCode, |
| char | propDir = 0 |
||
| ) |
◆ ~RKTrackRep() [2/2]
|
virtual |
Member Function Documentation
◆ calcForwardJacobianAndNoise() [1/2]
|
private |
Definition at line 774 of file RKTrackRep.cc.
775 {
776
778 std::cout << "RKTrackRep::calcForwardJacobianAndNoise " << std::endl;
779 }
780
782 Exception exc("RKTrackRep::calcForwardJacobianAndNoise ==> cache is empty. Extrapolation must run with a MeasuredStateOnPlane.",__LINE__,__FILE__);
783 throw exc;
784 }
785
786 // The Jacobians returned from RKutta are transposed.
787 TMatrixD jac(TMatrixD::kTransposed, TMatrixD(7, 7, ExtrapSteps_.back().jac7_));
788 TMatrixDSym noise(7, ExtrapSteps_.back().noise7_);
790 noise += TMatrixDSym(7, ExtrapSteps_[i].noise7_).Similarity(jac);
791 jac *= TMatrixD(TMatrixD::kTransposed, TMatrixD(7, 7, ExtrapSteps_[i].jac7_));
792 }
793
794 // Project into 5x5 space.
795 M1x3 pTilde = {startState7[3], startState7[4], startState7[5]};
796 const TVector3& normal = startPlane.getNormal();
797 double pTildeW = pTilde[0] * normal.X() + pTilde[1] * normal.Y() + pTilde[2] * normal.Z();
798 double spu = pTildeW > 0 ? 1 : -1;
801 }
802 M5x7 J_pM;
803 calcJ_pM_5x7(J_pM, startPlane.getU(), startPlane.getV(), pTilde, spu);
804 M7x5 J_Mp;
805 calcJ_Mp_7x5(J_Mp, destPlane.getU(), destPlane.getV(), destPlane.getNormal(), *((M1x3*) &destState7[3]));
806 jac.Transpose(jac); // Because the helper function wants transposed input.
811
815 }
816
817}
void calcJ_Mp_7x5(M7x5 &J_Mp, const TVector3 &U, const TVector3 &V, const TVector3 &W, const M1x3 &A) const
Definition RKTrackRep.cc:1582
void calcJ_pM_5x7(M5x7 &J_pM, const TVector3 &U, const TVector3 &V, const M1x3 &pTilde, double spu) const
Definition RKTrackRep.cc:1452
void J_pMTTxJ_MMTTxJ_MpTT(const M7x5 &J_pMT, const M7x7 &J_MMT, const M5x7 &J_MpT, M5x5 &J_pp)
Definition RKTools.cc:508
void J_MpTxcov7xJ_Mp(const M7x5 &J_Mp, const M7x7 &cov7, M5x5 &out5)
Definition RKTools.cc:196
◆ calcForwardJacobianAndNoise() [2/2]
|
private |
◆ calcJ_Mp_7x5() [1/2]
|
private |
Definition at line 1582 of file RKTrackRep.cc.
1582 {
1583
1584 /*if (debugLvl_ > 1) {
1585 std::cout << "RKTrackRep::calcJ_Mp_7x5 \n";
1586 std::cout << " U = "; U.Print();
1587 std::cout << " V = "; V.Print();
1588 std::cout << " W = "; W.Print();
1589 std::cout << " A = "; RKTools::printDim(A, 3,1);
1590 }*/
1591
1593
1597
1598 // J_Mp matrix is d(q/p,u',v',u,v) / d(x,y,z,ax,ay,az,q/p) (in is 7x7)
1599
1600 // d(u')/d(ax,ay,az)
1601 double fact = 1./(AtW*AtW);
1602 J_Mp[16] = fact * (U.X()*AtW - W.X()*AtU); // [3][1]
1603 J_Mp[21] = fact * (U.Y()*AtW - W.Y()*AtU); // [4][1]
1604 J_Mp[26] = fact * (U.Z()*AtW - W.Z()*AtU); // [5][1]
1605 // d(v')/d(ax,ay,az)
1609 // d(q/p)/d(q/p)
1610 J_Mp[30] = 1.; // [6][0] - not needed for array matrix multiplication
1611 //d(u)/d(x,y,z)
1612 J_Mp[3] = U.X(); // [0][3]
1613 J_Mp[8] = U.Y(); // [1][3]
1614 J_Mp[13] = U.Z(); // [2][3]
1615 //d(v)/d(x,y,z)
1619
1620 /*if (debugLvl_ > 1) {
1621 std::cout << " J_Mp_7x5_ = "; RKTools::printDim(J_Mp, 7,5);
1622 }*/
1623
1624}
◆ calcJ_Mp_7x5() [2/2]
|
private |
◆ calcJ_pM_5x7() [1/2]
|
private |
Definition at line 1452 of file RKTrackRep.cc.
1452 {
1453 /*if (debugLvl_ > 1) {
1454 std::cout << "RKTrackRep::calcJ_pM_5x7 \n";
1455 std::cout << " U = "; U.Print();
1456 std::cout << " V = "; V.Print();
1457 std::cout << " pTilde = "; RKTools::printDim(pTilde, 3,1);
1458 std::cout << " spu = " << spu << "\n";
1459 }*/
1460
1462
1463 const double pTildeMag = sqrt(pTilde[0]*pTilde[0] + pTilde[1]*pTilde[1] + pTilde[2]*pTilde[2]);
1464 const double pTildeMag2 = pTildeMag*pTildeMag;
1465
1466 const double utpTildeOverpTildeMag2 = (U.X()*pTilde[0] + U.Y()*pTilde[1] + U.Z()*pTilde[2]) / pTildeMag2;
1467 const double vtpTildeOverpTildeMag2 = (V.X()*pTilde[0] + V.Y()*pTilde[1] + V.Z()*pTilde[2]) / pTildeMag2;
1468
1469 //J_pM matrix is d(x,y,z,ax,ay,az,q/p) / d(q/p,u',v',u,v) (out is 7x7)
1470
1471 // d(x,y,z)/d(u)
1472 J_pM[21] = U.X(); // [3][0]
1473 J_pM[22] = U.Y(); // [3][1]
1474 J_pM[23] = U.Z(); // [3][2]
1475 // d(x,y,z)/d(v)
1479 // d(q/p)/d(q/p)
1480 J_pM[6] = 1.; // not needed for array matrix multiplication
1481 // d(ax,ay,az)/d(u')
1482 double fact = spu / pTildeMag;
1483 J_pM[10] = fact * ( U.X() - pTilde[0]*utpTildeOverpTildeMag2 ); // [1][3]
1484 J_pM[11] = fact * ( U.Y() - pTilde[1]*utpTildeOverpTildeMag2 ); // [1][4]
1485 J_pM[12] = fact * ( U.Z() - pTilde[2]*utpTildeOverpTildeMag2 ); // [1][5]
1486 // d(ax,ay,az)/d(v')
1490
1491 /*if (debugLvl_ > 1) {
1492 std::cout << " J_pM_5x7_ = "; RKTools::printDim(J_pM_5x7_, 5,7);
1493 }*/
1494}
◆ calcJ_pM_5x7() [2/2]
|
private |
◆ checkCache() [1/2]
|
private |
Definition at line 2412 of file RKTrackRep.cc.
2412 {
2413
2414 if (state.getRep() != this){
2415 Exception exc("RKTrackRep::checkCache ==> state is defined wrt. another TrackRep",__LINE__,__FILE__);
2416 exc.setFatal();
2417 throw exc;
2418 }
2419
2420 if (dynamic_cast<const MeasurementOnPlane*>(&state) != NULL) {
2421 Exception exc("RKTrackRep::checkCache - cannot extrapolate MeasurementOnPlane",__LINE__,__FILE__);
2422 exc.setFatal();
2423 throw exc;
2424 }
2425
2426 cachePos_ = 0;
2427 RKStepsFXStart_ = 0;
2428 RKStepsFXStop_ = 0;
2429 ExtrapSteps_.clear();
2430 initArrays();
2431
2432
2433 if (plane != NULL &&
2440
2441 // clean up cache. Only use steps with same sign.
2442 double firstStep(0);
2444 if (i == 0) {
2445 firstStep = RKSteps_.at(0).matStep_.stepSize_;
2446 continue;
2447 }
2449 if (RKSteps_.at(i-1).matStep_.materialProperties_ == RKSteps_.at(i).matStep_.materialProperties_) {
2451 }
2453 }
2454 }
2455
2457 std::cout << "RKTrackRep::checkCache: use cached material and step values.\n";
2458 }
2459 }
2460 else {
2461
2463 std::cout << "RKTrackRep::checkCache: can NOT use cached material and step values.\n";
2464
2465 if (plane != NULL) {
2467 std::cout << "state.getPlane() != lastStartState_.getPlane()\n";
2468 }
2469 else {
2471 std::cout << "state.getState() != lastStartState_.getState()\n";
2472 }
2475 }
2476 }
2477 }
2478 }
2479
2481 RKSteps_.clear();
2482
2484 }
2485}
std::vector< RKStep > RKSteps_
state where the last extrapolation has ended
Definition RKTrackRep.h:279
virtual TVector3 getPos(const StateOnPlane &state) const
Get the cartesian position of a state.
Definition RKTrackRep.cc:669
void setStatePlane(const TVectorD &state, const SharedPlanePtr &plane)
Definition StateOnPlane.h:68
◆ checkCache() [2/2]
|
private |
◆ clone() [1/2]
|
inlinevirtual |
◆ clone() [2/2]
|
inlinevirtual |
◆ estimateStep() [1/2]
|
private |
invalid: RKSteps_.back().state7_ = { state7[0], state7[1], state7[2], state7[3], state7[4], state7[5], state7[6] };
fNoMaterial
Definition at line 1963 of file RKTrackRep.cc.
1968 {
1969
1973 }
1974 else {
1976 // we need to step exactly to the plane, so don't use the cache!
1979 }
1980 else {
1982 std::cout << " RKTrackRep::estimateStep: use stepSize " << cachePos_ << " from cache: " << RKSteps_.at(cachePos_).matStep_.stepSize_ << "\n";
1983 }
1984 //for(int n = 0; n < 1*7; ++n) RKSteps_[cachePos_].state7_[n] = state7[n];
1985 ++RKStepsFXStop_;
1988 }
1989 }
1990 }
1991
1993
1995 std::cout << " RKTrackRep::estimateStep \n";
1996 std::cout << " position: "; TVector3(state7[0], state7[1], state7[2]).Print();
1997 std::cout << " direction: "; TVector3(state7[3], state7[4], state7[5]).Print();
1998 }
1999
2000 // calculate SL distance to surface
2001 double Dist = SU[3] - (state7[0]*SU[0] +
2002 state7[1]*SU[1] +
2003 state7[2]*SU[2]); // Distance between start coordinates and surface
2004 double An = state7[3]*SU[0] +
2005 state7[4]*SU[1] +
2006 state7[5]*SU[2]; // An = dir * N; component of dir normal to surface
2007
2008 double SLDist; // signed
2009 if (fabs(An) > 1.E-10)
2010 SLDist = Dist/An;
2011 else {
2012 SLDist = Dist*1.E10;
2013 if (An<0) SLDist *= -1.;
2014 }
2015
2016 limits.setLimit(stp_plane, SLDist);
2017 limits.setStepSign(SLDist);
2018
2020 std::cout << " Distance to plane: " << Dist << "\n";
2021 std::cout << " SL distance to plane: " << SLDist << "\n";
2022 if (limits.getStepSign()>0)
2023 std::cout << " Direction is pointing towards surface.\n";
2024 else
2025 std::cout << " Direction is pointing away from surface.\n";
2026 }
2027 // DONE calculate SL distance to surface
2028
2029 //
2030 // Limit according to curvature and magnetic field inhomogenities
2031 // and improve stepsize estimation to reach plane
2032 //
2033 double fieldCurvLimit(limits.getLowestLimitSignedVal()); // signed
2034 std::pair<double, double> distVsStep (9.E99, 9.E99); // first: smallest straight line distances to plane; second: RK steps
2035
2036 static const unsigned int maxNumIt = 10;
2038
2040
2042 // if max iterations are reached, take a safe value
2043 // (in previous iteration, fieldCurvLimit has been not more than doubled)
2044 // and break.
2045 fieldCurvLimit *= 0.5;
2046 break;
2047 }
2048
2049 M1x7 state7_temp = { state7[0], state7[1], state7[2], state7[3], state7[4], state7[5], state7[6] }; // invalid: M1x7 state7_temp(state7);
2050 M1x3 SA;
2051
2054 std::cout << " maxStepArg = " << fieldCurvLimit << "; q = " << q << " \n";
2055 }
2056
2057 // remember steps and resulting SL distances to plane for stepsize improvement
2058 // calculate distance to surface
2059 Dist = SU[3] - (state7_temp[0] * SU[0] +
2060 state7_temp[1] * SU[1] +
2061 state7_temp[2] * SU[2]); // Distance between position and surface
2062
2063 An = state7_temp[3] * SU[0] +
2064 state7_temp[4] * SU[1] +
2065 state7_temp[5] * SU[2]; // An = dir * N; component of dir normal to surface
2066
2067 if (fabs(Dist/An) < fabs(distVsStep.first)) {
2068 distVsStep.first = Dist/An;
2069 distVsStep.second = fieldCurvLimit;
2070 }
2071
2072 // resize limit according to q never grow step size more than
2073 // two-fold to avoid infinite grow-shrink loops with strongly
2074 // inhomogeneous fields.
2075 if (q>2) {
2076 fieldCurvLimit *= 2;
2077 break;
2078 }
2079
2080 fieldCurvLimit *= q * 0.95;
2081
2082 if (fabs(q-1) < 0.25 || // good enough!
2083 fabs(fieldCurvLimit) > limits.getLowestLimitVal()) // other limits are lower!
2084 break;
2085 }
2088 else
2089 limits.setLimit(stp_fieldCurv, fieldCurvLimit);
2090
2092 if (fabs(distVsStep.first) < 8.E99) {
2093 stepToPlane = distVsStep.first + distVsStep.second;
2094 }
2095 limits.setLimit(stp_plane, stepToPlane);
2096
2097
2098 //
2099 // Select direction
2100 //
2101 // auto select
2104 std::cout << " auto select direction";
2105 if (!plane.isFinite()) std::cout << ", plane is not finite";
2106 std::cout << ".\n";
2107 }
2108 }
2109 // see if straight line approximation is ok
2112 std::cout << " straight line approximation is fine.\n";
2113 }
2114
2115 // if direction is pointing to active part of surface
2116 if( plane.isInActive(state7[0], state7[1], state7[2], state7[3], state7[4], state7[5]) ) {
2118 std::cout << " direction is pointing to active part of surface. \n";
2119 }
2120 }
2121 // if we are near the plane, but not pointing to the active area, make a big step!
2122 else {
2123 limits.removeLimit(stp_plane);
2124 limits.setStepSign(propDir_);
2126 std::cout << " we are near the plane, but not pointing to the active area. make a big step! \n";
2127 }
2128 }
2129 }
2130 // propDir_ is set and we are not pointing to an active part of a plane -> propDir_ decides!
2131 else {
2133 limits.removeLimit(stp_plane);
2134 limits.setStepSign(propDir_);
2136 std::cout << " invert Step according to propDir_ and make a big step. \n";
2137 }
2138 }
2139 }
2140
2141
2142 // call stepper and reduce stepsize if step not too small
2143 static const RKStep defaultRKStep;
2144 RKSteps_.push_back( defaultRKStep );
2145 std::vector<RKStep>::iterator lastStep = RKSteps_.end() - 1;
2148 ++RKStepsFXStop_;
2149 if ( true){
2150
2152 M1x7 state7_temp = { state7[0], state7[1], state7[2], state7[3], state7[4], state7[5], state7[6] }; // invalid: ... = (state7);
2153
2155 state7_temp,
2156 charge/state7[6], // |p|
2157 relMomLoss,
2158 pdgCode_,
2159 lastStep->matStep_.materialProperties_,
2160 limits,
2161 true);
2162 }
2163 else { //assume material has not changed
2165 lastStep->matStep_.materialProperties_ = (lastStep - 1)->matStep_.materialProperties_;
2166 }
2167 }
2168 }
2169
2171 std::cout << " final limits:\n";
2172 limits.Print();
2173 }
2174
2175 double finalStep = limits.getLowestLimitSignedVal();
2176
2177 lastStep->matStep_.stepSize_ = finalStep;
2178 lastStep->limits_ = limits;
2179
2181 std::cout << " --> Step to be used: " << finalStep << "\n";
2182 }
2183
2184 return finalStep;
2185
2186}
char propDir_
propagation direction (-1, 0, 1) -> (backward, auto, forward)
Definition AbsTrackRep.h:347
static MaterialEffects * getInstance()
Definition MaterialEffects.cc:72
void stepper(const RKTrackRep *rep, M1x7 &state7, const double &mom, double &relMomLoss, const int &pdg, MaterialProperties ¤tMaterial, StepLimits &limits, bool varField=true)
Returns maximum length so that a specified momentum loss will not be exceeded.
Definition MaterialEffects.cc:196
double RKPropagate(M1x7 &state7, M7x7 *jacobian, M1x3 &SA, double S, bool varField=true, bool calcOnlyLastRowOfJ=false) const
The actual Runge Kutta propagation.
Definition RKTrackRep.cc:1143
◆ estimateStep() [2/2]
|
private |
◆ Extrap() [1/2]
|
private |
Handles propagation and material effects.
extrapolateToPlane(), extrapolateToPoint() and extrapolateToLine() etc. call this function. Extrap() needs a plane as an argument, hence extrapolateToPoint() and extrapolateToLine() create virtual detector planes. In this function, RKutta() is called and the resulting points and point paths are filtered so that the direction doesn't change and tiny steps are filtered out. After the propagation the material effects are called via the MaterialEffects singleton. Extrap() will loop until the plane is reached, unless the propagation fails or the maximum number of iterations is exceeded.
fNoMaterial &&
Definition at line 2201 of file RKTrackRep.cc.
2210 {
2211
2212 static const unsigned int maxNumIt(500);
2213 unsigned int numIt(0);
2214
2215 double coveredDistance(0.);
2216 double dqop(0.);
2217
2218 const TVector3 W(destPlane.getNormal());
2219 M1x4 SU = {W.X(), W.Y(), W.Z(), destPlane.distance(0., 0., 0.)};
2220
2221 // make SU vector point away from origin
2222 if (W*destPlane.getO() < 0) {
2223 SU[0] *= -1;
2224 SU[1] *= -1;
2225 SU[2] *= -1;
2226 }
2227
2228
2229 M1x7 startState7;
2231
2232 while(true){
2233
2235 std::cout << "\n============ RKTrackRep::Extrap loop nr. " << numIt << " ============\n";
2236 std::cout << "Start plane: "; startPlane.Print();
2237 std::cout << "fillExtrapSteps " << fillExtrapSteps << "\n";
2238 }
2239
2240 if(++numIt > maxNumIt){
2241 Exception exc("RKTrackRep::Extrap ==> maximum number of iterations exceeded",__LINE__,__FILE__);
2242 exc.setFatal();
2243 throw exc;
2244 }
2245
2246 // initialize jacobianT with unit matrix
2249
2250 M7x7* noise = NULL;
2251 isAtBoundary = false;
2252
2253 // Remember previous state
2254 M1x7 oldState7;
2256
2257 // propagation
2258 bool checkJacProj = false;
2261
2263 coveredDistance, checkJacProj, noiseProjection_,
2264 limits_, onlyOneStep, !fillExtrapSteps) ) {
2265 Exception exc("RKTrackRep::Extrap ==> Runge Kutta propagation failed",__LINE__,__FILE__);
2266 exc.setFatal();
2267 throw exc;
2268 }
2269
2272 isAtBoundary = true;
2273
2274
2276 std::cout<<"RKSteps \n";
2278 std::cout << "stepSize = " << it->matStep_.stepSize_ << "\t";
2279 it->matStep_.materialProperties_.Print();
2280 }
2281 std::cout<<"\n";
2282 }
2283
2284
2285
2286 // call MatFX
2287 if(fillExtrapSteps) {
2288 noise = &noiseArray_;
2290 }
2291
2293 if ( nPoints>0){
2294 // momLoss has a sign - negative loss means momentum gain
2296 RKStepsFXStart_,
2297 RKStepsFXStop_,
2298 fabs(charge/state7[6]), // momentum
2299 pdgCode_,
2300 noise);
2301
2303
2305 std::cout << "momLoss: " << momLoss << " GeV; relative: " << momLoss/fabs(charge/state7[6]) << "\n";
2307 std::cout << "7D noise: \n";
2308 RKTools::printDim(*noise, 7, 7);
2309 }
2310 }
2311
2312 // do momLoss only for defined 1/momentum .ne.0
2313 if(fabs(state7[6])>1.E-10) {
2314 double qop = charge/(fabs(charge/state7[6])-momLoss);
2315 dqop = qop - state7[6];
2316 state7[6] = qop;
2317
2318 // correct state7 with dx/dqop, dy/dqop ... Greatly improves extrapolation accuracy!
2320 std::cout << "correct state7 with dx/dqop, dy/dqop ...\n";
2321 }
2324 }
2325 }
2326 } // finished MatFX
2327
2328
2329 // fill ExtrapSteps_
2330 if (fillExtrapSteps) {
2331 static const ExtrapStep defaultExtrapStep;
2332 ExtrapSteps_.push_back(defaultExtrapStep);
2333 std::vector<ExtrapStep>::iterator lastStep = ExtrapSteps_.end() - 1;
2334
2335 // Store Jacobian of this step for final calculation.
2337
2338 if( checkJacProj == true ){
2339 //project the noise onto the destPlane
2341
2343 std::cout << "7D noise projected onto plane: \n";
2345 }
2346 }
2347
2348 // Store this step's noise for final calculation.
2350
2352 std::cout<<"ExtrapSteps \n";
2353 for (std::vector<ExtrapStep>::iterator it = ExtrapSteps_.begin(); it != ExtrapSteps_.end(); ++it){
2356 }
2357 std::cout<<"\n";
2358 }
2359 }
2360
2361
2362
2363 // check if at boundary
2364 if (stopAtBoundary) {
2366 std::cout << "stopAtBoundary -> break; \n ";
2367 }
2368 break;
2369 }
2370
2371 if (onlyOneStep) {
2373 std::cout << "onlyOneStep -> break; \n ";
2374 }
2375 break;
2376 }
2377
2378 //break if we arrived at destPlane
2379 if(atPlane) {
2381 std::cout << "arrived at destPlane with a distance of " << destPlane.distance(state7[0], state7[1], state7[2]) << " cm left. ";
2382 if (destPlane.isInActive(state7[0], state7[1], state7[2], state7[3], state7[4], state7[5]))
2383 std::cout << "In active area of destPlane. \n";
2384 else
2385 std::cout << "NOT in active area of plane. \n";
2386 }
2387 break;
2388 }
2389
2390 }
2391
2392 if (fillExtrapSteps) {
2393 // propagate cov and add noise
2394 calcForwardJacobianAndNoise(startState7, startPlane, state7, destPlane);
2395
2396 if (cov != NULL) {
2397 cov->Similarity(fJacobian_);
2398 *cov += fNoise_;
2399 }
2400
2402 if (cov != NULL) {
2403 std::cout << "final covariance matrix after Extrap: "; cov->Print();
2404 }
2405 }
2406 }
2407
2408 return coveredDistance;
2409}
double effects(const std::vector< genfit::RKStep > &steps, int materialsFXStart, int materialsFXStop, const double &mom, const int &pdg, M7x7 *noise=nullptr)
Calculates energy loss in the traveled path, optional calculation of noise matrix.
void calcForwardJacobianAndNoise(const M1x7 &startState7, const DetPlane &startPlane, const M1x7 &destState7, const DetPlane &destPlane) const
Definition RKTrackRep.cc:774
bool RKutta(const M1x4 &SU, const DetPlane &plane, double charge, M1x7 &state7, M7x7 *jacobianT, double &coveredDistance, bool &checkJacProj, M7x7 &noiseProjection, StepLimits &limits, bool onlyOneStep=false, bool calcOnlyLastRowOfJ=false) const
Propagates the particle through the magnetic field.
Definition RKTrackRep.cc:1705
void setLimit(StepLimitType type, double value)
absolute of value will be taken! If limit is already lower, it will be set to value anyway.
Definition StepLimits.h:89
std::pair< StepLimitType, double > getLowestLimit(double margin=1.E-3) const
Get the lowest limit.
Definition StepLimits.cc:44
void printDim(const double *mat, unsigned int dimX, unsigned int dimY)
Definition RKTools.cc:691
◆ Extrap() [2/2]
|
private |
Handles propagation and material effects.
extrapolateToPlane(), extrapolateToPoint() and extrapolateToLine() etc. call this function. Extrap() needs a plane as an argument, hence extrapolateToPoint() and extrapolateToLine() create virtual detector planes. In this function, RKutta() is called and the resulting points and point paths are filtered so that the direction doesn't change and tiny steps are filtered out. After the propagation the material effects are called via the MaterialEffects singleton. Extrap() will loop until the plane is reached, unless the propagation fails or the maximum number of iterations is exceeded.
◆ extrapolateBy() [1/2]
|
virtual |
Extrapolates the state by step (cm) and returns the extrapolation length and, via reference, the extrapolated state.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Implements genfit::AbsTrackRep.
Definition at line 579 of file RKTrackRep.cc.
582 {
583
585 std::cout << "RKTrackRep::extrapolateBy()\n";
586 }
587
588 checkCache(state, NULL);
589
590 static const unsigned int maxIt(1000);
591
592 // to 7D
593 M1x7 state7;
594 getState7(state, state7);
595
596 bool fillExtrapSteps(false);
597 if (dynamic_cast<MeasuredStateOnPlane*>(&state) != NULL) {
598 fillExtrapSteps = true;
599 }
600 else if (calcJacobianNoise)
601 fillExtrapSteps = true;
602
603 double tracklength(0.);
604
606
607 bool isAtBoundary(false);
608
609 DetPlane startPlane(*(state.getPlane()));
611 unsigned int iterations(0);
612
613 while(true){
614 if(++iterations == maxIt) {
615 Exception exc("RKTrackRep::extrapolateBy ==> maximum number of iterations reached",__LINE__,__FILE__);
616 exc.setFatal();
617 throw exc;
618 }
619
620 pos.SetXYZ(state7[0], state7[1], state7[2]);
621 dir.SetXYZ(state7[3], state7[4], state7[5]);
622
624
625 plane->setON(dest, dir);
626
627 tracklength += this->Extrap(startPlane, *plane, getCharge(state), isAtBoundary, state7, false, NULL, true, stopAtBoundary, (step-tracklength));
628
629 // check break conditions
630 if (stopAtBoundary && isAtBoundary) {
631 pos.SetXYZ(state7[0], state7[1], state7[2]);
632 dir.SetXYZ(state7[3], state7[4], state7[5]);
633 plane->setON(pos, dir);
634 break;
635 }
636
639 std::cout << "RKTrackRep::extrapolateBy(): reached after " << iterations << " iterations. \n";
640 }
641 pos.SetXYZ(state7[0], state7[1], state7[2]);
642 dir.SetXYZ(state7[3], state7[4], state7[5]);
643 plane->setON(pos, dir);
644 break;
645 }
646
647 startPlane = *plane;
648
649 }
650
651 if (fillExtrapSteps) { // now do the full extrapolation with covariance matrix
652 // make use of the cache
655
657 }
658 else {
659 state.setPlane(plane);
660 getState5(state, state7);
661 }
662
663 lastEndState_ = state;
664
665 return tracklength;
666}
virtual double getCharge(const StateOnPlane &state) const
Get the (fitted) charge of a state. This is not always equal the pdg charge (e.g. if the charge sign ...
Definition RKTrackRep.cc:712
double Extrap(const DetPlane &startPlane, const DetPlane &destPlane, double charge, bool &isAtBoundary, M1x7 &state7, bool fillExtrapSteps, TMatrixDSym *cov=nullptr, bool onlyOneStep=false, bool stopAtBoundary=false, double maxStep=1.E99) const
Handles propagation and material effects.
Definition RKTrackRep.cc:2201
void getState7(const StateOnPlane &state, M1x7 &state7) const
Definition RKTrackRep.cc:1351
void checkCache(const StateOnPlane &state, const SharedPlanePtr *plane) const
Definition RKTrackRep.cc:2412
void getState5(StateOnPlane &state, const M1x7 &state7) const
Definition RKTrackRep.cc:1385
virtual double extrapolateToPlane(StateOnPlane &state, const SharedPlanePtr &plane, bool stopAtBoundary=false, bool calcJacobianNoise=false) const
Extrapolates the state to plane, and returns the extrapolation length and, via reference,...
Definition RKTrackRep.cc:97
void setPlane(const SharedPlanePtr &plane)
Definition StateOnPlane.h:67
boost::shared_ptr< genfit::DetPlane > SharedPlanePtr
Shared Pointer to a DetPlane.
Definition SharedPlanePtr.h:43
◆ extrapolateBy() [2/2]
|
virtual |
Extrapolates the state by step (cm) and returns the extrapolation length and, via reference, the extrapolated state.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Implements genfit::AbsTrackRep.
◆ extrapolateToCylinder() [1/2]
|
virtual |
Extrapolates the state to the cylinder surface, and returns the extrapolation length and, via reference, the extrapolated state.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Implements genfit::AbsTrackRep.
Definition at line 353 of file RKTrackRep.cc.
358 {
359
361 std::cout << "RKTrackRep::extrapolateToCylinder()\n";
362 }
363
364 checkCache(state, NULL);
365
366 static const unsigned int maxIt(1000);
367
368 // to 7D
369 M1x7 state7;
370 getState7(state, state7);
371
372 bool fillExtrapSteps(false);
373 if (dynamic_cast<MeasuredStateOnPlane*>(&state) != NULL) {
374 fillExtrapSteps = true;
375 }
376 else if (calcJacobianNoise)
377 fillExtrapSteps = true;
378
379 double tracklength(0.), maxStep(1.E99);
380
382
383 bool isAtBoundary(false);
384
385 DetPlane startPlane(*(state.getPlane()));
387 unsigned int iterations(0);
388
389 while(true){
390 if(++iterations == maxIt) {
391 Exception exc("RKTrackRep::extrapolateToCylinder ==> maximum number of iterations reached",__LINE__,__FILE__);
392 exc.setFatal();
393 throw exc;
394 }
395
396 pos.SetXYZ(state7[0], state7[1], state7[2]);
397 dir.SetXYZ(state7[3], state7[4], state7[5]);
398
399 // solve quadratic equation
400 TVector3 AO = (pos - linePoint);
401 TVector3 AOxAB = (AO.Cross(lineDirection));
402 TVector3 VxAB = (dir.Cross(lineDirection));
403 float ab2 = (lineDirection * lineDirection);
408 if(arg < 0) {
409 Exception exc("RKTrackRep::extrapolateToCylinder ==> cannot solve",__LINE__,__FILE__);
410 exc.setFatal();
411 throw exc;
412 }
413 double term = sqrt(arg);
414 double k1, k2;
415 if (b<0) {
416 k1 = (-b + term)/(2.*a);
418 }
419 else {
421 k2 = (-b - term)/(2.*a);
422 }
423
424 // select smallest absolute solution -> closest cylinder surface
426 if (fabs(k2)<fabs(k))
427 k = k2;
428
431 }
432
434
435 plane->setO(dest);
436 plane->setUV((dest-linePoint).Cross(lineDirection), lineDirection);
437
438 tracklength += this->Extrap(startPlane, *plane, getCharge(state), isAtBoundary, state7, false, NULL, true, stopAtBoundary, maxStep);
439
440 // check break conditions
441 if (stopAtBoundary && isAtBoundary) {
442 pos.SetXYZ(state7[0], state7[1], state7[2]);
443 dir.SetXYZ(state7[3], state7[4], state7[5]);
444 plane->setO(pos);
445 plane->setUV((pos-linePoint).Cross(lineDirection), lineDirection);
446 break;
447 }
448
450
451 startPlane = *plane;
452
453 }
454
455 if (fillExtrapSteps) { // now do the full extrapolation with covariance matrix
456 // make use of the cache
459
461 }
462 else {
463 state.setPlane(plane);
464 getState5(state, state7);
465 }
466
467 lastEndState_ = state;
468
469 return tracklength;
470}
◆ extrapolateToCylinder() [2/2]
|
virtual |
Extrapolates the state to the cylinder surface, and returns the extrapolation length and, via reference, the extrapolated state.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Implements genfit::AbsTrackRep.
◆ extrapolateToLine() [1/6]
|
virtual |
Extrapolates the state to the POCA to a line, and returns the extrapolation length and, via reference, the extrapolated state.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Implements genfit::AbsTrackRep.
Definition at line 150 of file RKTrackRep.cc.
154 {
155
157 std::cout << "RKTrackRep::extrapolateToLine()\n";
158 }
159
160 checkCache(state, NULL);
161
162 static const unsigned int maxIt(1000);
163
164 // to 7D
165 M1x7 state7;
166 getState7(state, state7);
167
168 bool fillExtrapSteps(false);
169 if (dynamic_cast<MeasuredStateOnPlane*>(&state) != NULL) {
170 fillExtrapSteps = true;
171 }
172 else if (calcJacobianNoise)
173 fillExtrapSteps = true;
174
177 TVector3 dir(state7[3], state7[4], state7[5]);
178 TVector3 lastDir(0,0,0);
179 TVector3 poca, poca_onwire;
180 bool isAtBoundary(false);
181
182 DetPlane startPlane(*(state.getPlane()));
184 unsigned int iterations(0);
185
186 while(true){
187 if(++iterations == maxIt) {
188 Exception exc("RKTrackRep::extrapolateToLine ==> extrapolation to line failed, maximum number of iterations reached",__LINE__,__FILE__);
189 exc.setFatal();
190 throw exc;
191 }
192
193 lastStep = step;
194 lastDir = dir;
195
196 step = this->Extrap(startPlane, *plane, charge, isAtBoundary, state7, false, NULL, true, stopAtBoundary, maxStep);
197 tracklength += step;
198
199 dir.SetXYZ(state7[3], state7[4], state7[5]);
200 poca.SetXYZ(state7[0], state7[1], state7[2]);
201 poca_onwire = pocaOnLine(linePoint, lineDirection, poca);
202
203 // check break conditions
204 if (stopAtBoundary && isAtBoundary) {
205 plane->setON(dir, poca);
206 break;
207 }
208
209 angle = fabs(dir.Angle((poca_onwire-poca))-TMath::PiOver2()); // angle between direction and connection to point - 90 deg
210 distToPoca = (poca_onwire-poca).Mag();
212
213 // if lastStep and step have opposite sign, the real normal vector lies somewhere between the last two normal vectors (i.e. the directions)
214 // -> try mean value of the two (normalization not needed)
215 if (lastStep*step < 0){
216 dir += lastDir;
217 maxStep = 0.5*fabs(lastStep); // make it converge!
218 }
219
220 startPlane = *plane;
221 plane->setU(dir.Cross(lineDirection));
222 }
223
224 if (fillExtrapSteps) { // now do the full extrapolation with covariance matrix
225 // make use of the cache
228
230 }
231 else {
232 state.setPlane(plane);
233 getState5(state, state7);
234 }
235
237 std::cout << "RKTrackRep::extrapolateToLine(): Reached POCA after " << iterations+1 << " iterations. Distance: " << (poca_onwire-poca).Mag() << " cm. Angle deviation: " << dir.Angle((poca_onwire-poca))-TMath::PiOver2() << " rad \n";
238 }
239
240 lastEndState_ = state;
241
242 return tracklength;
243}
TVector3 pocaOnLine(const TVector3 &linePoint, const TVector3 &lineDirection, const TVector3 &point) const
Definition RKTrackRep.cc:2189
◆ extrapolateToLine() [2/6]
|
virtual |
Extrapolates the state to the POCA to a line, and returns the extrapolation length and, via reference, the extrapolated state.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Implements genfit::AbsTrackRep.
◆ extrapolateToLine() [3/6]
|
virtual |
Extrapolates the state to the POCA to a line, and returns the extrapolation length and, via reference, the extrapolated state.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Implements genfit::AbsTrackRep.
◆ extrapolateToLine() [4/6]
|
virtual |
Extrapolates the state to the POCA to a line, and returns the extrapolation length and, via reference, the extrapolated state.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Implements genfit::AbsTrackRep.
◆ extrapolateToLine() [5/6]
|
inlinevirtual |
Resembles the interface of GFAbsTrackRep in old versions of genfit.
This interface to extrapolateToLine is intended to resemble the interface of GFAbsTrackRep in old versions of genfit and is implemented by default via the preceding function.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Reimplemented from genfit::AbsTrackRep.
Definition at line 120 of file AbsTrackRep.h.
127 {
128 TVector3 wireDir(point2 - point1);
129 wireDir.Unit();
130 double retval = this->extrapolateToLine(state, point1, wireDir, stopAtBoundary, calcJacobianNoise);
131 poca = this->getPos(state);
132 dirInPoca = this->getMom(state);
133 dirInPoca.Unit();
134
135 poca_onwire = point1 + wireDir*((poca - point1)*wireDir);
136
137 return retval;
138 }
virtual TVector3 getMom(const StateOnPlane &state) const
Get the cartesian momentum vector of a state.
Definition RKTrackRep.cc:677
virtual double extrapolateToLine(StateOnPlane &state, const TVector3 &linePoint, const TVector3 &lineDirection, bool stopAtBoundary=false, bool calcJacobianNoise=false) const
Extrapolates the state to the POCA to a line, and returns the extrapolation length and,...
Definition RKTrackRep.cc:150
◆ extrapolateToLine() [6/6]
|
inlinevirtual |
Resembles the interface of GFAbsTrackRep in old versions of genfit.
This interface to extrapolateToLine is intended to resemble the interface of GFAbsTrackRep in old versions of genfit and is implemented by default via the preceding function.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Reimplemented from genfit::AbsTrackRep.
Definition at line 120 of file AbsTrackRep.h.
127 {
128 TVector3 wireDir(point2 - point1);
129 wireDir.Unit();
130 double retval = this->extrapolateToLine(state, point1, wireDir, stopAtBoundary, calcJacobianNoise);
131 poca = this->getPos(state);
132 dirInPoca = this->getMom(state);
133 dirInPoca.Unit();
134
135 poca_onwire = point1 + wireDir*((poca - point1)*wireDir);
136
137 return retval;
138 }
◆ extrapolateToPlane() [1/3]
|
virtual |
Extrapolates the state to plane, and returns the extrapolation length and, via reference, the extrapolated state.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Implements genfit::AbsTrackRep.
Definition at line 97 of file RKTrackRep.cc.
100 {
101
103 std::cout << "RKTrackRep::extrapolateToPlane()\n";
104 }
105
106
107 if (state.getPlane() == plane) {
109 std::cout << "state is already defined at plane. Do nothing! \n";
110 }
111 return 0;
112 }
113
114 checkCache(state, &plane);
115
116 // to 7D
117 M1x7 state7;
118 getState7(state, state7);
119
120 TMatrixDSym* covPtr(NULL);
121 bool fillExtrapSteps(false);
122 if (dynamic_cast<MeasuredStateOnPlane*>(&state) != NULL) {
123 covPtr = &(static_cast<MeasuredStateOnPlane*>(&state)->getCov());
124 fillExtrapSteps = true;
125 }
126 else if (calcJacobianNoise)
127 fillExtrapSteps = true;
128
129 // actual extrapolation
130 bool isAtBoundary(false);
131 double coveredDistance = Extrap(*(state.getPlane()), *plane, getCharge(state), isAtBoundary, state7, fillExtrapSteps, covPtr, false, stopAtBoundary);
132
133 if (stopAtBoundary && isAtBoundary) {
135 TVector3(state7[3], state7[4], state7[5]))));
136 }
137 else {
138 state.setPlane(plane);
139 }
140
141 // back to 5D
142 getState5(state, state7);
143
144 lastEndState_ = state;
145
146 return coveredDistance;
147}
◆ extrapolateToPlane() [2/3]
|
virtual |
Extrapolates the state to plane, and returns the extrapolation length and, via reference, the extrapolated state.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Implements genfit::AbsTrackRep.
◆ extrapolateToPlane() [3/3]
|
virtual |
Definition at line 84 of file RKTrackRep.cc.
87 {
88 // std::cout << "test "<<point.X()<<std::endl;
91 shared,
92 stopAtBoundary,
93 calcJacobianNoise);
94}
◆ extrapolateToPoint() [1/4]
|
inlinevirtual |
Extrapolates the state to the POCA to a point, and returns the extrapolation length and, via reference, the extrapolated state.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Implements genfit::AbsTrackRep.
Definition at line 95 of file RKTrackRep.h.
98 {
100 }
virtual double extrapToPoint(StateOnPlane &state, const TVector3 &point, const TMatrixDSym *G=NULL, bool stopAtBoundary=false, bool calcJacobianNoise=false) const
Definition RKTrackRep.cc:246
◆ extrapolateToPoint() [2/4]
|
inlinevirtual |
Extrapolates the state to the POCA to a point, and returns the extrapolation length and, via reference, the extrapolated state.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Implements genfit::AbsTrackRep.
Definition at line 100 of file RKTrackRep.h.
◆ extrapolateToPoint() [3/4]
|
inlinevirtual |
Extrapolates the state to the POCA to a point in the metric of G, and returns the extrapolation length and, via reference, the extrapolated state.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Implements genfit::AbsTrackRep.
Definition at line 102 of file RKTrackRep.h.
◆ extrapolateToPoint() [4/4]
|
inlinevirtual |
Extrapolates the state to the POCA to a point in the metric of G, and returns the extrapolation length and, via reference, the extrapolated state.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Implements genfit::AbsTrackRep.
Definition at line 107 of file RKTrackRep.h.
◆ extrapolateToSphere() [1/2]
|
virtual |
Extrapolates the state to the sphere surface, and returns the extrapolation length and, via reference, the extrapolated state.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Implements genfit::AbsTrackRep.
Definition at line 473 of file RKTrackRep.cc.
477 {
478
480 std::cout << "RKTrackRep::extrapolateToSphere()\n";
481 }
482
483 checkCache(state, NULL);
484
485 static const unsigned int maxIt(1000);
486
487 // to 7D
488 M1x7 state7;
489 getState7(state, state7);
490
491 bool fillExtrapSteps(false);
492 if (dynamic_cast<MeasuredStateOnPlane*>(&state) != NULL) {
493 fillExtrapSteps = true;
494 }
495 else if (calcJacobianNoise)
496 fillExtrapSteps = true;
497
498 double tracklength(0.), maxStep(1.E99);
499
501
502 bool isAtBoundary(false);
503
504 DetPlane startPlane(*(state.getPlane()));
506 unsigned int iterations(0);
507
508 while(true){
509 if(++iterations == maxIt) {
510 Exception exc("RKTrackRep::extrapolateToSphere ==> maximum number of iterations reached",__LINE__,__FILE__);
511 exc.setFatal();
512 throw exc;
513 }
514
515 pos.SetXYZ(state7[0], state7[1], state7[2]);
516 dir.SetXYZ(state7[3], state7[4], state7[5]);
517
518 // solve quadratic equation
519 TVector3 AO = (pos - point);
520 double dirAO = dir * AO;
521 double arg = dirAO*dirAO - AO*AO + radius*radius;
522 if(arg < 0) {
523 Exception exc("RKTrackRep::extrapolateToSphere ==> cannot solve",__LINE__,__FILE__);
524 exc.setFatal();
525 throw exc;
526 }
527 double term = sqrt(arg);
528 double k1, k2;
529 k1 = -dirAO + term;
530 k2 = -dirAO - term;
531
532 // select smallest absolute solution -> closest cylinder surface
534 if (fabs(k2)<fabs(k))
535 k = k2;
536
539 }
540
542
543 plane->setON(dest, dest-point);
544
545 tracklength += this->Extrap(startPlane, *plane, getCharge(state), isAtBoundary, state7, false, NULL, true, stopAtBoundary, maxStep);
546
547 // check break conditions
548 if (stopAtBoundary && isAtBoundary) {
549 pos.SetXYZ(state7[0], state7[1], state7[2]);
550 dir.SetXYZ(state7[3], state7[4], state7[5]);
551 plane->setON(pos, pos-point);
552 break;
553 }
554
556
557 startPlane = *plane;
558
559 }
560
561 if (fillExtrapSteps) { // now do the full extrapolation with covariance matrix
562 // make use of the cache
565
567 }
568 else {
569 state.setPlane(plane);
570 getState5(state, state7);
571 }
572
573 lastEndState_ = state;
574
575 return tracklength;
576}
◆ extrapolateToSphere() [2/2]
|
virtual |
Extrapolates the state to the sphere surface, and returns the extrapolation length and, via reference, the extrapolated state.
If stopAtBoundary is true, the extrapolation stops as soon as a material boundary is encountered.
If state has a covariance, jacobian and noise matrices will be calculated and the covariance will be propagated. If state has no covariance, jacobian and noise will only be calculated if calcJacobianNoise == true.
Implements genfit::AbsTrackRep.
◆ extrapToPoint() [1/2]
|
privatevirtual |
Definition at line 246 of file RKTrackRep.cc.
250 {
251
253 std::cout << "RKTrackRep::extrapolateToPoint()\n";
254 }
255
256 checkCache(state, NULL);
257
258 static const unsigned int maxIt(1000);
259
260 // to 7D
261 M1x7 state7;
262 getState7(state, state7);
263
264 bool fillExtrapSteps(false);
265 if (dynamic_cast<MeasuredStateOnPlane*>(&state) != NULL) {
266 fillExtrapSteps = true;
267 }
268 else if (calcJacobianNoise)
269 fillExtrapSteps = true;
270
272 TVector3 dir(state7[3], state7[4], state7[5]);
273 if (G != NULL) {
275 Exception exc("RKTrackRep::extrapolateToLine ==> G is not 3x3",__LINE__,__FILE__);
276 exc.setFatal();
277 throw exc;
278 }
279 dir = TMatrix(*G) * dir;
280 }
281 TVector3 lastDir(0,0,0);
282
283 TVector3 poca;
284 bool isAtBoundary(false);
285
286 DetPlane startPlane(*(state.getPlane()));
288 unsigned int iterations(0);
289
290 while(true){
291 if(++iterations == maxIt) {
292 Exception exc("RKTrackRep::extrapolateToPoint ==> extrapolation to point failed, maximum number of iterations reached",__LINE__,__FILE__);
293 exc.setFatal();
294 throw exc;
295 }
296
297 lastStep = step;
298 lastDir = dir;
299
300 step = this->Extrap(startPlane, *plane, getCharge(state), isAtBoundary, state7, false, NULL, true, stopAtBoundary, maxStep);
301 tracklength += step;
302
303 dir.SetXYZ(state7[3], state7[4], state7[5]);
304 if (G != NULL) {
305 dir = TMatrix(*G) * dir;
306 }
307 poca.SetXYZ(state7[0], state7[1], state7[2]);
308
309 // check break conditions
310 if (stopAtBoundary && isAtBoundary) {
311 plane->setON(dir, poca);
312 break;
313 }
314
315 angle = fabs(dir.Angle((point-poca))-TMath::PiOver2()); // angle between direction and connection to point - 90 deg
316 distToPoca = (point-poca).Mag();
318
319 // if lastStep and step have opposite sign, the real normal vector lies somewhere between the last two normal vectors (i.e. the directions)
320 // -> try mean value of the two (normalization not needed)
321 if (lastStep*step < 0){
322 dir += lastDir;
323 maxStep = 0.5*fabs(lastStep); // make it converge!
324 }
325
326 startPlane = *plane;
327 plane->setNormal(dir);
328 }
329
330 if (fillExtrapSteps) { // now do the full extrapolation with covariance matrix
331 // make use of the cache
334
336 }
337 else {
338 state.setPlane(plane);
339 getState5(state, state7);
340 }
341
342
344 std::cout << "RKTrackRep::extrapolateToPoint(): Reached POCA after " << iterations+1 << " iterations. Distance: " << (point-poca).Mag() << " cm. Angle deviation: " << dir.Angle((point-poca))-TMath::PiOver2() << " rad \n";
345 }
346
347 lastEndState_ = state;
348
349 return tracklength;
350}
◆ extrapToPoint() [2/2]
|
privatevirtual |
◆ get6DCov() [1/2]
|
virtual |
Get the 6D covariance.
Implements genfit::AbsTrackRep.
Definition at line 704 of file RKTrackRep.cc.
704 {
705 TMatrixDSym cov(6);
707
708 return cov;
709}
void transformPM6(const MeasuredStateOnPlane &state, M6x6 &out6x6) const
Definition RKTrackRep.cc:1497
◆ get6DCov() [2/2]
|
virtual |
Get the 6D covariance.
Implements genfit::AbsTrackRep.
◆ getBackwardJacobianAndNoise() [1/2]
|
virtual |
Get the jacobian and noise matrix of the last extrapolation if it would have been done in opposite direction.
Implements genfit::AbsTrackRep.
Definition at line 844 of file RKTrackRep.cc.
844 {
845
847 std::cout << "RKTrackRep::getBackwardJacobianAndNoise " << std::endl;
848 }
849
851 Exception exc("RKTrackRep::getBackwardJacobianAndNoise ==> cache is empty. Extrapolation must run with a MeasuredStateOnPlane.",__LINE__,__FILE__);
852 throw exc;
853 }
854
855 jacobian.ResizeTo(5,5);
856 jacobian = fJacobian_;
857 if (!useInvertFast) {
858 TDecompLU invertAlgo(jacobian);
860 if(status == 0){
861 Exception e("cannot invert matrix, status = 0", __LINE__,__FILE__);
862 e.setFatal();
863 throw e;
864 }
865 } else {
866 double det;
867 jacobian.InvertFast(&det);
868 if(det < 1e-80){
869 Exception e("cannot invert matrix, status = 0", __LINE__,__FILE__);
870 e.setFatal();
871 throw e;
872 }
873 }
874
875 noise.ResizeTo(5,5);
876 noise = fNoise_;
877 noise.Similarity(jacobian);
878
879 // lastStartState_ = jacobian * lastEndState_ + deltaState
880 deltaState.ResizeTo(5);
882}
◆ getBackwardJacobianAndNoise() [2/2]
|
virtual |
Get the jacobian and noise matrix of the last extrapolation if it would have been done in opposite direction.
Implements genfit::AbsTrackRep.
◆ getCharge() [1/2]
|
virtual |
Get the (fitted) charge of a state. This is not always equal the pdg charge (e.g. if the charge sign was flipped during the fit).
Implements genfit::AbsTrackRep.
Definition at line 712 of file RKTrackRep.cc.
712 {
713
714 if (dynamic_cast<const MeasurementOnPlane*>(&state) != NULL) {
715 Exception exc("RKTrackRep::getCharge - cannot get charge from MeasurementOnPlane",__LINE__,__FILE__);
716 exc.setFatal();
717 throw exc;
718 }
719
721
722 // return pdgCharge with sign of q/p
723 if (state.getState()(0) * pdgCharge < 0)
724 return -pdgCharge;
725 else
726 return pdgCharge;
727}
double getPDGCharge() const
Get the charge of the particle of the pdg code.
Definition AbsTrackRep.cc:93
◆ getCharge() [2/2]
|
virtual |
Get the (fitted) charge of a state. This is not always equal the pdg charge (e.g. if the charge sign was flipped during the fit).
Implements genfit::AbsTrackRep.
◆ getDim() [1/2]
|
inlinevirtual |
Get the dimension of the state vector used by the track representation.
Implements genfit::AbsTrackRep.
Definition at line 129 of file RKTrackRep.h.
129{return 5;}
◆ getDim() [2/2]
|
inlinevirtual |
Get the dimension of the state vector used by the track representation.
Implements genfit::AbsTrackRep.
Definition at line 134 of file RKTrackRep.h.
134{return 5;}
◆ getForwardJacobianAndNoise() [1/2]
|
virtual |
Get the jacobian and noise matrix of the last extrapolation.
Implements genfit::AbsTrackRep.
Definition at line 820 of file RKTrackRep.cc.
820 {
821
822 jacobian.ResizeTo(5,5);
823 jacobian = fJacobian_;
824
825 noise.ResizeTo(5,5);
826 noise = fNoise_;
827
828 // lastEndState_ = jacobian * lastStartState_ + deltaState
829 deltaState.ResizeTo(5);
830 // Calculate this without temporaries:
831 //deltaState = lastEndState_.getState() - jacobian * lastStartState_.getState()
833 deltaState *= jacobian;
835 deltaState *= -1;
836
837
839 std::cout << "delta state : "; deltaState.Print();
840 }
841}
◆ getForwardJacobianAndNoise() [2/2]
|
virtual |
Get the jacobian and noise matrix of the last extrapolation.
Implements genfit::AbsTrackRep.
◆ getMom() [1/2]
|
virtual |
Get the cartesian momentum vector of a state.
Implements genfit::AbsTrackRep.
Definition at line 677 of file RKTrackRep.cc.
◆ getMom() [2/2]
|
virtual |
Get the cartesian momentum vector of a state.
Implements genfit::AbsTrackRep.
◆ getMomMag() [1/2]
|
virtual |
get the magnitude of the momentum in GeV.
Implements genfit::AbsTrackRep.
Definition at line 730 of file RKTrackRep.cc.
730 {
731 // p = q / qop
733 assert (p>=0);
735}
◆ getMomMag() [2/2]
|
virtual |
get the magnitude of the momentum in GeV.
Implements genfit::AbsTrackRep.
◆ getMomVar() [1/2]
|
virtual |
get the variance of the absolute value of the momentum .
Implements genfit::AbsTrackRep.
Definition at line 738 of file RKTrackRep.cc.
738 {
739
740 if (dynamic_cast<const MeasurementOnPlane*>(&state) != NULL) {
741 Exception exc("RKTrackRep::getMomVar - cannot get momVar from MeasurementOnPlane",__LINE__,__FILE__);
742 exc.setFatal();
743 throw exc;
744 }
745
746 // p(qop) = q/qop
747 // dp/d(qop) = - q / (qop^2)
748 // (delta p) = (delta qop) * |dp/d(qop)| = delta qop * |q / (qop^2)|
749 // (var p) = (var qop) * q^2 / (qop^4)
750
751 // delta means sigma
752 // cov(0,0) is sigma^2
753
755}
◆ getMomVar() [2/2]
|
virtual |
get the variance of the absolute value of the momentum .
Implements genfit::AbsTrackRep.
◆ getPos() [1/2]
|
virtual |
Get the cartesian position of a state.
Implements genfit::AbsTrackRep.
Definition at line 669 of file RKTrackRep.cc.
◆ getPos() [2/2]
|
virtual |
Get the cartesian position of a state.
Implements genfit::AbsTrackRep.
◆ getPosMom() [1/2]
|
virtual |
Get cartesian position and momentum vector of a state.
Implements genfit::AbsTrackRep.
Definition at line 687 of file RKTrackRep.cc.
◆ getPosMom() [2/2]
|
virtual |
Get cartesian position and momentum vector of a state.
Implements genfit::AbsTrackRep.
◆ getPosMomCov() [1/2]
|
virtual |
Translates MeasuredStateOnPlane into 3D position, momentum and 6x6 covariance.
Implements genfit::AbsTrackRep.
Definition at line 697 of file RKTrackRep.cc.
697 {
698 getPosMom(state, pos, mom);
699 cov.ResizeTo(6,6);
701}
virtual void getPosMom(const StateOnPlane &state, TVector3 &pos, TVector3 &mom) const
Get cartesian position and momentum vector of a state.
Definition RKTrackRep.cc:687
◆ getPosMomCov() [2/2]
|
virtual |
Translates MeasuredStateOnPlane into 3D position, momentum and 6x6 covariance.
Implements genfit::AbsTrackRep.
◆ getQop() [1/2]
|
inlinevirtual |
Get charge over momentum.
Implements genfit::AbsTrackRep.
Definition at line 142 of file RKTrackRep.h.
142{return state.getState()(0);}
◆ getQop() [2/2]
|
inlinevirtual |
Get charge over momentum.
Implements genfit::AbsTrackRep.
Definition at line 147 of file RKTrackRep.h.
147{return state.getState()(0);}
◆ getRadiationLenght() [1/2]
|
virtual |
Get the accumulated X/X0 (path / radiation length) of the material crossed in the last extrapolation.
Implements genfit::AbsTrackRep.
Definition at line 905 of file RKTrackRep.cc.
905 {
906
907 // Todo: test
908
910 Exception exc("RKTrackRep::getRadiationLenght ==> cache is empty.",__LINE__,__FILE__);
911 throw exc;
912 }
913
914 double radLen(0);
915
917 radLen += RKSteps_.at(i).matStep_.stepSize_ / RKSteps_.at(i).matStep_.materialProperties_.getRadLen();
918 }
919
920 return radLen;
921}
◆ getRadiationLenght() [2/2]
|
virtual |
Get the accumulated X/X0 (path / radiation length) of the material crossed in the last extrapolation.
Implements genfit::AbsTrackRep.
◆ getSpu() [1/2]
| double genfit::RKTrackRep::getSpu | ( | const StateOnPlane & | state | ) | const |
Definition at line 758 of file RKTrackRep.cc.
758 {
759
760 if (dynamic_cast<const MeasurementOnPlane*>(&state) != NULL) {
761 Exception exc("RKTrackRep::getSpu - cannot get spu from MeasurementOnPlane",__LINE__,__FILE__);
762 exc.setFatal();
763 throw exc;
764 }
765
766 const TVectorD& auxInfo = state.getAuxInfo();
767 if (auxInfo.GetNrows() == 1)
768 return state.getAuxInfo()(0);
769 else
770 return 1.;
771}
◆ getSpu() [2/2]
| double genfit::RKTrackRep::getSpu | ( | const StateOnPlane & | state | ) | const |
◆ getState5() [1/2]
|
private |
Definition at line 1385 of file RKTrackRep.cc.
1385 {
1386
1387 // state5: (q/p, u', v'. u, v)
1388
1389 double spu(1.);
1390
1391 const TVector3& O(state.getPlane()->getO());
1392 const TVector3& U(state.getPlane()->getU());
1394 const TVector3& W(state.getPlane()->getNormal());
1395
1396 TVector3 posShift(state7[0], state7[1], state7[2]);
1397 posShift -= state.getPlane()->getO();
1398
1399 // force A to be in normal direction and set spu accordingly
1400 double AtW = state7[3]*W.X() + state7[4]*W.Y() + state7[5]*W.Z();
1401 if (AtW < 0) {
1402 //fDir *= -1.;
1403 //AtW *= -1.;
1404 spu = -1.;
1405 }
1406
1407 TVectorD& state5 = state.getState();
1408
1409 state5(0) = state7[6]; // q/p
1410 state5(1) = (state7[3]*U.X() + state7[4]*U.Y() + state7[5]*U.Z()) / AtW; // u' = (A * U) / (A * W)
1411 state5(2) = (state7[3]*V.X() + state7[4]*V.Y() + state7[5]*V.Z()) / AtW; // v' = (A * V) / (A * W)
1412 state5(3) = ((state7[0]-O.X())*U.X() +
1413 (state7[1]-O.Y())*U.Y() +
1414 (state7[2]-O.Z())*U.Z()); // u = (pos - O) * U
1415 state5(4) = ((state7[0]-O.X())*V.X() +
1416 (state7[1]-O.Y())*V.Y() +
1418
1419 setSpu(state, spu);
1420
1421}
void setSpu(StateOnPlane &state, double spu) const
Definition RKTrackRep.cc:1136
◆ getState5() [2/2]
|
private |
◆ getState7() [1/2]
|
private |
Definition at line 1351 of file RKTrackRep.cc.
1351 {
1352
1353 if (dynamic_cast<const MeasurementOnPlane*>(&state) != NULL) {
1354 Exception exc("RKTrackRep::getState7 - cannot get pos or mom from a MeasurementOnPlane",__LINE__,__FILE__);
1355 exc.setFatal();
1356 throw exc;
1357 }
1358
1359 const TVector3& U(state.getPlane()->getU());
1361 const TVector3& O(state.getPlane()->getO());
1362 const TVector3& W(state.getPlane()->getNormal());
1363
1364 assert(state.getState().GetNrows() == 5);
1365 const double* state5 = state.getState().GetMatrixArray();
1366
1368
1372
1376
1377 // normalize dir
1380
1381 state7[6] = state5[0]; // q/p
1382}
double getSpu(const StateOnPlane &state) const
Definition RKTrackRep.cc:758
◆ getState7() [2/2]
|
private |
◆ getSteps() [1/2]
|
virtual |
Get stepsizes and material properties of crossed materials of the last extrapolation.
Implements genfit::AbsTrackRep.
Definition at line 885 of file RKTrackRep.cc.
◆ getSteps() [2/2]
|
virtual |
Get stepsizes and material properties of crossed materials of the last extrapolation.
Implements genfit::AbsTrackRep.
◆ getTOF() [1/2]
|
virtual |
Get the time of flight [ns] of the last extrapolation.
Implements genfit::AbsTrackRep.
Definition at line 924 of file RKTrackRep.cc.
924 {
925
926 // Todo: test
927
929 Exception exc("RKTrackRep::getTOF ==> cache is empty.",__LINE__,__FILE__);
930 throw exc;
931 }
932
937
938 double tof(0);
939
941
945
946 if (fabs(p1-p2) < 1E-6) {
950 }
951 else {
952 // assume linear momentum loss
953 tof += 1.E9 / c / (p1 - p2) * trackLen *
954 (sqrt(m2 + p1*p1) - sqrt(m2 + p2*p2) +
956 }
957
958 p1 = p2;
959 }
960
961 return tof;
962}
double getMass(const StateOnPlane &state) const
Get tha particle mass in GeV/c^2.
Definition AbsTrackRep.cc:100
◆ getTOF() [2/2]
|
virtual |
Get the time of flight [ns] of the last extrapolation.
Implements genfit::AbsTrackRep.
◆ initArrays() [1/2]
|
private |
Definition at line 1332 of file RKTrackRep.cc.
1332 {
1336 noiseProjection_[i*8] = 1;
1338
1339 fJacobian_.UnitMatrix();
1340 fNoise_.Zero();
1342
1343 RKSteps_.reserve(100);
1344 ExtrapSteps_.reserve(100);
1345
1348}
◆ initArrays() [2/2]
|
private |
◆ isSame() [1/2]
|
virtual |
check if other is of same type (e.g. RKTrackRep) and has same pdg code.
Implements genfit::AbsTrackRep.
Definition at line 2502 of file RKTrackRep.cc.
2502 {
2504 return false;
2505
2507}
virtual bool isSameType(const AbsTrackRep *other)
check if other is of same type (e.g. RKTrackRep).
Definition RKTrackRep.cc:2494
◆ isSame() [2/2]
|
virtual |
check if other is of same type (e.g. RKTrackRep) and has same pdg code.
Implements genfit::AbsTrackRep.
◆ isSameType() [1/2]
|
virtual |
check if other is of same type (e.g. RKTrackRep).
Implements genfit::AbsTrackRep.
Definition at line 2494 of file RKTrackRep.cc.
2494 {
2496 return false;
2497
2498 return true;
2499}
◆ isSameType() [2/2]
|
virtual |
check if other is of same type (e.g. RKTrackRep).
Implements genfit::AbsTrackRep.
◆ makePlane()
| SharedPlanePtr genfit::RKTrackRep::makePlane | ( | const TVector3 & | o, |
| const TVector3 & | u, | ||
| const TVector3 & | v | ||
| ) |
Definition at line 79 of file RKTrackRep.cc.
◆ momMag() [1/2]
|
private |
Definition at line 2488 of file RKTrackRep.cc.
2488 {
2489 double momMag2 = state7[3]*state7[3] + state7[4]*state7[4] + state7[5]*state7[5];
2490 return sqrt(momMag2);
2491}
◆ momMag() [2/2]
|
private |
◆ pocaOnLine() [1/2]
|
private |
Definition at line 2189 of file RKTrackRep.cc.
2189 {
2190
2191 TVector3 retVal(lineDirection);
2192
2194 retVal *= t;
2195 retVal += linePoint;
2196 return retVal; // = linePoint + t*lineDirection
2197
2198}
◆ pocaOnLine() [2/2]
|
private |
◆ RKPropagate() [1/2]
| double genfit::RKTrackRep::RKPropagate | ( | M1x7 & | state7, |
| M7x7 * | jacobian, | ||
| M1x3 & | SA, | ||
| double | S, | ||
| bool | varField = true, |
||
| bool | calcOnlyLastRowOfJ = false |
||
| ) | const |
The actual Runge Kutta propagation.
propagate state7 with step S. Fills SA (Start directions derivatives dA/S). If jacobian is NULL, only the state is propagated, otherwise also the 7x7 jacobian is calculated. If varField is false, the magnetic field will only be evaluated at the starting position. The return value is an estimation on how good the extrapolation is, and it is usually fine if it is > 1. It gives a suggestion how you must scale S so that the quality will be sufficient.
Definition at line 1143 of file RKTrackRep.cc.
1148 {
1149
1150 // important fixed numbers
1151 static const double EC = 0.000149896229; // c/(2*10^12) resp. c/2Tera
1152 static const double P3 = 1./3.; // 1/3
1153 static const double DLT = .0002; // max. deviation for approximation-quality test
1155 // Aux parameters
1158 double S3(0), S4(0), PS2(0);
1159 M1x3 H0 = {0.,0.,0.}, H1 = {0.,0.,0.}, H2 = {0.,0.,0.};
1161 // Variables for Runge Kutta solver
1162 double A0(0), A1(0), A2(0), A3(0), A4(0), A5(0), A6(0);
1163 double B0(0), B1(0), B2(0), B3(0), B4(0), B5(0), B6(0);
1164 double C0(0), C1(0), C2(0), C3(0), C4(0), C5(0), C6(0);
1165
1166 //
1167 // Runge Kutta Extrapolation
1168 //
1169 S3 = P3*S;
1170 S4 = 0.25*S;
1171 PS2 = state7[6]*EC * S;
1172
1173 // First point
1175 FieldManager::getInstance()->getFieldVal(r[0], r[1], r[2], H0[0], H0[1], H0[2]); // magnetic field in 10^-1 T = kGauss
1176 H0[0] *= PS2; H0[1] *= PS2; H0[2] *= PS2; // H0 is PS2*(Hx, Hy, Hz) @ R0
1177 A0 = A[1]*H0[2]-A[2]*H0[1]; B0 = A[2]*H0[0]-A[0]*H0[2]; C0 = A[0]*H0[1]-A[1]*H0[0]; // (ax, ay, az) x H0
1180
1181 // Second point
1182 if (varField) {
1185 H1[0] *= PS2; H1[1] *= PS2; H1[2] *= PS2; // H1 is PS2*(Hx, Hy, Hz) @ (x, y, z) + 0.25*S * [(A0, B0, C0) + 2*(ax, ay, az)]
1186 }
1187 else { H1[0] = H0[0]; H1[1] = H0[1]; H1[2] = H0[2]; }; // invalid: H1 = H0; !!
1188 A3 = B2*H1[2]-C2*H1[1]+A[0]; B3 = C2*H1[0]-A2*H1[2]+A[1]; C3 = A2*H1[1]-B2*H1[0]+A[2]; // (A2, B2, C2) x H1 + (ax, ay, az)
1189 A4 = B3*H1[2]-C3*H1[1]+A[0]; B4 = C3*H1[0]-A3*H1[2]+A[1]; C4 = A3*H1[1]-B3*H1[0]+A[2]; // (A3, B3, C3) x H1 + (ax, ay, az)
1191
1192 // Last point
1193 if (varField) {
1196 H2[0] *= PS2; H2[1] *= PS2; H2[2] *= PS2; // H2 is PS2*(Hx, Hy, Hz) @ (x, y, z) + 0.25*S * (A4, B4, C4)
1197 }
1198 else { H2[0] = H0[0]; H2[1] = H0[1]; H2[2] = H0[2]; }; // invalid: H2 = H0; !!
1199 A6 = B5*H2[2]-C5*H2[1]; B6 = C5*H2[0]-A5*H2[2]; C6 = A5*H2[1]-B5*H2[0]; // (A5, B5, C5) x H2
1200
1201
1202 //
1203 // Derivatives of track parameters
1204 //
1205 if(jacobianT != NULL){
1206
1207 // jacobianT
1208 // 1 0 0 0 0 0 0
1209 // 0 1 0 0 0 0 0
1210 // 0 0 1 0 0 0 0
1211 // x x x x x x 0
1212 // x x x x x x 0
1213 // x x x x x x 0
1214 // x x x x x x 1
1215
1216 double dA0(0), dA2(0), dA3(0), dA4(0), dA5(0), dA6(0);
1217 double dB0(0), dB2(0), dB3(0), dB4(0), dB5(0), dB6(0);
1218 double dC0(0), dC2(0), dC3(0), dC4(0), dC5(0), dC6(0);
1219
1221
1222 if (!calcOnlyLastRowOfJ) {
1223
1224 if (!varField) {
1225 // d(x, y, z)/d(x, y, z) submatrix is unit matrix
1226 (*jacobianT)[0] = 1; (*jacobianT)[8] = 1; (*jacobianT)[16] = 1;
1227 // d(ax, ay, az)/d(ax, ay, az) submatrix is 0
1228 // start with d(x, y, z)/d(ax, ay, az)
1229 start = 3;
1230 }
1231
1233
1234 //first point
1238
1242
1243 //second point
1247
1251
1252 //last point
1256
1257 dA6 = dB5*H2[2]-dC5*H2[1]; // dA6/dp }
1258 dB6 = dC5*H2[0]-dA5*H2[2]; // dB6/dp } = (dA5, dB5, dC5) x H2
1259 dC6 = dA5*H2[1]-dB5*H2[0]; // dC6/dp }
1260
1261 // this gives the same results as multiplying the old with the new Jacobian
1262 (*jacobianT)[i] += (dA2+dA3+dA4)*S3; (*jacobianT)[i+3] = ((dA0+2.*dA3)+(dA5+dA6))*P3; // dR := dR + S3*[(dA2, dB2, dC2) + (dA3, dB3, dC3) + (dA4, dB4, dC4)]
1263 (*jacobianT)[i+1] += (dB2+dB3+dB4)*S3; (*jacobianT)[i+4] = ((dB0+2.*dB3)+(dB5+dB6))*P3; // dA := 1/3*[(dA0, dB0, dC0) + 2*(dA3, dB3, dC3) + (dA5, dB5, dC5) + (dA6, dB6, dC6)]
1265 }
1266
1267 } // end if (!calcOnlyLastRowOfJ)
1268
1270
1271 (*jacobianT)[i+3] *= state7[6]; (*jacobianT)[i+4] *= state7[6]; (*jacobianT)[i+5] *= state7[6];
1272
1273 //first point
1275 dB0 = H0[0]*(*jacobianT)[i+5]-H0[2]*(*jacobianT)[i+3] + B0; // dB0/dp } = dA x H0 + (A0, B0, C0)
1277
1281
1282 //second point
1284 dB3 = (*jacobianT)[i+4]+dC2*H1[0]-dA2*H1[2] + (B3-A[1]); // dB3/dp } = dA + (dA2, dB2, dC2) x H1
1286
1288 dB4 = (*jacobianT)[i+4]+dC3*H1[0]-dA3*H1[2] + (B4-A[1]); // dB4/dp } = dA + (dA3, dB3, dC3) x H1
1290
1291 //last point
1295
1296 dA6 = dB5*H2[2]-dC5*H2[1] + A6; // dA6/dp }
1297 dB6 = dC5*H2[0]-dA5*H2[2] + B6; // dB6/dp } = (dA5, dB5, dC5) x H2 + (A6, B6, C6)
1298 dC6 = dA5*H2[1]-dB5*H2[0] + C6; // dC6/dp }
1299
1300 // this gives the same results as multiplying the old with the new Jacobian
1301 (*jacobianT)[i] += (dA2+dA3+dA4)*S3/state7[6]; (*jacobianT)[i+3] = ((dA0+2.*dA3)+(dA5+dA6))*P3/state7[6]; // dR := dR + S3*[(dA2, dB2, dC2) + (dA3, dB3, dC3) + (dA4, dB4, dC4)]
1302 (*jacobianT)[i+1] += (dB2+dB3+dB4)*S3/state7[6]; (*jacobianT)[i+4] = ((dB0+2.*dB3)+(dB5+dB6))*P3/state7[6]; // dA := 1/3*[(dA0, dB0, dC0) + 2*(dA3, dB3, dC3) + (dA5, dB5, dC5) + (dA6, dB6, dC6)]
1303 (*jacobianT)[i+2] += (dC2+dC3+dC4)*S3/state7[6]; (*jacobianT)[i+5] = ((dC0+2.*dC3)+(dC5+dC6))*P3/state7[6];
1304
1305 }
1306
1307 //
1308 // Track parameters in last point
1309 //
1310 R[0] += (A2+A3+A4)*S3; A[0] += (SA[0]=((A0+2.*A3)+(A5+A6))*P3-A[0]); // R = R0 + S3*[(A2, B2, C2) + (A3, B3, C3) + (A4, B4, C4)]
1311 R[1] += (B2+B3+B4)*S3; A[1] += (SA[1]=((B0+2.*B3)+(B5+B6))*P3-A[1]); // A = 1/3*[(A0, B0, C0) + 2*(A3, B3, C3) + (A5, B5, C5) + (A6, B6, C6)]
1313
1314 // normalize A
1315 double CBA = 1./sqrt(A[0]*A[0]+A[1]*A[1]+A[2]*A[2]); // 1/|A|
1317
1318
1319 // Test approximation quality on given step
1320 double EST = fabs((A1+A6)-(A3+A4)) +
1321 fabs((B1+B6)-(B3+B4)) +
1322 fabs((C1+C6)-(C3+C4)); // EST = ||(ABC1+ABC6)-(ABC3+ABC4)||_1 = ||(axzy x H0 + ABC5 x H2) - (ABC2 x H1 + ABC3 x H1)||_1
1323 if (EST < 1.E-7) EST = 1.E-7; // prevent q from getting too large
1326 }
1327 return pow(DLT/EST, par);
1328}
TVector3 getFieldVal(const TVector3 &position)
This does NOT use the cache!
Definition FieldManager.h:63
◆ RKPropagate() [2/2]
| double genfit::RKTrackRep::RKPropagate | ( | M1x7 & | state7, |
| M7x7 * | jacobian, | ||
| M1x3 & | SA, | ||
| double | S, | ||
| bool | varField = true, |
||
| bool | calcOnlyLastRowOfJ = false |
||
| ) | const |
The actual Runge Kutta propagation.
propagate state7 with step S. Fills SA (Start directions derivatives dA/S). If jacobian is NULL, only the state is propagated, otherwise also the 7x7 jacobian is calculated. If varField is false, the magnetic field will only be evaluated at the starting position. The return value is an estimation on how good the extrapolation is, and it is usually fine if it is > 1. It gives a suggestion how you must scale S so that the quality will be sufficient.
◆ RKutta() [1/2]
|
private |
Propagates the particle through the magnetic field.
If the propagation is successful and the plane is reached, the function returns true. Propagated state and the jacobian of the extrapolation are written to state7 and jacobianT. The jacobian is only calculated if jacobianT != NULL. In the main loop of the Runge Kutta algorithm, the estimateStep() is called and may reduce the estimated stepsize so that a maximum momentum loss will not be exceeded, and stop at material boundaries. If this is the case, RKutta() will only propagate the reduced distance and then return. This is to ensure that material effects, which are calculated after the propagation, are taken into account properly.
Definition at line 1705 of file RKTrackRep.cc.
1715 {
1716
1717 // limits, check-values, etc. Can be tuned!
1718 static const double Wmax = 6000.; // max. way allowed [cm] TR 1/12/2014 3000->6000
1719 static const double AngleMax = 6.3; // max. total angle change of momentum. Prevents extrapolating a curler round and round if no active plane is found.
1720 static const double Pmin = 4.E-3; // minimum momentum for propagation [GeV]
1721 static const unsigned int maxNumIt = 1000; // maximum number of iterations in main loop
1722 // Aux parameters
1726 double Way = 0.; // Sum of absolute values of all extrapolation steps [cm]
1728 double relMomLoss = 0; // relative momentum loss in RKutta
1729 double deltaAngle = 0.; // total angle by which the momentum has changed during extrapolation
1731
1732
1734 std::cout << "RKTrackRep::RKutta \n";
1735 std::cout << "position: "; TVector3(R[0], R[1], R[2]).Print();
1736 std::cout << "direction: "; TVector3(A[0], A[1], A[2]).Print();
1738 std::cout << "destination: "; plane.Print();
1739 }
1740
1741 checkJacProj = false;
1742
1743 // check momentum
1744 if(momentum < Pmin){
1745 std::ostringstream sstream;
1747 Exception exc(sstream.str(),__LINE__,__FILE__);
1748 exc.setFatal();
1749 throw exc;
1750 }
1751
1753
1754 // Step estimation (signed)
1756
1757 //
1758 // Main loop of Runge-Kutta method
1759 //
1761
1763 Exception exc("RKTrackRep::RKutta ==> maximum number of iterations exceeded",__LINE__,__FILE__);
1764 exc.setFatal();
1765 throw exc;
1766 }
1767
1770 }
1771
1773 RKPropagate(state7, jacobianT, SA, S, true, calcOnlyLastRowOfJ); // the actual Runge Kutta propagation
1774
1775 // update paths
1777 Way += fabs(S);
1778
1779 // check way limit
1780 if(Way > Wmax){
1781 std::ostringstream sstream;
1782 sstream << "RKTrackRep::RKutta ==> Total extrapolation length is longer than length limit : " << Way << " cm !";
1783 Exception exc(sstream.str(),__LINE__,__FILE__);
1784 exc.setFatal();
1785 throw exc;
1786 }
1787
1788 if (onlyOneStep) return(true);
1789
1790 // if stepsize has been limited by material, break the loop and return. No linear extrapolation!
1793 std::cout<<" momLossExceeded -> return(true); \n";
1794 }
1795 return(true);
1796 }
1797
1798 // if stepsize has been limited by material boundary, break the loop and return. No linear extrapolation!
1801 std::cout<<" at boundary -> return(true); \n";
1802 }
1803 return(true);
1804 }
1805
1806
1807 // estimate Step for next loop or linear extrapolation
1809 limits.removeLimit(stp_fieldCurv);
1810 limits.removeLimit(stp_momLoss);
1811 limits.removeLimit(stp_boundary);
1812 limits.removeLimit(stp_plane);
1814
1816 fabs(S) < MINSTEP) {
1818 std::cout<<" (at Plane && fabs(S) < MINSTEP) -> break and do linear extrapolation \n";
1819 }
1820 break;
1821 }
1823 fabs(S) < MINSTEP) {
1825 std::cout<<" (momLossExceeded && fabs(S) < MINSTEP) -> return(true), no linear extrapolation; \n";
1826 }
1828 --RKStepsFXStop_;
1829 return(true); // no linear extrapolation!
1830 }
1831
1832 // check if total angle is bigger than AngleMax. Can happen if a curler should be fitted and it does not hit the active area of the next plane.
1833 deltaAngle += acos(ABefore[0]*A[0] + ABefore[1]*A[1] + ABefore[2]*A[2]);
1834 if (fabs(deltaAngle) > AngleMax){
1835 std::ostringstream sstream;
1836 sstream << "RKTrackRep::RKutta ==> Do not get to an active plane! Already extrapolated " << deltaAngle * 180 / TMath::Pi() << "°.";
1837 Exception exc(sstream.str(),__LINE__,__FILE__);
1838 exc.setFatal();
1839 throw exc;
1840 }
1841
1842 // check if we went back and forth multiple times -> we don't come closer to the plane!
1847 Exception exc("RKTrackRep::RKutta ==> Do not get closer to plane!",__LINE__,__FILE__);
1848 exc.setFatal();
1849 throw exc;
1850 }
1851 }
1852
1853 } //end of main loop
1854
1855
1856 //
1857 // linear extrapolation to plane
1858 //
1860
1863 std::cout << " RKutta - linear extrapolation to surface\n";
1864 }
1865 Sl = 1./Sl; // Sl = inverted last Stepsize Sl
1866
1867 // normalize SA
1868 SA[0]*=Sl; SA[1]*=Sl; SA[2]*=Sl; // SA/Sl = delta A / delta way; local derivative of A with respect to the length of the way
1869
1870 // calculate A
1874
1875 // normalize A
1876 CBA = 1./sqrt(A[0]*A[0]+A[1]*A[1]+A[2]*A[2]); // 1/|A|
1878
1879 R[0] += S*(A[0]-0.5*S*SA[0]); // R = R + S*(A - 0.5*S*SA); approximation for final point on surface
1882
1883
1884 coveredDistance += S;
1885 Way += fabs(S);
1886 }
1888 std::cout << " RKutta - last stepsize too small -> can't do linear extrapolation! \n";
1889 }
1890
1891 //
1892 // Project Jacobian of extrapolation onto destination plane
1893 //
1894 if (jacobianT != NULL) {
1895
1896 // projected jacobianT
1897 // x x x x x x 0
1898 // x x x x x x 0
1899 // x x x x x x 0
1900 // x x x x x x 0
1901 // x x x x x x 0
1902 // x x x x x x 0
1903 // x x x x x x 1
1904
1906 Exception exc("RKTrackRep::Extrap ==> covariance is projected onto destination plane again",__LINE__,__FILE__);
1907 throw exc;
1908 }
1909
1911 //std::cout << " Jacobian^T of extrapolation before Projection:\n";
1912 //RKTools::printDim(*jacobianT, 7,7);
1913 std::cout << " Project Jacobian of extrapolation onto destination plane\n";
1914 }
1916 An = (fabs(An) > 1.E-7 ? 1./An : 0); // 1/A_normal
1919 if (calcOnlyLastRowOfJ)
1920 i = 42;
1921
1923 double* jacPtr = *jacobianT;
1927 }
1928 checkJacProj = true;
1929
1930
1932 //std::cout << " Jacobian^T of extrapolation after Projection:\n";
1933 //RKTools::printDim(*jacobianT, 7,7);
1934 }
1935
1936 if (!calcOnlyLastRowOfJ) {
1937 for (int iRow = 0; iRow < 3; ++iRow) {
1938 for (int iCol = 0; iCol < 3; ++iCol) {
1939 noiseProjection[iRow*7 + iCol] = (iRow == iCol) - An * SU[iCol] * A[iRow];
1940 noiseProjection[(iRow + 3)*7 + iCol] = - An * SU[iCol] * SA[iRow];
1941
1942 }
1943 }
1944
1945 // noiseProjection will look like this:
1946 // x x x 0 0 0 0
1947 // x x x 0 0 0 0
1948 // x x x 0 0 0 0
1949 // x x x 1 0 0 0
1950 // x x x 0 1 0 0
1951 // x x x 0 0 1 0
1952 // 0 0 0 0 0 0 1
1953 }
1954
1955 }
1956 } // end of linear extrapolation to surface
1957
1958 return(true);
1959
1960}
double estimateStep(const M1x7 &state7, const M1x4 &SU, const DetPlane &plane, const double &charge, double &relMomLoss, StepLimits &limits) const
Definition RKTrackRep.cc:1963
◆ RKutta() [2/2]
|
private |
Propagates the particle through the magnetic field.
If the propagation is successful and the plane is reached, the function returns true. Propagated state and the jacobian of the extrapolation are written to state7 and jacobianT. The jacobian is only calculated if jacobianT != NULL. In the main loop of the Runge Kutta algorithm, the estimateStep() is called and may reduce the estimated stepsize so that a maximum momentum loss will not be exceeded, and stop at material boundaries. If this is the case, RKutta() will only propagate the reduced distance and then return. This is to ensure that material effects, which are calculated after the propagation, are taken into account properly.
◆ setChargeSign() [1/2]
|
virtual |
Set the sign of the charge according to charge.
Implements genfit::AbsTrackRep.
Definition at line 1122 of file RKTrackRep.cc.
1122 {
1123
1124 if (dynamic_cast<MeasurementOnPlane*>(&state) != NULL) {
1125 Exception exc("RKTrackRep::setChargeSign - cannot set charge of a MeasurementOnPlane",__LINE__,__FILE__);
1126 exc.setFatal();
1127 throw exc;
1128 }
1129
1130 if (state.getState()(0) * charge < 0) {
1131 state.getState()(0) *= -1.;
1132 }
1133}
◆ setChargeSign() [2/2]
|
virtual |
Set the sign of the charge according to charge.
Implements genfit::AbsTrackRep.
◆ setPosMom() [1/4]
|
virtual |
Set position and momentum of state.
Implements genfit::AbsTrackRep.
Definition at line 965 of file RKTrackRep.cc.
965 {
966
967 if (state.getRep() != this){
968 Exception exc("RKTrackRep::setPosMom ==> state is defined wrt. another TrackRep",__LINE__,__FILE__);
969 throw exc;
970 }
971
972 if (dynamic_cast<MeasurementOnPlane*>(&state) != NULL) {
973 Exception exc("RKTrackRep::setPosMom - cannot set pos/mom of a MeasurementOnPlane",__LINE__,__FILE__);
974 exc.setFatal();
975 throw exc;
976 }
977
978 // init auxInfo if that has not yet happened
979 TVectorD& auxInfo = state.getAuxInfo();
980 if (auxInfo.GetNrows() != 1) {
981 auxInfo.ResizeTo(1);
982 setSpu(state, 1.);
983 }
984
985 if (state.getPlane() != NULL && state.getPlane()->distance(pos) < MINSTEP) { // pos is on plane -> do not change plane!
986
987 M1x7 state7;
988
989 state7[0] = pos.X();
990 state7[1] = pos.Y();
991 state7[2] = pos.Z();
992
993 state7[3] = mom.X();
994 state7[4] = mom.Y();
995 state7[5] = mom.Z();
996
997 // normalize dir
1000 state7[i] *= norm;
1001
1003
1004 getState5(state, state7);
1005
1006 }
1007 else { // pos is not on plane -> create new plane!
1008
1009 // TODO: Raise Warning that a new plane has been created!
1011 state.setPlane(plane);
1012
1013 TVectorD& state5(state.getState());
1014
1016 state5(1) = 0.; // u'
1017 state5(2) = 0.; // v'
1018 state5(3) = 0.; // u
1019 state5(4) = 0.; // v
1020
1021 setSpu(state, 1.);
1022 }
1023
1024}
◆ setPosMom() [2/4]
|
virtual |
Set position and momentum of state.
Implements genfit::AbsTrackRep.
◆ setPosMom() [3/4]
|
virtual |
Set position and momentum of state.
Implements genfit::AbsTrackRep.
Definition at line 1027 of file RKTrackRep.cc.
1027 {
1028 if (state6.GetNrows()!=6){
1029 Exception exc("RKTrackRep::setPosMom ==> state has to be 6d (x, y, z, px, py, pz)",__LINE__,__FILE__);
1030 throw exc;
1031 }
1032 setPosMom(state, TVector3(state6(0), state6(1), state6(2)), TVector3(state6(3), state6(4), state6(5)));
1033}
virtual void setPosMom(StateOnPlane &state, const TVector3 &pos, const TVector3 &mom) const
Set position and momentum of state.
Definition RKTrackRep.cc:965
◆ setPosMom() [4/4]
|
virtual |
Set position and momentum of state.
Implements genfit::AbsTrackRep.
◆ setPosMomCov() [1/4]
|
virtual |
Set position, momentum and covariance of state.
Implements genfit::AbsTrackRep.
Definition at line 1078 of file RKTrackRep.cc.
1078 {
1079
1080 if (cov6x6.GetNcols()!=6 || cov6x6.GetNrows()!=6){
1081 Exception exc("RKTrackRep::setPosMomCov ==> cov has to be 6x6 (x, y, z, px, py, pz)",__LINE__,__FILE__);
1082 throw exc;
1083 }
1084
1086
1087 M1x7 state7;
1088 getState7(state, state7);
1089
1091
1092 transformM6P(cov6x6_, state7, state);
1093
1094}
void transformM6P(const M6x6 &in6x6, const M1x7 &state7, MeasuredStateOnPlane &state) const
Definition RKTrackRep.cc:1627
◆ setPosMomCov() [2/4]
|
virtual |
Set position, momentum and covariance of state.
Implements genfit::AbsTrackRep.
◆ setPosMomCov() [3/4]
|
virtual |
Set position, momentum and covariance of state.
Implements genfit::AbsTrackRep.
Definition at line 1096 of file RKTrackRep.cc.
1096 {
1097
1098 if (state6.GetNrows()!=6){
1099 Exception exc("RKTrackRep::setPosMomCov ==> state has to be 6d (x, y, z, px, py, pz)",__LINE__,__FILE__);
1100 throw exc;
1101 }
1102
1103 if (cov6x6.GetNcols()!=6 || cov6x6.GetNrows()!=6){
1104 Exception exc("RKTrackRep::setPosMomCov ==> cov has to be 6x6 (x, y, z, px, py, pz)",__LINE__,__FILE__);
1105 throw exc;
1106 }
1107
1108 TVector3 pos(state6(0), state6(1), state6(2));
1109 TVector3 mom(state6(3), state6(4), state6(5));
1111
1112 M1x7 state7;
1113 getState7(state, state7);
1114
1116
1117 transformM6P(cov6x6_, state7, state);
1118
1119}
◆ setPosMomCov() [4/4]
|
virtual |
Set position, momentum and covariance of state.
Implements genfit::AbsTrackRep.
◆ setPosMomErr() [1/2]
|
virtual |
Set position and momentum and error of state.
Implements genfit::AbsTrackRep.
Definition at line 1036 of file RKTrackRep.cc.
1036 {
1037
1038 // TODO: test!
1039
1040 setPosMom(state, pos, mom);
1041
1042 const TVector3& U(state.getPlane()->getU());
1044 TVector3 W(state.getPlane()->getNormal());
1045
1049
1050 TMatrixDSym& cov(state.getCov());
1051
1056
1057 cov(1,1) = pow((U.X()/pw - W.X()*pu/(pw*pw)),2.) * momErr.X()*momErr.X() +
1058 pow((U.Y()/pw - W.Y()*pu/(pw*pw)),2.) * momErr.Y()*momErr.Y() +
1059 pow((U.Z()/pw - W.Z()*pu/(pw*pw)),2.) * momErr.Z()*momErr.Z();
1060
1061 cov(2,2) = pow((V.X()/pw - W.X()*pv/(pw*pw)),2.) * momErr.X()*momErr.X() +
1062 pow((V.Y()/pw - W.Y()*pv/(pw*pw)),2.) * momErr.Y()*momErr.Y() +
1063 pow((V.Z()/pw - W.Z()*pv/(pw*pw)),2.) * momErr.Z()*momErr.Z();
1064
1065 cov(3,3) = posErr.X()*posErr.X() * U.X()*U.X() +
1066 posErr.Y()*posErr.Y() * U.Y()*U.Y() +
1067 posErr.Z()*posErr.Z() * U.Z()*U.Z();
1068
1072
1073}
◆ setPosMomErr() [2/2]
|
virtual |
Set position and momentum and error of state.
Implements genfit::AbsTrackRep.
◆ setQop() [1/2]
|
inlinevirtual |
Set charge/momentum.
Implements genfit::AbsTrackRep.
Definition at line 163 of file RKTrackRep.h.
163{state.getState()(0) = qop;}
◆ setQop() [2/2]
|
inlinevirtual |
Set charge/momentum.
Implements genfit::AbsTrackRep.
Definition at line 168 of file RKTrackRep.h.
168{state.getState()(0) = qop;}
◆ setSpu() [1/2]
| void genfit::RKTrackRep::setSpu | ( | StateOnPlane & | state, |
| double | spu | ||
| ) | const |
Definition at line 1136 of file RKTrackRep.cc.
1136 {
1137 state.getAuxInfo().ResizeTo(1);
1138 (state.getAuxInfo())(0) = spu;
1139}
◆ setSpu() [2/2]
| void genfit::RKTrackRep::setSpu | ( | StateOnPlane & | state, |
| double | spu | ||
| ) | const |
◆ transformM6P() [1/2]
|
private |
Definition at line 1627 of file RKTrackRep.cc.
1629 { // plane and charge must already be set!
1630
1631 // get vectors and aux variables
1632 const TVector3& U(state.getPlane()->getU());
1634 const TVector3& W(state.getPlane()->getNormal());
1635
1636 const double AtU = state7[3]*U.X() + state7[4]*U.Y() + state7[5]*U.Z();
1638 const double AtW = state7[3]*W.X() + state7[4]*W.Y() + state7[5]*W.Z();
1639
1640 // J_Mp matrix is d(q/p,u',v',u,v) / d(x,y,z,px,py,pz) (in is 6x6)
1641
1642 const double qop = state7[6];
1644
1645 M6x5 J_Mp_6x5;
1647
1648 //d(u)/d(x,y,z)
1649 J_Mp_6x5[3] = U.X(); // [0][3]
1650 J_Mp_6x5[8] = U.Y(); // [1][3]
1651 J_Mp_6x5[13] = U.Z(); // [2][3]
1652 //d(v)/d(x,y,z)
1656 // d(q/p)/d(px,py,pz)
1657 double fact = (-1.) * qop / p;
1658 J_Mp_6x5[15] = fact * state7[3]; // [3][0]
1659 J_Mp_6x5[20] = fact * state7[4]; // [4][0]
1660 J_Mp_6x5[25] = fact * state7[5]; // [5][0]
1661 // d(u')/d(px,py,pz)
1662 fact = 1./(p*AtW*AtW);
1663 J_Mp_6x5[16] = fact * (U.X()*AtW - W.X()*AtU); // [3][1]
1664 J_Mp_6x5[21] = fact * (U.Y()*AtW - W.Y()*AtU); // [4][1]
1665 J_Mp_6x5[26] = fact * (U.Z()*AtW - W.Z()*AtU); // [5][1]
1666 // d(v')/d(px,py,pz)
1670
1671 // do the transformation
1672 // out5x5 = J_Mp^T * in * J_Mp
1674 RKTools::J_MpTxcov6xJ_Mp(J_Mp_6x5, in6x6, out5x5_);
1675
1676}
void J_MpTxcov6xJ_Mp(const M6x5 &J_Mp, const M6x6 &cov6, M5x5 &out5)
Definition RKTools.cc:277
◆ transformM6P() [2/2]
|
private |
◆ transformM7P() [1/2]
|
private |
Definition at line 1559 of file RKTrackRep.cc.
1561 { // plane must already be set!
1562
1563 // get vectors and aux variables
1564 const TVector3& U(state.getPlane()->getU());
1566 const TVector3& W(state.getPlane()->getNormal());
1567
1569
1570 M5x7 J_Mp;
1571 calcJ_Mp_7x5(J_Mp, U, V, W, A);
1572
1573 // since the Jacobian contains a lot of zeros, and the resulting cov has to be symmetric,
1574 // the multiplication can be done much faster directly on array level
1575 // out5x5 = J_Mp^T * in * J_Mp
1577 RKTools::J_MpTxcov7xJ_Mp(J_Mp, in7x7, out5x5_);
1578
1579}
◆ transformM7P() [2/2]
|
private |
◆ transformPM6() [1/2]
|
private |
Definition at line 1497 of file RKTrackRep.cc.
1498 {
1499
1500 // get vectors and aux variables
1501 const TVector3& U(state.getPlane()->getU());
1503 const TVector3& W(state.getPlane()->getNormal());
1504
1505 const TVectorD& state5(state.getState());
1507
1508 M1x3 pTilde;
1512
1513 const double pTildeMag = sqrt(pTilde[0]*pTilde[0] + pTilde[1]*pTilde[1] + pTilde[2]*pTilde[2]);
1514 const double pTildeMag2 = pTildeMag*pTildeMag;
1515
1516 const double utpTildeOverpTildeMag2 = (U.X()*pTilde[0] + U.Y()*pTilde[1] + U.Z()*pTilde[2]) / pTildeMag2;
1517 const double vtpTildeOverpTildeMag2 = (V.X()*pTilde[0] + V.Y()*pTilde[1] + V.Z()*pTilde[2]) / pTildeMag2;
1518
1519 //J_pM matrix is d(x,y,z,px,py,pz) / d(q/p,u',v',u,v) (out is 6x6)
1520
1521 const double qop = state5(0);
1523
1524 M5x6 J_pM_5x6;
1526
1527 // d(px,py,pz)/d(q/p)
1529 J_pM_5x6[3] = fact * pTilde[0]; // [0][3]
1530 J_pM_5x6[4] = fact * pTilde[1]; // [0][4]
1531 J_pM_5x6[5] = fact * pTilde[2]; // [0][5]
1532 // d(px,py,pz)/d(u')
1533 fact = p * spu / pTildeMag;
1534 J_pM_5x6[9] = fact * ( U.X() - pTilde[0]*utpTildeOverpTildeMag2 ); // [1][3]
1535 J_pM_5x6[10] = fact * ( U.Y() - pTilde[1]*utpTildeOverpTildeMag2 ); // [1][4]
1536 J_pM_5x6[11] = fact * ( U.Z() - pTilde[2]*utpTildeOverpTildeMag2 ); // [1][5]
1537 // d(px,py,pz)/d(v')
1541 // d(x,y,z)/d(u)
1542 J_pM_5x6[18] = U.X(); // [3][0]
1543 J_pM_5x6[19] = U.Y(); // [3][1]
1544 J_pM_5x6[20] = U.Z(); // [3][2]
1545 // d(x,y,z)/d(v)
1549
1550
1551 // do the transformation
1552 // out = J_pM^T * in5x5 * J_pM
1554 RKTools::J_pMTxcov5xJ_pM(J_pM_5x6, in5x5_, out6x6);
1555
1556}
void J_pMTxcov5xJ_pM(const M5x7 &J_pM, const M5x5 &cov5, M7x7 &out7)
Definition RKTools.cc:36
◆ transformPM6() [2/2]
|
private |
◆ transformPM7() [1/2]
|
private |
Definition at line 1425 of file RKTrackRep.cc.
1426 {
1427
1428 // get vectors and aux variables
1429 const TVector3& U(state.getPlane()->getU());
1431 const TVector3& W(state.getPlane()->getNormal());
1432
1433 const TVectorD& state5(state.getState());
1435
1436 M1x3 pTilde;
1440
1441 M5x7 J_pM;
1442 calcJ_pM_5x7(J_pM, U, V, pTilde, spu);
1443
1444 // since the Jacobian contains a lot of zeros, and the resulting cov has to be symmetric,
1445 // the multiplication can be done much faster directly on array level
1446 // out = J_pM^T * in5x5 * J_pM
1448 RKTools::J_pMTxcov5xJ_pM(J_pM, in5x5_, out7x7);
1449}
◆ transformPM7() [2/2]
|
private |
Member Data Documentation
◆ cachePos_
|
mutableprivate |
use cached RKSteps_ for extrapolation
Definition at line 288 of file RKTrackRep.h.
◆ ExtrapSteps_
|
mutableprivate |
Definition at line 282 of file RKTrackRep.h.
◆ fJacobian_
|
mutableprivate |
steps made in Extrap during last extrapolation
Definition at line 284 of file RKTrackRep.h.
◆ fNoise_
|
mutableprivate |
Definition at line 285 of file RKTrackRep.h.
◆ J_MMT_
|
mutableprivate |
Definition at line 295 of file RKTrackRep.h.
◆ lastEndState_
|
mutableprivate |
state where the last extrapolation has started
Definition at line 278 of file RKTrackRep.h.
◆ lastStartState_
|
mutableprivate |
Definition at line 277 of file RKTrackRep.h.
◆ limits_
|
mutableprivate |
Definition at line 292 of file RKTrackRep.h.
◆ noiseArray_
|
mutableprivate |
Definition at line 293 of file RKTrackRep.h.
◆ noiseProjection_
|
mutableprivate |
noise matrix of the last extrapolation
Definition at line 294 of file RKTrackRep.h.
◆ RKSteps_
|
mutableprivate |
state where the last extrapolation has ended
Definition at line 279 of file RKTrackRep.h.
◆ RKStepsFXStart_
|
mutableprivate |
RungeKutta steps made in the last extrapolation.
Definition at line 280 of file RKTrackRep.h.
◆ RKStepsFXStop_
|
mutableprivate |
Definition at line 281 of file RKTrackRep.h.
◆ useCache_
|
mutableprivate |
Definition at line 287 of file RKTrackRep.h.
The documentation for this class was generated from the following files:
- genfit/include/genfit/RKTrackRep.h
- genfit/trackReps/include/RKTrackRep.h
- genfit/trackReps/src/RKTrackRep.cc
