9#include <analysis/modules/DistanceCalculator/DistanceCalculatorModule.h>
11#include <analysis/dataobjects/Btube.h>
12#include <analysis/dataobjects/Particle.h>
19const float eps = 1
E-6;
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");
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"));
50 m_distance = Eigen::Matrix<double, Eigen::Dynamic, 1>::Zero(3);
62 B2FATAL(
"Please provide a decay string for DistanceCalculatorModule");
67Eigen::Vector3d getDocaTrackVertex(
const Particle* pTrack,
const Particle* pVertex)
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;
77TMatrixFSym getDocaTrackVertexError(
const Particle* pTrack,
const Particle* pVertex)
82 Eigen::Vector3d v = V.normalized();
87 TMatrixFSym Jacobian(3);
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++)
96 return err_r.Similarity(Jacobian);
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);
111 if (n.norm() < eps) {
114 return getDocaTrackVertex(p1, p2);
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;
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);
133 Eigen::Vector3d n_dir = n.normalized();
138 TMatrixFSym Jacobian(3);
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);
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;
155TMatrixFSym getDistanceVerticesErrors(
const Particle* p1,
const Particle* p2)
157 TMatrixFSym err_r = p1->getVertexErrorMatrix() + p2->getVertexErrorMatrix();
163 if (
m_mode ==
"2tracks") {
166 if (
m_mode ==
"2vertices") {
169 if (
m_mode ==
"vertextrack") {
172 if (
m_mode ==
"trackvertex") {
179 if (
m_mode ==
"2tracks") {
182 if (
m_mode ==
"2vertices") {
185 if (
m_mode ==
"vertextrack") {
188 if (
m_mode ==
"trackvertex") {
193Eigen::Vector3d getDocaBtubeVertex(
const Particle* pVertex,
const Btube* tube)
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;
203TMatrixFSym getDocaBtubeVertexError(
const Particle* pVertex,
const Btube* tube)
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++)
215 return err_r.Similarity(Jacobian);
218Eigen::Vector3d getDocaBtubeTrack(
const Particle* pTrack,
const Btube* tube)
222 Eigen::Vector3d n =
v1.cross(v2);
223 if (n.norm() < eps) {
226 return getDocaBtubeVertex(pTrack, tube);
228 Eigen::Vector3d n_dir = n.normalized();
229 Eigen::Vector3d pv(pTrack->
getX(), pTrack->
getY(), pTrack->
getZ());
231 Eigen::Vector3d r = pv - p2v;
232 return r.dot(n_dir) * n_dir;
235TMatrixFSym getDocaBtubeTrackError(
const Particle* pTrack,
const Btube* tube)
238 Eigen::Vector3d p1p(pTrack->
getPx(), pTrack->
getPy(), pTrack->
getPz());
240 Eigen::Vector3d n = p1p.cross(p2p);
241 if (n.norm() < eps) {
242 return getDocaBtubeVertexError(pTrack, tube);
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);
254 if (
m_mode ==
"vertexbtube") {
257 if (
m_mode ==
"trackbtube") {
264 if (
m_mode ==
"vertexbtube") {
267 if (
m_mode ==
"trackbtube") {
274 unsigned int n =
m_plist->getListSize();
275 for (
unsigned int i = 0; i < n; 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>();
282 B2FATAL(
"No associated Btube found");
284 const Particle* p0 = selectedParticles[0];
289 if (selectedParticles.size() == 2) {
290 const Particle* p1 = selectedParticles[0];
291 const Particle* p2 = selectedParticles[1];
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++) {
305 if (DistanceMag != 0) {
306 Eigen::Vector3d dist_dir =
m_distance.normalized();
307 dist_err =
sqrt(((dist_dir.transpose()) * (DistanceCovMatrix_eigen * dist_dir)).norm());
309 particle->writeExtraInfo(
"CalculatedDistance", DistanceMag);
310 particle->writeExtraInfo(
"CalculatedDistanceError", dist_err);
312 particle->writeExtraInfo(
"CalculatedDistanceVector_X",
m_distance(0));
313 particle->writeExtraInfo(
"CalculatedDistanceVector_Y",
m_distance(1));
314 particle->writeExtraInfo(
"CalculatedDistanceVector_Z",
m_distance(2));
For each MCParticle with hits in the CDC, this class stores some summarising information on those hit...
TMatrixFSym getTubeCenterErrorMatrix() const
Returns Btube Center Error Matrix.
Eigen::Matrix< double, 3, 1 > getTubeDirection() const
Returns Btube direction.
Eigen::Matrix< double, 3, 1 > getTubeCenter() const
Returns Btube center.
TMatrixFSym m_distanceCovMatrix
covariance matrix of distance
virtual void initialize() override
declare data store elements
virtual void event() override
process event
std::string m_decayString
decay string
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
virtual ~DistanceCalculatorModule()
Destructor.
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
std::string m_mode
option string
void setDescription(const std::string &description)
Sets the description of the module.
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Class to store reconstructed particles.
double getPx() const
Returns x component of momentum.
TMatrixFSym getVertexErrorMatrix() const
Returns the 3x3 position error sub-matrix.
double getPz() const
Returns z component of momentum.
double getX() const
Returns x component of vertex position.
double getPy() const
Returns y component of momentum.
double getZ() const
Returns z component of vertex position.
double getY() const
Returns y component of vertex position.
void addParam(const std::string &name, T ¶mVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
double sqrt(double a)
sqrt for double
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.