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