JShowerBrightPointRegressor.hh

Go to the documentation of this file.

#ifndef __JFIT__JSHOWERBRIGHTPOINTREGRESSOR__
#define __JFIT__JSHOWERBRIGHTPOINTREGRESSOR__
 
#include <array>
 
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JNPETable.hh"
#include "JPhysics/JConstants.hh"
 
#include "JTools/JResult.hh"
 
#include "JMath/JZero.hh"
 
#include "JFit/JGandalf.hh"
#include "JFit/JPoint4E.hh"
#include "JFit/JMEstimator.hh"
#include "JFit/JRegressor.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JTimeRange.hh"
 
#include "Jeep/JMessage.hh"
 
/**
 * \file
 * Data regression method for JFIT::JPoint4E from a bright point isoptropic emission PDF.
 *
 * \author adomi, vcarretero
 */
 
namespace JFIT {}
namespace JPP { using namespace JFIT; }
 
namespace JFIT {
 
  using JPHYSICS::JPDFType_t;
  using JPHYSICS::DIRECT_LIGHT_FROM_BRIGHT_POINT;                                   
  using JPHYSICS::SCATTERED_LIGHT_FROM_BRIGHT_POINT;
  using JTOOLS::get_value;
 
  /**
   * Function to constrain the energy during the fit, to prevent unphysical values.
   *
   * \param  value                       model (I/O)
   */
  void model(JPoint4E& value)
  {
    using namespace std;
 
    double E = max(0.0,value.getE());
    value = JPoint4E(static_cast<const JPoint4D&>(value), E);
  }
 
  /**
   * Template specialisation for storage of PDF tables.
   */
  template<>
  struct JRegressorStorage<JPoint4E, JGandalf>
  {
    typedef JTOOLS::JSplineFunction1S_t                                      JFunction1D_t;
    typedef JTOOLS::JMAPLIST<JTOOLS::JPolint2FunctionalMap,
                             JTOOLS::JPolint1FunctionalGridMap>::maplist     JPDFMapList_t;
    typedef JPHYSICS::JPDFTable<JFunction1D_t, JPDFMapList_t>                JPDF_t;
 
    static const int    NUMBER_OF_PDFS = 2;
 
    typedef std::array<JPDF_t, NUMBER_OF_PDFS>                               JPDFs_t;   //!< PDFs
    
    /**
     * Default constructor.
     */
    JRegressorStorage()
    {}
 
    /**
     * Parameterized constructor
     *
     * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which
     * will be replaced by the corresponding PDF types listed in JRegressorStorage<JPoint4E, JGandalf>::pdf_t.
     *
     * The <tt>TTS</tt> corresponds to the additional time smearing applied to the PDFs.
     *
     * \param  fileDescriptor     PDF file descriptor
     * \param  T_ns               time range [ns]
     * \param  TTS                TTS [ns]
     * \param  numberOfPoints     number of points for Gauss-Hermite integration of TTS
     * \param  epsilon            precision        for Gauss-Hermite integration of TTS
     */
    JRegressorStorage(const std::string& fileDescriptor,
                      const JTimeRange&  T_ns, 
                      const double       TTS,
                      const int          numberOfPoints = 25,
                      const double       epsilon        = 1.0e-10) :
      T_ns(T_ns)
 
    {
      using namespace std;
      using namespace JPP;
      
      const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
      
      for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
        
        const string file_name = getFilename(fileDescriptor, pdf_t[i]);
          
        _pdf[i].load(file_name.c_str());
 
        _pdf[i].setExceptionHandler(supervisor);
      }
 
      // Add PDFs
      for (int i = 1; i < NUMBER_OF_PDFS; i += 2) {
  
        _pdf[ i ].add(_pdf[i-1]); 
 
        JPDF_t buffer;
 
        _pdf[i-1].swap(buffer);
        
        if (TTS > 0.0) {
          _pdf[i].blur(TTS, numberOfPoints, epsilon);
        } 
      }    
    }
 
