1#ifndef __JFIT__JLINE1ZESTIMATOR__
2#define __JFIT__JLINE1ZESTIMATOR__
19namespace JPP {
using namespace JFIT; }
111 (*this)(__begin, __end);
134 const int N =
distance(__begin, __end);
136 if (N >= NUMBER_OF_PARAMETERS) {
145 for (T i = __begin; i != __end; ++i) {
152 const double W = 1.0/N;
161 t0 *= getSpeedOfLight();
169 double xi =
j->getX() - getX();
170 double yi =
j->getY() - getY();
171 double ti = (
j->getT() * getSpeedOfLight() - t0 -
j->getZ() + getZ()) /
getKappaC();
173 for (
bool done =
false; !done; ) {
175 if ((done = (++j == __end))) {
179 double xj =
j->getX() - getX();
180 double yj =
j->getY() - getY();
181 double tj = (
j->getT() * getSpeedOfLight() - t0 -
j->getZ() + getZ()) /
getKappaC();
187 const double y = ((xj + xi) * dx +
211 t0 *= getInverseSpeedOfLight();
219 if (fabs(svd.S.a11) < MINIMAL_SVD_WEIGHT * fabs(svd.S.a00)) {
220 THROW(
JValueOutOfRange,
"JEstimator<JLine1Z>::JEstimator(): singular value " << svd.S.a11 <<
' ' << svd.S.a00);
223 V = svd.invert(MINIMAL_SVD_WEIGHT);
225 __x += V.a00 * y0 + V.a01 * y1 + V.a02 * y2;
226 __y += V.a10 * y0 + V.a11 * y1 + V.a12 * y2;
227 __t = V.a20 * y0 + V.a21 * y1 + V.a22 * y2;
229 __t *=
getKappaC() * getInverseSpeedOfLight();
233 throw JValueOutOfRange(
"JEstimator<JLine1Z>::JEstimator(): Not enough data points.");
257 void update(T __begin, T __end,
const JMatrixNZ& A)
262 const int N =
distance(__begin, __end);
264 if (N != (
int) A.size()) {
265 THROW(
JValueOutOfRange,
"JEstimator<JLine1Z>::update(): Wrong number of points " << N <<
" != " << A.size());
268 if (N >= NUMBER_OF_PARAMETERS) {
278 for (
size_t row = 0; row != A.size(); ++row, ++i) {
280 const double dx = i->getX() - getX();
281 const double dy = i->getY() - getY();
283 const double rt = sqrt(dx*dx + dy*dy);
285 double xr =
getKappaC() * getInverseSpeedOfLight();
286 double yr =
getKappaC() * getInverseSpeedOfLight();
296 for (
size_t col = 0; col != A.size(); ++col, ++
j) {
298 const double dx =
j->getX() - getX();
299 const double dy =
j->getY() - getY();
300 const double dz =
j->getZ() - getZ();
302 const double rt = sqrt(dx*dx + dy*dy);
304 double xc =
getKappaC() * getInverseSpeedOfLight();
305 double yc =
getKappaC() * getInverseSpeedOfLight();
313 const double ts =
j->getT() - (dz + rt *
getKappaC()) * getInverseSpeedOfLight();
315 const double vs = A(row,col);
321 V.a00 += xr * vs * xc;
322 V.a01 += xr * vs * yc;
323 V.a02 += xr * vs * tc;
324 V.a11 += yr * vs * yc;
325 V.a12 += yr * vs * tc;
326 V.a22 += tr * vs * tc;
336 if (fabs(svd.S.a11) < MINIMAL_SVD_WEIGHT * fabs(svd.S.a00)) {
337 THROW(
JValueOutOfRange,
"JEstimator<JLine1Z>::update(): singular value " << svd.S.a11 <<
' ' << svd.S.a00);
340 V = svd.invert(MINIMAL_SVD_WEIGHT);
342 __x += V.a00 * x1 + V.a01 * y1 + V.a02 * t1;
343 __y += V.a10 * x1 + V.a11 * y1 + V.a12 * t1;
344 __t = V.a20 * x1 + V.a21 * y1 + V.a22 * t1;
352 static const int NUMBER_OF_PARAMETERS = 3;
355 static double MINIMAL_SVD_WEIGHT;
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Template definition of linear fit.
Data structure for fit of straight line paralel to z-axis.
Determination of the co-variance matrix of hits for a track along z-axis (JFIT::JLine1Z).
Exception for accessing a value in a collection that is outside of its range.
Singular value decomposition.
Auxiliary classes and methods for linear and iterative data regression.
double getKappaC()
Get average R-dependence of arrival time of Cherenkov light (a.k.a.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).