JShower3EZRegressor.hh

Go to the documentation of this file.

#ifndef __JFIT__JSHOWER3EZREGRESSOR__
#define __JFIT__JSHOWER3EZREGRESSOR__
 
#include <array>
 
#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 "JPhysics/JConstants.hh"
#include "JTools/JRange.hh"
#include "JTools/JResult.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JAbstractHistogram.hh"
 
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JVersor3D.hh"
#include "JGeometry3D/JDirection3D.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/JShower3EZ.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JLine3Z.hh"
 
#include "Jeep/JMessage.hh"
 
/**
 * \file
 * Data regression method for JFIT::JShower3EZ.
 * \author mdejong, vcarretero
 */
 
 
namespace JFIT {
 
  using JPHYSICS::JPDFType_t;
  using JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER;
  using JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWER;
  using JTOOLS::get_value;
 
  /**
   * Constrain PMT angle to [0,pi].
   *
   * \param  angle        angle [rad]
   * \return              angle [rad]
   */
  inline double getPMTAngle(const double angle)
  {
    const double epsilon = 1.0e-6;
    const JTOOLS::JRange<double> range(epsilon, JMATH::PI - epsilon);
 
    return range.constrain(fabs(angle));
  }
 
   
  /**
   * Function to constrain the versor and energy during the fit, to prevent unphysical values.
   *
   * \param  value                       model (I/O)
   */
  void model(JShower3EZ& value)
  {
    using namespace std;
    
    
    double Tx = value.getDX();
    double Ty = value.getDY();
    double E = max(0.0,value.getE());
    const double u = hypot(Tx, Ty);
 
    if (u > 1.0) {
      Tx /= u;
      Ty /= u;
    }
 
    value = JShower3EZ(static_cast<const JPoint4D&>(value), JVersor3Z(Tx,Ty), E, value.getBy());
    
  }
 
  /**
   * Template specialisation for storage of PDF tables.
   */
  template<>
  struct JRegressorStorage<JShower3EZ, JSimplex>
  {
    typedef JTOOLS::JPolint1Function1D_t                                     JFunction1D_t;
    typedef JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap,
                             JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap,
                                              JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap,
                                                               JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap> > > >        JPDFMaplist_t;
    typedef JPHYSICS::JPDFTable<JFunction1D_t, JPDFMaplist_t>                JPDF_t;
 
    typedef JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap,
                             JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap,
                                              JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap,
                                                               JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap> > > >        JNPEMaplist_t;
    typedef JPHYSICS::JNPETable<double, double, JNPEMaplist_t>               JNPE_t;
    
    
    static const int    NUMBER_OF_PDFS = 2;
 
    typedef std::array<JNPE_t, NUMBER_OF_PDFS>                            JNPEs_t;   //!< NPEs
 
    /**
     * Default constructor.
     */
    JRegressorStorage()
    {}
 
    /**
     * 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<JShower3Z, JSimplex>::pdf_t.
     *
     * \param  T_ns               time range [ns]
     * \param  fileDescriptor     PDF file descriptor
     */
 
