#ifndef __JAANET__JHONDAFLUXINTERPOLATOR__ #define __JAANET__JHONDAFLUXINTERPOLATOR__ #include #include #include #include "Jeep/JPrint.hh" #include "Jeep/JMessage.hh" #include "JLang/JClonable.hh" #include "JLang/JException.hh" #include "JTools/JRange.hh" #include "JTools/JPolint.hh" #include "JTools/JMapList.hh" #include "JTools/JCollection.hh" #include "JTools/JGridCollection.hh" #include "JTools/JFunction1D_t.hh" #include "JTools/JMultiFunction.hh" #include "JTools/JFunctionalMap_t.hh" #include "JAAnet/JDiffuseFlux.hh" /** * \author bjung */ namespace JAANET {} namespace JPP { using namespace JAANET; } namespace JAANET { using JEEP::JMessage; using JLANG::JClonable; using JTOOLS::JRange; using JTOOLS::JArray; using JTOOLS::JMAPLIST; using JTOOLS::JMultiFunction; using JTOOLS::JConstantFunction1D; using JTOOLS::JPolint1FunctionalMap; static const char* const HONDA_ENERGY_KEYWORD = "Enu(GeV)"; static const char* const HONDA_COSINE_ZENITH_ANGLE_KEYWORD = "cosZ"; static const char* const HONDA_AZIMUTH_ANGLE_KEYWORD = "phi_Az"; static const char* const HONDA_MUON_NEUTRINO_FLUX_KEYWORD = "NuMu"; static const char* const HONDA_MUON_ANTINEUTRINO_FLUX_KEYWORD = "NuMubar"; static const char* const HONDA_ELECTRON_NEUTRINO_FLUX_KEYWORD = "NuE"; static const char* const HONDA_ELECTRON_ANTINEUTRINO_FLUX_KEYWORD = "NuEbar"; static const char* const HONDA_BINSPECS_SEPARATOR = "--"; static const char HONDA_BINSPECS_BEGIN = '['; static const char HONDA_BINSPECS_END = ']'; /** * Auxiliary data structure for reading Honda bin ranges. */ struct JHondaBinRange : JRange { typedef JRange JRange_t; /** * Default constructor. */ JHondaBinRange() {} /** * Constructor. * * \param minimum minimum value * \param maximum maximum value */ JHondaBinRange(const double minimum, const double maximum) : JRange_t(minimum, maximum) {} /** * Get central value of this range. * * \return central value */ double getCentralValue() const { return (getLowerLimit() + getUpperLimit()) / 2.0; } /** * Read bin range from input. * * \param in input stream * \param object bin range * \return input stream */ friend inline std::istream& operator>>(std::istream& in, JHondaBinRange& object) { using namespace std; double min = JRange_t::getMinimum(); double max = JRange_t::getMaximum(); string buffer; if (in >> min >> buffer >> max && min <= max && buffer == HONDA_BINSPECS_SEPARATOR) { static_cast(object) = JRange_t(min, max); } return in; } }; /** * Auxiliary data structure for reading angular binning specifications of Honda flux table. */ struct JHondaAngularBinSpecs { /** * Default constructor. */ JHondaAngularBinSpecs() : coszRange(-1.0, 1.0), phiRange ( 0.0, 360.0) {} /** * Constructor. * * \param minCosz minimum cosine zenith angle * \param maxCosz maximum cosine zenith angle * \param minPhi minimum azimuth angle * \param maxPhi maximum azimuth angle */ JHondaAngularBinSpecs(const double minCosz, const double maxCosz, const double minPhi, const double maxPhi) : coszRange(minCosz, maxCosz), phiRange (minPhi, maxPhi) {} /** * Read bin specifications from input. * * \param in input stream * \param object bin specifications * \return input stream */ friend inline std::istream& operator>>(std::istream& in, JHondaAngularBinSpecs& binspecs) { static const JEquationParameters eqpars("=", ",", "./", "#"); JProperties properties(eqpars, 1); properties[HONDA_COSINE_ZENITH_ANGLE_KEYWORD] = binspecs.coszRange; properties[HONDA_AZIMUTH_ANGLE_KEYWORD] = binspecs.phiRange; in >> properties; return in >> properties; } JHondaBinRange coszRange; //!< Cosine zenith-angle range JHondaBinRange phiRange; //!< Azimuthal angle range [deg] }; /** * Template definition for Honda flux table interpolator. */ template >, template class JCoszFunctionalMap_t = JPolint1FunctionalMap, template class JEnergyFunctionalMap_t = JPolint1FunctionalMap> class JHondaFluxInterpolator final : public JMultiFunction ::maplist>, public JClonable >, public JMessage > { public: typedef typename JMAPLIST::maplist JFunctionalMaplist_t; typedef JHondaFluxInterpolator interpolator_type; typedef JMultiFunction multifunction_type; enum { NUMBER_OF_DIMENSIONS = multifunction_type::NUMBER_OF_DIMENSIONS }; typedef typename multifunction_type::abscissa_type abscissa_type; typedef typename multifunction_type::argument_type argument_type; typedef typename multifunction_type::result_type result_type; typedef typename multifunction_type::value_type value_type; typedef typename multifunction_type::multimap_type multimap_type; typedef typename multifunction_type::super_const_iterator super_const_iterator; typedef typename multifunction_type::super_iterator super_iterator; typedef typename multifunction_type::function_type function_type; typedef JRange JRange_t; /** * Default constructor. */ JHondaFluxInterpolator() : table(), coszRange(0.0, -1.0), phiRange (0.0, -1.0) {} /** * Constructor. * * \param fileName oscillation probability table filename */ JHondaFluxInterpolator(const char* fileName) : JHondaFluxInterpolator() { this->load(fileName); this->check_validity(); } /** * Load oscillation probability table from file. * * \param fileName oscillation probability table filename */ void load(const char* const fileName) { using namespace std; using namespace JPP; try { NOTICE("Loading Honda flux table from file " << fileName << "... " << flush); multimap_type& zmap = static_cast(*this); ifstream ifs(fileName); JHondaAngularBinSpecs binspecs; for (string line; getline(ifs, line); ) { double energy = 0.0; result_type values; istringstream iss1(line); const double cosz = binspecs.coszRange.getCentralValue(); const double phi = binspecs.phiRange .getCentralValue(); coszRange.include(cosz); phiRange .include(phi); if (iss1 >> energy >> values) { const double log10E = log10(energy); for (typename result_type::iterator i = values.begin(); i != values.end(); ++i) { *i = log(*i); } zmap[log10E][cosz][phi] = values; } else { const size_t p1 = line.find (HONDA_BINSPECS_BEGIN); const size_t p2 = line.rfind(HONDA_BINSPECS_END); if (p1 != string::npos && p2 != string::npos) { const string binspecstr = line.substr(p1 + 1, p2 - p1 - 1); istringstream iss2(binspecstr); iss2 >> binspecs; } } } } catch(const exception& error) { THROW(JFileReadException, "JHondaFluxInterpolator::load(): Error reading file " << fileName); } } /** * Check whether this interpolator is valid. * * \return true if valid; else false */ bool is_valid() const override final { return this->size() > 0 && coszRange.getLength() > 0.0 && phiRange.getLength() > 0.0; } /** * Get diffuse flux for given particle PDG-identifier, energy, cosine zenith angle and azimuth angle. * * \param type PDG particle type * \param log10E logarithmic neutrino energy [GeV] * \param costh cosine zenith angle * \param phi azimuth angle (defined clockwise w.r.t. the North) [deg] * \return diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1] */ double getFactor(const int type, const double log10E, const double costh, const double phi) const { using namespace std; using namespace JPP; static const double epsilon_costh = 1.0e-3; //!< [--] static const double epsilon_phi = 1.0; //!< [deg] static const JRange_t mod(0.0, 360.0); const multifunction_type& function = static_cast(*this); const double _costh = min(max(costh, -1.0), 1.0); const double _phi = mod.mod(180.0 - phi); //!< Honda azimuth angle is defined counterclockwise w.r.t. the South\n //!< (c.f. arXiv:1102.2688), so a conversion is necessary. result_type lnFluxes; if (coszRange(_costh) && phiRange(_phi)) { lnFluxes = function(log10E, _costh, _phi); } else { // Perform linear extrapolation for edges of angular ranges if (_costh > coszRange.getUpperLimit() && phiRange(_phi)) { const abscissa_type x1 = coszRange.getUpperLimit(); const abscissa_type x2 = coszRange.getUpperLimit() - epsilon_costh; const result_type& y1 = function(log10E, x1, _phi); const result_type& y2 = function(log10E, x2, _phi); lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_costh - x1); } else if (_costh < coszRange.getLowerLimit() && phiRange(_phi)) { const abscissa_type x1 = coszRange.getLowerLimit(); const abscissa_type x2 = coszRange.getLowerLimit() + epsilon_costh; const result_type& y1 = function(log10E, x1, _phi); const result_type& y2 = function(log10E, x2, _phi); lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_costh - x1); } else if (coszRange(_costh) && _phi > phiRange.getUpperLimit()) { const abscissa_type x1 = phiRange.getUpperLimit(); const abscissa_type x2 = phiRange.getUpperLimit() - epsilon_costh; const result_type& y1 = function(log10E, _costh, x1); const result_type& y2 = function(log10E, _costh, x2); lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_phi - x1); } else if (coszRange(_costh) && _phi < phiRange.getLowerLimit()) { const abscissa_type x1 = phiRange.getLowerLimit(); const abscissa_type x2 = phiRange.getLowerLimit() + epsilon_costh; const result_type& y1 = function(log10E, _costh, x1); const result_type& y2 = function(log10E, _costh, x2); lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_phi - x1); } else if (_costh > coszRange.getUpperLimit() && _phi > phiRange.getUpperLimit()) { const abscissa_type c1 = coszRange.getUpperLimit(); const abscissa_type c2 = coszRange.getUpperLimit() - epsilon_costh; const abscissa_type p1 = phiRange .getUpperLimit(); const abscissa_type p2 = phiRange .getUpperLimit() - epsilon_phi; const result_type& y1 = function(log10E, c1, p1); const result_type& y2 = function(log10E, c2, p2); lnFluxes = y1 + ((y2 - y1) / (c2 - c1) / (p2 - p1)) * (_costh - p1) * (_phi - p1); } else if (_costh > coszRange.getUpperLimit() && _phi < phiRange.getUpperLimit()) { const abscissa_type c1 = coszRange.getUpperLimit(); const abscissa_type c2 = coszRange.getUpperLimit() - epsilon_costh; const abscissa_type p1 = phiRange .getLowerLimit(); const abscissa_type p2 = phiRange .getLowerLimit() + epsilon_phi; const result_type& y1 = function(log10E, c1, p1); const result_type& y2 = function(log10E, c2, p2); lnFluxes = y1 + ((y2 - y1) / (c2 - c1) / (p2 - p1)) * (_costh - p1) * (_phi - p1); } else if (_costh < coszRange.getLowerLimit() && phi > phiRange.getUpperLimit()) { const abscissa_type c1 = coszRange.getLowerLimit(); const abscissa_type c2 = coszRange.getLowerLimit() + epsilon_costh; const abscissa_type p1 = phiRange .getUpperLimit(); const abscissa_type p2 = phiRange .getUpperLimit() - epsilon_phi; const result_type& y1 = function(log10E, c1, p1); const result_type& y2 = function(log10E, c2, p2); lnFluxes = y1 + ((y2 - y1) / (c2 - c1) / (p2 - p1)) * (_costh - p1) * (_phi - p1); } else { const abscissa_type c1 = coszRange.getLowerLimit(); const abscissa_type c2 = coszRange.getLowerLimit() + epsilon_costh; const abscissa_type p1 = phiRange .getLowerLimit(); const abscissa_type p2 = phiRange .getLowerLimit() + epsilon_phi; const result_type& y1 = function(log10E, c1, p1); const result_type& y2 = function(log10E, c2, p2); lnFluxes = y1 + ((y2 - y1) / (c2 - c1) / (p2 - p1)) * (_costh - p1) * (_phi - p1); } } switch(type) { case TRACK_TYPE_NUMU: return exp(lnFluxes[0]); case TRACK_TYPE_ANTINUMU: return exp(lnFluxes[1]); case TRACK_TYPE_NUE: return exp(lnFluxes[2]); case TRACK_TYPE_ANTINUE: return exp(lnFluxes[3]); default: return 0.0; } } /** * Get flux of given event. * * \param type PDG particle type * \param log10E logarithmic neutrino energy [GeV] * \param costh cosine zenith angle * \param phi azimuth angle (defined clockwise w.r.t. the North) [deg] * \return diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1] */ double getFlux(const int type, const double log10E, const double costh, const double phi) const { return getFactor(type, log10E, costh, phi); } /** * Get diffuse flux for given particle PDG-identifier, energy, cosine zenith angle and azimuth angle. * * \param type PDG particle type * \param log10E logarithmic neutrino energy [GeV] * \param costh cosine zenith angle * \param phi azimuth angle (defined clockwise w.r.t. the North) [deg] * \return diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1] */ double operator()(const int type, const double log10E, const double costh, const double phi) const { return getFactor(type, log10E, costh, phi); } /** * Get flux of given event. * * \param evt event * \return flux [GeV^-1 * m^-2 * sr^-1 * s^-1] */ double getFactor(const Evt& evt) const override final { const Trk& neutrino = get_neutrino(evt); const double log10E = log10(neutrino.E); const double costh = neutrino.dir.z; //!< The gSeaGen azimuth angle is defined clockwise w.r.t. the North (c.f. arXiv:2003.14040) const double phi = atan(neutrino.dir.y / neutrino.dir.x) * 180.0 / M_PI; return getFactor(neutrino.type, log10E, costh, phi); } /** * Get flux of given event. * * \param event event * \return flux [GeV^-1 * m^-2 * sr^-1 * s^-1] */ double getFlux(const Evt& event) const { return getFactor(event); } /** * Get properties of this class. * * \param eqpars equation parameters */ JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) override final { return JHondaFluxInterpolatorHelper(*this, eqpars); } /** * Get properties of this class. * * \param eqpars equation parameters */ JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) const override final { return JHondaFluxInterpolatorHelper(*this, eqpars); } /** * Stream input. * * \param in input stream * \return input stream */ std::istream& read(std::istream& in) override final { using namespace std; streampos pos = in.tellg(); if (!(in >> table)) { in.clear(); in.seekg(pos); JProperties properties = getProperties(); in >> properties; } load(table.c_str()); this->check_validity(); return in; } private: /** * Auxiliary class for I/O of Honda flux interpolator. */ struct JHondaFluxInterpolatorHelper : public JProperties { /** * Constructor. * * \param interpolator Honda flux interpolator * \param eqpars equation parameters */ template JHondaFluxInterpolatorHelper(const JHondaFluxInterpolator_t& interpolator, const JEquationParameters& eqpars) : JProperties(eqpars, 1) { (*this)[JEvtWeightFactor::getTypeKey()] = "Honda flux interpolator"; this->insert(gmake_property(interpolator.table)); } }; std::string table; //!< Table filename JRange_t coszRange; //!< Zenith angle range JRange_t phiRange; //!< Azimuth angle range [deg] using JMessage::debug; //!< Debug level }; /** * Template specialisation for interpolator of azimuth-averaged Honda flux table. */ template class JCoszFunctionalMap_t, template class JEnergyFunctionalMap_t> class JHondaFluxInterpolator >, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> final : public JMultiFunction >, typename JMAPLIST::maplist>, public JClonable >, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >, public JMessage >, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> > { public: typedef typename JMAPLIST::maplist JFunctionalMaplist_t; typedef JConstantFunction1D > JPhiFunction_t; typedef JHondaFluxInterpolator interpolator_type; typedef JMultiFunction multifunction_type; enum { NUMBER_OF_DIMENSIONS = multifunction_type::NUMBER_OF_DIMENSIONS }; typedef typename multifunction_type::abscissa_type abscissa_type; typedef typename multifunction_type::argument_type argument_type; typedef typename multifunction_type::result_type result_type; typedef typename multifunction_type::value_type value_type; typedef typename multifunction_type::multimap_type multimap_type; typedef typename multifunction_type::super_const_iterator super_const_iterator; typedef typename multifunction_type::super_iterator super_iterator; typedef typename multifunction_type::function_type function_type; typedef JRange JRange_t; /** * Default constructor. */ JHondaFluxInterpolator() : table(), coszRange(0.0, -1.0) {} /** * Constructor. * * \param fileName oscillation probability table filename */ JHondaFluxInterpolator(const char* fileName) : JHondaFluxInterpolator() { this->load(fileName); this->check_validity(); } /** * Load oscillation probability table from file. * * \param fileName oscillation probability table filename */ void load(const char* const fileName) { using namespace std; using namespace JPP; try { NOTICE("Loading Honda flux table from file " << fileName << "... " << flush); multimap_type& zmap = static_cast(*this); ifstream ifs(fileName); JHondaAngularBinSpecs binspecs; for (string line; getline(ifs, line); ) { double energy = 0.0; result_type values; istringstream iss1(line); const double cosz = binspecs.coszRange.getCentralValue(); coszRange.include(cosz); if (iss1 >> energy >> values) { const double log10E = log10(energy); for (typename result_type::iterator i = values.begin(); i != values.end(); ++i) { *i = log(*i); } zmap[log10E][cosz] = values; } else { const size_t p1 = line.find (HONDA_BINSPECS_BEGIN); const size_t p2 = line.rfind(HONDA_BINSPECS_END); if (p1 != string::npos && p2 != string::npos) { const string binspecstr = line.substr(p1 + 1, p2 - p1 - 1); istringstream iss2(binspecstr); iss2 >> binspecs; } } } NOTICE("OK" << endl); } catch(const exception& error) { THROW(JFileReadException, "JHondaFluxInterpolator::load(): Error reading file " << fileName); } } /** * Check whether this interpolator is valid. * * \return true if valid; else false */ bool is_valid() const override final { return this->size() > 0 && coszRange.getLength(); } /** * Get diffuse flux for given particle PDG-identifier, energy and zenith-angle. * * \param type PDG particle type * \param log10E logarithmic neutrino energy [GeV] * \param costh cosine zenith angle * \return diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1] */ double dNdEdOmega(const int type, const double log10E, const double costh) const override final { using namespace std; using namespace JPP; static const double epsilon_costh = 1.0e-3; const multifunction_type& function = static_cast(*this); const double _costh = min(max(costh, -1.0), 1.0); result_type lnFluxes; if (coszRange(_costh)) { lnFluxes = function(log10E, _costh); } else { // Perform linear extrapolation for edges of cosine zenith angle range if (_costh > coszRange.getUpperLimit()) { const abscissa_type x1 = coszRange.getUpperLimit(); const abscissa_type x2 = coszRange.getUpperLimit() - epsilon_costh; const result_type& y1 = function(log10E, x1); const result_type& y2 = function(log10E, x2); lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_costh - x1); } else { const abscissa_type x1 = coszRange.getLowerLimit(); const abscissa_type x2 = coszRange.getLowerLimit() + epsilon_costh; const result_type& y1 = function(log10E, x1); const result_type& y2 = function(log10E, x2); lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_costh - x1); } } switch(type) { case TRACK_TYPE_NUMU: return exp(lnFluxes[0]); case TRACK_TYPE_ANTINUMU: return exp(lnFluxes[1]); case TRACK_TYPE_NUE: return exp(lnFluxes[2]); case TRACK_TYPE_ANTINUE: return exp(lnFluxes[3]); default: return 0.0; } } /** * Get properties of this class. * * \param eqpars equation parameters */ JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) override final { return JHondaFluxInterpolatorHelper(*this, eqpars); } /** * Get properties of this class. * * \param eqpars equation parameters */ JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) const override final { return JHondaFluxInterpolatorHelper(*this, eqpars); } /** * Stream input. * * \param in input stream * \return input stream */ std::istream& read(std::istream& in) override final { using namespace std; streampos pos = in.tellg(); if (!(in >> table)) { in.clear(); in.seekg(pos); JProperties properties = getProperties(); in >> properties; } load(table.c_str()); this->check_validity(); return in; } private: /** * Auxiliary class for I/O of Honda flux interpolator. */ struct JHondaFluxInterpolatorHelper : public JProperties { /** * Constructor. * * \param interpolator Honda flux interpolator * \param eqpars equation parameters */ template JHondaFluxInterpolatorHelper(const JHondaFluxInterpolator_t& interpolator, const JEquationParameters& eqpars) : JProperties(eqpars, 1) { (*this)[JEvtWeightFactor::getTypeKey()] = "Honda flux interpolator"; this->insert(gmake_property(interpolator.table)); } }; std::string table; //!< Table filename JRange_t coszRange; //!< Cosine zenith angle range using JMessage::debug; //!< Debug level }; /** * Alias template definition for 2-dimensional Honda flux interpolator. */ template class JCoszFunctionalMap_t = JPolint1FunctionalMap, template class JEnergyFunctionalMap_t = JPolint1FunctionalMap> using JHondaFluxInterpolator2D = JHondaFluxInterpolator >, JCoszFunctionalMap_t, JEnergyFunctionalMap_t>; } namespace JEEP { /** * JMessage template specialization for oscillation probability interpolators. */ template class JCoszFunctionalMap_t, template class JEnergyFunctionalMap_t> struct JMessage > { static int debug; }; /** * Default verbosity for oscillation probability interpolators. */ template class JCoszFunctionalMap_t, template class JEnergyFunctionalMap_t> int JMessage >::debug = (int) notice_t; } #endif