11 #include <analysis/modules/DistanceCalculator/DistanceCalculatorModule.h>
12 #include <analysis/dataobjects/Particle.h>
13 #include <analysis/dataobjects/ParticleList.h>
14 #include <analysis/dataobjects/Btube.h>
18 #include <TMatrixFSym.h>
22 #include <Eigen/Geometry>
27 const float eps = 1E-6;
41 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");
44 addParam(
"listName", m_listName,
"", std::string(
""));
45 addParam(
"decayString", m_decayString,
"", std::string(
""));
46 addParam(
"mode", m_mode,
47 "Specifies how the distance is calculated:\n"
48 "vertextrack: calculate the distance of closest approach between a track and a vertex, taking the first candidate as vertex, default\n"
49 "trackvertex: calculate the distance of closest approach between a track and a vertex, taking the first candidate as track,\n"
50 "2tracks: calculates the distance of closest approach between two tracks,\n"
51 "2vertices: calculates the distance between two vertices,\n"
52 "vertexbtube: calculates the distance of closest approach between a vertex and a Btube,\n"
53 "trackbtube: calculates the distance of closest approach between a track and a Btube",
54 std::string(
"vertextrack"));
57 DistanceCalculatorModule::~DistanceCalculatorModule()
61 void DistanceCalculatorModule::initialize()
63 if (m_decayString !=
"") {
64 m_decayDescriptor.init(m_decayString);
66 B2FATAL(
"Please provide a decay string for DistanceCalculatorModule");
70 Eigen::Vector3d getDocaTrackVertex(
const Particle* pTrack,
const Particle* pVertex)
73 Eigen::Vector3d p_ray(pTrack->
getX(), pTrack->
getY(), pTrack->
getZ());
74 Eigen::Vector3d p(pVertex->
getX(), pVertex->
getY(), pVertex->
getZ());
75 Eigen::Vector3d v_dir = v.normalized();
76 Eigen::Vector3d r = p - p_ray;
77 return r - r.dot(v_dir) * v_dir;
80 TMatrixFSym getDocaTrackVertexError(
const Particle* pTrack,
const Particle* pVertex)
85 Eigen::Vector3d v = V.normalized();
90 TMatrixFSym Jacobian(3);
93 for (
int i = 0; i < 3; i++)
94 for (
int j = 0; j < 3; j++)
95 Jacobian(i, j) = -v(i) * v(j);
96 for (
int i = 0; i < 3; i++)
99 return err_r.Similarity(Jacobian);
111 Eigen::Vector3d
v1(p1->getPx(), p1->getPy(), p1->getPz());
112 Eigen::Vector3d
v2(p2->getPx(), p2->getPy(), p2->getPz());
113 Eigen::Vector3d n =
v1.cross(v2);
114 if (n.norm() < eps) {
117 return getDocaTrackVertex(p1, p2);
119 Eigen::Vector3d n_dir = n.normalized();
120 Eigen::Vector3d p1v(p1->getX(), p1->getY(), p1->getZ());
121 Eigen::Vector3d p2v(p2->getX(), p2->getY(), p2->getZ());
122 Eigen::Vector3d r = p1v - p2v;
123 return r.dot(n_dir) * n_dir;
129 TMatrixFSym err_r = p1->getVertexErrorMatrix() + p2->getVertexErrorMatrix();
130 Eigen::Vector3d p1p(p1->getPx(), p1->getPy(), p1->getPz());
131 Eigen::Vector3d p2p(p2->getPx(), p2->getPy(), p2->getPz());
132 Eigen::Vector3d n = p1p.cross(p2p);
133 if (n.norm() < eps) {
134 return getDocaTrackVertexError(p1, p2);
136 Eigen::Vector3d n_dir = n.normalized();
141 TMatrixFSym Jacobian(3);
144 for (
int i = 0; i < 3; i++)
145 for (
int j = 0; j < 3; j++)
146 Jacobian(i, j) = n_dir(i) * n_dir(j);
147 return err_r.Similarity(Jacobian);
152 Eigen::Vector3d p1v(p1->getX(), p1->getY(), p1->getZ());
153 Eigen::Vector3d p2v(p2->getX(), p2->getY(), p2->getZ());
154 Eigen::Vector3d r = p1v - p2v;
157 TMatrixFSym getDistanceVerticesErrors(
const Particle* p1,
const Particle* p2)
159 TMatrixFSym err_r = p1->getVertexErrorMatrix() + p2->getVertexErrorMatrix();
163 Eigen::Vector3d getDistance(
const Particle* p1,
const Particle* p2,
const std::string& mode)
165 if (mode ==
"2tracks") {
166 return getDocaTracks(p1, p2);
168 if (mode ==
"2vertices") {
169 return getDistanceVertices(p1, p2);
171 if (mode ==
"trackvertex") {
172 return getDocaTrackVertex(p2, p1);
175 return getDocaTrackVertex(p1, p2);
177 TMatrixFSym getDistanceErrors(
const Particle* p1,
const Particle* p2,
const std::string& mode)
179 if (mode ==
"2tracks") {
180 return getDocaTracksError(p1, p2);
182 if (mode ==
"2vertices") {
183 return getDistanceVerticesErrors(p1, p2);
185 if (mode ==
"trackvertex") {
186 return getDocaTrackVertexError(p2, p1);
189 return getDocaTrackVertexError(p1, p2);
192 Eigen::Vector3d getDocaBtubeVertex(
const Particle* pVertex,
const Btube* tube)
196 Eigen::Vector3d p(pVertex->
getX(), pVertex->
getY(), pVertex->
getZ());
197 Eigen::Vector3d v_dir = v.normalized();
198 Eigen::Vector3d r = p - p_ray;
199 return r - r.dot(v_dir) * v_dir;
202 TMatrixFSym getDocaBtubeVertexError(
const Particle* pVertex,
const Btube* tube)
206 Eigen::Vector3d v = V.normalized();
207 TMatrixFSym Jacobian(3);
208 for (
int i = 0; i < 3; i++)
209 for (
int j = 0; j < 3; j++)
210 Jacobian(i, j) = -v(i) * v(j);
211 for (
int i = 0; i < 3; i++)
214 return err_r.Similarity(Jacobian);
217 Eigen::Vector3d getDocaBtubeTrack(
const Particle* pTrack,
const Btube* tube)
221 Eigen::Vector3d n =
v1.cross(v2);
222 if (n.norm() < eps) {
225 return getDocaBtubeVertex(pTrack, tube);
227 Eigen::Vector3d n_dir = n.normalized();
228 Eigen::Vector3d pv(pTrack->
getX(), pTrack->
getY(), pTrack->
getZ());
230 Eigen::Vector3d r = pv - p2v;
231 return r.dot(n_dir) * n_dir;
234 TMatrixFSym getDocaBtubeTrackError(
const Particle* pTrack,
const Btube* tube)
237 Eigen::Vector3d p1p(pTrack->
getPx(), pTrack->
getPy(), pTrack->
getPz());
239 Eigen::Vector3d n = p1p.cross(p2p);
240 if (n.norm() < eps) {
241 return getDocaBtubeVertexError(pTrack, tube);
243 Eigen::Vector3d n_dir = n.normalized();
244 TMatrixFSym Jacobian(3);
245 for (
int i = 0; i < 3; i++)
246 for (
int j = 0; j < 3; j++)
247 Jacobian(i, j) = n_dir(i) * n_dir(j);
248 return err_r.Similarity(Jacobian);
250 Eigen::Vector3d getBtubeDistance(
const Particle* p,
const Btube* t,
const std::string& mode)
252 if (mode ==
"vertexbtube") {
253 return getDocaBtubeVertex(p, t);
256 return getDocaBtubeTrack(p, t);
259 TMatrixFSym getBtubeDistanceErrors(
const Particle* p,
const Btube* t,
const std::string& mode)
261 if (mode ==
"vertexbtube") {
262 return getDocaBtubeVertexError(p, t);
265 return getDocaBtubeTrackError(p, t);
269 void DistanceCalculatorModule::event()
272 unsigned int n = plist->getListSize();
273 for (
unsigned int i = 0; i < n; i++) {
274 Particle* particle = plist->getParticle(i);
275 std::vector<const Particle*> selectedParticles = m_decayDescriptor.getSelectionParticles(particle);
276 Eigen::Vector3d Distance;
277 TMatrixFSym DistanceCovMatrix(3);
278 if ((m_mode ==
"vertexbtube") || (m_mode ==
"trackbtube")) {
279 if (selectedParticles.size() == 1) {
280 const Btube* t = particle->getRelatedTo<
Btube>();
282 B2FATAL(
"No associated Btube found");
284 const Particle* p0 = selectedParticles[0];
285 Distance = getBtubeDistance(p0, t, m_mode);
286 DistanceCovMatrix = getBtubeDistanceErrors(particle, t, m_mode);
289 if (selectedParticles.size() == 2) {
290 const Particle* p1 = selectedParticles[0];
291 const Particle* p2 = selectedParticles[1];
292 Distance = getDistance(p1, p2, m_mode);
293 DistanceCovMatrix = getDistanceErrors(p1, p2, m_mode);
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) = DistanceCovMatrix[row][column];
303 float DistanceMag = Distance.norm();
305 if (DistanceMag != 0) {
306 Eigen::Vector3d dist_dir = Distance.normalized();
307 dist_err = TMath::Sqrt(((dist_dir.transpose()) * (DistanceCovMatrix_eigen * dist_dir)).norm());
309 particle->writeExtraInfo(
"CalculatedDistance", DistanceMag);
310 particle->writeExtraInfo(
"CalculatedDistanceError", dist_err);
311 particle->writeExtraInfo(
"CalculatedDistanceVector_X", Distance(0));
312 particle->writeExtraInfo(
"CalculatedDistanceVector_Y", Distance(1));
313 particle->writeExtraInfo(
"CalculatedDistanceVector_Z", Distance(2));
315 particle->writeExtraInfo(
"CalculatedDistanceCovMatrixXX", DistanceCovMatrix[0][0]);
316 particle->writeExtraInfo(
"CalculatedDistanceCovMatrixXY", DistanceCovMatrix[0][1]);
317 particle->writeExtraInfo(
"CalculatedDistanceCovMatrixXZ", DistanceCovMatrix[0][2]);
318 particle->writeExtraInfo(
"CalculatedDistanceCovMatrixYX", DistanceCovMatrix[1][0]);
319 particle->writeExtraInfo(
"CalculatedDistanceCovMatrixYY", DistanceCovMatrix[1][1]);
320 particle->writeExtraInfo(
"CalculatedDistanceCovMatrixYZ", DistanceCovMatrix[1][2]);
321 particle->writeExtraInfo(
"CalculatedDistanceCovMatrixZX", DistanceCovMatrix[2][0]);
322 particle->writeExtraInfo(
"CalculatedDistanceCovMatrixZY", DistanceCovMatrix[2][1]);
323 particle->writeExtraInfo(
"CalculatedDistanceCovMatrixZZ", DistanceCovMatrix[2][2]);