    JRegressorStorage(const std::string& fileDescriptor, const JTimeRange&  T_ns):
      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]);
            
        JPDF_t _pdf;
        _pdf.load(file_name.c_str());
          
        _pdf.setExceptionHandler(supervisor);
        
        _npe[ i ] = JNPE_t(_pdf, T_ns);
      }
      // Add NPEs
      for (int i = 1; i < NUMBER_OF_PDFS; i += 2) {
 
        _npe[ i ].add(_npe[i-1]);
 
        JNPE_t buffer;
 
        _npe[i-1].swap(buffer);
      }    
    }
 
    
 
    /**
     * Get NPEs.
     *
     * \return                    PDFs
     */
    const JNPEs_t& getNPE() const
    {
      return _npe;
    }
 
    /**
     * PDF types.
     */
    static const JPDFType_t pdf_t[NUMBER_OF_PDFS];
    JTimeRange T_ns;                    //!< Time window with respect to Cherenkov hypothesis [ns]
 
  private:
    JNPEs_t _npe;                       //!< PDFs
  };
 
  /**
   * PDF types.
   */
  const JPDFType_t JRegressorStorage<JShower3EZ, JSimplex>::pdf_t[] = {
    DIRECT_LIGHT_FROM_EMSHOWER,
    SCATTERED_LIGHT_FROM_EMSHOWER
  };
  
  /**
   * Regressor function object for JShower3EZ fit using JGandalf minimiser.
   */
  template<>
  struct JRegressor<JShower3EZ, JSimplex> :
    public JAbstractRegressor<JShower3EZ, JSimplex>,
    public JRegressorStorage <JShower3EZ, JSimplex>
  {
    using JAbstractRegressor<JShower3EZ, JSimplex>::operator();
 
    typedef JRegressorStorage<JShower3EZ, JSimplex>                          storage_type;
 
    /**
     * Default constructor
     */
    JRegressor() :
      storage_type(),
      npe(getNPE()),
      estimator() 
    {}
 
    /**
     * Constructor.
     *
     * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which
     * will be replaced by the PDF types listed in JRegressorStorage<JShower3EZ, JSimplex>::pdf_t.
     *
     * \param  fileDescriptor     PDF file descriptor
     * \param  T_ns               time range [ns]
     */
    JRegressor(const std::string& fileDescriptor, const JTimeRange& T_ns) :
      storage_type(fileDescriptor, T_ns),
      npe(getNPE()),
      estimator(new JMEstimatorNull())
    {}
 
    /**
     * Constructor.
     *
     * \param  storage            PDF storage
     */
    JRegressor(const storage_type& storage) :
      npe(storage.getNPE()),
      estimator(new JMEstimatorNull())
    {
      T_ns = storage.T_ns;
    }
 
 
    /**
     * Fit function.
     * This method is used to determine the chi2 of given PMT with respect to shower hypothesis.
     *
     * \param  shower             shower
     * \param  pmt                pmt
     * \return                    chi2 
     */
    double operator()(const JShower3EZ& shower, const JPMTW0& pmt) const
    {
      using namespace JPP;
      
      JPosition3D  D(pmt.getPosition());
      JDirection3D U(pmt.getDirection());
 
      D.sub(shower.getPosition());
 
      const double z  = D.getDot(shower.getDirection());
      const double x  = D.getX();
      const double y  = D.getY();
      const double cd = z/D.getLength();                   // cosine angle between shower direction and PMT position
 
      U.rotate(JRotation3Z(-atan2(y,x)));                  // rotate PMT axis to x-z plane
 
      const double theta = getPMTAngle(U.getTheta());
      const double phi   = getPMTAngle(U.getPhi());
 
      JNPE_t::result_type H0 = getH0(pmt.getR());  // background hypothesis value for time integrated PDF. 
      JNPE_t::result_type H1 = getH1(D.getLength(), cd, theta, phi, shower.getE());  // signal hypothesis value for time integrated PDF. 
 
      if (get_value(H1) >= Vmax_npe) {
        H1 *= Vmax_npe / get_value(H1);
      }
 
      H1 += H0; // now H1 is signal + background
 
      const bool hit = pmt.getN() != 0;
      const double u = getChi2(get_value(H1), hit) - getChi2(get_value(H0), hit); // - log likelihood ratio
 
      return estimator->getRho(u);
    }
 
    /**
     * Get background hypothesis value for time integrated PDF.
     *
     * \param  R_Hz               rate [Hz]
     * \return                    hypothesis value
     */
    JNPE_t::result_type getH0(const double R_Hz) const
    {
      return JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength());
    }
 
    /**
     * Get signal hypothesis value for time integrated PDF.
     *
     * \param  D                  PMT distance from shower [m]
     * \param  cd                 cosine angle between shower direction and PMT position
     * \param  theta              PMT zenith  angle [deg]
     * \param  phi                PMT azimuth angle [deg]
     * \param  E                  shower energy [GeV]
     * \return                    hypothesis value
     */
    JNPE_t::result_type getH1(const double D,
                              const double cd,
                              const double theta,
                              const double phi,
                              const double E) const
    {
      JNPE_t::result_type h1 = JMATH::zero;
        
      for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
 
        if (!npe[i].empty() && D <= npe[i].getXmax()) {
 
          try {
 
            JNPE_t::result_type y1 = E * npe[i](std::max(D, npe[i].getXmin()), cd, theta, phi);
 
            // safety measure
 
            if(y1 < 0){
              y1 = 0;  
            }
            
            h1 += y1;
          
          }
          catch(JLANG::JException& error) {
            ERROR(error << std::endl);
          }
        }
      }
 
      return h1;
    }
 
 
    static double       Vmax_npe;                  //!< Maximal integral of PDF [npe]
 
    const JNPEs_t&      npe;
 
    std::shared_ptr<JMEstimator>       estimator;  //!< M-Estimator function
  };  
 
  /**
   * Template specialisation for storage of PDF tables.
   */
  template<>
  struct JRegressorStorage<JShower3EZ, JGandalf>
  {
    typedef JTOOLS::JPolint1Function1D_t                                     JFunction1D_t;
    typedef JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap,
                             JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap,
                                              JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap,
                                                               JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap> > > >        JPDFMaplist_t;
    typedef JPHYSICS::JPDFTable<JFunction1D_t, JPDFMaplist_t>                JPDF_t;
 
    typedef JTOOLS::JMapList<JTOOLS::JPolint1FunctionalMap,
                             JTOOLS::JMapList<JTOOLS::JPolint1FunctionalMapH,
                                              JTOOLS::JMapList<JTOOLS::JPolint1FunctionalGridMap,
                                                               JTOOLS::JMapList<JTOOLS::JPolint1FunctionalGridMap> > > >        JNPEMaplist_t;
    typedef JPHYSICS::JNPETable<double, double, JNPEMaplist_t>               JNPE_t;
 
 
    static const int    NUMBER_OF_PDFS = 2;
 
    typedef std::array<JNPE_t, NUMBER_OF_PDFS>                               JNPEs_t;   //!< NPEs
 
    /**
     * 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<JShower3Z, JGandalf>::pdf_t.
     *
     * \param  fileDescriptor     PDF file descriptor
     * \param  T_ns               time range [ns]
 
     */
 
    JRegressorStorage(const std::string& fileDescriptor, const JTimeRange&  T_ns):
      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]);
 
        JPDF_t _pdf;
        _pdf.load(file_name.c_str());
 
        _pdf.setExceptionHandler(supervisor);
 
        _npe[i] = JNPE_t(_pdf, T_ns);
 
      }
 
      // Add PDFs
      for (int i = 1; i < NUMBER_OF_PDFS; i += 2) {
 
        _npe[ i ].add(_npe[i-1]);
 
        JNPE_t buffer;
 
        _npe[i-1].swap(buffer);
      }
    }
 
 
    /**
     * Get NPEs.
     *
     * \return                    PDFs
     */
    const JNPEs_t& getNPE() const
    {
      return _npe;
    }
 
    /**
     * PDF types.
     */
    static const JPDFType_t pdf_t[NUMBER_OF_PDFS];
    JTimeRange T_ns;                    //!< Time window with respect to Cherenkov hypothesis [ns]
 
  private:
    JNPEs_t _npe;                       //!< PDFs
  };
 
  /**
   * PDF types.
   */
  const JPDFType_t JRegressorStorage<JShower3EZ, JGandalf>::pdf_t[] = {
    DIRECT_LIGHT_FROM_EMSHOWER,
    SCATTERED_LIGHT_FROM_EMSHOWER
  };
 
  /**
   * Regressor function object for JShower3EZ fit using JGandalf minimiser.
   */
  template<>
  struct JRegressor<JShower3EZ, JGandalf> :
    public JAbstractRegressor<JShower3EZ, JGandalf>,
    public JRegressorStorage <JShower3EZ, JGandalf>
  {
    using JAbstractRegressor<JShower3EZ, JGandalf>::operator();
 
    typedef JRegressorStorage<JShower3EZ, JGandalf>                          storage_type;
 
    /**
     * Default constructor
     */
    JRegressor() :
      storage_type(),
      npe(getNPE()),
      estimator()
    {}
 
    /**
     * Constructor.
     *
     * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which
     * will be replaced by the PDF types listed in JRegressorStorage<JShower3EZ, JSimplex>::pdf_t.
     *
     * \param  fileDescriptor     PDF file descriptor
     * \param  T_ns               time range [ns]
     */
    JRegressor(const std::string& fileDescriptor, const JTimeRange&  T_ns) :
      storage_type(fileDescriptor, T_ns),
      npe(getNPE()),
      estimator(new JMEstimatorNull())
    {}
 
    /**
     * Constructor.
     *
     * \param  storage            PDF storage
     */
    JRegressor(const storage_type& storage) :
      npe(storage.getNPE()),
      estimator(new JMEstimatorNull())
    {
      T_ns = storage.T_ns;
    }
 
 
    /**
     * Fit function.
     * This method is used to determine the chi2 of given PMT with respect to shower hypothesis.
     *
     * \param  shower             shower
     * \param  pmt                pmt
     * \return                    chi2 
     */
    result_type operator()(const JShower3EZ& shower, const JPMTW0& pmt) const
    {
      using namespace JPP;
      using namespace std;
 
      JPosition3D  D(pmt.getPosition());
      JDirection3D U(pmt.getDirection());
 
      D.sub(shower.getPosition());
 
      const double x  = D.getX();
      const double y  = D.getY();
      const double d  = D.getLength();  
      const double cd = D.getDot(shower.getDirection())/d;  // cosine angle between shower direction and PMT position
 
      U.rotate(JRotation3Z(-atan2(y,x)));                  // rotate PMT axis to x-z plane
 
      const double theta = getPMTAngle(U.getTheta());
      const double phi   = getPMTAngle(U.getPhi());
 
      JNPE_t::result_type H0 = getH0(pmt.getR());                        // background hypothesis value for time integrated PDF. 
      JNPE_t::result_type H1 = getH1(d, cd, theta, phi, shower.getE());  // signal hypothesis value for time integrated PDF. 
 
      if (get_value(H1) >= Vmax_npe) {
        H1 *= Vmax_npe / get_value(H1);
      }
 
      double signal_npe = get_value(H1);
 
      H1 += H0; // now H1 is signal + background
 
      double expectation = get_value(H1);
 
      const bool hit = pmt.getN() != 0;
      const double u = H1.getChi2(hit) - H0.getChi2(hit);
 
      result_type result;
 
      result.chi2     = estimator->getRho(u);
 
      double energy_gradient = signal_npe/shower.getE(); // dP/dE
      if(hit) energy_gradient *= -exp(-expectation)/(1-exp(-expectation)); //dchi2/d(H1), if !hit is 1
 
      result.gradient = JShower3EZ(JPoint4D(JVector3D(0,       // d(cos_th0)/d(x)
                                                      0,       // d(cos_th0)/d(y)
                                                      0),          // d(cos_th0)/d(z) 
                                            0.0),                                 // d(cos_th0)/d(t) 
                                   JVersor3Z(x/d,                // d(cos_th0)/d(dx)
                                             y/d),               // d(cos_th0)/d(dy)  
                                   energy_gradient);                                          // d(chi2)/d(E)
     
      result.gradient.mul(estimator->getPsi(u));
      static_cast<JShower3Z&>(result.gradient).mul(H1.getDerivativeOfChi2(hit) - H0.getDerivativeOfChi2(hit));                    // x d(chi2)/d(cos_th0)
      
      return result;
    }
 
    /**
     * Get background hypothesis value for time integrated PDF.
     *
     * \param  R_Hz               rate [Hz]
     * \return                    hypothesis value
     */
    JNPE_t::result_type getH0(const double R_Hz) const
    {
      return JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength(), 0.0);
    }
 
    /**
     * Get signal hypothesis value for time integrated PDF.
     *
     * \param  D                  PMT distance from shower [m]
     * \param  cd                 cosine angle between shower direction and PMT position
     * \param  theta              PMT zenith  angle [deg]
     * \param  phi                PMT azimuth angle [deg]
     * \param  E                  shower energy [GeV]
     * \return                    hypothesis value
     */
    JNPE_t::result_type getH1(const double D,
                              const double cd,
                              const double theta,
                              const double phi,
                              const double E) const
    {
      JNPE_t::result_type h1 = JMATH::zero;
        
      for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
 
        if (!npe[i].empty() && D <= npe[i].getXmax()) {
 
          try {
 
            JNPE_t::result_type y1 = E * npe[i](std::max(D, npe[i].getXmin()), cd, theta, phi);
          
            if (get_value(y1) > 0.0) {
              h1 += y1;
            }
          
          }
          catch(JLANG::JException& error) {
            ERROR(error << std::endl);
          }
        }
      }
 
      return h1;
    }
 
 
    static double       Vmax_npe;                  //!< Maximal integral of PDF [npe]
    
    const JNPEs_t&      npe;
 
    std::shared_ptr<JMEstimator>       estimator;  //!< M-Estimator function
  };  
 
 
 
  /**
   * Template specialisation for storage of PDF tables.
   */
  template<>
  struct JRegressorStorage<JShower3EZ, JAbstractMinimiser>
  {
    typedef JTOOLS::JPolint1Function1D_t                                     JFunction1D_t;
    typedef JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap,
                             JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap,
                                              JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap,
                                                               JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap> > > >        JPDFMaplist_t;
    typedef JPHYSICS::JPDFTable<JFunction1D_t, JPDFMaplist_t>                JPDF_t;
 
    typedef JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap,
                             JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap,
                                              JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap,
                                                               JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap> > > >        JNPEMaplist_t;
    typedef JPHYSICS::JNPETable<double, double, JNPEMaplist_t>               JNPE_t;
 
 
    static const int    NUMBER_OF_PDFS = 2;
 
    typedef std::array<JNPE_t, NUMBER_OF_PDFS>                               JNPEs_t;   //!< NPEs
 
    /**
     * 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<JShower3Z, JAbstractMinimiser>::pdf_t.
     *
     * \param  fileDescriptor     PDF file descriptor
     * \param  T_ns               time range [ns]
 
     */
 
    JRegressorStorage(const std::string& fileDescriptor, const JTimeRange&  T_ns):
      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]);
 
        JPDF_t _pdf;
        _pdf.load(file_name.c_str());
 
        _pdf.setExceptionHandler(supervisor);
 
        _npe[i] = JNPE_t(_pdf, T_ns);
 
      }
 
      // Add PDFs
      for (int i = 1; i < NUMBER_OF_PDFS; i += 2) {
 
        _npe[i].add(_npe[i-1]);
 
        JNPE_t buffer;
 
        _npe[i-1].swap(buffer);
      }
    }
 
 
    
    /**
     * Get NPEs.
     *
     * \return                    PDFs
     */
    const JNPEs_t& getNPE() const
    {
      return _npe;
    }
 
    /**
     * PDF types.
     */
    static const JPDFType_t pdf_t[NUMBER_OF_PDFS];
    JTimeRange T_ns;                    //!< Time window with respect to Cherenkov hypothesis [ns]
 
  private:
    JNPEs_t _npe;                       //!< PDFs
  };
  
  /**
   * PDF types.
   */
  const JPDFType_t JRegressorStorage<JShower3EZ, JAbstractMinimiser>::pdf_t[] = {
    DIRECT_LIGHT_FROM_EMSHOWER,
    SCATTERED_LIGHT_FROM_EMSHOWER
  };
 
  /**
   * Regressor function object for JShower3EZ fit using JGandalf minimiser.
   */
  template<>
  struct JRegressor<JShower3EZ, JAbstractMinimiser> :
    public JAbstractRegressor<JShower3EZ, JAbstractMinimiser>,
    public JRegressorStorage <JShower3EZ, JAbstractMinimiser>
  {
    using JAbstractRegressor<JShower3EZ, JAbstractMinimiser>::operator();
 
    typedef JRegressorStorage<JShower3EZ, JAbstractMinimiser>                          storage_type;
 
    /**
     * Default constructor
     */
    JRegressor() :
      storage_type(),
      npe(getNPE()),
      estimator()
    {}
 
    /**
     * Constructor.
     *
     * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which
     * will be replaced by the PDF types listed in JRegressorStorage<JShower3EZ, JSimplex>::pdf_t.
     *
     * \param  fileDescriptor     PDF file descriptor
     * \param  T_ns               time range [ns]
 
     */
    JRegressor(const std::string& fileDescriptor, const JTimeRange&  T_ns) :
      storage_type(fileDescriptor, T_ns),
      npe(getNPE()),
      estimator(new JMEstimatorNull())
    {}
 
    /**
     * Constructor.
     *
     * \param  storage            PDF storage
     */
    JRegressor(const storage_type& storage) :
      npe(storage.getNPE()),
      estimator(new JMEstimatorNull())
    {
      T_ns = storage.T_ns;
    }
    
    /**
     * Fit function.
     * This method is used to determine the chi2 of given PMT with respect to shower hypothesis.
     *
     * \param  shower             shower
     * \param  pmt                pmt
     * \return                    chi2 
     */
    double operator()(const JShower3EZ& shower, const JPMTW0& pmt) const
    {
      using namespace JPP;
      
      JPosition3D  D(pmt.getPosition());
      JDirection3D U(pmt.getDirection());
 
      D.sub(shower.getPosition());
 
      const double z  = D.getDot(shower.getDirection());
      const double x  = D.getX();
      const double y  = D.getY();
      const double cd = z/D.getLength();                   // cosine angle between shower direction and PMT position
 
      U.rotate(JRotation3Z(-atan2(y,x)));                  // rotate PMT axis to x-z plane
 
      const double theta = getPMTAngle(U.getTheta());
      const double phi   = getPMTAngle(U.getPhi());
 
      JNPE_t::result_type H0 = getH0(pmt.getR());                                    // background hypothesis value for time integrated PDF. 
      JNPE_t::result_type H1 = getH1(D.getLength(), cd, theta, phi, shower.getE());  // signal hypothesis value for time integrated PDF. 
 
      if (get_value(H1) >= Vmax_npe) {
        H1 *= Vmax_npe / get_value(H1);
      }
 
      H1 += H0; // now H1 is signal + background
 
      const bool hit = pmt.getN() != 0;
      const double u = getChi2(get_value(H1), hit) - getChi2(get_value(H0), hit); // -log likelihood ratio
 
      return estimator->getRho(u);
    }
 
    /**
     * Get background hypothesis value for time integrated PDF.
     *
     * \param  R_Hz               rate [Hz]
     * \return                    hypothesis value
     */
    JNPE_t::result_type getH0(const double R_Hz) const
    {
      return JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength());
    }
 
    /**
     * Get signal hypothesis value for time integrated PDF.
     *
     * \param  D                  PMT distance from shower [m]
     * \param  cd                 cosine angle between shower direction and PMT position
     * \param  theta              PMT zenith  angle [deg]
     * \param  phi                PMT azimuth angle [deg]
     * \param  E                  shower energy [GeV]
     * \return                    hypothesis value
     */
    JNPE_t::result_type getH1(const double D,
                              const double cd,
                              const double theta,
                              const double phi,
                              const double E) const
    {
      JNPE_t::result_type h1 = JMATH::zero;
        
      for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
 
        if (!npe[i].empty() && D <= npe[i].getXmax()) {
 
          try {
 
            JNPE_t::result_type y1 = E * npe[i](std::max(D, npe[i].getXmin()), cd, theta, phi);
 
            // safety measure
 
            if(y1 < 0){
              y1 = 0;  
            }
            
            h1 += y1;
          
          }
          catch(JLANG::JException& error) {
            ERROR(error << std::endl);
          }
        }
      }
 
      return h1;
    }
 
    /**
     * Get signal hypothesis value for time integrated PDF.
     *
     * \param  shower             shower
     * \param  pmt                pmt
     * \return                    hypothesis value
     */
    JNPE_t::result_type getH1(const JShower3EZ& shower, const JPMTW0& pmt) const
    {
      using namespace JPP; 
    
      JPosition3D  D(pmt.getPosition());
      JDirection3D U(pmt.getDirection());
 
 
      const double z  = D.getDot(shower.getDirection());
      const double x  = D.getX();
      const double y  = D.getY();
      const double cd = z/D.getLength();                   // cosine angle between shower direction and PMT position
 
      U.rotate(JRotation3Z(-atan2(y,x)));                  // rotate PMT axis to x-z plane
 
      const double theta = getPMTAngle(U.getTheta());
      const double phi   = getPMTAngle(U.getPhi());
 
      JNPE_t::result_type H1 = getH1(D.getLength(), cd, theta, phi, 1.0);  // signal hypothesis value for time integrated PDF. E=1 because it is linear with E.
 
      if (get_value(H1) >= Vmax_npe) {
        H1 *= Vmax_npe / get_value(H1);
      }
 
      return H1;
    }
    static double       Vmax_npe;                  //!< Maximal integral of PDF [npe]
 
    const JNPEs_t&      npe;
 
    std::shared_ptr<JMEstimator>       estimator;  //!< M-Estimator function   
  };
 
  /**
   * Default values.
   */
  double     JRegressor<JShower3EZ, JSimplex>::Vmax_npe = std::numeric_limits<double>::max();
 
  double     JRegressor<JShower3EZ, JGandalf>::Vmax_npe = std::numeric_limits<double>::max();
 
  double     JRegressor<JShower3EZ, JAbstractMinimiser>::Vmax_npe = std::numeric_limits<double>::max();
 
 
 
}
 
#endif