Jpp
JPoint4DEstimator.hh
Go to the documentation of this file.
1 #ifndef __JFIT__JPOINT4DESTIMATOR__
2 #define __JFIT__JPOINT4DESTIMATOR__
3 
4 #include "JTools/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} \\
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<>
69  class JEstimator<JPoint4D> :
70  public JPoint4D
71  {
72  public:
73 
74  /**
75  * Fit constructor.
76  */
77  JEstimator() :
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  * \return chi2
109  */
110  template<class T>
111  void operator ()(T __begin,T __end)
112  {
113  using namespace std;
114  using namespace JTOOLS;
115 
116  const int N = distance(__begin, __end);
117 
118  if (N >= NUMBER_OF_PARAMETERS) {
119 
120  double t0 = 0.0;
121 
122  __x = 0.0;
123  __y = 0.0;
124  __z = 0.0;
125 
126  for (T i = __begin; i != __end; ++i) {
127  __x += i->getX();
128  __y += i->getY();
129  __z += i->getZ();
130  t0 += i->getT();
131  }
132 
133  div(N);
134  t0 /= N;
135 
136  V.reset();
137 
138  t0 *= getSpeedOfLight();
139 
140  double y0 = 0.0;
141  double y1 = 0.0;
142  double y2 = 0.0;
143  double y3 = 0.0;
144 
145  T j = __begin;
146 
147  double xi = j->getX() - getX();
148  double yi = j->getY() - getY();
149  double zi = j->getZ() - getZ();
150  double ti = (j->getT() * getSpeedOfLight() - t0) / getIndexOfRefraction();
151 
152  for (bool done = false; !done; ) {
153 
154  if ((done = (++j == __end))) {
155  j = __begin;
156  }
157 
158  double xj = j->getX() - getX();
159  double yj = j->getY() - getY();
160  double zj = j->getZ() - getZ();
161  double tj = (j->getT() * getSpeedOfLight() - t0) / getIndexOfRefraction();
162 
163  double dx = xj - xi;
164  double dy = yj - yi;
165  double dz = zj - zi;
166  double dt = ti - tj; // opposite sign!
167 
168  const double y = ((xj + xi) * dx +
169  (yj + yi) * dy +
170  (zj + zi) * dz +
171  (tj + ti) * dt);
172 
173  dx *= 2;
174  dy *= 2;
175  dz *= 2;
176  dt *= 2;
177 
178  V.a00 += dx * dx;
179  V.a01 += dx * dy;
180  V.a02 += dx * dz;
181  V.a03 += dx * dt;
182  V.a11 += dy * dy;
183  V.a12 += dy * dz;
184  V.a13 += dy * dt;
185  V.a22 += dz * dz;
186  V.a23 += dz * dt;
187  V.a33 += dt * dt;
188 
189  y0 += dx * y;
190  y1 += dy * y;
191  y2 += dz * y;
192  y3 += dt * y;
193 
194  xi = xj;
195  yi = yj;
196  zi = zj;
197  ti = tj;
198  }
199 
200  t0 *= getInverseSpeedOfLight();
201 
202  V.a10 = V.a01;
203  V.a20 = V.a02;
204  V.a30 = V.a03;
205  V.a21 = V.a12;
206  V.a31 = V.a13;
207  V.a32 = V.a23;
208 
209  V.invert();
210 
211  __x += V.a00 * y0 + V.a01 * y1 + V.a02 * y2 + V.a03 * y3;
212  __y += V.a10 * y0 + V.a11 * y1 + V.a12 * y2 + V.a13 * y3;
213  __z += V.a20 * y0 + V.a21 * y1 + V.a22 * y2 + V.a23 * y3;
214  __t = V.a30 * y0 + V.a31 * y1 + V.a32 * y2 + V.a33 * y3;
215 
216  __t *= getInverseSpeedOfLight() * getIndexOfRefraction();
217  __t += t0;
218 
219  } else {
220  throw JValueOutOfRange("JEstimator<JPoint4D>::operator(): Not enough data points.");
221  }
222  }
223 
224  static const int NUMBER_OF_PARAMETERS = 4; //!< number of parameters of fit
225  JMATH::JMatrix4S V; //!< co-variance matrix of fit parameters
226 
227  };
228 }
229 
230 #endif
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
JPoint4D.hh
JTOOLS::j
int j
Definition: JPolint.hh:634
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
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
JTOOLS::getIndexOfRefraction
double getIndexOfRefraction()
Get average index of refraction of water.
Definition: JConstants.hh:111
JMATH::JMatrix4S
4 x 4 symmetric matrix
Definition: JMatrix4S.hh:26
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
JFIT::JEstimator< JPoint4D >::V
JMATH::JMatrix4S V
co-variance matrix of fit parameters
Definition: JPoint4DEstimator.hh:225
JFIT::JPoint4D
Data structure for vertex fit.
Definition: JPoint4D.hh:22
JFIT::JEstimator< JPoint4D >::JEstimator
JEstimator(T __begin, T __end)
Fit constructor.
Definition: JPoint4DEstimator.hh:95
std
Definition: jaanetDictionary.h:36
JTOOLS
Auxiliary classes and methods for multi-dimensional interpolations and histograms.
Definition: JAbstractCollection.hh:9
JMatrix4S.hh
JFIT::JEstimator< JPoint4D >::JEstimator
JEstimator()
Fit constructor.
Definition: JPoint4DEstimator.hh:77
Generated by   doxygen 1.8.16