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