Jpp  19.1.0-rc.1
the software that should make you happy
JPoint4DEstimator.hh
Go to the documentation of this file.
1 #ifndef __JFIT__JPOINT4DESTIMATOR__
2 #define __JFIT__JPOINT4DESTIMATOR__
3 
4 #include "JPhysics/JConstants.hh"
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  */
14 namespace JFIT {}
15 namespace JPP { using namespace JFIT; }
16 
17 namespace 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 
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.
Definition: JException.hh:180
4 x 4 symmetric matrix
Definition: JMatrix4S.hh:28
Auxiliary classes and methods for linear and iterative data regression.
Definition: JEnergy.hh:15
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
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).
int j
Definition: JPolint.hh:792
Definition: JSTDTypes.hh:14