#ifndef __JFIT__JPOINT4DESTIMATOR__ #define __JFIT__JPOINT4DESTIMATOR__ #include "JPhysics/JConstants.hh" #include "JMath/JMatrix4S.hh" #include "JFit/JPoint4D.hh" #include "JFit/JEstimator.hh" /** * \file * Linear fit of JFIT::JPoint4D. * \author mdejong */ namespace JFIT {} namespace JPP { using namespace JFIT; } namespace JFIT { /** * Linear fit of bright point (position and time) between hits (objects with position and time). * \f{center}\setlength{\unitlength}{0.6cm}\begin{picture}(12,7) \put( 6.0, 1.0){\circle*{0.3}} \put( 6.0, 0.0){\makebox(0,0)[b]{$(x_{0},y_{0},z_{0})$}} \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)} \put( 1.0, 6.0){\circle*{0.2}} \put( 1.0, 6.5){\makebox(0,0)[b]{$(x_i,y_i,z_i,t_i)$}} \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)} \put(11.0, 6.0){\circle*{0.2}} \put(11.0, 6.5){\makebox(0,0)[b]{$(x_j,y_j,z_j,t_j)$}} \end{picture} \f} * \f[ t_j = t_0 + \frac{c}{n} \times \sqrt((x_j - x_0)^2 + (y_j - y_0)^2 + (z_j - z_0)^2) \f] * * where: * \f{eqnarray*}{ x_0 & = & \textrm{x position of vertex (fit parameter)} \\ y_0 & = & \textrm{y position of vertex (fit parameter)} \\ z_0 & = & \textrm{z position of vertex (fit parameter)} \\ t_0 & = & \textrm{time at vertex (fit parameter)} \\ \\ c & = & \textrm{speed of light (in vacuum)} \\ n & = & \textrm{index of refraction corresponding to the group velocity of light} \\ \f} * * Defining: * \f{eqnarray*}{ t_j' & \equiv & nct_j \\ t_0' & \equiv & nct_0 \\ \f} * \f[ \Rightarrow (t_j' - t_0')^2 = (x_j - x_0)^2 + (y_j - y_0)^2 + (z_j - z_0)^2 \f] * * The parameters \f$ \{x_0, y_0, z_0, t_0\} \f$ are estimated in the constructor of this class based on * 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. */ template<> class JEstimator : public JPoint4D { public: /** * Fit constructor. */ JEstimator() : JPoint4D() {} /** * Fit constructor. * * The template argument T refers to an iterator of a data structure which should have the following member methods: * - double %getX(); // [m] * - double %getY(); // [m] * - double %getZ(); // [m] * - double %getT(); // [ns] * * \param __begin begin of data * \param __end end of data */ template JEstimator(T __begin, T __end) : JPoint4D() { (*this)(__begin, __end); } /** * Fit function. * This method is used to find the vertex of a given set of hits * * \param __begin begin of data * \param __end end of data */ template void operator ()(T __begin,T __end) { using namespace std; using namespace JPP; const int N = distance(__begin, __end); if (N >= NUMBER_OF_PARAMETERS) { double t0 = 0.0; __x = 0.0; __y = 0.0; __z = 0.0; for (T i = __begin; i != __end; ++i) { __x += i->getX(); __y += i->getY(); __z += i->getZ(); t0 += i->getT(); } div(N); t0 /= N; V.reset(); t0 *= getSpeedOfLight(); double y0 = 0.0; double y1 = 0.0; double y2 = 0.0; double y3 = 0.0; T j = __begin; double xi = j->getX() - getX(); double yi = j->getY() - getY(); double zi = j->getZ() - getZ(); double ti = (j->getT() * getSpeedOfLight() - t0) / getIndexOfRefraction(); for (bool done = false; !done; ) { if ((done = (++j == __end))) { j = __begin; } double xj = j->getX() - getX(); double yj = j->getY() - getY(); double zj = j->getZ() - getZ(); double tj = (j->getT() * getSpeedOfLight() - t0) / getIndexOfRefraction(); double dx = xj - xi; double dy = yj - yi; double dz = zj - zi; double dt = ti - tj; // opposite sign! const double y = ((xj + xi) * dx + (yj + yi) * dy + (zj + zi) * dz + (tj + ti) * dt); dx *= 2; dy *= 2; dz *= 2; dt *= 2; V.a00 += dx * dx; V.a01 += dx * dy; V.a02 += dx * dz; V.a03 += dx * dt; V.a11 += dy * dy; V.a12 += dy * dz; V.a13 += dy * dt; V.a22 += dz * dz; V.a23 += dz * dt; V.a33 += dt * dt; y0 += dx * y; y1 += dy * y; y2 += dz * y; y3 += dt * y; xi = xj; yi = yj; zi = zj; ti = tj; } t0 *= getInverseSpeedOfLight(); V.a10 = V.a01; V.a20 = V.a02; V.a30 = V.a03; V.a21 = V.a12; V.a31 = V.a13; V.a32 = V.a23; V.invert(); __x += V.a00 * y0 + V.a01 * y1 + V.a02 * y2 + V.a03 * y3; __y += V.a10 * y0 + V.a11 * y1 + V.a12 * y2 + V.a13 * y3; __z += V.a20 * y0 + V.a21 * y1 + V.a22 * y2 + V.a23 * y3; __t = V.a30 * y0 + V.a31 * y1 + V.a32 * y2 + V.a33 * y3; __t *= getInverseSpeedOfLight() * getIndexOfRefraction(); __t += t0; } else { throw JValueOutOfRange("JEstimator::operator(): Not enough data points."); } } static const int NUMBER_OF_PARAMETERS = 4; //!< number of parameters of fit JMATH::JMatrix4S V; //!< co-variance matrix of fit parameters }; } #endif