Belle II Software light-2406-ragdoll
RecoKlong.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 <framework/logging/Logger.h>
11#include <framework/geometry/B2Vector3.h>
12
13#include <analysis/dataobjects/Particle.h>
14#include <mdst/dataobjects/KLMCluster.h>
15
16#include <analysis/VertexFitting/TreeFitter/RecoKlong.h>
17#include <analysis/VertexFitting/TreeFitter/FitParams.h>
18#include <analysis/VertexFitting/TreeFitter/ErrCode.h>
19
20namespace TreeFitter {
21
23 RecoParticle(particle, mother),
24 m_dim(3),
25 m_init(false),
26 m_useEnergy(false),
27 m_clusterPars(),
28 m_covariance()
29 {
30 initParams() ;
31 }
32
34 {
35 const int posindexmother = mother()->posIndex();
36
37 Eigen::Matrix<double, 1, 3> vertexToCluster = Eigen::Matrix<double, 1, 3>::Zero(1, 3);
38 for (unsigned int i = 0; i < 3; ++i) {
39 vertexToCluster(i) = m_clusterPars(i) - fitparams.getStateVector()(posindexmother + i);
40 }
41
42 const double distanceToMother = vertexToCluster.norm();
43 const double energy = m_clusterPars(3);
44 const double pdgMass = particle()->getPDGMass();
45 const double absMom = std::sqrt(energy * energy - pdgMass * pdgMass);
46
47 const int momindex = momIndex();
48
49 for (unsigned int i = 0; i < 3; ++i) {
50 //px = |p| dx/|dx|
51 fitparams.getStateVector()(momindex + i) = absMom * vertexToCluster(i) / distanceToMother;
52 }
53
54 return ErrCode(ErrCode::Status::success);
55 }
56
58 {
59 return ErrCode(ErrCode::Status::success);
60 }
61
63 {
64 const int momindex = momIndex();
65 const int posindex = mother()->posIndex();
66
67 const double factorE = 1000 * m_covariance(3, 3);
68 const double factorX = 1000; // ~ 10cm error on initial vertex
69
70 fitparams.getCovariance().block<4, 4>(momindex, momindex) =
71 Eigen::Matrix<double, 4, 4>::Identity(4, 4) * factorE;
72
73 fitparams.getCovariance().block<3, 3>(posindex, posindex) =
74 Eigen::Matrix<double, 3, 3>::Identity(3, 3) * factorX;
75
76 return ErrCode(ErrCode::Status::success);
77 }
78
80 {
81 const Belle2::KLMCluster* cluster = particle()->getKLMCluster();
82
83 const Belle2::B2Vector3D cluster_pos = cluster->getClusterPosition();
84
85 m_init = true;
86 m_covariance = Eigen::Matrix<double, 4, 4>::Zero(4, 4);
87
88 TMatrixDSym cov7 = cluster->getError7x7();
89
90 for (int row = 0; row < 3; ++row) {
91 for (int col = 0; col < 3; ++col) {
92 m_covariance(row, col) = cov7[row + 4][col + 4] ;
93 }
94 }
95
99 if (0 == m_covariance(3, 3)) {m_covariance(3, 3) = .214;}
100
101 m_clusterPars(0) = cluster_pos.X();
102 m_clusterPars(1) = cluster_pos.Y();
103 m_clusterPars(2) = cluster_pos.Z();
104 m_clusterPars(3) = sqrt(particle()->getPDGMass() * particle()->getPDGMass() +
105 cluster->getMomentumMag() * cluster->getMomentumMag());
106
107 auto p_vec = particle()->getMomentum();
108 // find highest momentum, eliminate dim with highest mom
109 if ((std::abs(p_vec.X()) >= std::abs(p_vec.Y())) && (std::abs(p_vec.X()) >= std::abs(p_vec.Z()))) {
110 m_i1 = 0; m_i2 = 1; m_i3 = 2;
111 } else if ((std::abs(p_vec.Y()) >= std::abs(p_vec.X())) && (std::abs(p_vec.Y()) >= std::abs(p_vec.Z()))) {
112 m_i1 = 1; m_i2 = 0; m_i3 = 2;
113 } else if ((std::abs(p_vec.Z()) >= std::abs(p_vec.Y())) && (std::abs(p_vec.Z()) >= std::abs(p_vec.X()))) {
114 m_i1 = 2; m_i2 = 1; m_i3 = 0;
115 } else {
116 B2ERROR("Could not estimate highest momentum for Klong constraint. Aborting this fit.\n px: "
117 << p_vec.X() << " py: " << p_vec.Y() << " pz: " << p_vec.Z() << " calculated from Ec: " << m_clusterPars(3));
118 return ErrCode(ErrCode::Status::photondimerror);
119 }
120
121 return ErrCode(ErrCode::Status::success);
122 }
123
125 {
126 const int momindex = momIndex() ;
127 const int posindex = mother()->posIndex();
128
129 const Eigen::Matrix<double, 1, 3> x_vertex = fitparams.getStateVector().segment(posindex, 3);
130 const Eigen::Matrix<double, 1, 3> p_vec = fitparams.getStateVector().segment(momindex, 3);
131
132 if (0 == p_vec[m_i1]) { return ErrCode(ErrCode::Status::klongdimerror); }
133
134 // p_vec[m_i1] must not be 0
135 const double elim = (m_clusterPars[m_i1] - x_vertex[m_i1]) / p_vec[m_i1];
136 const double mom = p_vec.norm();
137 const double pdgMass = particle()->getPDGMass();
138 const double energy = std::sqrt(mom * mom + pdgMass * pdgMass);
139
140 // r'
141 Eigen::Matrix<double, 3, 1> residual3 = Eigen::Matrix<double, 3, 1>::Zero(3, 1);
142 residual3(0) = m_clusterPars[m_i2] - x_vertex[m_i2] - p_vec[m_i2] * elim;
143 residual3(1) = m_clusterPars[m_i3] - x_vertex[m_i3] - p_vec[m_i3] * elim;
144 residual3(2) = m_clusterPars[3] - energy;
145
146 Eigen::Matrix<double, 3, 4> P = Eigen::Matrix<double, 3, 4>::Zero(3, 4);
147
148 // dr'/dm | m:={xc,yc,zc,Ec} the measured quantities
149 P(0, m_i2) = 1;
150 P(0, m_i1) = - p_vec[m_i2] / p_vec[m_i1];
151
152 P(1, m_i3) = 1;
153 P(1, m_i1) = - p_vec[m_i3] / p_vec[m_i1];
154 P(2, 3) = 1; // dE/dEc
155
156 p.getResiduals().segment(0, 3) = residual3;
157
158 p.getV() = P * m_covariance.selfadjointView<Eigen::Lower>() * P.transpose();
159
160 // dr'/dm | m:={x,y,z,px,py,pz,E}
161 // x := x_vertex (decay vertex of mother)
162 p.getH()(0, posindex + m_i1) = p_vec[m_i2] / p_vec[m_i1];
163 p.getH()(0, posindex + m_i2) = -1.0;
164 p.getH()(0, posindex + m_i3) = 0;
165
166 p.getH()(1, posindex + m_i1) = p_vec[m_i3] / p_vec[m_i1];
167 p.getH()(1, posindex + m_i2) = 0;
168 p.getH()(1, posindex + m_i3) = -1.0;
169
170 p.getH()(0, momindex + m_i1) = p_vec[m_i2] * elim / p_vec[m_i1];
171 p.getH()(0, momindex + m_i2) = -1. * elim;
172 p.getH()(0, momindex + m_i3) = 0;
173
174 p.getH()(1, momindex + m_i1) = p_vec[m_i3] * elim / p_vec[m_i1];;
175 p.getH()(1, momindex + m_i2) = 0;
176 p.getH()(1, momindex + m_i3) = -1. * elim;
177
178 p.getH()(2, momindex + m_i1) = -1. * p_vec[m_i1] / mom;
179 p.getH()(2, momindex + m_i2) = -1. * p_vec[m_i2] / mom;
180 p.getH()(2, momindex + m_i3) = -1. * p_vec[m_i3] / mom;
181
182 return ErrCode(ErrCode::Status::success);
183 }
184
185}
186
187
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:435
DataType X() const
access variable X (= .at(0) without boundary check)
Definition: B2Vector3.h:431
DataType Y() const
access variable Y (= .at(1) without boundary check)
Definition: B2Vector3.h:433
KLM cluster data.
Definition: KLMCluster.h:28
Class to store reconstructed particles.
Definition: Particle.h:75
const KLMCluster * getKLMCluster() const
Returns the pointer to the KLMCluster object that was used to create this Particle (ParticleType == c...
Definition: Particle.cc:926
double getPDGMass(void) const
Returns uncertainty on the invariant mass (requires valid momentum error matrix)
Definition: Particle.cc:604
ROOT::Math::XYZVector getMomentum() const
Returns momentum vector.
Definition: Particle.h:560
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
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
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
Eigen::Matrix< double, 4, 4 > m_covariance
covariance (x_c,y_c,z_c,E_c) of measured pars
Definition: RecoKlong.h:76
ErrCode initCovariance(FitParams &fitparams) const override
init covariance
Definition: RecoKlong.cc:62
ErrCode initParams()
update or init params
Definition: RecoKlong.cc:79
int m_i3
another random index
Definition: RecoKlong.h:83
int m_i1
index with the highest momentum.
Definition: RecoKlong.h:79
virtual ErrCode initParticleWithMother(FitParams &fitparams) override
init particle with mother
Definition: RecoKlong.cc:33
int m_i2
random index
Definition: RecoKlong.h:81
RecoKlong(Belle2::Particle *bc, const ParticleBase *mother)
constructor
Definition: RecoKlong.cc:22
virtual ErrCode initMotherlessParticle(FitParams &fitparams) override
init particle without mother
Definition: RecoKlong.cc:57
bool m_init
was initialized*
Definition: RecoKlong.h:67
ErrCode projectRecoConstraint(const FitParams &fitparams, Projection &p) const override
project klong constraint
Definition: RecoKlong.cc:124
Eigen::Matrix< double, 1, 4 > m_clusterPars
constrains measured params (x_c, y_c, z_c, E_c)
Definition: RecoKlong.h:73
base for RecoPhoton RecoTrack
Definition: RecoParticle.h:16
virtual int momIndex() const override
get momentum index
Definition: RecoParticle.h:42