JPoint4DRegressor.hh

Go to the documentation of this file.

#ifndef __JPOINT4DREGRESSOR__
#define __JPOINT4DREGRESSOR__ 
 
#include <memory>
#include <cmath>
 
#include "JPhysics/JConstants.hh"
 
#include "JGeometry3D/JVector3D.hh"
 
#include "JFit/JTimeRange.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JGandalf.hh"
#include "JFit/JMEstimator.hh"
#include "JFit/JRegressor.hh"
#include "JFit/JPoint4D.hh"
 
#include "Jeep/JMessage.hh"
 
/**
 * \author adomi
 */
 
namespace JFIT {
 
  /**
   * Regressor function object for JPoint4D fit using JSimplex minimiser.
   */
  template<>
  struct JRegressor<JPoint4D, JSimplex> :
    public JAbstractRegressor<JPoint4D, JSimplex>
  {
    using JAbstractRegressor<JPoint4D, JSimplex>::operator();
 
 
    /**
     * Constructor.
     */
    JRegressor(double sigma)
    {
      this->sigma = sigma;
    }
 
    /*  Fit Function
     *  This method is used to determine the chi2 of given hit with respect
     *  to Shower's vertex (brightest point) 
     * 
     *  \param vx           Shower's vertex      
     *  \param hit          hit
     *  \return             chi2
     */
    template<class JHit_t>
    double operator()(const JPoint4D& vx, const JHit_t& hit) const
    {
      const double dt = hit.getT() - vx.getT(hit.getPosition());
      
      const double u  = dt / sigma;
      
      return estimator->getRho(u) * hit.getW();
    }
    
    std::shared_ptr<JMEstimator>       estimator;  //!< M-Estimator function    
    double                             sigma;      //!< Time resolution [ns]
  };
 
 
  template<>
  struct JRegressor<JPoint4D, JGandalf> :
    public JAbstractRegressor<JPoint4D, JGandalf>
  {
    using JAbstractRegressor<JPoint4D, JGandalf>::operator();
 
    constexpr static const double DMIN = 0.1;    //!< minimal distance [m]
 
    /**
     * Constructor.
     */
    JRegressor(double sigma)
    {
      this->sigma = sigma;
    }
 
    /*  Fit Function
     *  This method is used to determine the chi2 of given hit with respect
     *  to Shower's vertex (brightest point)
     *
     *  \param vx           Shower's vertex
     *  \param hit          hit
     *  \return             chi2
     */
    template<class JHit_t>
    result_type operator()(const JPoint4D& vx, const JHit_t& hit) const
    {
      using namespace JPP;
 
      result_type result;
 
      const double dt = hit.getT() - vx.getT(hit.getPosition());
 
      const double u  = dt / sigma;
 
      const JVector3D pos = hit.getPosition() - vx.getPosition();
      const double    D   = pos.getLength();  
 
      const double    weight = getIndexOfRefraction() * getInverseSpeedOfLight() / ((D < DMIN ? DMIN : D) * sigma);
      
      result.chi2     = estimator->getRho(u) * hit.getW();
 
      result.gradient = hit.getW() * JPoint4D(pos * weight, -1.0 / sigma);
 
      result.gradient.mul(0.5 * estimator->getPsi(u));
      
      return result;
    }
 
    std::shared_ptr<JMEstimator>       estimator;  //!< M-Estimator function
    double                             sigma;      //!< Time resolution [ns]
  };
}
 
#endif