#ifndef __JFIT__JSHOWER3EZREGRESSOR__ #define __JFIT__JSHOWER3EZREGRESSOR__ #include "JPhysics/JPDFTypes.hh" #include "JPhysics/JPDFTable.hh" #include "JPhysics/JPDFToolkit.hh" #include "JPhysics/JNPETable.hh" #include "JTools/JFunction1D_t.hh" #include "JTools/JFunctionalMap_t.hh" #include "JPhysics/JConstants.hh" #include "JTools/JRange.hh" #include "JTools/JResult.hh" #include "JTools/JFunction1D_t.hh" #include "JTools/JFunctionalMap_t.hh" #include "JTools/JAbstractHistogram.hh" #include "JGeometry3D/JVector3D.hh" #include "JGeometry3D/JVersor3D.hh" #include "JGeometry3D/JDirection3D.hh" #include "JMath/JZero.hh" #include "JFit/JTimeRange.hh" #include "JFit/JPMTW0.hh" #include "JFit/JSimplex.hh" #include "JFit/JGandalf.hh" #include "JFit/JMEstimator.hh" #include "JFit/JRegressor.hh" #include "JFit/JShower3EZ.hh" #include "JFit/JFitToolkit.hh" #include "JFit/JLine3Z.hh" #include "Jeep/JMessage.hh" /** * \file * Data regression method for JFIT::JShower3EZ. * \author mdejong, vcarretero */ namespace JFIT { using JPHYSICS::JPDFType_t; using JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER; using JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWER; using JTOOLS::get_value; /** * Constrain PMT angle to [0,pi]. * * \param angle angle [rad] * \return angle [rad] */ inline double getPMTAngle(const double angle) { const double epsilon = 1.0e-6; const JTOOLS::JRange range(epsilon, JMATH::PI - epsilon); return range.constrain(fabs(angle)); } /** * Function to constrain the versor and energy during the fit, to prevent unphysical values. * * \param value model (I/O) */ void model(JShower3EZ& value) { using namespace std; double Tx = value.getDX(); double Ty = value.getDY(); double E = max(0.0,value.getE()); const double u = hypot(Tx, Ty); if (u > 1.0) { Tx /= u; Ty /= u; } value = JShower3EZ(static_cast(value), JVersor3Z(Tx,Ty), E, value.getBy()); } /** * 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 > > > JPDFMaplist_t; typedef JPHYSICS::JPDFTable JPDF_t; typedef JTOOLS::JMapList > > > JNPEMaplist_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); npe[ i ] = JNPE_t(pdf); } catch(const JException& error) { FATAL(error.what() << endl); } } // Add PDFs for (int i = 1; i < NUMBER_OF_PDFS; i += 2) { npe[ i ].add(npe[i-1]); JNPE_t buffer; npe[i-1].swap(buffer); } } /** * Fit function. * This method is used to determine the chi2 of given PMT with respect to shower hypothesis. * * \param shower shower * \param pmt pmt * \return chi2 */ double operator()(const JShower3EZ& shower, const JPMTW0& pmt) const { using namespace JPP; JPosition3D D(pmt.getPosition()); JDirection3D U(pmt.getDirection()); D.sub(shower.getPosition()); const double z = D.getDot(shower.getDirection()); const double x = D.getX(); const double y = D.getY(); const double cd = z/D.getLength(); // cosine angle between shower direction and PMT position U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane const double theta = getPMTAngle(U.getTheta()); const double phi = getPMTAngle(U.getPhi()); JNPE_t::result_type H0 = getH0(pmt.getR()); // background hypothesis value for time integrated PDF. JNPE_t::result_type H1 = getH1(D.getLength(), cd, theta, phi, shower.getE()); // signal hypothesis value for time integrated PDF. if (get_value(H1) >= Vmax_npe) { H1 *= Vmax_npe / get_value(H1); } H1 += H0; // now H1 is signal + background const bool hit = pmt.getN() != 0; const double u = getChi2(get_value(H1), hit); // -log(lik) return estimator->getRho(u); } /** * Get background hypothesis value for time integrated PDF. * * \param R_Hz rate [Hz] * \return hypothesis value */ JNPE_t::result_type getH0(const double R_Hz) const { return JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength()); } /** * Get signal hypothesis value for time integrated PDF. * * \param D PMT distance from shower [m] * \param cd cosine angle between shower direction and PMT position * \param theta PMT zenith angle [deg] * \param phi PMT azimuth angle [deg] * \param E shower energy [GeV] * \return hypothesis value */ JNPE_t::result_type getH1(const double D, const double cd, const double theta, const double phi, const double E) const { JNPE_t::result_type h1 = JMATH::zero; for (int i = 0; i != NUMBER_OF_PDFS; ++i) { if (!npe[i].empty() && D <= npe[i].getXmax()) { try { JNPE_t::result_type y1 = E * npe[i](std::max(D, npe[i].getXmin()), cd, theta, phi); // safety measure if(y1 < 0){ y1 = 0; } h1 += y1; } catch(JLANG::JException& error) { ERROR(error << std::endl); } } } return h1; } 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]; JNPE_t npe[NUMBER_OF_PDFS]; //!< PDF JLANG::JSharedPointer estimator; //!< M-Estimator function }; /** * Regressor function object for JShower3EZ fit using JGandalf minimiser. */ template<> struct JRegressor : public JAbstractRegressor { using JAbstractRegressor::operator(); typedef JTOOLS::JSplineFunction1S_t JFunction1D_t; typedef JTOOLS::JMapList > > > JPDFMaplist_t; typedef JPHYSICS::JPDFTable JPDF_t; typedef JTOOLS::JMapList > > > JNPEMaplist_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); npe[i] = JNPE_t(pdf); } catch(const JException& error) { FATAL(error.what() << endl); } } // Add PDFs for (int i = 1; i < NUMBER_OF_PDFS; i += 2) { npe[ i ].add(npe[i-1]); JNPE_t buffer; npe[i-1].swap(buffer); } } /** * Fit function. * This method is used to determine the chi2 of given PMT with respect to shower hypothesis. * * \param shower shower * \param pmt pmt * \return chi2 */ result_type operator()(const JShower3EZ& shower, const JPMTW0& pmt) const { using namespace JPP; using namespace std; JPosition3D D(pmt.getPosition()); JDirection3D U(pmt.getDirection()); D.sub(shower.getPosition()); const double x = D.getX(); const double y = D.getY(); const double d = D.getLength(); const double cd = D.getDot(shower.getDirection())/d; // cosine angle between shower direction and PMT position U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane const double theta = getPMTAngle(U.getTheta()); const double phi = getPMTAngle(U.getPhi()); JNPE_t::result_type H0 = getH0(pmt.getR()); // background hypothesis value for time integrated PDF. JNPE_t::result_type H1 = getH1(d, cd, theta, phi, shower.getE()); // signal hypothesis value for time integrated PDF. if (get_value(H1) >= Vmax_npe) { H1 *= Vmax_npe / get_value(H1); } double signal_npe = get_value(H1); H1 += H0; // now H1 is signal + background double expectation = get_value(H1); const bool hit = pmt.getN() != 0; const double u = H1.getChi2(hit); result_type result; result.chi2 = estimator->getRho(u); double energy_gradient = signal_npe/shower.getE(); // dP/dE if(hit) energy_gradient *= -exp(-expectation)/(1-exp(-expectation)); //dchi2/d(H1), if !hit is 1 result.gradient = JShower3EZ(JPoint4D(JVector3D(0, // d(cos_th0)/d(x) 0, // d(cos_th0)/d(y) 0), // d(cos_th0)/d(z) 0.0), // d(cos_th0)/d(t) JVersor3Z(x/d, // d(cos_th0)/d(dx) y/d), // d(cos_th0)/d(dy) energy_gradient); // d(chi2)/d(E) result.gradient.mul(estimator->getPsi(u)); static_cast(result.gradient).mul(H1.getDerivativeOfChi2(hit)); // x d(chi2)/d(cos_th0) return result; } /** * Get background hypothesis value for time integrated PDF. * * \param R_Hz rate [Hz] * \return hypothesis value */ JNPE_t::result_type getH0(const double R_Hz) const { return JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength(), 0.0); } /** * Get signal hypothesis value for time integrated PDF. * * \param D PMT distance from shower [m] * \param cd cosine angle between shower direction and PMT position * \param theta PMT zenith angle [deg] * \param phi PMT azimuth angle [deg] * \param E shower energy [GeV] * \return hypothesis value */ JNPE_t::result_type getH1(const double D, const double cd, const double theta, const double phi, const double E) const { JNPE_t::result_type h1 = JMATH::zero; for (int i = 0; i != NUMBER_OF_PDFS; ++i) { if (!npe[i].empty() && D <= npe[i].getXmax()) { try { JNPE_t::result_type y1 = E * npe[i](std::max(D, npe[i].getXmin()), cd, theta, phi); if (get_value(y1) > 0.0) { h1 += y1; } } catch(JLANG::JException& error) { ERROR(error << std::endl); } } } return h1; } 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 JFIT::JPDFType_t pdf_t[NUMBER_OF_PDFS]; JNPE_t npe[NUMBER_OF_PDFS]; //!< PDF JLANG::JSharedPointer estimator; //!< M-Estimator function }; /** * Regressor function object for JShower3EZ fit using Abstract minimiser, that just computes the chi2 without a fit. */ template<> struct JRegressor : public JAbstractRegressor { using JAbstractRegressor::operator(); typedef JTOOLS::JSplineFunction1S_t JFunction1D_t; typedef JTOOLS::JMapList > > > JPDFMaplist_t; typedef JPHYSICS::JPDFTable JPDF_t; typedef JTOOLS::JMapList > > > JNPEMaplist_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); npe[ i ] = JNPE_t(pdf); } catch(const JException& error) { FATAL(error.what() << endl); } } // Add PDFs for (int i = 1; i < NUMBER_OF_PDFS; i += 2) { npe[ i ].add(npe[i-1]); JNPE_t buffer; npe[i-1].swap(buffer); } } /** * Fit function. * This method is used to determine the chi2 of given PMT with respect to shower hypothesis. * * \param shower shower * \param pmt pmt * \return chi2 */ double operator()(const JShower3EZ& shower, const JPMTW0& pmt) const { using namespace JPP; JPosition3D D(pmt.getPosition()); JDirection3D U(pmt.getDirection()); D.sub(shower.getPosition()); const double z = D.getDot(shower.getDirection()); const double x = D.getX(); const double y = D.getY(); const double cd = z/D.getLength(); // cosine angle between shower direction and PMT position U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane const double theta = getPMTAngle(U.getTheta()); const double phi = getPMTAngle(U.getPhi()); JNPE_t::result_type H0 = getH0(pmt.getR()); // background hypothesis value for time integrated PDF. JNPE_t::result_type H1 = getH1(D.getLength(), cd, theta, phi, shower.getE()); // signal hypothesis value for time integrated PDF. if (get_value(H1) >= Vmax_npe) { H1 *= Vmax_npe / get_value(H1); } H1 += H0; // now H1 is signal + background const bool hit = pmt.getN() != 0; const double u = getChi2(get_value(H1), hit); // -log(lik) return estimator->getRho(u); } /** * Get background hypothesis value for time integrated PDF. * * \param R_Hz rate [Hz] * \return hypothesis value */ JNPE_t::result_type getH0(const double R_Hz) const { return JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength()); } /** * Get signal hypothesis value for time integrated PDF. * * \param D PMT distance from shower [m] * \param cd cosine angle between shower direction and PMT position * \param theta PMT zenith angle [deg] * \param phi PMT azimuth angle [deg] * \param E shower energy [GeV] * \return hypothesis value */ JNPE_t::result_type getH1(const double D, const double cd, const double theta, const double phi, const double E) const { JNPE_t::result_type h1 = JMATH::zero; for (int i = 0; i != NUMBER_OF_PDFS; ++i) { if (!npe[i].empty() && D <= npe[i].getXmax()) { try { JNPE_t::result_type y1 = E * npe[i](std::max(D, npe[i].getXmin()), cd, theta, phi); // safety measure if(y1 < 0){ y1 = 0; } h1 += y1; } catch(JLANG::JException& error) { ERROR(error << std::endl); } } } return h1; } /** * Get signal hypothesis value for time integrated PDF. * * \param shower shower * \param pmt pmt * \return hypothesis value */ JNPE_t::result_type getH1(const JShower3EZ& shower, const JPMTW0& pmt) const { using namespace JPP; JPosition3D D(pmt.getPosition()); JDirection3D U(pmt.getDirection()); const double z = D.getDot(shower.getDirection()); const double x = D.getX(); const double y = D.getY(); const double cd = z/D.getLength(); // cosine angle between shower direction and PMT position U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane const double theta = getPMTAngle(U.getTheta()); const double phi = getPMTAngle(U.getPhi()); JNPE_t::result_type H1 = getH1(D.getLength(), cd, theta, phi, 1.0); // signal hypothesis value for time integrated PDF. E=1 because it is linear with E. if (get_value(H1) >= Vmax_npe) { H1 *= Vmax_npe / get_value(H1); } return H1; } 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]; JNPE_t npe[NUMBER_OF_PDFS]; //!< PDF JLANG::JSharedPointer estimator; //!< M-Estimator function }; /** * PDF types. */ const JPDFType_t JRegressor::pdf_t[] = { DIRECT_LIGHT_FROM_EMSHOWER, SCATTERED_LIGHT_FROM_EMSHOWER }; const JPDFType_t JRegressor::pdf_t[] = { DIRECT_LIGHT_FROM_EMSHOWER, SCATTERED_LIGHT_FROM_EMSHOWER }; const JPDFType_t JRegressor::pdf_t[] = { DIRECT_LIGHT_FROM_EMSHOWER, SCATTERED_LIGHT_FROM_EMSHOWER }; /** * Default values. */ JTimeRange JRegressor::T_ns; double JRegressor::Vmax_npe = std::numeric_limits::max(); JTimeRange JRegressor::T_ns; double JRegressor::Vmax_npe = std::numeric_limits::max(); JTimeRange JRegressor::T_ns; double JRegressor::Vmax_npe = std::numeric_limits::max(); } #endif