    /**
     * Get PDFs.
     *
     * \return                    PDFs
     */
    const JPDFs_t& getPDF() const
    {
      return _pdf;
    }
 
 
    /**
     * PDF types.
     */
    static const JPDFType_t pdf_t[NUMBER_OF_PDFS];
    JTimeRange T_ns;                    //!< Time window with respect to Cherenkov hypothesis [ns]
 
  private:
    JPDFs_t _pdf;                       //!< PDFs
  };
   
  /**
   * PDF types.
   */
  const JPDFType_t JRegressorStorage<JPoint4E, JGandalf>::pdf_t[] = {
    DIRECT_LIGHT_FROM_BRIGHT_POINT,
    SCATTERED_LIGHT_FROM_BRIGHT_POINT
  };    
 
 
  /**
   * Regressor function object for JLine3Z fit using JGandalf minimiser.
   */
  template<>
  struct JRegressor<JPoint4E, JGandalf> :
    public JAbstractRegressor<JPoint4E, JGandalf>,
    public JRegressorStorage <JPoint4E, JGandalf>
  {
    using JAbstractRegressor<JPoint4E, JGandalf>::operator();
 
    typedef JRegressorStorage<JPoint4E, JGandalf>                          storage_type;
 
    /**
     * Default constructor
     */
    JRegressor() :
      storage_type(),
      pdf(getPDF())
    {}
 
    /**
     * Constructor.
     *
     * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which
     * will be replaced by the PDF types listed in JRegressorStorage<JPoint4E, JGandalf>::pdf_t.
     *
     * The <tt>TTS</tt> corresponds to the additional time smearing applied to the PDFs.
     *
     * \param  fileDescriptor     PDF file descriptor
     * \param  T_ns               time range [ns]
     * \param  TTS                TTS [ns]
     * \param  numberOfPoints     number of points for Gauss-Hermite integration of TTS
     * \param  epsilon            precision        for Gauss-Hermite integration of TTS
     */
    JRegressor(const std::string& fileDescriptor,
               const JTimeRange&  T_ns, 
               const double       TTS,
               const int          numberOfPoints = 25,
               const double       epsilon        = 1.0e-10) :
      storage_type(fileDescriptor, T_ns, TTS, numberOfPoints, epsilon),
      pdf(getPDF())
    {}
 
    /**
     * Constructor.
     *
     * \param  storage            PDF storage
     */
    JRegressor(const storage_type& storage) :
      pdf(storage.getPDF())
    {
      T_ns = storage.T_ns;
    }
 
    /**
     * Fit function.
     * This method is used to determine the chi2 and gradient of given hit with respect a bright point emitting isotropically
     *
     * JHit_t refers to a data structure which should have the following member methods:
     *   - double getX();       // [m]
     *   - double getY();       // [m]
     *   - double getZ();       // [m]
     *   - double getDX();      // [u]
     *   - double getDY();      // [u]
     *   - double getDZ();      // [u]
     *   - double getT();       // [ns]
     *
     * \param  vx                 shower vertex
     * \param  hit                hit
     * \return                    chi2 and gradient
     */
    template<class JHit_t>
    result_type operator()(const JPoint4E& vx, const JHit_t& hit) const
    {      
      using namespace JPP;
 
      JPosition3D  D(hit.getPosition());
      JDirection3D U(hit.getDirection());
 
      D.sub(vx.getPosition());
      double length = D.getLength();
      double ct = U.getDot(D) / length;
 
      if (ct > +1.0) { ct = +1.0; }
      if (ct < -1.0) { ct = -1.0; }
 
      const double t  = vx.getT() + (length * getIndexOfRefraction() * getInverseSpeedOfLight());
 
      const double dt = T_ns.constrain(hit.getT() - t);
 
      JPDF_t::result_type H0 = getH0(hit.getR(), dt);                    // getH0 = Get background hypothesis value 
      JPDF_t::result_type H1 = getH1(length, ct, dt);  // getH1 = Get signal hypothesis value / 1 GeV
 
      if (get_value(H1) >= Vmax_npe) {
        H1 *= Vmax_npe / get_value(H1);
      }
 
      double H1_value = get_value(H1);
      double v_H1 = H1.v; //Integral from tmin to t of just H1
      double V_H1 = H1.V; //Integral from tmin to tmax of just H1
      H1 *= vx.getE();
      
      JPDF_t::result_type HT = H1+H0; //now H1 is signal + background
      double HT_value = get_value(HT);
      result_type result;
      result.chi2 = HT.getChi2() - H0.getChi2();  // Likelihood ratio
      
      double exp_V_HT = exp(-HT.V); //V is the integral from tmin to tmax of EH1+H0
 
      double energy_gradient = -1 / HT_value; //dPdE
      energy_gradient *= (H1_value -  HT_value * v_H1) * (1-exp_V_HT)  -  HT_value * exp_V_HT * V_H1; //Numerator
      energy_gradient /= (1-exp_V_HT); // Denominator
 
      /*
       * Here it is evaluated: d(chi2)/d(ct) * d(ct)/d(x0,y0,z0,t0)  + d(chi2)/dE
       */
      result.gradient = JPoint4E(JPoint4D(JVector3D(-getIndexOfRefraction() * D.getX() / length,    // d(ct)/d(x0) 
                                                    -getIndexOfRefraction() * D.getY() / length,    // d(ct)/d(y0) 
                                                    -getIndexOfRefraction() * D.getZ() / length),   // d(ct)/d(z0) 
                                          getSpeedOfLight()),                                       // d(ct)/d(t0)
                                 energy_gradient);                                                  // d(chi2)/d(E)
      
      static_cast<JPoint4D&>(result.gradient).mul(getInverseSpeedOfLight() * (HT.getDerivativeOfChi2() -
                                                                              H0.getDerivativeOfChi2()));  // x d(chi2)/d(ct1)
 
      return result;
 
    }
 
    /**
     * Get background hypothesis value for time differentiated PDF.
     *
     * \param  R_Hz               rate [Hz]
     * \param  t1                 time [ns]
     * \return                    hypothesis value
     */
    JPDF_t::result_type getH0(const double R_Hz,
                              const double t1) const
    {
      using namespace JPP;
 
      return JPDF_t::result_type(R_Hz * 1e-9, t1, T_ns);
    }
    
    /**
     * Get signal hypothesis value per 1 GeV for bright point emission PDF.
     *
     * \param  D                hit distance from shower vertex [m]
     * \param  ct               cosine of the HIT angle
     * \param  t                arrival time of the light
     * \return                  hypothesis value / GeV
     */
    JPDF_t::result_type getH1(const double D,
                              const double ct,
                              const double t) const
    {
      using namespace JPP;
 
      JPDF_t::result_type h1 = JMATH::zero;
        
      for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
 
        if (!pdf[i].empty() && D <= pdf[i].getXmax()) {
 
          try {
 
            JPDF_t::result_type y1 = pdf[i](std::max(D, pdf[i].getXmin()), ct, t);
            
            // safety measures
            
            if (y1.f <= 0.0) {
              y1.f  = 0.0;
              y1.fp = 0.0;
            }
            
            if (y1.v <= 0.0) {
              y1.v  = 0.0;
            }
            
            h1 += y1;     
          }
          catch(JLANG::JException& error) {
            ERROR(error << std::endl);
          }
        }
      }
 
      return h1;
    }
 
    /**
     * Get maximal road width of PDF.
     *
     * \return                    road width [m]
     */
    inline double getRmax() const
    {
      using namespace JPP;
 
      double xmax = 0.0;
 
      for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
 
        if (!pdf[i].empty() && pdf[i].getXmax() > xmax) {
          xmax = pdf[i].getXmax();
        }
 
      }
 
      return xmax;
    }
 
    static double       Vmax_npe;                  //!< Maximal integral of PDF [npe]
 
    const JPDFs_t&      pdf;                       //!< PDF
  };  
 
    
  /**
   * Default values.
   */
  double     JRegressor<JPoint4E, JGandalf>::Vmax_npe = std::numeric_limits<double>::max();
 
}
 
#endif