JEnergyRegressor.hh

Go to the documentation of this file.

#ifndef __JFIT__JENERGYREGRESSOR__
#define __JFIT__JENERGYREGRESSOR__
 
#include <string>
#include <vector>
#include <memory>
 
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JNPETable.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JTimeRange.hh"
#include "JPhysics/JConstants.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JGeometry3D/JAxis3D.hh"
#include "JFit/JRegressor.hh"
#include "JFit/JEnergy.hh"
#include "JFit/JNPE.hh"
#include "JFit/JNPEHit.hh"
#include "JFit/JMEstimator.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JTimeRange.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 * Data regression method for JFIT::JEnergy.
 * \author mdejong
 */
namespace JFIT {}
namespace JPP { using namespace JFIT; }
 
namespace JFIT {
 
  using JPHYSICS::JPDFType_t;
  using JPHYSICS::DIRECT_LIGHT_FROM_MUON;
  using JPHYSICS::SCATTERED_LIGHT_FROM_MUON;
  using JPHYSICS::DIRECT_LIGHT_FROM_DELTARAYS;
  using JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS;
  using JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWERS;
  using JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWERS;
  using JGEOMETRY3D::JAxis3D;
  
 
  /**
   * Template specialisation for storage of PDF tables.
   */
  template<>
  struct JRegressorStorage<JEnergy>
  {
    typedef JTOOLS::JMAPLIST<JTOOLS::JPolint1FunctionalMap,
                             JTOOLS::JPolint1FunctionalGridMap,
                             JTOOLS::JPolint1FunctionalGridMap>::maplist  JNPEMaplist_t;
    typedef JPHYSICS::JNPETable<double, double, JNPEMaplist_t>            JNPE_t;    //!< time integrated PDF
 
    static const int           NUMBER_OF_PDFS  =  6;                                 //!< Number of PDFs
 
    typedef std::vector<JNPE_t>                                           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 PDF types listed in JRegressorStorage<JEnergy>::pdf_t.
     *
     * \param  fileDescriptor     PDF file descriptor
     * \param  range              time range [ns]
     */
    JRegressorStorage(const std::string& fileDescriptor, const JTimeRange& range = JTimeRange())
    {
      using namespace std;
      using namespace JPP;
 
      typedef JSplineFunction1D_t                                     JFunction1D_t;
      typedef JPDFTable<JFunction1D_t, JNPEMaplist_t>                 JPDF_t;
 
      const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
 
      for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
 
        JPDF_t pdf;
 
        const JPDFType_t type      = pdf_t[i];
        const string     file_name = getFilename(fileDescriptor, type);
 
        pdf.load(file_name.c_str());
 
        pdf.setExceptionHandler(supervisor);
  
        if      (is_bremsstrahlung(type))
          _YB.push_back(JNPE_t(pdf, range));
        else if (is_deltarays(type))
          _YA.push_back(JNPE_t(pdf, range));
        else
          _Y1.push_back(JNPE_t(pdf, range));
      }
 
      // Add PDFs
 
      if (_Y1.size() == 2) { _Y1[1].add(_Y1[0]); _Y1.erase(_Y1.begin()); }
      if (_YA.size() == 2) { _YA[1].add(_YA[0]); _YA.erase(_YA.begin()); }
      if (_YB.size() == 2) { _YB[1].add(_YB[0]); _YB.erase(_YB.begin()); }
    }
 
 
    /**
     * Get light from muon.
     *
     * \return                     NPEs
     */
    const JNPEs_t& getY1() const
    {
      return _Y1;
    }
 
 
    /**
     * Get light from delta-rays.
     *
     * \return                     NPEs
     */
    const JNPEs_t& getYA() const
    {
      return _YA;
    }
 
 
    /**
     * Get light from EM showers.
     *
     * \return                     NPEs
     */
    const JNPEs_t& getYB() const
    {
      return _YB;
    }
 
 
    /**
     * PDF types.
     */
    static const JPDFType_t pdf_t[NUMBER_OF_PDFS];
 
