Belle II Software light-2406-ragdoll
DistanceCalculatorModule.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <analysis/modules/DistanceCalculator/DistanceCalculatorModule.h>
10
11#include <Eigen/Dense>
12#include <iostream>
13
14using namespace std;
15using namespace Belle2;
16
17const float eps = 1E-6;
18
19//-----------------------------------------------------------------
20// Register the Module
21//-----------------------------------------------------------------
22REG_MODULE(DistanceCalculator);
23
24//-----------------------------------------------------------------
25// Implementation
26//-----------------------------------------------------------------
27
29 m_distanceCovMatrix(3)
30{
31 // Set module properties
32 setDescription("Calculates distance between two vertices, distance of closest approach between a vertex and a track, distance of closest approach between two tracks, distance of closest approach between a vertex/track and Btube");
33
34 // Parameter definitions
35 addParam("listName", m_listName, "", std::string(""));
36 addParam("decayString", m_decayString, "", std::string(""));
37 addParam("mode", m_mode,
38 "Specifies how the distance is calculated:\n"
39 "vertextrack: calculate the distance of closest approach between a track and a vertex, taking the first candidate as vertex, default\n"
40 "trackvertex: calculate the distance of closest approach between a track and a vertex, taking the first candidate as track,\n"
41 "2tracks: calculates the distance of closest approach between two tracks,\n"
42 "2vertices: calculates the distance between two vertices,\n"
43 "vertexbtube: calculates the distance of closest approach between a vertex and a Btube,\n"
44 "trackbtube: calculates the distance of closest approach between a track and a Btube",
45 std::string("vertextrack"));
46
47 m_distance = Eigen::Matrix<double, Eigen::Dynamic, 1>::Zero(3);
48}
49
51{
52}
53
55{
56 if (m_decayString != "") {
58 } else {
59 B2FATAL("Please provide a decay string for DistanceCalculatorModule");
60 }
61 m_plist.isRequired(m_listName);
62}
63
64Eigen::Vector3d getDocaTrackVertex(const Particle* pTrack, const Particle* pVertex)
65{
66 Eigen::Vector3d v(pTrack->getPx(), pTrack->getPy(), pTrack->getPz());
67 Eigen::Vector3d p_ray(pTrack->getX(), pTrack->getY(), pTrack->getZ());
68 Eigen::Vector3d p(pVertex->getX(), pVertex->getY(), pVertex->getZ());
69 Eigen::Vector3d v_dir = v.normalized();
70 Eigen::Vector3d r = p - p_ray;
71 return r.dot(v_dir) * v_dir - r;
72}
73
74TMatrixFSym getDocaTrackVertexError(const Particle* pTrack, const Particle* pVertex)
75{
76 //covariance matrix of r
77 TMatrixFSym err_r = pTrack->getVertexErrorMatrix() + pVertex->getVertexErrorMatrix();
78 Eigen::Vector3d V(pTrack->getPx(), pTrack->getPy(), pTrack->getPz());
79 Eigen::Vector3d v = V.normalized();
80 //d_j = r_j - v_j * v_k r_k
81 //Jij = del_i d_j = delta_ij - v_i * v_j
82 //Since the vector of closest approach is a linear function of r, its
83 //propagation of errors is exact
84 TMatrixFSym Jacobian(3);
85 //Fill the jacobian matrix according to the equation:
86 // J_ij = delta_ij -v(i)v(j)
87 for (int i = 0; i < 3; i++)
88 for (int j = 0; j < 3; j++)
89 Jacobian(i, j) = -v(i) * v(j);
90 for (int i = 0; i < 3; i++)
91 Jacobian(i, i) += 1;
92
93 return err_r.Similarity(Jacobian); //calculates J * err_r * J^T, and returns it
94}
95
96/*
97 * This method calculates the Doca vector between two tracks (i.e. lines)
98 * This is done using the following equation:
99 * d = (n*(x1-x2))n
100 * where d is the doca vector, n is a unit vector that is orthogonal to both lines,
101 * and x1 and x2 are points on the first and second line respectively
102 */
103Eigen::Vector3d getDocaTracks(const Particle* p1, const Particle* p2)
104{
105 Eigen::Vector3d v1(p1->getPx(), p1->getPy(), p1->getPz());
106 Eigen::Vector3d v2(p2->getPx(), p2->getPy(), p2->getPz());
107 Eigen::Vector3d n = v1.cross(v2); //The doca vector must be orthogonal to both lines.
108 if (n.norm() < eps) {
109 //if tracks are parallel then the distance is independent of the point on the line
110 //and we can use the old way to calculate the DOCA
111 return getDocaTrackVertex(p1, p2);
112 }
113 Eigen::Vector3d n_dir = n.normalized();
114 Eigen::Vector3d p1v(p1->getX(), p1->getY(), p1->getZ());
115 Eigen::Vector3d p2v(p2->getX(), p2->getY(), p2->getZ());
116 Eigen::Vector3d r = p1v - p2v;
117 return r.dot(n_dir) * n_dir;
118}
119
120TMatrixFSym getDocaTracksError(const Particle* p1, const Particle* p2)
121{
122 //covariance matrix of r
123 TMatrixFSym err_r = p1->getVertexErrorMatrix() + p2->getVertexErrorMatrix();
124 Eigen::Vector3d p1p(p1->getPx(), p1->getPy(), p1->getPz());
125 Eigen::Vector3d p2p(p2->getPx(), p2->getPy(), p2->getPz());
126 Eigen::Vector3d n = p1p.cross(p2p);
127 if (n.norm() < eps) {
128 return getDocaTrackVertexError(p1, p2);
129 }
130 Eigen::Vector3d n_dir = n.normalized();
131 //d_j = n_dir_j * n_dir_k r_k
132 //Jij = del_i d_j = n_dir_i * n_dir_j
133 //Since the vector of closest approach is a linear function of r, it's
134 //propagation of errors is exact
135 TMatrixFSym Jacobian(3);
136 //Fill the jacobian matrix according to the equation:
137 // J_ij = n(i)n(j)
138 for (int i = 0; i < 3; i++)
139 for (int j = 0; j < 3; j++)
140 Jacobian(i, j) = n_dir(i) * n_dir(j);
141 return err_r.Similarity(Jacobian); //calculates J * err_r * J^T, and returns it
142}
143
144Eigen::Vector3d getDistanceVertices(const Particle* p1, const Particle* p2)
145{
146 Eigen::Vector3d p1v(p1->getX(), p1->getY(), p1->getZ());
147 Eigen::Vector3d p2v(p2->getX(), p2->getY(), p2->getZ());
148 Eigen::Vector3d r = p2v - p1v;
149 return r;
150}
151
152TMatrixFSym getDistanceVerticesErrors(const Particle* p1, const Particle* p2)
153{
154 TMatrixFSym err_r = p1->getVertexErrorMatrix() + p2->getVertexErrorMatrix();
155 return err_r;
156}
157
159{
160 if (m_mode == "2tracks") {
161 m_distance = getDocaTracks(p1, p2);
162 }
163 if (m_mode == "2vertices") {
164 m_distance = getDistanceVertices(p1, p2);
165 }
166 if (m_mode == "vertextrack") {
167 m_distance = getDocaTrackVertex(p2, p1);
168 }
169 if (m_mode == "trackvertex") {
170 m_distance = getDocaTrackVertex(p1, p2);
171 }
172}
173
175{
176 if (m_mode == "2tracks") {
177 m_distanceCovMatrix = getDocaTracksError(p1, p2);
178 }
179 if (m_mode == "2vertices") {
180 m_distanceCovMatrix = getDistanceVerticesErrors(p1, p2);
181 }
182 if (m_mode == "vertextrack") {
183 m_distanceCovMatrix = getDocaTrackVertexError(p2, p1);
184 }
185 if (m_mode == "trackvertex") {
186 m_distanceCovMatrix = getDocaTrackVertexError(p1, p2);
187 }
188}
189
190Eigen::Vector3d getDocaBtubeVertex(const Particle* pVertex, const Btube* tube)
191{
192 Eigen::Vector3d v(tube->getTubeDirection()(0, 0), tube->getTubeDirection()(1, 0), tube->getTubeDirection()(2, 0));
193 Eigen::Vector3d p_ray(tube->getTubeCenter()(0, 0), tube->getTubeCenter()(1, 0), tube->getTubeCenter()(2, 0));
194 Eigen::Vector3d p(pVertex->getX(), pVertex->getY(), pVertex->getZ());
195 Eigen::Vector3d v_dir = v.normalized();
196 Eigen::Vector3d r = p - p_ray;
197 return r.dot(v_dir) * v_dir - r;
198}
199
200TMatrixFSym getDocaBtubeVertexError(const Particle* pVertex, const Btube* tube)
201{
202 TMatrixFSym err_r = tube->getTubeCenterErrorMatrix() + pVertex->getVertexErrorMatrix();
203 Eigen::Vector3d V(tube->getTubeDirection()(0, 0), tube->getTubeDirection()(1, 0), tube->getTubeDirection()(2, 0));
204 Eigen::Vector3d v = V.normalized();
205 TMatrixFSym Jacobian(3);
206 for (int i = 0; i < 3; i++)
207 for (int j = 0; j < 3; j++)
208 Jacobian(i, j) = -v(i) * v(j);
209 for (int i = 0; i < 3; i++)
210 Jacobian(i, i) += 1;
211
212 return err_r.Similarity(Jacobian);
213}
214
215Eigen::Vector3d getDocaBtubeTrack(const Particle* pTrack, const Btube* tube)
216{
217 Eigen::Vector3d v1(pTrack->getPx(), pTrack->getPy(), pTrack->getPz());
218 Eigen::Vector3d v2 = tube->getTubeDirection();
219 Eigen::Vector3d n = v1.cross(v2); //The doca vector must be orthogonal to both lines.
220 if (n.norm() < eps) {
221 //if tracks are parallel then the distance is independent of the point on the line
222 //and we can use the old way to calculate the DOCA
223 return getDocaBtubeVertex(pTrack, tube);
224 }
225 Eigen::Vector3d n_dir = n.normalized();
226 Eigen::Vector3d pv(pTrack->getX(), pTrack->getY(), pTrack->getZ());
227 Eigen::Vector3d p2v(tube->getTubeCenter()(0, 0), tube->getTubeCenter()(1, 0), tube->getTubeCenter()(2, 0));
228 Eigen::Vector3d r = pv - p2v;
229 return r.dot(n_dir) * n_dir;
230}
231
232TMatrixFSym getDocaBtubeTrackError(const Particle* pTrack, const Btube* tube)
233{
234 TMatrixFSym err_r = pTrack->getVertexErrorMatrix() + tube->getTubeCenterErrorMatrix();
235 Eigen::Vector3d p1p(pTrack->getPx(), pTrack->getPy(), pTrack->getPz());
236 Eigen::Vector3d p2p(tube->getTubeDirection()(0, 0), tube->getTubeDirection()(1, 0), tube->getTubeDirection()(2, 0));
237 Eigen::Vector3d n = p1p.cross(p2p);
238 if (n.norm() < eps) {
239 return getDocaBtubeVertexError(pTrack, tube);
240 }
241 Eigen::Vector3d n_dir = n.normalized();
242 TMatrixFSym Jacobian(3);
243 for (int i = 0; i < 3; i++)
244 for (int j = 0; j < 3; j++)
245 Jacobian(i, j) = n_dir(i) * n_dir(j);
246 return err_r.Similarity(Jacobian);
247}
248
250{
251 if (m_mode == "vertexbtube") {
252 m_distance = getDocaBtubeVertex(p, t);
253 }
254 if (m_mode == "trackbtube") {
255 m_distance = getDocaBtubeTrack(p, t);
256 }
257}
258
260{
261 if (m_mode == "vertexbtube") {
262 m_distanceCovMatrix = getDocaBtubeVertexError(p, t);
263 }
264 if (m_mode == "trackbtube") {
265 m_distanceCovMatrix = getDocaBtubeTrackError(p, t);
266 }
267}
268
270{
271 unsigned int n = m_plist->getListSize();
272 for (unsigned int i = 0; i < n; i++) {
273 Particle* particle = m_plist->getParticle(i);
274 std::vector<const Particle*> selectedParticles = m_decayDescriptor.getSelectionParticles(particle);
275 if ((m_mode == "vertexbtube") || (m_mode == "trackbtube")) {
276 if (selectedParticles.size() == 1) {
277 const Btube* t = particle->getRelatedTo<Btube>();
278 if (!t) {
279 B2FATAL("No associated Btube found");
280 }
281 const Particle* p0 = selectedParticles[0];
282 getBtubeDistance(p0, t);
283 getBtubeDistanceErrors(particle, t);
284 }
285 } else {
286 if (selectedParticles.size() == 2) {
287 const Particle* p1 = selectedParticles[0];
288 const Particle* p2 = selectedParticles[1];
289 getDistance(p1, p2);
290 getDistanceErrors(p1, p2);
291 }
292 }
293 Eigen::Matrix3d DistanceCovMatrix_eigen;
294 DistanceCovMatrix_eigen.resize(3, 3);
295 for (int row = 0; row < 3; row++) {
296 for (int column = 0; column < 3; column++) {
297 DistanceCovMatrix_eigen(row, column) = m_distanceCovMatrix[row][column]; //converting from TMatrixFSym to Eigen
298 }
299 }
300 float DistanceMag = m_distance.norm();
301 double dist_err = 0;
302 if (DistanceMag != 0) {
303 Eigen::Vector3d dist_dir = m_distance.normalized();
304 dist_err = sqrt(((dist_dir.transpose()) * (DistanceCovMatrix_eigen * dist_dir)).norm());
305 }
306 particle->writeExtraInfo("CalculatedDistance", DistanceMag);
307 particle->writeExtraInfo("CalculatedDistanceError", dist_err);
308 //The Distance vector always points from the "vertex" to the "track/btube", if a single vertex is involved. In case of two vertices, it points from the "first" vertex to the "second" vertex. For two tracks or one track and btube, the direction is along "first track x second track" or "track x btube"
309 particle->writeExtraInfo("CalculatedDistanceVector_X", m_distance(0));
310 particle->writeExtraInfo("CalculatedDistanceVector_Y", m_distance(1));
311 particle->writeExtraInfo("CalculatedDistanceVector_Z", m_distance(2));
312 //Save the Covariance Matrix
313 particle->writeExtraInfo("CalculatedDistanceCovMatrixXX", m_distanceCovMatrix[0][0]);
314 particle->writeExtraInfo("CalculatedDistanceCovMatrixXY", m_distanceCovMatrix[0][1]);
315 particle->writeExtraInfo("CalculatedDistanceCovMatrixXZ", m_distanceCovMatrix[0][2]);
316 particle->writeExtraInfo("CalculatedDistanceCovMatrixYX", m_distanceCovMatrix[1][0]);
317 particle->writeExtraInfo("CalculatedDistanceCovMatrixYY", m_distanceCovMatrix[1][1]);
318 particle->writeExtraInfo("CalculatedDistanceCovMatrixYZ", m_distanceCovMatrix[1][2]);
319 particle->writeExtraInfo("CalculatedDistanceCovMatrixZX", m_distanceCovMatrix[2][0]);
320 particle->writeExtraInfo("CalculatedDistanceCovMatrixZY", m_distanceCovMatrix[2][1]);
321 particle->writeExtraInfo("CalculatedDistanceCovMatrixZZ", m_distanceCovMatrix[2][2]);
322 }
323}
For each MCParticle with hits in the CDC, this class stores some summarising information on those hit...
Definition: Btube.h:27
TMatrixFSym getTubeCenterErrorMatrix() const
Returns Btube Center Error Matrix.
Definition: Btube.h:85
Eigen::Matrix< double, 3, 1 > getTubeDirection() const
Returns Btube direction.
Definition: Btube.h:71
Eigen::Matrix< double, 3, 1 > getTubeCenter() const
Returns Btube center.
Definition: Btube.h:64
bool init(const std::string &str)
Initialise the DecayDescriptor from given string.
std::vector< const Particle * > getSelectionParticles(const Particle *particle)
Get a vector of pointers with selected daughters in the decay tree.
TMatrixFSym m_distanceCovMatrix
covariance matrix of distance
virtual void initialize() override
declare data store elements
virtual void event() override
process event
std::string m_listName
name of particle list
void getBtubeDistanceErrors(const Particle *p, const Btube *t)
returns the error on the distance between the btube and the object (depending on m_mode)
DistanceCalculatorModule()
Constructor: Sets the description, the properties and the parameters of the module.
void getDistanceErrors(const Particle *p1, const Particle *p2)
returns the error on the distance between tracks or vertex objects (depending on the m_mode)
Eigen::Vector3d m_distance
distance between two objects (track/vertex/Btube and vertex/track)
DecayDescriptor m_decayDescriptor
decay descriptor which specifies which particles are used to calculate the distance
void getDistance(const Particle *p1, const Particle *p2)
returns the distance between tracks or vertex objects (depending on the m_mode)
void getBtubeDistance(const Particle *p, const Btube *t)
returns the distance between the btube and the object (depending on m_mode)
StoreObjPtr< ParticleList > m_plist
input particle list
Base class for Modules.
Definition: Module.h:72
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
Class to store reconstructed particles.
Definition: Particle.h:75
double getPx() const
Returns x component of momentum.
Definition: Particle.h:587
TMatrixFSym getVertexErrorMatrix() const
Returns the 3x3 position error sub-matrix.
Definition: Particle.cc:447
double getPz() const
Returns z component of momentum.
Definition: Particle.h:605
double getX() const
Returns x component of vertex position.
Definition: Particle.h:640
void writeExtraInfo(const std::string &name, const double value)
Sets the user defined extraInfo.
Definition: Particle.cc:1308
double getPy() const
Returns y component of momentum.
Definition: Particle.h:596
double getZ() const
Returns z component of vertex position.
Definition: Particle.h:658
double getY() const
Returns y component of vertex position.
Definition: Particle.h:649
TO * getRelatedTo(const std::string &name="", const std::string &namedRelation="") const
Get the object to which this object has a relation.
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition: Module.h:560
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.