JPoint4DEstimator.hh

Go to the documentation of this file.

#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<JPoint4D> :
    public JPoint4D
  {
  public:
  
    /**
     * Fit constructor.
     */
    JEstimator() :
      JPoint4D()
    {}


    /**
      * Fit constructor.
      *
      * The template argument <tt>T</tt> 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<class T>
    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<class T>
    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<JPoint4D>::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