Jpp  18.0.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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.
JMATH::JMatrix4S V
co-variance matrix of fit parameters
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
Data structure for vertex fit.
Definition: JPoint4D.hh:22
Template definition of linear fit.
Definition: JEstimator.hh:25
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
do set_variable OUTPUT_DIRECTORY $WORKDIR T
Physics constants.
JEstimator(T __begin, T __end)
Fit constructor.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:36
const double getSpeedOfLight()
Get speed of light.
const double getInverseSpeedOfLight()
Get inverse speed of light.
int j
Definition: JPolint.hh:703
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:162
4 x 4 symmetric matrix
Definition: JMatrix4S.hh:26
esac done
Definition: JAddHDE.sh:21