#ifndef __JFIT__JSHOWERBJORKENYREGRESSOR__ #define __JFIT__JSHOWERBJORKENYREGRESSOR__ #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/JMEstimator.hh" #include "JFit/JRegressor.hh" #include "JFit/JShowerEH.hh" #include "JFit/JFitToolkit.hh" #include "Jeep/JMessage.hh" /** * \file * Data regression method for JFIT::JShowerEH. * \author adomi */ namespace JFIT { using JPHYSICS::JPDFType_t; using JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER; using JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWER; using JPHYSICS::DIRECT_LIGHT_FROM_BRIGHT_POINT; using JPHYSICS::SCATTERED_LIGHT_FROM_BRIGHT_POINT; using JTOOLS::get_value; /** * Regressor function object for JShowerEH 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; // Part for the Isotropic PDF typedef JTOOLS::JMAPLIST::maplist JMapList_t2; typedef JPHYSICS::JPDFTable JPDF_t2; typedef JTOOLS::JMAPLIST::maplist JNPEMapList_t2; typedef JPHYSICS::JNPETable JNPE_t2; /** * 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::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)); const JPDF_t2::JSupervisor supervisor2(new JPDF_t2::JDefaultResult(JMATH::zero)); for (int i = 0; i != NUMBER_OF_PDFS; ++i) { try { JPDF_t pdf; JPDF_t2 pdf2; const string file_name = getFilename(fileDescriptor, pdf_t[i]); NOTICE("loading PDF from file " << file_name << "... " << flush); if(i < 2){ pdf.load(file_name.c_str()); pdf.setExceptionHandler(supervisor); npe[ i ] = JNPE_t(pdf); } else { pdf2.load(file_name.c_str()); NOTICE("OK" << endl); pdf2.setExceptionHandler(supervisor2); npe2[ i-2 ] = JNPE_t2(pdf2); } } catch(const JException& error) { FATAL(error.what() << endl); } } // Add PDFs for (int i = 1; i < (NUMBER_OF_PDFS-2); i += 2) { npe[ i ].add(npe[i-1]); JNPE_t buffer; npe[i-1].swap(buffer); npe2[ i ].add(npe2[i-1]); JNPE_t2 buffer2; npe2[i-1].swap(buffer2); } } /** * 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 JShowerEH& shower, const JPMTW0& pmt) const { using namespace JPP; using namespace std; JPosition3D D(pmt.getPosition()); JDirection3D U(pmt.getDirection()); D.sub(shower.getPosition()); double ct = U.getDot(D) / D.getLength(); JVersor3D shower_dir(shower.getDirection()); const double z = D.getDot(shower_dir); const double x = D.getX() - z * shower.getDX(); const double y = D.getY() - z * shower.getDY(); const double cosDelta = z/D.getLength(); // Delta = angle between shower direction and PMT position U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane const double theta = U.getTheta(); const double phi = fabs(U.getPhi()); double H0 = getH0(pmt.getR()); // background double H1 = getH1(D.getLength(), ct, cosDelta, theta, phi, shower.getEem(), shower.getEh(), shower.getBy()); // signal if (H1 >= Vmax_npe) { H1 *= Vmax_npe / H1; } H1 += H0; // now H1 is signal + background const bool hit = pmt.getN() != 0; const double u = getChi2(H1, hit); // -log(lik) return estimator->getRho(u); } /** * Get background hypothesis value for time integrated PDF. * * \param R_Hz rate [Hz] * \return hypothesis value */ double getH0(const double R_Hz) const { return get_value(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 ct angle between shower direction and PMT position * \param cosDelta angle between shower direction and PMT position * \param theta PMT zenith angle [deg] * \param phi PMT azimuth angle [deg] * \param Eem EM shower energy [GeV] * \param Eh H shower energy [GeV] * \param Y Bjorken Y * \return hypothesis value */ double getH1(const double D, const double ct, const double cosDelta, const double theta, const double phi, const double Eem, const double Eh, const double Y) const { double h1 = 0; for (int i = 0; i != (NUMBER_OF_PDFS-1); ++i) { if (!npe[i].empty() && D <= npe[i].getXmax() && !npe2[i].empty() && D <= npe2[i].getXmax()) { try { JNPE_t::result_type P_em; JNPE_t2::result_type P_h; P_em = fabs(Eem) * npe[i](std::max(D, npe[i].getXmin()), cosDelta, theta, phi); P_h = fabs(Eh) * npe2[i](std::max(D, npe2[i].getXmin()), ct); double y1 = get_value(P_em) + get_value(P_h); if(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 = 4; static const JPDFType_t pdf_t[NUMBER_OF_PDFS]; JNPE_t npe[NUMBER_OF_PDFS-2]; //!< PDF JNPE_t2 npe2[NUMBER_OF_PDFS-2]; //!< PDF JLANG::JSharedPointer estimator; //!< M-Estimator function }; /** * PDF types. */ const JPDFType_t JRegressor::pdf_t[] = { DIRECT_LIGHT_FROM_EMSHOWER, SCATTERED_LIGHT_FROM_EMSHOWER, 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