JEnergyRegressor.hh

Go to the documentation of this file.

#ifndef __JFIT__JENERGYREGRESSOR__
#define __JFIT__JENERGYREGRESSOR__

#include <string>

#include "JLang/JSinglePointer.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JNPETable.hh"
#include "JPhysics/JPDFToolkit.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;
  

  /**
   * Regressor function object for JEnergy fit.
   */
  template<>
  struct JRegressor<JEnergy> :
    public JAbstractRegressor<JEnergy>
  {
    using JAbstractRegressor<JEnergy>::operator();
    using JMessage< JAbstractMinimiser<JEnergy> >::debug;

    typedef JTOOLS::JMAPLIST<JTOOLS::JPolint1FunctionalMap,
                             JTOOLS::JPolint1FunctionalGridMap,
                             JTOOLS::JPolint1FunctionalGridMap>::maplist    JNPEMaplist_t;
    typedef JPHYSICS::JNPETable<double, double, JNPEMaplist_t>              JNPE_t;

    
    /**
     * Constructor.
     *
     * The PDF file descriptor should contain the wild card character JPHYSICS::WILD_CARD
     * which will be replaced by the PDF types listed in JRegressor<JEnergy, JSimplex>::pdf_t.
     *
     * \param  fileDescriptor     PDF file descriptor
     */
    JRegressor(const std::string& fileDescriptor) :
      estimator(new JMEstimatorNormal())
    {
      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) {

        try {
          
          JPDF_t pdf;

          const JPDFType_t type      = pdf_t[i];
          const string     file_name = getFilename(fileDescriptor, type);
          
          NOTICE("loading PDF from file " << file_name << "... " << flush);
          
          pdf.load(file_name.c_str());
          
          NOTICE("OK" << endl);

          pdf.setExceptionHandler(supervisor);
          
          if      (is_bremsstrahlung(type))
            YB.push_back(JNPE_t(pdf));
          else if (is_deltarays(type))
            YA.push_back(JNPE_t(pdf));
          else
            Y1.push_back(JNPE_t(pdf));
        }
        catch(const JLANG::JException& error) {
          FATAL(error.what() << endl);
        }
      }

      // 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()); }
    }


    /**
     * 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);
    }

    
    /**
     * 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(JK40(T_ns, R_Hz), y1, yA, yB, z);
    }


    /**
     * Create data structure for handling light yields for set of PMTs.
     * Note that the PMT geometry should be relative to the muon trajectory,
     * conform method JGEOMETRY3D::JAxis3D::transform.
     *
     * \param  __begin            begin of PMTs
     * \param  __end              end   of PMTs
     * \param  rates_Hz           K40 singles and multiples rates [Hz]
     * \return                    light yields
     */
    template<class T>
    inline JNPE getNPE(T __begin,
                       T __end,
                       const JK40Rates& rates_Hz) const
    {
      using namespace std;
      using namespace JPP;

      double y1 = 0.0;
      double yA = 0.0;
      double yB = 0.0;
      double z  = 0.0;

      for (T i = __begin; i != __end; ++i) {

        const JNPE   npe = this->getNPE(*i, rates_Hz.getSinglesRate());
        const double x   = i->getX();
        const double y   = i->getY();
        const double R   = sqrt(x*x + y*y);

        y1 += npe.getY1();
        yA += npe.getYA();
        yB += npe.getYB();
        z  += i->getZ()  -  R / getTanThetaC();
      }

      const int N = distance(__begin, __end);

      return JNPE(JK40(T_ns, N, rates_Hz), y1, yA, yB, z / N);
    }


    /**
     * 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)));
    }


    std::vector<JNPE_t> Y1;     //!< light from muon
    std::vector<JNPE_t> YA;     //!< light from delta-rays
    std::vector<JNPE_t> YB;     //!< light from EM showers

    static JTimeRange   T_ns;   //!< Time window with respect to Cherenkov hypothesis [ns]


    static const int NUMBER_OF_PDFS = 6;

    static const JPDFType_t pdf_t[NUMBER_OF_PDFS];

    JLANG::JSinglePointer<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 std::vector<JNPE_t>& NPE)
    {
      double xmax = 0.0;

      for (std::vector<JNPE_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 std::vector<JNPE_t>& NPE,
                                const double R,
                                const double theta,
                                const double phi)
    {
      using JTOOLS::get_value;

      double npe = 0.0;

      for (std::vector<JNPE_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;
    }    
  };


  /**
   * PDF types.
   */
  const JPDFType_t JRegressor<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 };


  /**
   * Time range.
   */
  JTimeRange JRegressor<JEnergy>::T_ns;
}

#endif