Jpp  master_rocky-40-g5f0272dcd
the software that should make you happy
JLine1ZEstimator.hh
Go to the documentation of this file.
1 #ifndef __JFIT__JLINE1ZESTIMATOR__
2 #define __JFIT__JLINE1ZESTIMATOR__
3 
4 #include "JPhysics/JConstants.hh"
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  */
18 namespace JFIT {}
19 namespace JPP { using namespace JFIT; }
20 
21 namespace 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  */
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
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.
Definition: JException.hh:712
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.
Linear fit of straight line parallel to z-axis to set of hits (objects with position and time).
const JEstimator< JLine1Z > & operator()(T __begin, T __end)
Fit.
JEstimator()
Default constructor.
JEstimator(T __begin, T __end)
Constructor.
void update(T __begin, T __end, const JMatrixNZ &A)
Update track parameters using updated co-variance matrix (e.g. matrix with one hit switched off).
static double MINIMAL_SVD_WEIGHT
minimal SVD weight.
JMATH::JMatrix3S V
co-variance matrix of fit parameters
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.
Definition: JException.hh:180
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.
const double getInverseSpeedOfLight()
Get inverse speed of light.
const double getSpeedOfLight()
Get speed of light.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
std::vector< double > vs
int j
Definition: JPolint.hh:792
Definition: JSTDTypes.hh:14
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:202