Jpp  15.0.1-rc.1-highQE
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 "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  * Default constructor.
94  */
96  JLine1Z()
97  {}
98 
99 
100  /**
101  * Constructor.
102  *
103  * \param __begin begin of data
104  * \param __end end of data
105  */
106  template<class T>
107  JEstimator(T __begin, T __end) :
108  JLine1Z()
109  {
110  (*this)(__begin, __end);
111  }
112 
113 
114  /**
115  * Fit.
116  *
117  * The template argument <tt>T</tt> refers to an iterator of a data structure which should have the following member methods:
118  * - double %getX(); // [m]
119  * - double %getY(); // [m]
120  * - double %getZ(); // [m]
121  * - double %getT(); // [ns]
122  *
123  * \param __begin begin of data
124  * \param __end end of data
125  */
126  template<class T>
127  const JEstimator<JLine1Z>& operator()(T __begin, T __end)
128  {
129  using namespace std;
130  using namespace JPP;
131 
132  const int N = distance(__begin, __end);
133 
134  if (N >= NUMBER_OF_PARAMETERS) {
135 
136  __x = 0.0;
137  __y = 0.0;
138  __z = 0.0;
139  __t = 0.0;
140 
141  double t0 = 0.0;
142 
143  for (T i = __begin; i != __end; ++i) {
144  __x += i->getX();
145  __y += i->getY();
146  __z += i->getZ();
147  t0 += i->getT();
148  }
149 
150  const double W = 1.0/N;
151 
152  __x *= W;
153  __y *= W;
154  __z *= W;
155  t0 *= W;
156 
157  V.reset();
158 
159  t0 *= getSpeedOfLight();
160 
161  double y0 = 0.0;
162  double y1 = 0.0;
163  double y2 = 0.0;
164 
165  T j = __begin;
166 
167  double xi = j->getX() - getX();
168  double yi = j->getY() - getY();
169  double ti = (j->getT() * getSpeedOfLight() - t0 - j->getZ() + getZ()) / getKappaC();
170 
171  for (bool done = false; !done; ) {
172 
173  if ((done = (++j == __end))) {
174  j = __begin;
175  }
176 
177  double xj = j->getX() - getX();
178  double yj = j->getY() - getY();
179  double tj = (j->getT() * getSpeedOfLight() - t0 - j->getZ() + getZ()) / getKappaC();
180 
181  double dx = xj - xi;
182  double dy = yj - yi;
183  double dt = ti - tj; // opposite sign!
184 
185  const double y = ((xj + xi) * dx +
186  (yj + yi) * dy +
187  (tj + ti) * dt);
188 
189  dx *= 2;
190  dy *= 2;
191  dt *= 2;
192 
193  V.a00 += dx * dx;
194  V.a01 += dx * dy;
195  V.a02 += dx * dt;
196  V.a11 += dy * dy;
197  V.a12 += dy * dt;
198  V.a22 += dt * dt;
199 
200  y0 += dx * y;
201  y1 += dy * y;
202  y2 += dt * y;
203 
204  xi = xj;
205  yi = yj;
206  ti = tj;
207  }
208 
209  t0 *= getInverseSpeedOfLight();
210 
211  V.a10 = V.a01;
212  V.a20 = V.a02;
213  V.a21 = V.a12;
214 
215  V.invert();
216 
217  __x += V.a00 * y0 + V.a01 * y1 + V.a02 * y2;
218  __y += V.a10 * y0 + V.a11 * y1 + V.a12 * y2;
219  __t = V.a20 * y0 + V.a21 * y1 + V.a22 * y2;
220 
221  __t *= getKappaC() * getInverseSpeedOfLight();
222  __t += t0;
223 
224  } else {
225  throw JValueOutOfRange("JEstimator<JLine1Z>::JEstimator(): Not enough data points.");
226  }
227 
228  return *this;
229  }
230 
231 
232  /**
233  * Update track parameters using updated co-variance matrix (e.g.\ matrix with one hit switched off).
234  *
235  * In this, it is assumed that the changes in x and y positions are small
236  * compared to the actual distances between the track and the hits.
237  *
238  * The template argument <tt>T</tt> refers to an iterator of a data structure which should have the following member methods:
239  * - double %getX(); // [m]
240  * - double %getY(); // [m]
241  * - double %getZ(); // [m]
242  * - double %getT(); // [ns]
243  *
244  * \param __begin begin of data
245  * \param __end end of data
246  * \param A co-variance matrix of hits
247  */
248  template<class T>
249  void update(T __begin, T __end, const JMatrixNZ& A)
250  {
251  using namespace std;
252  using namespace JPP;
253 
254  const int N = distance(__begin, __end);
255 
256  if (N != (int) A.size()) {
257  THROW(JValueOutOfRange, "JEstimator<JLine1Z>::update(): Wrong number of points " << N << " != " << A.size());
258  }
259 
260  if (N >= NUMBER_OF_PARAMETERS) {
261 
262  double x1 = 0.0;
263  double y1 = 0.0;
264  double t1 = 0.0;
265 
266  V.reset();
267 
268  T i = __begin;
269 
270  for (size_t row = 0; row != A.size(); ++row, ++i) {
271 
272  const double dx = i->getX() - getX();
273  const double dy = i->getY() - getY();
274 
275  const double rt = sqrt(dx*dx + dy*dy);
276 
277  double xr = getKappaC() * getInverseSpeedOfLight();
278  double yr = getKappaC() * getInverseSpeedOfLight();
279  double tr = 1.0;
280 
281  if (rt != 0.0) {
282  xr *= -dx / rt;
283  yr *= -dy / rt;
284  }
285 
286  T j = __begin;
287 
288  for (size_t col = 0; col != A.size(); ++col, ++j) {
289 
290  const double dx = j->getX() - getX();
291  const double dy = j->getY() - getY();
292  const double dz = j->getZ() - getZ();
293 
294  const double rt = sqrt(dx*dx + dy*dy);
295 
296  double xc = getKappaC() * getInverseSpeedOfLight();
297  double yc = getKappaC() * getInverseSpeedOfLight();
298  double tc = 1.0;
299 
300  if (rt != 0.0) {
301  xc *= -dx / rt;
302  yc *= -dy / rt;
303  }
304 
305  const double ts = j->getT() - (dz + rt * getKappaC()) * getInverseSpeedOfLight();
306 
307  const double vs = A(row,col);
308 
309  x1 += xr * vs * ts;
310  y1 += yr * vs * ts;
311  t1 += tr * vs * ts;
312 
313  V.a00 += xr * vs * xc;
314  V.a01 += xr * vs * yc;
315  V.a02 += xr * vs * tc;
316  V.a11 += yr * vs * yc;
317  V.a12 += yr * vs * tc;
318  V.a22 += tr * vs * tc;
319  }
320  }
321 
322  V.a10 = V.a01;
323  V.a20 = V.a02;
324  V.a21 = V.a12;
325 
326  V.invert();
327 
328  __x += V.a00 * x1 + V.a01 * y1 + V.a02 * t1;
329  __y += V.a10 * x1 + V.a11 * y1 + V.a12 * t1;
330  __t = V.a20 * x1 + V.a21 * y1 + V.a22 * t1;
331 
332  } else {
333  THROW(JValueOutOfRange, "JEstimator<JLine1Z>::update(): Not enough data points " << N);
334  }
335  }
336 
337 
338  static const int NUMBER_OF_PARAMETERS = 3; //!< number of parameters of fit
339  JMATH::JMatrix3S V; //!< co-variance matrix of fit parameters
340  };
341 }
342 
343 #endif
Linear fit methods.
Linear fit of straight line parallel to z-axis to set of hits (objects with position and time)...
double getKappaC()
Get average R-dependence of arrival time of Cherenkov light (a.k.a.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:670
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
JEstimator()
Default constructor.
Template definition of linear fit.
Definition: JEstimator.hh:25
JEstimator(T __begin, T __end)
Constructor.
3 x 3 symmetric matrix
Definition: JMatrix3S.hh:26
Determination of the co-variance matrix of hits for a track along z-axis (JFIT::JLine1Z).
Definition: JMatrixNZ.hh:28
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)...
do set_variable OUTPUT_DIRECTORY $WORKDIR T
const JEstimator< JLine1Z > & operator()(T __begin, T __end)
Fit.
Physics constants.
JMATH::JMatrix3S V
co-variance matrix of fit parameters
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:157
const double getSpeedOfLight()
Get speed of light.
const double getInverseSpeedOfLight()
Get inverse speed of light.
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
int j
Definition: JPolint.hh:666
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:162
set_variable NUMBER_OF_ITERATIONS set_variable EPSILON cat acoustics_fit_parameters txt<< EOF $CONFIGURATION[*]Nmin=3;sigma_s=100.0e-6;stdev=10.0;mestimator=0;fixStrings=0;EOF for STRING in $STRINGS[*];do#fit stretching and(z) position of given string set_variable DETECTOR_TMP/tmp/detector_A.datx JEditDetector-a $DETECTOR-o $DETECTOR_TMP-r $STRING JEditDetector-a $DETECTOR-o $DETECTOR-k $STRING for MUL in 0.005 0.001;do DX_M=0.2 for((N=0;$N< $NUMBER_OF_ITERATIONS;++N));do CHI2[3]=$CHI2[1] fitPositionOfString $STRING Z $DX_M fitStretchingOfString $STRING $MUL if(($CHI2[3]-$CHI2[1]< $EPSILON));then break fi done if(($N >=$NUMBER_OF_ITERATIONS));then printf"warning: reached maximum number of iterations %d - converenge %7.3f\n"$N $(($CHI2[3]-$CHI2[1])) fi done JMergeDetector-a $DETECTOR-a $DETECTOR_TMP-o $DETECTOR rm-f $DETECTOR_TMP JConvertDetectorFormat-a $DETECTOR-o $DETECTOR-r-d 0 > &dev null done
source $JPP_DIR setenv csh $JPP_DIR eval JShellParser o a A