#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 "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 : public JAbstractRegressor { using JAbstractRegressor::operator(); typedef JTOOLS::JSplineFunction1S_t JFunction1D_t; typedef JTOOLS::JMAPLIST::maplist JMapList_t; typedef JPHYSICS::JPDFTable JPDF_t; typedef JPHYSICS::JNPETable 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::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& NPE, const double D, const double cd) { double npe = 0.0; for (std::vector::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 Y; //!< light from EM showers JLANG::JSharedPointer estimator; //!< M-Estimator function }; /** * PDF types. */ const JPDFType_t JRegressor::pdf_t[] = { DIRECT_LIGHT_FROM_BRIGHT_POINT, SCATTERED_LIGHT_FROM_BRIGHT_POINT }; /** * Default values. */ JTimeRange JRegressor::T_ns; double JRegressor::Vmax_npe = std::numeric_limits::max(); } #endif