  private:
    JNPEs_t _Y1;                        //!< light from muon
    JNPEs_t _YA;                        //!< light from delta-rays
    JNPEs_t _YB;                        //!< light from EM showers
  };
 
 
  /**
   * PDF types.
   */
  const JPDFType_t JRegressorStorage<JEnergy>::pdf_t[] = { DIRECT_LIGHT_FROM_MUON,
                                                           SCATTERED_LIGHT_FROM_MUON,
                                                           DIRECT_LIGHT_FROM_DELTARAYS,
                                                           SCATTERED_LIGHT_FROM_DELTARAYS,
                                                           DIRECT_LIGHT_FROM_EMSHOWERS,
                                                           SCATTERED_LIGHT_FROM_EMSHOWERS };
 
 
 
  /**
   * Regressor function object for fit of muon energy.
   */
  template<>
  struct JRegressor<JEnergy> :
    public JAbstractRegressor<JEnergy>,
    public JRegressorStorage <JEnergy>
  {
    using JAbstractRegressor<JEnergy>::operator();
 
    typedef JRegressorStorage<JEnergy>                                    storage_type;
 
    using storage_type::getY1;
    using storage_type::getYA;
    using storage_type::getYB;
    
 
    /**
     * Constructor.
     *
     * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which
     * will be replaced by the PDF types listed in JRegressor<JEnergy>::pdf_t.
     *
     * \param  fileDescriptor     PDF file descriptor
     * \param  range              time range [ns]
     */
    JRegressor(const std::string& fileDescriptor, const JTimeRange& range = JTimeRange()) :
      storage_type(fileDescriptor, range),
      Y1(getY1()),
      YA(getYA()),
      YB(getYB()),
      estimator(new JMEstimatorNormal()) 
    {}
 
 
    /**
     * Constructor.
     *
     * \param  storage            PDF storage
     */
    JRegressor(const storage_type& storage) :
      Y1(storage.getY1()),
      YA(storage.getYA()),
      YB(storage.getYB()),
      estimator(new JMEstimatorNormal()) 
    {}
 
 
    /**
     * Fit function.
     * This method is used to determine chi2 of given number of photo-electrons for given energy of muon.
     *
     * \param  x                  energy
     * \param  npe                npe 
     * \return                    chi2
     */
    inline double operator()(const JEnergy& x, const JNPEHit& npe) const
    {
      const double E = x.getE();
      const double u = npe.getChi2(E);
        
      return estimator->getRho(u);
    }
 
    
    /**
     * Get number of photo-electrons from muon.
     * Note that the PMT geometry should be relative to the muon trajectory,
     * conform method JGEOMETRY3D::JAxis3D::transform.
     *
     * \param  axis               PMT axis
     * \return                    npe
     */
    inline double getY1(const JAxis3D& axis) const
    {
      return getNPE(Y1, hypot(axis.getX(), axis.getY()), axis.getTheta(), fabs(axis.getPhi()));
    }
 
    
    /**
     * Get number of photo-electrons from delta-rays.
     * Note that the PMT geometry should be relative to the muon trajectory,
     * conform method JGEOMETRY3D::JAxis3D::transform.
     *
     * \param  axis               PMT axis
     * \return                    npe
     */
    inline double getYA(const JAxis3D& axis) const
    {
      return getNPE(YA, hypot(axis.getX(), axis.getY()), axis.getTheta(), fabs(axis.getPhi()));
    }
 
    
    /**
     * Get number of photo-electrons from Bremsstrahlung.
     * Note that the PMT geometry should be relative to the muon trajectory,
     * conform method JGEOMETRY3D::JAxis3D::transform.
     *
     * \param  axis               PMT axis
     * \return                    npe
     */
    inline double getYB(const JAxis3D& axis) const
    {
      return getNPE(YB, hypot(axis.getX(), axis.getY()), axis.getTheta(), fabs(axis.getPhi()));
    }
    
 
    /**
     * Create data structure for handling light yields for PMT.
     * Note that the PMT geometry should be relative to the muon trajectory,
     * conform method JGEOMETRY3D::JAxis3D::transform.
     *
     * \param  axis               PMT axis
     * \param  R_Hz               singles rate [Hz]
     * \return                    light yields
     */
    inline JNPE getNPE(const JAxis3D& axis,
                       const double   R_Hz) const
    {
      using namespace std;
      using namespace JPP;
 
      const double x  = axis.getX();
      const double y  = axis.getY();
      const double R  = sqrt(x*x + y*y);
      const double z  = axis.getZ()  -  R / getTanThetaC();
 
      const double theta = axis.getTheta();
      const double phi   = fabs(axis.getPhi());
 
      const double y1 = getNPE(Y1, R, theta, phi);
      const double yA = getNPE(YA, R, theta, phi);
      const double yB = getNPE(YB, R, theta, phi);
 
      return JNPE(getN(T_ns, R_Hz * 1.0e-9), y1, yA, yB, z);
    }
 
 
    /**
     * Get maximal road width of NPE.
     *
     * \return                    road width [m]
     */
    inline double getRmax() const
    {
      return std::max(getRmax(Y1), std::max(getRmax(YA), getRmax(YB)));
    }
 
 
    static JTimeRange   T_ns;                      //!< Time window with respect to Cherenkov hypothesis [ns]
 
    const JNPEs_t&      Y1;                        //!< PDF
    const JNPEs_t&      YA;                        //!< PDF
    const JNPEs_t&      YB;                        //!< PDF
 
    std::shared_ptr<JMEstimator> estimator;        //!< M-Estimator function
    
  protected:
    /**
     * Get maximal road width of PDF.
     *
     * \param  NPE                NPE tables
     * \return                    road width [m]
     */
    static inline double getRmax(const JNPEs_t& NPE)
    {
      double xmax = 0.0;
 
      for (JNPEs_t::const_iterator i = NPE.begin(); i != NPE.end(); ++i) {
        if (!i->empty() && i->getXmax() > xmax) {
          xmax = i->getXmax();
        }
      }
 
      return xmax;
    }
 
 
    /**
     * Get number of photo-electrons.
     *
     * \param  NPE                NPE tables
     * \param  R                  distance between muon and PMT [m]
     * \param  theta              zenith  angle orientation PMT [rad]
     * \param  phi                azimuth angle orientation PMT [rad]
     * \return                    number of photo-electrons
     */
    static inline double getNPE(const JNPEs_t& NPE,
                                const double R,
                                const double theta,
                                const double phi)
    {
      using JTOOLS::get_value;
 
      double npe = 0.0;
 
      for (JNPEs_t::const_iterator i = NPE.begin(); i != NPE.end(); ++i) {
 
        if (R <= i->getXmax()) {
 
          try {
 
            const double y = get_value((*i)(std::max(R, i->getXmin()), theta, phi));
 
            if (y > 0.0) {
              npe += y;
            }
          }
          catch(const JLANG::JException& error) {
            ERROR(error << std::endl);
          }
        }
      }
 
      return npe;
    }    
  };
 
 
  /**
   * Time range.
   */
  JTimeRange JRegressor<JEnergy>::T_ns;
}
 
#endif