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 "JTools/JConstants.hh"
#include "JTools/JResult.hh"

#include "JMath/JZero.hh"

#include "JFit/JSimplex.hh"
#include "JFit/JShower3EZ.hh"
#include "JFit/JMEstimator.hh"
#include "JFit/JRegressor.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JShowerNPE.hh"
#include "JFit/JShowerNPEHit.hh"

#include "JDetector/JPMT.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 JDETECTOR::JPMT;

  /**
   * 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::WILD_CARD 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  pmt                PMT
     * \param  R_Hz               singles rate [Hz]
     * \return                    light yields
     */
    inline JShowerNPE getNPE(const JPMT&  pmt,
                             const double R_Hz) const
    {
      using namespace JPP;

      const JPosition3D D(pmt.getPosition());
      const JDirection3D U(pmt.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(JK40(T_ns, R_Hz), 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