#ifndef __JFIT__JLINE3ZREGRESSOR__ #define __JFIT__JLINE3ZREGRESSOR__ #include #include "JPhysics/JPDFTypes.hh" #include "JPhysics/JPDFTable.hh" #include "JPhysics/JNPETable.hh" #include "JPhysics/JPDFToolkit.hh" #include "JPhysics/JGeane.hh" #include "JTools/JFunction1D_t.hh" #include "JTools/JFunctionalMap_t.hh" #include "JPhysics/JConstants.hh" #include "JTools/JRange.hh" #include "JGeometry3D/JPosition3D.hh" #include "JGeometry3D/JDirection3D.hh" #include "JMath/JZero.hh" #include "JLang/JSharedPointer.hh" #include "JFit/JTimeRange.hh" #include "JFit/JLine3Z.hh" #include "JFit/JPMTW0.hh" #include "JFit/JSimplex.hh" #include "JFit/JGandalf.hh" #include "JFit/JMEstimator.hh" #include "JFit/JRegressor.hh" #include "JFit/JRegressorHelper.hh" #include "Jeep/JMessage.hh" /** * \file * Data regression method for JFIT::JLine3Z. * \author mdejong */ namespace JFIT {} namespace JPP { using namespace JFIT; } namespace JFIT { /** * Regressor function object for JLine3Z fit using JSimplex minimiser. */ template<> struct JRegressor : public JAbstractRegressor { using JAbstractRegressor::operator(); /** * Constructor. * * \param sigma time resolution of hit [ns] */ JRegressor(double sigma) : estimator(new JMEstimatorNormal()) { this->sigma = sigma; } /** * Fit function. * This method is used to determine the chi2 of given hit with respect to trajectory of muon. * * The template argument JHit_t refers to a data structure which should have the following member methods: * - double getX(); // [m] * - double getY(); // [m] * - double getZ(); // [m] * - double getT(); // [ns] * * \param track track * \param hit hit * \return chi2 */ template double operator()(const JLine3Z& track, const JHit_t& hit) const { using namespace JPP; JPosition3D D(hit.getX(), hit.getY(), hit.getZ()); D.sub(track.getPosition()); const double z = D.getDot(track.getDirection()); const double R = sqrt(D.getLengthSquared() - z*z); const double t1 = track.getT() + (z + R * getKappaC()) * getInverseSpeedOfLight(); const double u = (t1 - hit.getT()) / sigma; return estimator->getRho(u) * hit.getW(); } JLANG::JSharedPointer estimator; //!< M-Estimator function double sigma; //!< Time resolution [ns] }; /** * Regressor function object for JLine3Z fit using JGandalf minimiser. */ template<> struct JRegressor : public JAbstractRegressor, public JRegressorStorage { using JAbstractRegressor::operator(); typedef JRegressorStorage JRegressorStorage_t; /** * Default constructor */ JRegressor() : JRegressorStorage_t(), pdf(getPDF()), npe(getNPE()) {} /** * 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. * * The TTS corresponds to the additional time smearing applied to the PDFs. * * \param fileDescriptor PDF file descriptor * \param TTS TTS [ns] * \param numberOfPoints number of points for Gauss-Hermite integration of TTS * \param epsilon precision for Gauss-Hermite integration of TTS */ JRegressor(const std::string& fileDescriptor, const double TTS, const int numberOfPoints = 25, const double epsilon = 1.0e-10) : JRegressorStorage_t(fileDescriptor, TTS, numberOfPoints, epsilon), pdf(getPDF()), npe(getNPE()) {} /** * Constructor. * * \param storage PDF storage */ JRegressor(const JRegressorStorage& storage) : pdf(storage.getPDF()), npe(storage.getNPE()) {} /** * Fit function. * This method is used to determine the chi2 and gradient of given hit with respect to trajectory of muon. * * The template argument JHit_t refers to a data structure which should have the following member methods: * - double getX(); // [m] * - double getY(); // [m] * - double getZ(); // [m] * - double getDX(); // [u] * - double getDY(); // [u] * - double getDZ(); // [u] * - double getT(); // [ns] * * \param track track * \param hit hit * \return chi2 and gradient */ template result_type operator()(const JLine3Z& track, const JHit_t& hit) const { using namespace JPP; JPosition3D D(hit.getX(), hit.getY(), hit.getZ()); JDirection3D U(hit.getDX(), hit.getDY(), hit.getDZ()); D.sub(track.getPosition()); const double z = D.getDot(track.getDirection()); const double x = D.getX() - z * track.getDX(); const double y = D.getY() - z * track.getDY(); const double R2 = D.getLengthSquared() - z*z; const double R = (R2 > Rmin_m*Rmin_m ? sqrt(R2) : Rmin_m); const double t1 = track.getT() + (z + R * getTanThetaC()) * getInverseSpeedOfLight(); U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane const double theta = U.getTheta(); const double phi = fabs(U.getPhi()); const double E = gWater.getE(E_GeV, z); const double dt = T_ns.constrain(hit.getT() - t1); JPDF_t::result_type H0 = getH0(hit.getR(), dt); JPDF_t::result_type H1 = getH1(E, R, theta, phi, dt); if (H1.V >= Vmax_npe) { H1 *= Vmax_npe / H1.V; } H1 += H0; // signal + background result_type result; result.chi2 = H1.getChi2() - H0.getChi2(); // Likelihood ratio const double wc = 1.0 - getTanThetaC() * z / R; // d(ct1)/d(z) result.gradient = JLine3Z(JLine1Z(JPosition3D(-getTanThetaC() * D.getX() / R, // d(ct1)/d(x) -getTanThetaC() * D.getY() / R, // d(ct1)/d(y) 0.0), getSpeedOfLight()), // d(ct1)/d(t) JVersor3Z(wc * (D.getX() - D.getZ()*track.getDX()/track.getDZ()), // d(ct1)/d(dx) wc * (D.getY() - D.getZ()*track.getDY()/track.getDZ()))); // d(ct1)/d(dy) result.gradient.mul(getInverseSpeedOfLight() * (H1.getDerivativeOfChi2() - H0.getDerivativeOfChi2())); // x d(chi2)/d(ct1) return result; } /** * Fit function. * This method is used to determine the chi2 and gradient of given PMT with respect to trajectory of muon. * * \param track track * \param pmt pmt * \return chi2 and gradient */ result_type operator()(const JLine3Z& track, const JPMTW0& pmt) const { using namespace JGEOMETRY3D; using namespace JPHYSICS; JPosition3D D(pmt.getPosition()); JDirection3D U(pmt.getDirection()); D.sub(track.getPosition()); const double z = D.getDot(track.getDirection()); const double x = D.getX() - z * track.getDX(); const double y = D.getY() - z * track.getDY(); const double R2 = D.getLengthSquared() - z*z; const double R = (R2 > Rmin_m*Rmin_m ? sqrt(R2) : Rmin_m); U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane const double theta = U.getTheta(); const double phi = fabs(U.getPhi()); const double E = gWater.getE(E_GeV, z); JNPE_t::result_type H0 = getH0(pmt.getR()); JNPE_t::result_type H1 = getH1(E, R, theta, phi); if (H1.f >= Vmax_npe) { H1 *= Vmax_npe / H1.f; } H1 += H0; const bool hit = pmt.getN() != 0; const double u = H1.getChi2(hit); result_type result; result.chi2 = estimator->getRho(u); result.gradient = JLine3Z(JLine1Z(JPosition3D(-D.getX() / R, // d(R)/d(x) -D.getY() / R, // d(R)/d(y) 0.0), 0.0), JVersor3Z(-z * D.getX() / R, // d(R)/d(dx) -z * D.getY() / R)); // d(R)/d(dy) result.gradient.mul(estimator->getPsi(u)); result.gradient.mul(H1.getDerivativeOfChi2(hit)); // x d(chi2)/d(R) return result; } /** * Get background hypothesis value for time differentiated PDF. * * \param R_Hz rate [Hz] * \param t1 time [ns] * \return hypothesis value */ JPDF_t::result_type getH0(const double R_Hz, const double t1) const { return JPDF_t::result_type(R_Hz * 1e-9, t1, T_ns); } /** * Get signal hypothesis value for time differentiated PDF. * * \param E muon energy at minimum distance of approach [GeV] * \param R minimum distance of approach [m] * \param theta PMT zenith angle [rad] * \param phi PMT azimuth angle [rad] * \param t1 arrival time relative to Cherenkov hypothesis [ns] * \return hypothesis value */ JPDF_t::result_type getH1(const double E, const double R, const double theta, const double phi, const double t1) const { using namespace std; using namespace JPP; JPDF_t::result_type h1 = zero; for (int i = 0; i != NUMBER_OF_PDFS; ++i) { if (!pdf[i].empty() && R <= pdf[i].getXmax()) { JPDF_t::result_type y1 = pdf[i](max(R, pdf[i].getXmin()), theta, phi, t1); // safety measures if (y1.f <= 0.0) { y1.f = 0.0; y1.fp = 0.0; } if (y1.v <= 0.0) { y1.v = 0.0; } // energy dependence if (is_deltarays(pdf_t[i])) { y1 *= getDeltaRaysFromMuon(E); } else if (is_bremsstrahlung(pdf_t[i])) { y1 *= E; } h1 += y1; } } return h1; } /** * 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 E muon energy at minimum distance of approach [GeV] * \param R minimum distance of approach [m] * \param theta PMT zenith angle [rad] * \param phi PMT azimuth angle [rad] * \return hypothesis value */ JNPE_t::result_type getH1(const double E, const double R, const double theta, const double phi) const { using namespace std; using namespace JPP; JNPE_t::result_type h1 = JMATH::zero; for (int i = 0; i != NUMBER_OF_PDFS; ++i) { if (!npe[i].empty() && R <= npe[i].getXmax()) { try { JNPE_t::result_type y1 = npe[i](max(R, npe[i].getXmin()), theta, phi); // safety measures if (y1.f > 0.0) { // energy dependence if (is_deltarays(pdf_t[i])) { y1 *= getDeltaRaysFromMuon(E); } else if (is_bremsstrahlung(pdf_t[i])) { y1 *= E; } h1 += y1; } } catch(JException& error) { ERROR(error << endl); } } } return h1; } /** * Get maximal road width of PDF. * * \return road width [m] */ inline double getRmax() const { double xmax = 0.0; for (int i = 0; i != NUMBER_OF_PDFS; ++i) { if (!pdf[i].empty() && pdf[i].getXmax() > xmax) { xmax = pdf[i].getXmax(); } } return xmax; } static JTimeRange T_ns; //!< Time window with respect to Cherenkov hypothesis [ns] static double Vmax_npe; //!< Maximal integral of PDF [npe] static double Rmin_m; //!< Minimal distance of [m] const JPDFs_t& pdf; //!< PDF const JNPEs_t& npe; //!< PDF double E_GeV = 0.0; //!< Energy of muon at vertex [GeV] JLANG::JSharedPointer estimator = new JMEstimatorNormal(); //!< M-Estimator function }; /** * Default values. */ JTimeRange JRegressor::T_ns; double JRegressor::Vmax_npe = std::numeric_limits::max(); double JRegressor::Rmin_m = 0.1; } #endif