JShowerBrightPointRegressor.hh

Go to the documentation of this file.

#ifndef __JFIT__JSHOWERBRIGHTPOINTREGRESSOR__
#define __JFIT__JSHOWERBRIGHTPOINTREGRESSOR__

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

#include "JTools/JConstants.hh"
#include "JTools/JResult.hh"

#include "JMath/JZero.hh"

#include "JFit/JHitW0.hh"
#include "JFit/JGandalf.hh"
#include "JFit/JPoint4D.hh"
#include "JFit/JMEstimator.hh"
#include "JFit/JRegressor.hh"
#include "JFit/JFitToolkit.hh"

#include "Jeep/JMessage.hh"

/**
 * \file
 * Data regression method for JFIT::JPoint4D from a bright point isoptropic emission PDF.
 *
 * \author adomi
 */

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;

  /**
   * Regressor function object for JPoint4D fit using JGandalf minimiser.
   */
  template<>
  struct JRegressor<JPoint4D, JGandalf> :
    public JAbstractRegressor<JPoint4D, JGandalf>
  {
    using JAbstractRegressor<JPoint4D, JGandalf>::operator();

    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;


    /**
     * Parameterized constructor
     *
     * The PDF file descriptor should contain the wild card character JPHYSICS::WILD_CARD which
     * will be replaced by the corresponding PDF types.
     *
     * \param  fileDescriptor     PDF file descriptor
     * \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 double       TTS,
               const int          numberOfPoints = 25,
               const double       epsilon        = 1.0e-10):
      E_GeV(0.0)
    {
      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) {
        
        try {

          const string file_name = getFilename(fileDescriptor, pdf_t[i]);
          
          NOTICE("loading PDF from file " << file_name << "... " << flush);

          pdf[i].load(file_name.c_str());

          NOTICE("OK" << endl);

          pdf[i].setExceptionHandler(supervisor);
        }
        catch(const JException& error) {
          FATAL(error.what() << endl);
        }
      }

      // 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);
        } else if (TTS < 0.0) {
          ERROR("Illegal value of TTS [ns]: " << TTS << endl);
        }
      }    
    }

    /**
     * Fit function.
     * This method is used to determine the chi2 and gradient of given hit with respect a bright point emitting isotropically
     *
     * JHitW0 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
     */
    result_type operator()(const JPoint4D& vx, const JHitW0& hit) const
    {      
      using namespace JPP;

      JPosition3D D(hit.getPosition());
      JDirection3D U(hit.getDirection());

      D.sub(vx.getPosition());

      double ct = U.getDot(D) / D.getLength();
      if (ct > +1.0) { ct = +1.0; }
      if (ct < -1.0) { ct = -1.0; }

      const double t = vx.getT() + (D.getLength() * 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(D.getLength(), ct, dt, E_GeV);  // getH1 = Get signal hypothesis value

      if (get_value(H1) >= Vmax_npe) {
        H1 *= Vmax_npe / get_value(H1);
      }

      H1 += H0;  // now H1 is signal + background

      result_type result;
      
      result.chi2 = H1.getChi2() - H0.getChi2();  // Likelihood ratio

      /*
       * Here it is evaluated: d(chi2)/d(ct) * d(ct)/d(x0,y0,z0,t0)
       */
      result.gradient = JPoint4D(JVector3D(-getIndexOfRefraction() * D.getX() / D.getLength(),    // d(ct)/d(x0) 
                                           -getIndexOfRefraction() * D.getY() / D.getLength(),    // d(ct)/d(y0) 
                                           -getIndexOfRefraction() * D.getZ() / D.getLength()),   // d(ct)/d(z0) 
                                 getSpeedOfLight());                                              // d(ct)/d(t0)
      
      result.gradient.mul(getInverseSpeedOfLight() * (H1.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 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
     * \param  E                shower energy [GeV]
     * \return                  hypothesis value
     */
    JPDF_t::result_type getH1(const double D,
                              const double ct,
                              const double t,
                              const double E) 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 = E * pdf[i](std::max(D, pdf[i].getXmin()), ct, t);
          
            if (get_value(y1) > 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 JTimeRange   T_ns;                      //!< Time window with respect to Cherenkov hypothesis [ns]
    static double       Vmax_npe;                  //!< Maximal integral of PDF [npe]

    static const int    NUMBER_OF_PDFS = 2;

    static const JPDFType_t pdf_t[NUMBER_OF_PDFS];

    JPDF_t              pdf[NUMBER_OF_PDFS];       //!< PDF

    double              E_GeV;                     //!< Energy of the shower [GeV]
  };  


  /**
   * PDF types.
   */
  const JPDFType_t JRegressor<JPoint4D, JGandalf>::pdf_t[] = { DIRECT_LIGHT_FROM_BRIGHT_POINT, 
                                                               SCATTERED_LIGHT_FROM_BRIGHT_POINT };

  /**
   * Default values.
   */
  JTimeRange JRegressor<JPoint4D, JGandalf>::T_ns;
  double     JRegressor<JPoint4D, JGandalf>::Vmax_npe = std::numeric_limits<double>::max();

}

#endif