JPoint4DRegressor.hh

Go to the documentation of this file.

#ifndef __JPOINT4DREGRESSOR__
#define __JPOINT4DREGRESSOR__ 

#include <string>
#include <iostream>
#include <cmath>

#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JNPETable.hh"

#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JConstants.hh"

#include "JTrigger/JHitL0.hh"

#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JVersor3D.hh"

#include "JMath/JZero.hh"

#include "JFit/JTimeRange.hh"
#include "JFit/JPMTW0.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JGandalf.hh"
#include "JFit/JMEstimator.hh"
#include "JFit/JRegressor.hh"
#include <JFit/JPoint4D.hh>
#include <JFit/JShower3Z.hh>

#include "Jeep/JMessage.hh"

#include <TH2D.h>

/**
 * \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(const bool __prefit):
      estimator(new JMEstimatorTukey(30.0)), 
      prefit(__prefit)
    {}


    /**
     * Constructor for a 2D histogram as PDF      
     */
    JRegressor(TH2D *hPDF2D, bool __prefit):
      estimator(new JMEstimatorTukey(30.0)),
      prefit(__prefit)
    {
      hpdf2d[0] = (TH2D*)hPDF2D->Clone();
      hpdf2d[0]->Smooth();
    }
    

    /*  Fit Function
     *  This method is used to determine the chi2 of given hit with respect
     *  to Shower's vertex (brightest point) 
     * 
     *  \param vertex       Shower's vertex      
     *  \param hit          L0 hit
     *  \return             -log(likelihood)
     */
    template<class JHit_t>
    double operator()(const JPoint4D& vertex, const JHit_t& hit) const
    {
      using namespace std;
      using namespace JGEOMETRY3D;
      using namespace JTOOLS;
      
      const JPosition3D hit_pos(hit.getPosition());
      const double D = hit_pos.getDistance(vertex);
      const double t_res = hit.getT() - vertex.getT(hit_pos);
      const JDirection3D photonDir(hit_pos - vertex.getPosition());
      const double cosT = photonDir.getDot(hit.getDirection());
      const double cosT_th = d_ref / sqrt(d_ref*d_ref + D*D);

      double LogLik = 0;
      
      if(prefit){
        
        LogLik = evaluatePrefit(t_res, cosT, cosT_th);
        
      } else {
        
        LogLik = evaluateFit(t_res, D);
        
      }

      return LogLik;
    }


    double evaluatePrefit(const double t_res, const double cosT, const double cosT_th) const
    {

      double p_T = 0;

      if(cosT < 0){
        p_T = 0;
      } else if(cosT > cosT_th){
        p_T = p_max;
      } else {
        p_T = p_max * cosT/cosT_th;
      }
      
      double lik = 1 / (sqrt(a*a + t_res*t_res) + p_T);
      double u   = -log(lik);
      
      return estimator->getRho(u);      
    }


    double evaluateFit(const double t_res, const double D) const
    {
     
      int x_bin = hpdf2d[0]->GetXaxis()->FindBin(t_res);
      int y_bin = hpdf2d[0]->GetYaxis()->FindBin(D);
      
      int nbinsX = hpdf2d[0]->GetXaxis()->GetNbins();
      int nbinsY = hpdf2d[0]->GetYaxis()->GetNbins();
      
      double lik = 0;
      
      if(x_bin > 1 && x_bin < nbinsX){
        if(y_bin > 1 && y_bin < nbinsY){
          lik = hpdf2d[0]->Interpolate(t_res, D);
        }
      }
      
      if(lik == 0){
        lik = 0.000004;  // value of K40 bg                               
      }

      double u   = -log(lik);
      return estimator->getRho(u);
      
    }

    static const double a;  //!< constant to prevent the likelihood going to infinity
    static const double d_ref;
    static const double p_max;
    
    static const int    NUMBER_OF_PDFS = 1;
    TH2D               *hpdf2d[NUMBER_OF_PDFS];    //!< PDF 2D

    JLANG::JSharedPointer<JMEstimator> estimator;  //!< M-Estimator function    
    const bool prefit;
  };
  
  const double JRegressor<JPoint4D, JSimplex>::a = 2.0;  //!< constant to prevent the likelihood going to infinity
  const double JRegressor<JPoint4D, JSimplex>::d_ref = 10.0;
  const double JRegressor<JPoint4D, JSimplex>::p_max = 100;

}

#endif