Jpp 19.3.0-rc.3
the software that should make you happy
Loading...
Searching...
No Matches
JPoint4DEstimator.hh
Go to the documentation of this file.
1#ifndef __JFIT__JPOINT4DESTIMATOR__
2#define __JFIT__JPOINT4DESTIMATOR__
3
5#include "JMath/JMatrix4S.hh"
6#include "JFit/JPoint4D.hh"
7#include "JFit/JEstimator.hh"
8
9/**
10 * \file
11 * Linear fit of JFIT::JPoint4D.
12 * \author mdejong
13 */
14namespace JFIT {}
15namespace JPP { using namespace JFIT; }
16
17namespace JFIT {
18
19 /**
20 * Linear fit of bright point (position and time) between hits (objects with position and time).
21 *
22 \f{center}\setlength{\unitlength}{0.6cm}\begin{picture}(12,7)
23
24 \put( 6.0, 1.0){\circle*{0.3}}
25 \put( 6.0, 0.0){\makebox(0,0)[b]{$(x_{0},y_{0},z_{0})$}}
26
27 \multiput(6.0, 1.0)(-0.5, 0.5){10}{\qbezier(0.0,0.0)(+0.1,0.35)(-0.25,0.25)\qbezier(-0.25,0.25)(-0.6,0.15)(-0.5,0.5)}
28 \put( 1.0, 6.0){\circle*{0.2}}
29 \put( 1.0, 6.5){\makebox(0,0)[b]{$(x_i,y_i,z_i,t_i)$}}
30
31 \multiput(6.0, 1.0)( 0.5, 0.5){10}{\qbezier(0.0,0.0)(0.35,-0.1)( 0.25,0.25)\qbezier( 0.25,0.25)(0.15, 0.6)( 0.5,0.5)}
32 \put(11.0, 6.0){\circle*{0.2}}
33 \put(11.0, 6.5){\makebox(0,0)[b]{$(x_j,y_j,z_j,t_j)$}}
34
35 \end{picture}
36 \f}
37 *
38 \f[
39 t_j = t_0 + \frac{c}{n} \times \sqrt((x_j - x_0)^2 + (y_j - y_0)^2 + (z_j - z_0)^2)
40 \f]
41 *
42 * where:
43 *
44 \f{eqnarray*}{
45 x_0 & = & \textrm{x position of vertex (fit parameter)} \\
46 y_0 & = & \textrm{y position of vertex (fit parameter)} \\
47 z_0 & = & \textrm{z position of vertex (fit parameter)} \\
48 t_0 & = & \textrm{time at vertex (fit parameter)} \\
49 \\
50 c & = & \textrm{speed of light (in vacuum)} \\
51 n & = & \textrm{index of refraction corresponding to the group velocity of light} \\
52 \f}
53 *
54 * Defining:
55 *
56 \f{eqnarray*}{
57 t_j' & \equiv & nct_j \\
58 t_0' & \equiv & nct_0 \\
59 \f}
60 *
61 \f[
62 \Rightarrow (t_j' - t_0')^2 = (x_j - x_0)^2 + (y_j - y_0)^2 + (z_j - z_0)^2
63 \f]
64 *
65 * The parameters \f$ \{x_0, y_0, z_0, t_0\} \f$ are estimated in the constructor of this class based on
66 * consecutive pairs of equations by which the quadratic terms in \f$ x_0 \f$, \f$ y_0 \f$, \f$ z_0 \f$ and \f$ t_0 \f$ are eliminated.
67 */
68 template<>
70 public JPoint4D
71 {
72 public:
73
74 /**
75 * Fit constructor.
76 */
78 JPoint4D()
79 {}
80
81
82 /**
83 * Fit constructor.
84 *
85 * The template argument <tt>T</tt> refers to an iterator of a data structure which should have the following member methods:
86 * - double %getX(); // [m]
87 * - double %getY(); // [m]
88 * - double %getZ(); // [m]
89 * - double %getT(); // [ns]
90 *
91 * \param __begin begin of data
92 * \param __end end of data
93 */
94 template<class T>
95 JEstimator(T __begin, T __end) :
96 JPoint4D()
97 {
98 (*this)(__begin, __end);
99 }
100
101
102 /**
103 * Fit function.
104 * This method is used to find the vertex of a given set of hits
105 *
106 * \param __begin begin of data
107 * \param __end end of data
108 */
109 template<class T>
110 void operator ()(T __begin,T __end)
111 {
112 using namespace std;
113 using namespace JPP;
114
115 const int N = distance(__begin, __end);
116
117 if (N >= NUMBER_OF_PARAMETERS) {
118
119 double t0 = 0.0;
120
121 __x = 0.0;
122 __y = 0.0;
123 __z = 0.0;
124
125 for (T i = __begin; i != __end; ++i) {
126 __x += i->getX();
127 __y += i->getY();
128 __z += i->getZ();
129 t0 += i->getT();
130 }
131
132 div(N);
133 t0 /= N;
134
135 V.reset();
136
137 t0 *= getSpeedOfLight();
138
139 double y0 = 0.0;
140 double y1 = 0.0;
141 double y2 = 0.0;
142 double y3 = 0.0;
143
144 T j = __begin;
145
146 double xi = j->getX() - getX();
147 double yi = j->getY() - getY();
148 double zi = j->getZ() - getZ();
149 double ti = (j->getT() * getSpeedOfLight() - t0) / getIndexOfRefraction();
150
151 for (bool done = false; !done; ) {
152
153 if ((done = (++j == __end))) {
154 j = __begin;
155 }
156
157 double xj = j->getX() - getX();
158 double yj = j->getY() - getY();
159 double zj = j->getZ() - getZ();
160 double tj = (j->getT() * getSpeedOfLight() - t0) / getIndexOfRefraction();
161
162 double dx = xj - xi;
163 double dy = yj - yi;
164 double dz = zj - zi;
165 double dt = ti - tj; // opposite sign!
166
167 const double y = ((xj + xi) * dx +
168 (yj + yi) * dy +
169 (zj + zi) * dz +
170 (tj + ti) * dt);
171
172 dx *= 2;
173 dy *= 2;
174 dz *= 2;
175 dt *= 2;
176
177 V.a00 += dx * dx;
178 V.a01 += dx * dy;
179 V.a02 += dx * dz;
180 V.a03 += dx * dt;
181 V.a11 += dy * dy;
182 V.a12 += dy * dz;
183 V.a13 += dy * dt;
184 V.a22 += dz * dz;
185 V.a23 += dz * dt;
186 V.a33 += dt * dt;
187
188 y0 += dx * y;
189 y1 += dy * y;
190 y2 += dz * y;
191 y3 += dt * y;
192
193 xi = xj;
194 yi = yj;
195 zi = zj;
196 ti = tj;
197 }
198
199 t0 *= getInverseSpeedOfLight();
200
201 V.a10 = V.a01;
202 V.a20 = V.a02;
203 V.a30 = V.a03;
204 V.a21 = V.a12;
205 V.a31 = V.a13;
206 V.a32 = V.a23;
207
208 V.invert();
209
210 __x += V.a00 * y0 + V.a01 * y1 + V.a02 * y2 + V.a03 * y3;
211 __y += V.a10 * y0 + V.a11 * y1 + V.a12 * y2 + V.a13 * y3;
212 __z += V.a20 * y0 + V.a21 * y1 + V.a22 * y2 + V.a23 * y3;
213 __t = V.a30 * y0 + V.a31 * y1 + V.a32 * y2 + V.a33 * y3;
214
215 __t *= getInverseSpeedOfLight() * getIndexOfRefraction();
216 __t += t0;
217
218 } else {
219 throw JValueOutOfRange("JEstimator<JPoint4D>::operator(): Not enough data points.");
220 }
221 }
222
223 static const int NUMBER_OF_PARAMETERS = 4; //!< number of parameters of fit
224 JMATH::JMatrix4S V; //!< co-variance matrix of fit parameters
225
226 };
227}
228
229#endif
Linear fit methods.
Physics constants.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
JMATH::JMatrix4S V
co-variance matrix of fit parameters
JEstimator(T __begin, T __end)
Fit constructor.
Template definition of linear fit.
Definition JEstimator.hh:25
Data structure for vertex fit.
Definition JPoint4D.hh:24
Exception for accessing a value in a collection that is outside of its range.
4 x 4 symmetric matrix
Definition JMatrix4S.hh:28
Auxiliary classes and methods for linear and iterative data regression.
Definition JEnergy.hh:15
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).