Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
JLine1ZEstimator.hh
Go to the documentation of this file.
1#ifndef __JFIT__JLINE1ZESTIMATOR__
2#define __JFIT__JLINE1ZESTIMATOR__
3
5#include "JMath/JMatrix3S.hh"
6#include "JMath/JSVD3D.hh"
7#include "JFit/JMatrixNZ.hh"
8#include "JFit/JLine1Z.hh"
9#include "JFit/JEstimator.hh"
10
11
12
13/**
14 * \file
15 * Linear fit of JFIT::JLine1Z.
16 * \author mdejong
17 */
18namespace JFIT {}
19namespace JPP { using namespace JFIT; }
20
21namespace JFIT {
22
23 /**
24 * Linear fit of straight line parallel to z-axis to set of hits (objects with position and time).
25 *
26 \f{center}\setlength{\unitlength}{0.7cm}\begin{picture}(8,12)
27
28 \put( 1.0, 0.5){\vector(0,1){9}}
29 \put( 1.0,10.0){\makebox(0,0){$z$}}
30 \put( 1.0, 0.0){\makebox(0,0){muon}}
31
32 \put( 1.0, 8.0){\line(1,0){6}}
33 \put( 4.0, 8.5){\makebox(0,0)[c]{$R$}}
34
35 \multiput( 1.0, 2.0)(0.5, 0.5){12}{\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)}
36 \put( 4.5, 4.5){\makebox(0,0)[l]{photon}}
37
38 \put( 1.0, 2.0){\circle*{0.2}}
39 \put( 0.5, 2.0){\makebox(0,0)[r]{$(x_{0},y_{0},z_{0},t_{0})$}}
40
41 \put( 1.0, 8.0){\circle*{0.2}}
42 \put( 0.5, 8.0){\makebox(0,0)[r]{$(x_{0},y_{0},z_{j})$}}
43
44 \put( 7.0, 8.0){\circle*{0.2}}
45 \put( 7.0, 9.0){\makebox(0,0)[c]{PMT}}
46 \put( 7.5, 8.0){\makebox(0,0)[l]{$(x_{j},y_{j},z_{j},t_{j})$}}
47
48 \put( 1.1, 2.1){
49 \put(0.0,1.00){\vector(-1,0){0.1}}
50 \qbezier(0.0,1.0)(0.25,1.0)(0.5,0.75)
51 \put(0.5,0.75){\vector(1,-1){0.1}}
52 \put(0.4,1.5){\makebox(0,0){$\theta_{c}$}}
53 }
54
55 \end{picture}
56 \f}
57
58 \f[
59 ct_j = ct_0 + (z_j - z_0) + \tan(\theta_{c}) \sqrt((x_j - x_0)^2 + (y_j - y_0)^2)
60 \f]
61 *
62 * where:
63 *
64 \f{eqnarray*}{
65 x_0 & = & \textrm{x position of track (fit parameter)} \\
66 y_0 & = & \textrm{y position of track (fit parameter)} \\
67 z_0 & = & \textrm{z position of track} \\
68 t_0 & = & \textrm{time of track at } z = z_0 \textrm{ (fit parameter)} \\
69 \\
70 c & = & \textrm{speed of light in vacuum} \\
71 \theta_{c} & = & \textrm{Cherenkov angle} \\
72 \f}
73 *
74 * Defining:
75 *
76 \f{eqnarray*}{
77 t_j' & \equiv & ct_j / \tan(\theta_{c}) - (z_j - z_0)/\tan(\theta_{c}) \\
78 t_0' & \equiv & ct_0 / \tan(\theta_{c}) \\
79 \f}
80 *
81 \f[
82 \Rightarrow (t_j' - t_0')^2 = (x_j - x_0)^2 + (y_j - y_0)^2
83 \f]
84 *
85 * The parameters \f$ \{x_0, y_0, t_0\} \f$ are estimated in the constructor of this class based on
86 * consecutive pairs of equations by which the quadratic terms in \f$ x_0 \f$, \f$ y_0 \f$ and \f$ t_0 \f$ are eliminated.
87 */
88 template<>
89 class JEstimator<JLine1Z> :
90 public JLine1Z
91 {
92 public:
93 /**
94 * Default constructor.
95 */
96 JEstimator() :
97 JLine1Z()
98 {}
99
100
101 /**
102 * Constructor.
103 *
104 * \param __begin begin of data
105 * \param __end end of data
106 */
107 template<class T>
108 JEstimator(T __begin, T __end) :
109 JLine1Z()
110 {
111 (*this)(__begin, __end);
112 }
113
114
115 /**
116 * Fit.
117 *
118 * The template argument <tt>T</tt> refers to an iterator of a data structure which should have the following member methods:
119 * - double %getX(); // [m]
120 * - double %getY(); // [m]
121 * - double %getZ(); // [m]
122 * - double %getT(); // [ns]
123 *
124 * \param __begin begin of data
125 * \param __end end of data
126 * \return this fit
127 */
128 template<class T>
129 const JEstimator<JLine1Z>& operator()(T __begin, T __end)
130 {
131 using namespace std;
132 using namespace JPP;
133
134 const int N = distance(__begin, __end);
135
136 if (N >= NUMBER_OF_PARAMETERS) {
137
138 __x = 0.0;
139 __y = 0.0;
140 __z = 0.0;
141 __t = 0.0;
142
143 double t0 = 0.0;
144
145 for (T i = __begin; i != __end; ++i) {
146 __x += i->getX();
147 __y += i->getY();
148 __z += i->getZ();
149 t0 += i->getT();
150 }
151
152 const double W = 1.0/N;
153
154 __x *= W;
155 __y *= W;
156 __z *= W;
157 t0 *= W;
158
159 V.reset();
160
161 t0 *= getSpeedOfLight();
162
163 double y0 = 0.0;
164 double y1 = 0.0;
165 double y2 = 0.0;
166
167 T j = __begin;
168
169 double xi = j->getX() - getX();
170 double yi = j->getY() - getY();
171 double ti = (j->getT() * getSpeedOfLight() - t0 - j->getZ() + getZ()) / getKappaC();
172
173 for (bool done = false; !done; ) {
174
175 if ((done = (++j == __end))) {
176 j = __begin;
177 }
178
179 double xj = j->getX() - getX();
180 double yj = j->getY() - getY();
181 double tj = (j->getT() * getSpeedOfLight() - t0 - j->getZ() + getZ()) / getKappaC();
182
183 double dx = xj - xi;
184 double dy = yj - yi;
185 double dt = ti - tj; // opposite sign!
186
187 const double y = ((xj + xi) * dx +
188 (yj + yi) * dy +
189 (tj + ti) * dt);
190
191 dx *= 2;
192 dy *= 2;
193 dt *= 2;
194
195 V.a00 += dx * dx;
196 V.a01 += dx * dy;
197 V.a02 += dx * dt;
198 V.a11 += dy * dy;
199 V.a12 += dy * dt;
200 V.a22 += dt * dt;
201
202 y0 += dx * y;
203 y1 += dy * y;
204 y2 += dt * y;
205
206 xi = xj;
207 yi = yj;
208 ti = tj;
209 }
210
211 t0 *= getInverseSpeedOfLight();
212
213 V.a10 = V.a01;
214 V.a20 = V.a02;
215 V.a21 = V.a12;
216
217 svd.decompose(V);
218
219 if (fabs(svd.S.a11) < MINIMAL_SVD_WEIGHT * fabs(svd.S.a00)) {
220 THROW(JValueOutOfRange, "JEstimator<JLine1Z>::JEstimator(): singular value " << svd.S.a11 << ' ' << svd.S.a00);
221 }
222
223 V = svd.invert(MINIMAL_SVD_WEIGHT);
224
225 __x += V.a00 * y0 + V.a01 * y1 + V.a02 * y2;
226 __y += V.a10 * y0 + V.a11 * y1 + V.a12 * y2;
227 __t = V.a20 * y0 + V.a21 * y1 + V.a22 * y2;
228
229 __t *= getKappaC() * getInverseSpeedOfLight();
230 __t += t0;
231
232 } else {
233 throw JValueOutOfRange("JEstimator<JLine1Z>::JEstimator(): Not enough data points.");
234 }
235
236 return *this;
237 }
238
239
240 /**
241 * Update track parameters using updated co-variance matrix (e.g.\ matrix with one hit switched off).
242 *
243 * In this, it is assumed that the changes in x and y positions are small
244 * compared to the actual distances between the track and the hits.
245 *
246 * The template argument <tt>T</tt> refers to an iterator of a data structure which should have the following member methods:
247 * - double %getX(); // [m]
248 * - double %getY(); // [m]
249 * - double %getZ(); // [m]
250 * - double %getT(); // [ns]
251 *
252 * \param __begin begin of data
253 * \param __end end of data
254 * \param A co-variance matrix of hits
255 */
256 template<class T>
257 void update(T __begin, T __end, const JMatrixNZ& A)
258 {
259 using namespace std;
260 using namespace JPP;
261
262 const int N = distance(__begin, __end);
263
264 if (N != (int) A.size()) {
265 THROW(JValueOutOfRange, "JEstimator<JLine1Z>::update(): Wrong number of points " << N << " != " << A.size());
266 }
267
268 if (N >= NUMBER_OF_PARAMETERS) {
269
270 double x1 = 0.0;
271 double y1 = 0.0;
272 double t1 = 0.0;
273
274 V.reset();
275
276 T i = __begin;
277
278 for (size_t row = 0; row != A.size(); ++row, ++i) {
279
280 const double dx = i->getX() - getX();
281 const double dy = i->getY() - getY();
282
283 const double rt = sqrt(dx*dx + dy*dy);
284
285 double xr = getKappaC() * getInverseSpeedOfLight();
286 double yr = getKappaC() * getInverseSpeedOfLight();
287 double tr = 1.0;
288
289 if (rt != 0.0) {
290 xr *= -dx / rt;
291 yr *= -dy / rt;
292 }
293
294 T j = __begin;
295
296 for (size_t col = 0; col != A.size(); ++col, ++j) {
297
298 const double dx = j->getX() - getX();
299 const double dy = j->getY() - getY();
300 const double dz = j->getZ() - getZ();
301
302 const double rt = sqrt(dx*dx + dy*dy);
303
304 double xc = getKappaC() * getInverseSpeedOfLight();
305 double yc = getKappaC() * getInverseSpeedOfLight();
306 double tc = 1.0;
307
308 if (rt != 0.0) {
309 xc *= -dx / rt;
310 yc *= -dy / rt;
311 }
312
313 const double ts = j->getT() - (dz + rt * getKappaC()) * getInverseSpeedOfLight();
314
315 const double vs = A(row,col);
316
317 x1 += xr * vs * ts;
318 y1 += yr * vs * ts;
319 t1 += tr * vs * ts;
320
321 V.a00 += xr * vs * xc;
322 V.a01 += xr * vs * yc;
323 V.a02 += xr * vs * tc;
324 V.a11 += yr * vs * yc;
325 V.a12 += yr * vs * tc;
326 V.a22 += tr * vs * tc;
327 }
328 }
329
330 V.a10 = V.a01;
331 V.a20 = V.a02;
332 V.a21 = V.a12;
333
334 svd.decompose(V);
335
336 if (fabs(svd.S.a11) < MINIMAL_SVD_WEIGHT * fabs(svd.S.a00)) {
337 THROW(JValueOutOfRange, "JEstimator<JLine1Z>::update(): singular value " << svd.S.a11 << ' ' << svd.S.a00);
338 }
339
340 V = svd.invert(MINIMAL_SVD_WEIGHT);
341
342 __x += V.a00 * x1 + V.a01 * y1 + V.a02 * t1;
343 __y += V.a10 * x1 + V.a11 * y1 + V.a12 * t1;
344 __t = V.a20 * x1 + V.a21 * y1 + V.a22 * t1;
345
346 } else {
347 THROW(JValueOutOfRange, "JEstimator<JLine1Z>::update(): Not enough data points " << N);
348 }
349 }
350
351
352 static const int NUMBER_OF_PARAMETERS = 3; //!< number of parameters of fit
353 JMATH::JMatrix3S V; //!< co-variance matrix of fit parameters
354 JMATH::JSVD3D svd;
355 static double MINIMAL_SVD_WEIGHT; //!< minimal SVD weight.
356 };
357
358 /**
359 * Set default minimal SVD weight.
360 */
362}
363
364#endif
Linear fit methods.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
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.
Template definition of linear fit.
Definition JEstimator.hh:25
Data structure for fit of straight line paralel to z-axis.
Definition JLine1Z.hh:29
Determination of the co-variance matrix of hits for a track along z-axis (JFIT::JLine1Z).
Definition JMatrixNZ.hh:30
Exception for accessing a value in a collection that is outside of its range.
3 x 3 symmetric matrix
Definition JMatrix3S.hh:31
Singular value decomposition.
Definition JSVD3D.hh:27
Auxiliary classes and methods for linear and iterative data regression.
Definition JEnergy.hh:15
double getKappaC()
Get average R-dependence of arrival time of Cherenkov light (a.k.a.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
int j
Definition JPolint.hh:801