JShowerEnergyRegressor.hh

Go to the documentation of this file.

#ifndef __JFIT__JSHOWERENERGYREGRESSOR__
#define __JFIT__JSHOWERENERGYREGRESSOR__
 
#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 "JLang/JSharedPointer.hh"
#include "JMath/JZero.hh"
 
#include "JGeometry3D/JAxis3D.hh"
 
#include "JFit/JEnergy.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JMEstimator.hh"
#include "JFit/JRegressor.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JShowerNPE.hh"
#include "JFit/JShowerNPEHit.hh"
#include "JFit/JTimeRange.hh"
 
#include "Jeep/JMessage.hh"
 
/**
 * \file
 * Data regression method for JFIT::JShower3EZ only focused on the energy estimation from 
 * a bright point emission PDF and considering the hit/non hit information for each PMT.
 *
 * \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;
  using JGEOMETRY3D::JAxis3D;
 
  /**
   * Regressor function object for JShower3EZ fit using JSimplex minimiser.
   */
  template<>
  struct JRegressor<JEnergy, JSimplex> :
    public JAbstractRegressor<JEnergy, JSimplex>
  {
    using JAbstractRegressor<JEnergy, JSimplex>::operator();
 
    typedef JTOOLS::JSplineFunction1S_t                                      JFunction1D_t;
    typedef JTOOLS::JMAPLIST<JTOOLS::JPolint1FunctionalMap,
                             JTOOLS::JPolint1FunctionalGridMap>::maplist     JMapList_t;
    typedef JPHYSICS::JPDFTable<JFunction1D_t, JMapList_t>                   JPDF_t;
 
    typedef JPHYSICS::JNPETable<double, double, JMapList_t>                  JNPE_t;
    
    /**
     * 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 JRegressor<JEnergy, JSimplex>::pdf_t.
     *
     * \param  fileDescriptor     PDF file descriptor
     */
 
    JRegressor(const std::string& fileDescriptor):
      estimator(new JMEstimatorNull())
    {
      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 {
 
          JPDF_t pdf;
 
          const string file_name = getFilename(fileDescriptor, pdf_t[i]);
          
          NOTICE("loading PDF from file " << file_name << "... " << flush);
 
          pdf.load(file_name.c_str());
 
          NOTICE("OK" << endl);
 
          pdf.setExceptionHandler(supervisor);
                  
          Y.push_back(JNPE_t(pdf));
        }
        catch(const JException& error) {
          FATAL(error.what() << endl);
        }
      }
 
      // Add PDFs
 
      for (int i = 1; i < NUMBER_OF_PDFS; i += 2) {
        Y[i].add(Y[i-1]);
      }    
 
      Y.erase(Y.begin());
    }
 
    /**
     * Fit function.
     * This method is used to determine the chi2 of given PMT with respect to shower hypothesis.
     *
     * \param  x                  energy
     * \param  npe                number of photoelectrons
     * \return                    chi2 
     */
    double operator()(const JEnergy& x, const JShowerNPEHit& npe) const
    {      
      using namespace JPP;
 
      const double E = x.getE();
      const double u = npe.getChi2(E);
 
      return estimator->getRho(u);
    }
 
    /**
     * Create data structure for handling light yields for PMT.
     *
     * \param  axis               PMT axis
     * \param  R_Hz               singles rate [Hz]
     * \return                    light yields
     */
    inline JShowerNPE getNPE(const JAxis3D& axis,
                             const double   R_Hz) const
    {
      using namespace JPP;
 
      const JPosition3D  D(axis.getPosition());
      const JDirection3D U(axis.getDirection());      
 
      const double U_length = std::sqrt(U.getDX()*U.getDX() + 
                                        U.getDY()*U.getDY() + 
                                        U.getDZ()*U.getDZ());
 
      const double ct = U.getDot(D) / (D.getLength()*U_length);
      
      const double y = getNPE(Y, D.getLength(), ct);
 
      return JShowerNPE(getN(T_ns, R_Hz * 1.0e-9), y);      
    }
 
 
    /**
     * Get number of photo-electrons.
     *
     * \param  NPE              NPE tables
     * \param  D                PMT distance from shower [m]
     * \param  cd               cosine of the PMT angle wrt the photon direction
     * \return                  number of photo-electrons
     */
    static inline double getNPE(const std::vector<JNPE_t>& NPE,
                                const double D,
                                const double cd)
    {
      double npe = 0.0;
       
      for (std::vector<JNPE_t>::const_iterator i = NPE.begin(); i != NPE.end(); ++i) {
         
        if (D <= i->getXmax()) {
           
          try {
             
            const double y = get_value((*i)(std::max(D, i->getXmin()), cd));
             
            if (y > 0.0) {
              npe += y;
            }
          }
          catch(const JLANG::JException& error) {
            ERROR(error << std::endl);
          }
        }
      }
       
      return npe;
    }    
      
 
    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];
 
    std::vector<JNPE_t> Y;     //!< light from EM showers
 
    JLANG::JSharedPointer<JMEstimator> estimator;  //!< M-Estimator function
  };  
 
  /**
   * PDF types.
   */
  const JPDFType_t JRegressor<JEnergy, JSimplex>::pdf_t[] = { DIRECT_LIGHT_FROM_BRIGHT_POINT, 
                                                              SCATTERED_LIGHT_FROM_BRIGHT_POINT };
 
  /**
   * Default values.
   */
  JTimeRange JRegressor<JEnergy, JSimplex>::T_ns;
  double     JRegressor<JEnergy, JSimplex>::Vmax_npe = std::numeric_limits<double>::max();
}
 
#endif