diff --git a/cmake/NEUTSetup.cmake b/cmake/NEUTSetup.cmake index 1e9b0a4..97c9034 100644 --- a/cmake/NEUTSetup.cmake +++ b/cmake/NEUTSetup.cmake @@ -1,147 +1,146 @@ # Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret ################################################################################ # This file is part of NUISANCE. # # NUISANCE is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # NUISANCE is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with NUISANCE. If not, see . ################################################################################ if(NEUT_ROOT STREQUAL "") cmessage(FATAL_ERROR "Variable NEUT_ROOT is not defined. Please export environment variable NEUT_ROOT or configure with -DNEUT_ROOT=/path/to/NEUT. This must be set to point to a prebuilt NEUT instance.") endif() if(CERN STREQUAL "") cmessage(FATAL_ERROR "Variable CERN is not defined. Please export environment variable CERN or configure with -DCERN=/path/to/CERNLIB. This must be set to point to a prebuilt CERNLIB instance.") endif() if(CERN_LEVEL STREQUAL "") cmessage(FATAL_ERROR "Variable CERN_LEVEL is not defined. Please export environment variable CERN_LEVEL or configure with -DCERN_LEVEL=XXXX (likely to be 2005).") endif() if(${NEUT_VERSION} VERSION_LESS 5.4.0) set(NEUT_LIB_DIR ${NEUT_ROOT}/lib/Linux_pc) else() set(NEUT_LIB_DIR ${NEUT_ROOT}/lib) endif() set(NEUT_CLASS ${NEUT_ROOT}/src/neutclass) LIST(APPEND NEUT_CXX_FLAGS -DUSE_NEUT -DNEUT_VERSION=${NEUT_VERSION}) if(${NEUT_VERSION} VERSION_GREATER 5.4.0) LIST(APPEND NEUT_CXX_FLAGS -DNEUT_COMMON_QEAV) endif() LIST(APPEND NEUT_INCLUDE_DIRS ${NEUT_ROOT}/include ${NEUT_ROOT}/src/neutclass) LIST(APPEND NEUT_LINK_DIRS ${NEUT_LIB_DIR} ${CERN}/${CERN_LEVEL}/lib) -if(USE_GENERATOR_REWEIGHT) - LIST(APPEND NEUT_INCLUDE_DIRS - ${NEUT_ROOT}/src/reweight) - LIST(APPEND NEUT_LINK_DIRS - ${NEUT_ROOT}/src/reweight) - LIST(APPEND NEUT_LIBS NReWeight) -endif() - if(${NEUT_VERSION} VERSION_GREATER 5.4.1.999) LIST(APPEND NEUT_LIBS neutcore_5.4.2 nuccorspl_5.4.2 #typo in NEUT, may hopefully disappear nuceff_5.4.2 partnuck_5.4.2 skmcsvc_5.4.2 tauola_5.4.2 HT2p2h_5.4.0 N1p1h_5.4.0) LIST(APPEND NEUT_CXX_FLAGS -DNEUT_COMMON_QEAV) elseif(${NEUT_VERSION} VERSION_GREATER 5.3.999) LIST(APPEND NEUT_LIBS neutcore_5.4.0 nuccorspl_5.4.0 #typo in NEUT, may hopefully disappear nuceff_5.4.0 partnuck_5.4.0 skmcsvc_5.4.0 tauola_5.4.0 HT2p2h_5.4.0 N1p1h_5.4.0) else() LIST(APPEND NEUT_LIBS - NReWeight neutcore nuccorrspl nuceff partnuck skmcsvc tauola) endif() LIST(APPEND NEUT_LIBS jetset74 pdflib804 mathlib packlib pawlib gfortran) +if(USE_GENERATOR_REWEIGHT) + LIST(APPEND NEUT_INCLUDE_DIRS + ${NEUT_ROOT}/src/reweight) + LIST(APPEND NEUT_LINK_DIRS + ${NEUT_ROOT}/src/reweight) + LIST(APPEND NEUT_LIBS NReWeight) +endif() + set(NEUT_ROOT_SHAREDOBJS) LIST(APPEND NEUT_ROOT_SHAREDOBJS ${NEUT_CLASS}/neutctrl.so ${NEUT_CLASS}/neutfsivert.so) # Check for new versions of NEUT with NUCLEON FSI if(EXISTS "${NEUT_CLASS}/neutnucfsistep.so") set(NEUT_NUCFSI 1) LIST(APPEND NEUT_CXX_FLAGS -DNEUT_NUCFSI_ENABLED) LIST(APPEND NEUT_ROOT_SHAREDOBJS ${NEUT_CLASS}/neutnucfsistep.so ${NEUT_CLASS}/neutnucfsivert.so ) endif() if(${NEUT_VERSION} VERSION_LESS 5.4.0) LIST(APPEND NEUT_ROOT_SHAREDOBJS ${NEUT_CLASS}/neutrootTreeSingleton.so) endif() LIST(APPEND NEUT_ROOT_SHAREDOBJS ${NEUT_CLASS}/neutvtx.so ${NEUT_CLASS}/neutfsipart.so ${NEUT_CLASS}/neutpart.so ${NEUT_CLASS}/neutvect.so ) foreach(so ${NEUT_ROOT_SHAREDOBJS}) get_filename_component(SONAME ${so} NAME_WE) add_library(${SONAME} SHARED IMPORTED) set_property(TARGET ${SONAME} PROPERTY IMPORTED_LOCATION ${so}) LIST(APPEND NEUT_IMPORTED_TARGETS ${SONAME}) endforeach() PrefixList(NEUT_LINK_DIRS "-L" ${NEUT_LINK_DIRS}) cmessage(STATUS "NEUT") cmessage(STATUS " Version : ${NEUT_VERSION}") cmessage(STATUS " Flags : ${NEUT_CXX_FLAGS}") cmessage(STATUS " Includes : ${NEUT_INCLUDE_DIRS}") cmessage(STATUS " Link Dirs : ${NEUT_LINK_DIRS}") cmessage(STATUS " Libs : ${NEUT_LIBS}") cmessage(STATUS " SOs : ${NEUT_IMPORTED_TARGETS}") diff --git a/cmake/NuWroSetup.cmake b/cmake/NuWroSetup.cmake index daa9c4c..13dd274 100644 --- a/cmake/NuWroSetup.cmake +++ b/cmake/NuWroSetup.cmake @@ -1,39 +1,41 @@ # Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret ################################################################################ # This file is part of NUISANCE. # # NUISANCE is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # NUISANCE is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with NUISANCE. If not, see . ################################################################################ if(NUWRO STREQUAL "") cmessage(FATAL_ERROR "Variable NUWRO is not defined. " "This must be set to point to a prebuilt NuWro instance.") endif() LIST(APPEND NUWRO_CXX_FLAGS -DUSE_NUWRO -Wno-sign-compare -Wno-unused-variable -Wno-reorder) LIST(APPEND NUWRO_INCLUDE_DIRS ${NUWRO}/src) add_library(NuWro_event1 SHARED IMPORTED) set_property(TARGET NuWro_event1 PROPERTY IMPORTED_LOCATION ${NUWRO}/bin/event1.so) LIST(APPEND NUWRO_IMPORTED_TARGETS NuWro_event1) set(NEED_PYTHIA6 TRUE) set(NEED_ROOTPYTHIA6 TRUE) cmessage(STATUS "NuWro") cmessage(STATUS " Flags : ${NUWRO_CXX_FLAGS}") cmessage(STATUS " Includes : ${NUWRO_INCLUDE_DIRS}") cmessage(STATUS " SOs : ${NUWRO_IMPORTED_TARGETS}") + +INSTALL(FILES ${NUWRO}/src/event1dict_rdict.pcm DESTINATION lib) diff --git a/src/generator/CMakeLists.txt b/src/generator/CMakeLists.txt index 324965c..fdb7fa5 100644 --- a/src/generator/CMakeLists.txt +++ b/src/generator/CMakeLists.txt @@ -1,17 +1,20 @@ add_subdirectory(input) add_subdirectory(utility) -add_subdirectory(variation) + +if(DEFINED USE_GENERATOR_REWEIGHT AND USE_GENERATOR_REWEIGHT) + add_subdirectory(variation) +endif() SET(input_implementation_files GeneratorManager.cxx) SET(input_header_files GeneratorManager.hxx) add_library(nuis_generator SHARED ${input_implementation_files}) install(TARGETS nuis_generator DESTINATION lib) install(FILES ${input_header_files} DESTINATION include/generator) LIST(APPEND GENERATOR_DEPENDENT_TARGETS nuis_generator) SET(GENERATOR_DEPENDENT_TARGETS ${GENERATOR_DEPENDENT_TARGETS} PARENT_SCOPE) diff --git a/src/utility/HistogramUtility.hxx b/src/utility/HistogramUtility.hxx index d9bb5d8..9e038df 100644 --- a/src/utility/HistogramUtility.hxx +++ b/src/utility/HistogramUtility.hxx @@ -1,471 +1,480 @@ // Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #pragma once #include "utility/ROOTUtility.hxx" #include "exception/exception.hxx" #include "string_parsers/from_string.hxx" #include "fhiclcpp/ParameterSet.h" #include "TAxis.h" #include "TH1D.h" #include "TH1F.h" #include "TH2D.h" #include "TH2F.h" #include "TH2Poly.h" #include #include #include #include namespace nuis { namespace utility { NEW_NUIS_EXCEPT(unimplemented_GetHistogram_method); NEW_NUIS_EXCEPT(invalid_histogram_descriptor); NEW_NUIS_EXCEPT(invalid_histogram_name); NEW_NUIS_EXCEPT(failed_to_clone); bool IsFlowBin(TAxis const &ax, Int_t bin_it); bool IsInHistogramRange(TAxis const &ax, double v); template struct HType_traits {}; template <> struct HType_traits { using type = TH1; static size_t const NDim = 1; using NumericT = double; static std::string name() { return "TH1"; } }; template <> struct HType_traits { using type = TH1D; static size_t const NDim = 1; using NumericT = double; static std::string name() { return "TH1D"; } }; template <> struct HType_traits { using type = TH1F; static size_t const NDim = 1; using NumericT = float; static std::string name() { return "TH1F"; } }; template <> struct HType_traits { using type = TH2; static size_t const NDim = 2; using NumericT = double; static std::string name() { return "TH2"; } }; template <> struct HType_traits { using type = TH2D; static size_t const NDim = 2; using NumericT = double; static std::string name() { return "TH2D"; } }; template <> struct HType_traits { using type = TH2F; static size_t const NDim = 2; using NumericT = float; static std::string name() { return "TH2F"; } }; template <> struct HType_traits { using type = TH2Poly; static size_t const NDim = 0; using NumericT = double; static std::string name() { return "TH2Poly"; } }; template struct HType_Helper {}; template <> struct HType_Helper<1, void> { using type = TH1; static size_t const NDim = HType_traits::NDim; using NumericT = HType_traits::NumericT; static std::string name() { return HType_traits::name(); } }; template <> struct HType_Helper<1, double> { using type = TH1D; static size_t const NDim = HType_traits::NDim; using NumericT = HType_traits::NumericT; static std::string name() { return HType_traits::name(); } }; template <> struct HType_Helper<1, float> { using type = TH1F; static size_t const NDim = HType_traits::NDim; using NumericT = HType_traits::NumericT; static std::string name() { return HType_traits::name(); } }; template <> struct HType_Helper<2, void> { using type = TH2; static size_t const NDim = HType_traits::NDim; using NumericT = HType_traits::NumericT; static std::string name() { return HType_traits::name(); } }; template <> struct HType_Helper<2, double> { using type = TH2D; static size_t const NDim = HType_traits::NDim; using NumericT = HType_traits::NumericT; static std::string name() { return HType_traits::name(); } }; template <> struct HType_Helper<2, float> { using type = TH2F; static size_t const NDim = HType_traits::NDim; using NumericT = HType_traits::NumericT; static std::string name() { return HType_traits::name(); } }; template struct TH_Helper {}; template struct TH_Helper::NDim == 1>::type> { static size_t const NDim = HType_traits::NDim; using NumericT = typename HType_traits::NumericT; static std::string name() { return HType_traits::name(); } static Int_t NbinsIncludeFlow(HT const &h) { return h.GetXaxis()->GetNbins() + 2; } static Int_t Nbins(HT const &h) { return h.GetXaxis()->GetNbins(); } static bool IsFlowBin(HT const &h, Int_t bin_it) { return nuis::utility::IsFlowBin(*h.GetXaxis(), bin_it); } static void Fill(HT &h, std::array const &v, double w = 1) { h.Fill(v[0], w); } static void Scale(HT &h, NumericT SF, char const *opt = "") { h.Scale(SF, opt); } static double Integral(HT const &h, char const *opt = "") { return h.Integral(opt); } static Int_t NbinsIncludeFlow(std::unique_ptr const &h) { return NbinsIncludeFlow(*h); } static Int_t Nbins(std::unique_ptr const &h) { return Nbins(*h); } static bool IsFlowBin(std::unique_ptr const &h, Int_t bin_it) { return IsFlowBin(*h, bin_it); } static void Fill(std::unique_ptr &h, std::array const &v, double w = 1) { Fill(*h, v, w); } static void Scale(std::unique_ptr &h, NumericT SF, char const *opt = "") { Scale(*h, SF, opt); } static double Integral(std::unique_ptr const &h, char const *opt = "") { return Integral(*h, opt); } }; template struct TH_Helper::NDim == 2>::type> { static size_t const NDim = HType_traits::NDim; using NumericT = typename HType_traits::NumericT; static std::string name() { return HType_traits::name(); } // TH2 *************************************************************** static Int_t NbinsIncludeFlow(HT const &h) { return (h.GetXaxis()->GetNbins() + 2) * (h.GetYaxis()->GetNbins() + 2); } static Int_t Nbins(HT const &h) { return (h.GetXaxis()->GetNbins()) * (h.GetYaxis()->GetNbins()); } static bool IsFlowBin(HT const &h, Int_t xbin_it, Int_t ybin_it) { return nuis::utility::IsFlowBin(*h.GetXaxis(), xbin_it) || nuis::utility::IsFlowBin(*h.GetYaxis(), ybin_it); } static void Fill(HT &h, std::array const &v, double w = 1) { h.Fill(v[0], v[1], w); } static void Scale(HT &h, NumericT SF, char const *opt = "") { h.Scale(SF, opt); } static double Integral(HT const &h, char const *opt = "") { return h.Integral(opt); } static Int_t NbinsIncludeFlow(std::unique_ptr const &h) { return NbinsIncludeFlow(*h); } static Int_t Nbins(std::unique_ptr const &h) { return Nbins(*h); } static bool IsFlowBin(std::unique_ptr const &h, Int_t xbin_it, Int_t ybin_it) { return IsFlowBin(*h, xbin_it, ybin_it); } static void Fill(std::unique_ptr &h, std::array const &v, double w = 1) { Fill(*h, v, w); } static void Scale(std::unique_ptr &h, NumericT SF, char const *opt = "") { Scale(*h, SF, opt); } static double Integral(std::unique_ptr const &h, char const *opt = "") { return Integral(*h, opt); } }; template struct TH_Helper< HT, typename std::enable_if::value>::type> { static size_t const NDim = 2; using NumericT = typename HType_traits::NumericT; static std::string name() { return HType_traits::name(); } // TH2Poly *************************************************************** static Int_t NbinsIncludeFlow(HT const &h) { return h.GetNumberOfBins() + 9; } static Int_t Nbins(HT const &h) { return h.GetNumberOfBins(); } static bool IsFlowBin(HT const &h, Int_t bin_it) { return (bin_it < 0); } static void Fill(HT &h, std::array const &v, double w = 1) { h.Fill(v[0], v[1], w); } static void Scale(HT &h, NumericT SF, char const *opt = "") { bool width = (std::string(opt).find("width") != std::string::npos); size_t nbins = Nbins(h); for (size_t bin_it = 0; bin_it < nbins; ++bin_it) { double bin_area = 1; if (width) { TH2PolyBin *poly_bin = dynamic_cast(h.GetBins()->At(bin_it)); bin_area = poly_bin->GetArea(); } h.SetBinContent(bin_it + 1, h.GetBinContent(bin_it + 1) * (SF / bin_area)); h.SetBinError(bin_it + 1, h.GetBinError(bin_it + 1) * (SF / bin_area)); } } static double Integral(HT &h, char const *opt = "") { bool width = (std::string(opt).find("width") != std::string::npos); size_t nbins = Nbins(h); double integral = 0; for (size_t bin_it = 0; bin_it < nbins; ++bin_it) { double bin_area = 1; if (width) { TH2PolyBin *poly_bin = dynamic_cast(h.GetBins()->At(bin_it)); bin_area = poly_bin->GetArea(); } integral += h.GetBinContent(bin_it + 1) * bin_area; } return integral; } static Int_t NbinsIncludeFlow(std::unique_ptr const &h) { return NbinsIncludeFlow(*h); } static Int_t Nbins(std::unique_ptr const &h) { return Nbins(*h); } static bool IsFlowBin(std::unique_ptr const &h, Int_t bin_it) { return IsFlowBin(*h, bin_it); } static void Fill(std::unique_ptr &h, std::array const &v, double w = 1) { Fill(*h, v, w); } static void Scale(std::unique_ptr &h, NumericT SF, char const *opt = "") { Scale(*h, SF, opt); } static double Integral(std::unique_ptr &h, char const *opt = "") { return Integral(*h, opt); } }; template void Clear(typename std::enable_if::NDim != 0, HT>::type &h) { for (Int_t bin_it = 0; bin_it < TH_Helper::NbinsIncludeFlow(h); ++bin_it) { h.SetBinContent(bin_it, 0); h.SetBinError(bin_it, 0); } } template void Clear( typename std::enable_if::value, HT>::type &h) { h.ClearBinContents(); } template inline std::unique_ptr Clone(HT const &source, bool clear = false, std::string const &clone_name = "") { std::unique_ptr target(dynamic_cast( source.Clone(clone_name.size() ? clone_name.c_str() : ""))); if (!target) { throw failed_to_clone() << "[ERROR]: Failed to clone a " << TH_Helper::name() << "."; } target->SetDirectory(nullptr); if (clear) { Clear(*target); } return target; } template inline std::unique_ptr Clone(std::unique_ptr const &source, bool clear = false, std::string const &clone_name = "") { return Clone(*source, clear, clone_name); } template inline std::unique_ptr GetHistogramFromROOTFile(TFile_ptr &f, std::string const &hname, bool ThrowIfMissing = true) { TObject *obj = f->Get(hname.c_str()); - if(std::string(obj->ClassName()) == "TList"){ + if (!obj) { + if (ThrowIfMissing) { + throw invalid_histogram_name() + << "[ERROR]: Failed to get " << TH_Helper::name() << " named " + << std::quoted(hname) << " from input file " + << std::quoted(f->GetName()); + } else { + return nullptr; + } + } + if (std::string(obj->ClassName()) == "TList") { obj = static_cast(obj)->At(0); } - HT *h = dynamic_cast(obj); if (!h) { if (ThrowIfMissing) { throw invalid_histogram_name() - << "[ERROR]: Failed to get " << TH_Helper::name() << " named " - << std::quoted(hname) << " from input file " + << "[ERROR]: Failed to cast " << std::quoted(hname) << " to " + << TH_Helper::name() << " from input file " << std::quoted(f->GetName()); } else { return nullptr; } } std::unique_ptr clone = Clone(*h); return clone; } template inline std::unique_ptr GetHistogramFromROOTFile(std::string const &fname, std::string const &hname, bool ThrowIfMissing = true) { TFile_ptr temp = CheckOpenTFile(fname, "READ"); return GetHistogramFromROOTFile(temp, hname, ThrowIfMissing); } template std::unique_ptr GetHistogram(std::string const &input_descriptor) { std::vector split_descriptor = fhicl::string_parsers::ParseToVect(input_descriptor, ";", true, true); if (split_descriptor.size() == 1) { // assume text throw unimplemented_GetHistogram_method(); } else if (split_descriptor.size() == 2) { return GetHistogramFromROOTFile(split_descriptor[0], split_descriptor[1]); } else { throw invalid_histogram_descriptor() << "[ERROR]: Failed to parse histogram descriptor: " << std::quoted(input_descriptor) << " as an input histogram (Text/ROOT)."; } } struct PolyBinSpecifier { double X, Y; bool UseXAxis; }; constexpr PolyBinSpecifier XPolyBinSpec(double X, double Y) { return PolyBinSpecifier{X + std::numeric_limits::epsilon() * 1E2, Y + std::numeric_limits::epsilon() * 1E2, true}; } constexpr PolyBinSpecifier YPolyBinSpec(double X, double Y) { return PolyBinSpecifier{X + std::numeric_limits::epsilon() * 1E2, Y + std::numeric_limits::epsilon() * 1E2, false}; } std::vector> GetTH2PolySlices( std::unique_ptr &hinp, std::vector> const &BinsSpecifiers); template inline typename std::enable_if::NDim == 1, std::unique_ptr>::type BuildHistFromFHiCL(fhicl::ParameterSet const &ps) { std::array xaxis = ps.get>("xaxis_descriptor"); std::unique_ptr rtn = std::make_unique( ps.get("name", "").c_str(), ps.get("title", "").c_str(), xaxis[0], xaxis[1], xaxis[2]); rtn->SetDirectory(nullptr); return rtn; } template inline typename std::enable_if::NDim == 2, std::unique_ptr>::type BuildHistFromFHiCL(fhicl::ParameterSet const &ps) { std::array xaxis = ps.get>("xaxis_descriptor"); std::array yaxis = ps.get>("yaxis_descriptor"); std::unique_ptr rtn = std::make_unique(ps.get("name", "").c_str(), ps.get("title", "").c_str(), xaxis[0], xaxis[1], xaxis[2], yaxis[0], yaxis[1], yaxis[2]); rtn->SetDirectory(nullptr); return rtn; } static bool const kYSlice = true; static bool const kXSlice = false; void SliceNorm(std::unique_ptr &hist, bool AlongY = kYSlice, char const *opt = ""); } // namespace utility } // namespace nuis