Belle II Software light-2406-ragdoll
RecoTrack.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * External Contributor: Wouter Hulsbergen *
5 * *
6 * See git log for contributors and copyright holders. *
7 * This file is licensed under LGPL-3.0, see LICENSE.md. *
8 **************************************************************************/
9
10#include <analysis/VertexFitting/TreeFitter/RecoTrack.h>
11#include <analysis/VertexFitting/TreeFitter/FitParams.h>
12#include <analysis/VertexFitting/TreeFitter/HelixUtils.h>
13
14#include <framework/gearbox/Const.h>
15#include <mdst/dataobjects/Track.h>
16
17#include <TMath.h>
18
19namespace TreeFitter {
20 constexpr double pi = TMath::Pi();
21 constexpr double twoPi = TMath::TwoPi();
22
24 RecoParticle(particle, mother),
25 m_bfield(0),
26 m_trackfit(particle->getTrackFitResult()),
27 m_cached(false),
28 m_flt(0),
29 m_params(5),
30 m_covariance(5, 5),
31 m_momentumScalingFactor(particle->getEffectiveMomentumScale())
32 {
33 m_bfield = Belle2::BFieldManager::getFieldInTesla(ROOT::Math::XYZVector(0, 0, 0)).Z(); //Bz in Tesla
34 m_covariance = Eigen::Matrix<double, 5, 5>::Zero(5, 5);
35 }
36
38 {
39 //initPar2
40 if (m_flt == 0) {
41 const_cast<RecoTrack*>(this)->updFltToMother(fitparams);
42 }
44 const int momindex = momIndex();
45 fitparams.getStateVector()(momindex) = recoP.X();
46 fitparams.getStateVector()(momindex + 1) = recoP.Y();
47 fitparams.getStateVector()(momindex + 2) = recoP.Z();
48 return ErrCode(ErrCode::Status::success);
49 }
50
52 {
53 return ErrCode(ErrCode::Status::success);
54 }
55
57 {
58 // we only need a rough estimate of the covariance
59 TMatrixFSym p4Err = particle()->getMomentumErrorMatrix();
60 const int momindex = momIndex();
61
62 for (int row = 0; row < 3; ++row) {
63 fitparams.getCovariance()(momindex + row, momindex + row) = 1000 * p4Err[row][row];
64 }
65
66 return ErrCode(ErrCode::Status::success);
67 }
68
69
71 {
72 const int posindexmother = mother()->posIndex();
74 fitparams.getStateVector()(posindexmother),
75 fitparams.getStateVector()(posindexmother + 1));
76 return ErrCode(ErrCode::Status::success);
77 } ;
78
80 {
81 m_flt = flt;
82 std::vector<float> trackParameter = m_trackfit->getTau();
83 for (unsigned int i = 0; i < trackParameter.size(); ++i) {
84 m_params(i) = trackParameter[i];
85 if (i == 2) m_params(2) /= m_momentumScalingFactor; // index 2 is the curvature of the track
86 }
87 TMatrixDSym cov = m_trackfit->getCovariance5();
88 for (int row = 0; row < 5; ++row) {
89 for (int col = 0; col <= row; ++col) {
90 m_covariance(row, col) = cov[row][col];
91 }
92 }
93
94 m_cached = true;
95 return ErrCode(ErrCode::Status::success);
96 }
97
99 {
100 ErrCode status;
101 const int posindexmother = mother()->posIndex();
102 const int momindex = momIndex();
103
104 Eigen::Matrix<double, 1, 6> positionAndMom = Eigen::Matrix<double, 1, 6>::Zero(1, 6);
105 positionAndMom.segment(0, 3) = fitparams.getStateVector().segment(posindexmother, 3);
106 positionAndMom.segment(3, 3) = fitparams.getStateVector().segment(momindex, 3);
107 Eigen::Matrix<double, 5, 6> jacobian = Eigen::Matrix<double, 5, 6>::Zero(5, 6);
108
109 const Belle2::Helix helix = Belle2::Helix(ROOT::Math::XYZVector(positionAndMom(0), positionAndMom(1), positionAndMom(2)),
110 ROOT::Math::XYZVector(positionAndMom(3), positionAndMom(4), positionAndMom(5)),
111 charge(), m_bfield);
112
114 positionAndMom(0),
115 positionAndMom(1),
116 positionAndMom(2),
117 positionAndMom(3),
118 positionAndMom(4),
119 positionAndMom(5),
120 m_bfield,
121 charge()
122 );
123 if (!m_cached) {
124 auto* nonconst = const_cast<RecoTrack*>(this);
125 if (m_flt == 0) { nonconst->updFltToMother(fitparams); }
126 nonconst->updateParams(m_flt);
127 }
128
129 Eigen::Matrix<double, 1, 5> helixpars(5);
130 helixpars(0) = helix.getD0();
131 helixpars(1) = helix.getPhi0();
132 helixpars(2) = helix.getOmega();
133 helixpars(3) = helix.getZ0();
134 helixpars(4) = helix.getTanLambda();
135
136 p.getResiduals().segment(0, 5) = m_params - helixpars;
137
138 //account for periodic boundary in phi residual
139 double phiResidual = p.getResiduals().segment(0, 5)(1);
140 phiResidual = std::fmod(phiResidual + pi, twoPi);
141 if (phiResidual < 0) phiResidual += twoPi;
142 phiResidual -= pi;
143 p.getResiduals().segment(0, 5)(1) = phiResidual;
144
145 p.getV().triangularView<Eigen::Lower>() = m_covariance.triangularView<Eigen::Lower>();
146
147 p.getH().block<5, 3>(0, posindexmother) = -1.0 * jacobian.block<5, 3>(0, 0);
148 p.getH().block<5, 3>(0, momindex) = -1.0 * jacobian.block<5, 3>(0, 3);
149
150 return status;
151 }
152
153}//end recotrack
static ROOT::Math::XYZVector getFieldInTesla(const ROOT::Math::XYZVector &pos)
return the magnetic field at a given position in Tesla.
Definition: BFieldManager.h:61
This class represents an ideal helix in perigee parameterization.
Definition: Helix.h:59
double getArcLength2DAtXY(const double &x, const double &y) const
Calculates the two dimensional arc length at which the circle in the xy projection is closest to the ...
Definition: Helix.cc:141
double getOmega() const
Getter for omega, which is a signed curvature measure of the track.
Definition: Helix.h:387
double getD0() const
Getter for d0, which is the signed distance to the perigee in the r-phi plane.
Definition: Helix.h:372
double getTanLambda() const
Getter for tan lambda, which is the z over two dimensional arc length slope of the track.
Definition: Helix.h:393
double getZ0() const
Getter for z0, which is the z coordinate of the perigee.
Definition: Helix.h:390
ROOT::Math::XYZVector getMomentumAtArcLength2D(const double &arcLength2D, const double &bz) const
Calculates the momentum vector at the given two dimensional arc length.
Definition: Helix.cc:270
double getPhi0() const
Getter for phi0, which is the azimuth angle of the transverse momentum at the perigee.
Definition: Helix.h:378
Class to store reconstructed particles.
Definition: Particle.h:75
TMatrixFSym getMomentumErrorMatrix() const
Returns the 4x4 momentum error matrix.
Definition: Particle.cc:435
Helix getHelix() const
Conversion to framework Helix (without covariance).
TMatrixDSym getCovariance5() const
Getter for covariance matrix of perigee parameters in matrix form.
std::vector< float > getTau() const
Getter for all perigee parameters.
abstract errorocode be aware that the default is success
Definition: ErrCode.h:14
Class to store and manage fitparams (statevector)
Definition: FitParams.h:20
Eigen::Matrix< double, -1, 1, 0, MAX_MATRIX_SIZE, 1 > & getStateVector()
getter for the fit parameters/statevector
Definition: FitParams.h:65
Eigen::Matrix< double, -1, -1, 0, MAX_MATRIX_SIZE, MAX_MATRIX_SIZE > & getCovariance()
getter for the states covariance
Definition: FitParams.h:53
static void getJacobianToCartesianFrameworkHelix(Eigen::Matrix< double, 5, 6 > &jacobian, const double x, const double y, const double z, const double px, const double py, const double pz, const double bfield, const double charge)
get the jacobian dh={helix pars}/dx={x,y,z,px,py,pz} for the implementation of the framework helix.
Definition: HelixUtils.cc:407
base class for all particles
Definition: ParticleBase.h:25
Belle2::Particle * particle() const
get basf2 particle
Definition: ParticleBase.h:92
virtual int posIndex() const
get vertex index (in statevector!)
Definition: ParticleBase.h:122
int charge() const
get charge
Definition: ParticleBase.h:144
const ParticleBase * mother() const
getMother() / hasMother()
Definition: ParticleBase.h:98
class to store the projected residuals and the corresponding jacobian as well as the covariance matri...
Definition: Projection.h:18
base for RecoPhoton RecoTrack
Definition: RecoParticle.h:16
virtual int momIndex() const override
get momentum index
Definition: RecoParticle.h:42
representation of all charged final states FIXME rename since this name is taken in tracking
Definition: RecoTrack.h:18
ErrCode initCovariance(FitParams &fitparams) const override
init covariance matrix of this particle constraint
Definition: RecoTrack.cc:56
const Belle2::TrackFitResult * m_trackfit
trackfit result from reconstruction
Definition: RecoTrack.h:72
Eigen::Matrix< double, 1, 5 > m_params
column vector to store the measurement
Definition: RecoTrack.h:81
const float m_momentumScalingFactor
scale the momenta by this correction factor
Definition: RecoTrack.h:87
virtual ErrCode projectRecoConstraint(const FitParams &, Projection &) const override
project the constraint (calculate residuals)
Definition: RecoTrack.cc:98
ErrCode updateParams(double flt)
updated the cached parameters
Definition: RecoTrack.cc:79
virtual ErrCode initParticleWithMother(FitParams &fitparams) override
init with mother particle (replacing initPar2)
Definition: RecoTrack.cc:37
double m_bfield
B field along z
Definition: RecoTrack.h:69
virtual ErrCode initMotherlessParticle(FitParams &fitparams) override
init without mother particle
Definition: RecoTrack.cc:51
ErrCode updFltToMother(const FitParams &fitparams)
update flight length to mother
Definition: RecoTrack.cc:70
double m_flt
helix arc length at vertex
Definition: RecoTrack.h:78
Eigen::Matrix< double, 5, 5 > m_covariance
only lower triangle filled!
Definition: RecoTrack.h:84
bool m_cached
flag to mark the particle as initialised
Definition: RecoTrack.h:75
RecoTrack(Belle2::Particle *bc, const ParticleBase *mother)
constructor
Definition: RecoTrack.cc:23