diff --git a/cmake/NEUTSetup.cmake b/cmake/NEUTSetup.cmake index 098c889..acf03aa 100644 --- a/cmake/NEUTSetup.cmake +++ b/cmake/NEUTSetup.cmake @@ -1,214 +1,217 @@ # 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 . ################################################################################ -find_program(neut-config NEUTCONFIGFOUND) +find_program(NEUTCONFIGFOUND NAMES neut-config) LIST(APPEND EXTRA_CXX_FLAGS -DNEED_FILL_NEUT_COMMONS) SET(HAVENEUTCONFIG FALSE) # We are dealing with shiny NEUT if(NOT "${NEUTCONFIGFOUND} " STREQUAL " ") SET(HAVENEUTCONFIG TRUE) + cmessage(STATUS "Found neut-config, using it to determine configuration.") +else() + cmessage(STATUS "Failed to find neut-config, assuming older NEUT build.") endif() if(HAVENEUTCONFIG) execute_process (COMMAND neut-config --version OUTPUT_VARIABLE NEUT_VERSION OUTPUT_STRIP_TRAILING_WHITESPACE) execute_process (COMMAND neut-config --incdir OUTPUT_VARIABLE NEUT_INCLUDE_DIRS OUTPUT_STRIP_TRAILING_WHITESPACE) execute_process (COMMAND neut-config --libdir OUTPUT_VARIABLE NEUT_LINK_DIRS OUTPUT_STRIP_TRAILING_WHITESPACE) GETLIBDIRS(neut-config --cernflags CERN_LIB_DIR) LIST(APPEND NEUT_LINK_DIRS ${CERN_LIB_DIR}) GETLIBS(neut-config --cernflags CERN_LIBS) if(USE_GENERATOR_REWEIGHT) execute_process (COMMAND neut-config --rwlibflags OUTPUT_VARIABLE NEUT_RWLIBS OUTPUT_STRIP_TRAILING_WHITESPACE) GETLIBS(neut-config --rwlibflags NEUT_RWLIBS) LIST(APPEND NEUT_LIBS ${NEUT_RWLIBS}) else() GETLIBS(neut-config --libflags NEUT_GENLIBS) GETLIBS(neut-config --iolibflags NEUT_LIBS) LIST(APPEND NEUT_LIBS ${NEUT_IOLIBS}) LIST(APPEND NEUT_LIBS ${NEUT_GENLIBS}) endif() LIST(APPEND NEUT_LIBS ${CERN_LIBS};gfortran) PrefixList(NEUT_LINK_DIRS "-L" ${NEUT_LINK_DIRS}) LIST(APPEND NEUT_CXX_FLAGS -DUSE_NEUT -DNEUT_VERSION=${NEUT_VERSION}) 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}") else() # Everything better be set up already 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_INCLUDE_DIRS ${NEUT_ROOT}/include ${NEUT_ROOT}/src/neutclass) LIST(APPEND NEUT_LINK_DIRS ${NEUT_LIB_DIR} ${CERN}/${CERN_LEVEL}/lib) 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_EQUAL 5.4.1) LIST(APPEND NEUT_LIBS neutcore_5.4.1 nuccorspl_5.4.1 #typo in NEUT, may hopefully disappear nuceff_5.4.1 partnuck_5.4.1 skmcsvc_5.4.1 tauola_5.4.1 HT2p2h_5.4.1 N1p1h_5.4.1) 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 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(REVERSE NEUT_LIBS) LIST(APPEND NEUT_LIBS NReWeight) LIST(REVERSE NEUT_LIBS) 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}) LIST(APPEND NEUT_CXX_FLAGS -DUSE_NEUT -DNEUT_VERSION=${NEUT_VERSION}) 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}") endif() SET(USE_NEUT TRUE CACHE BOOL "Whether to enable NEUT support. Requires external libraries. " FORCE) diff --git a/cmake/c++CompilerSetup.cmake b/cmake/c++CompilerSetup.cmake index 5083311..20ffe38 100644 --- a/cmake/c++CompilerSetup.cmake +++ b/cmake/c++CompilerSetup.cmake @@ -1,57 +1,58 @@ # 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 . ################################################################################ set(CXX_WARNINGS -Wall -Wextra) +set(CXX_IGNORE_WARNINGS -Wno-delete-non-virtual-dtor -Wno-unused -Wno-misleading-indentation) -LIST(APPEND EXTRA_CXX_FLAGS ${CXX_WARNINGS} -Werror -Wno-delete-non-virtual-dtor -Wno-unused "-D__FILENAME__=\"$(subst ${CMAKE_SOURCE_DIR}/,,$(abspath $<))\"") +LIST(APPEND EXTRA_CXX_FLAGS ${CXX_WARNINGS} -Werror ${CXX_IGNORE_WARNINGS} "-D__FILENAME__=\"$(subst ${CMAKE_SOURCE_DIR}/,,$(abspath $<))\"") BuildFlagString(NUISANCE_LINK_DIRS "-L" ${EXTRA_LINK_DIRS}) LIST(APPEND EXTRA_LIBS dl) BuildLibraryFlagString(STR_EXTRA_LIBS ${EXTRA_LIBS}) BuildFlagString(STR_EXTRA_SHAREDOBJS " " ${EXTRA_SHAREDOBJS}) #This ends up holding all of the libraries and search paths for extenal dependencies CatStringsIfNotEmpty(NUISANCE_DEPEND_LIBS ${STR_EXTRA_SHAREDOBJS} ${STR_EXTRA_LIBS}) BuildFlagString(STR_EXTRA_LINK_FLAGS " " ${EXTRA_LINK_FLAGS}) CatStringsIfNotEmpty(CMAKE_LINK_FLAGS ${CMAKE_LINK_FLAGS} ${STR_EXTRA_LINK_FLAGS}) get_directory_property(NUISANCE_INCLUDE_DIRS INCLUDE_DIRECTORIES) BuildFlagString(NUISANCE_CXX_FLAGS " " ${EXTRA_CXX_FLAGS}) CatStringsIfNotEmpty(NUISANCE_CXX_FLAGS ${CMAKE_CXX_FLAGS} ${NUISANCE_CXX_FLAGS} ) if (VERBOSE) cmessage (STATUS "C++ Compiler : ${CXX_COMPILER_NAME}") cmessage (STATUS " Flags : ${NUISANCE_CXX_FLAGS}") cmessage (STATUS " Release Flags : ${CMAKE_CXX_FLAGS_RELEASE}") cmessage (STATUS " Debug Flags : ${CMAKE_CXX_FLAGS_DEBUG}") cmessage (STATUS " Include Dirs : ${NUISANCE_INCLUDE_DIRS}") cmessage (STATUS " Linker Flags : ${CMAKE_LINK_FLAGS}") cmessage (STATUS " Link Dirs : ${NUISANCE_LINK_DIRS}") cmessage (STATUS " Lib Flags : ${NUISANCE_DEPEND_LIBS}") endif() add_compile_options(${EXTRA_CXX_FLAGS}) diff --git a/src/app/nuisevsum.cxx b/src/app/nuisevsum.cxx index 5461e22..c50f41b 100644 --- a/src/app/nuisevsum.cxx +++ b/src/app/nuisevsum.cxx @@ -1,79 +1,88 @@ #include "config/GlobalConfiguration.hxx" #include "input/IInputHandler.hxx" #include "event/MinimalEvent.hxx" #include "samples/IEventProcessor.hxx" #include "plugins/Instantiate.hxx" #include "exception/exception.hxx" #include "fhiclcpp/make_ParameterSet.h" #include "string_parsers/from_string.hxx" #include NEW_NUIS_EXCEPT(invalid_cli_arguments); size_t NMax = std::numeric_limits::max(); std::string input_file; std::string input_type; +bool show_intermediates = false; void SayUsage(char const *argv[]) { std::cout << "[USAGE]: " << argv[0] << "\n" "\t-i : Input file passed to named " "IInputHandler instance \n" "\t-H : Name of IInputHandler subclass " "capable of reading NUISANCE events from the argument of -i.\n" "\t-n : Maximum number of events to " "read. Will read entire input file by default.\n" + "\t--show-intermediates : Show pre-FSI particles where" + " relevant.\n" << std::endl; } void handleOpts(int argc, char const *argv[]) { int opt = 1; while (opt < argc) { if ((std::string(argv[opt]) == "-?") || (std::string(argv[opt]) == "--help")) { SayUsage(argv); exit(0); } else if (std::string(argv[opt]) == "-i") { input_file = argv[++opt]; } else if (std::string(argv[opt]) == "-H") { input_type = argv[++opt]; } else if (std::string(argv[opt]) == "-n") { NMax = fhicl::string_parsers::str2T(argv[++opt]); + } else if (std::string(argv[opt]) == "--show-intermediates") { + show_intermediates = true; } else { std::cout << "[ERROR]: Unknown option: " << argv[opt] << std::endl; SayUsage(argv); exit(1); } opt++; } } int main(int argc, char const *argv[]) { nuis::config::EnsureConfigurationRead("nuis.global.config.fcl"); handleOpts(argc, argv); if (!input_type.length() || !input_file.length()) { SayUsage(argv); throw invalid_cli_arguments() << "[ERROR]: Require both -i and -H cli options to be passed."; } fhicl::ParameterSet sample_config; sample_config.put("input_type", input_type); sample_config.put("file", input_file); + if (show_intermediates) { + sample_config.put("keep_intermediates", true); + } - nuis::plugins::plugin_traits::unique_ptr_t VerboseEventSummary = - nuis::plugins::Instantiate("VerboseEventSummary"); + nuis::plugins::plugin_traits::unique_ptr_t + VerboseEventSummary = + nuis::plugins::Instantiate("VerboseEventSummary"); VerboseEventSummary->Initialize(sample_config); VerboseEventSummary->ProcessSample(NMax); } diff --git a/src/generator/input/GENIEInputHandler.cxx b/src/generator/input/GENIEInputHandler.cxx index 1bfb67d..f820fd9 100644 --- a/src/generator/input/GENIEInputHandler.cxx +++ b/src/generator/input/GENIEInputHandler.cxx @@ -1,215 +1,222 @@ #include "generator/input/GENIEInputHandler.hxx" #include "generator/utility/GENIEUtility.hxx" #include "utility/PDGCodeUtility.hxx" #ifdef GENIE_V3_INTERFACE +#include "Framework/Conventions/Units.h" #include "Framework/EventGen/EventRecord.h" #include "Framework/GHEP/GHepParticle.h" #include "Framework/GHEP/GHepRecord.h" #else +#include "Conventions/Units.h" #include "EVGCore/EventRecord.h" #include "GHEP/GHepParticle.h" #include "GHEP/GHepRecord.h" #endif #include "fhiclcpp/ParameterSet.h" using namespace nuis; using namespace nuis::event; using namespace nuis::utility; using namespace nuis::genietools; GENIEInputHandler::GENIEInputHandler() : fInputTree(), fGenieNtpl(nullptr), fFileWeight(1) {} GENIEInputHandler::GENIEInputHandler(GENIEInputHandler &&other) : fInputTree(std::move(other.fInputTree)), fReaderEvent(std::move(other.fReaderEvent)), fGenieNtpl(nullptr), fFileWeight(1) {} std::set GENIEInputHandler::GetSplineList() { std::set splinenames; - size_t NEvents = 1000; // GetNEvents(); + size_t NEvents = std::min(size_t(10000), GetNEvents()); size_t ShoutEvery = NEvents / 100; + std::cout << "[INFO]: Determining list of required GENIE splines." + << std::endl; std::cout << "[INFO]: Read " << 0 << "/" << NEvents << " GENIE events." << std::flush; for (size_t ev_it = 0; ev_it < NEvents; ++ev_it) { if (fGenieNtpl) { fGenieNtpl->Clear(); } fInputTree.tree->GetEntry(ev_it); if (ShoutEvery && !(ev_it % ShoutEvery)) { std::cout << "\r[INFO]: Read " << ev_it << "/" << NEvents << " GENIE events." << std::flush; } genie::GHepRecord *GHep = static_cast(fGenieNtpl->event); if (!GHep) { throw invalid_GENIE_event() << "[ERROR]: GENIE event " << ev_it << " failed to contain a GHepRecord"; } splinenames.insert(GHep->Summary()->AsString()); } std::cout << "[INFO]: Having read GENIE events, found " << splinenames.size() << " splines: " << std::endl; for (auto const &sn : splinenames) { std::cout << "\t" << sn << std::endl; } return splinenames; } void GENIEInputHandler::Initialize(fhicl::ParameterSet const &ps) { fInputTree = CheckGetTTree(ps.get("file"), "gtree"); fInputTree.tree->SetBranchAddress("gmcrec", &fGenieNtpl); fKeepIntermediates = ps.get("keep_intermediates", false); fKeepNuclearParticles = ps.get("keep_nuclear_particles", false); fhicl::ParameterSet XSecInfo = ps.get("xsec_info", {}); if (XSecInfo.has_key("weight")) { fFileWeight = XSecInfo.get("weight"); std::cout << "[INFO]: Average GENIE XSecWeight = " << fFileWeight << "" << std::endl; - } else if (XSecInfo.has_key("flux") || XSecInfo.has_key("target") || - XSecInfo.has_key("spline_file") || - XSecInfo.has_key("filter_splines_by_event_tree")) { + } else if (XSecInfo.has_key("spline_file")) { std::cout << "[INFO]: Attempting to build GENIE file weight with input splines." << std::endl; if (XSecInfo.get("filter_splines_by_event_tree", false)) { std::cout << "[INFO]: Ensuring we only read relevant splines... this " "will take a while." << std::endl; fFileWeight = genietools::GetFileWeight(XSecInfo, GetSplineList()) / double(GetNEvents()); } else { fFileWeight = genietools::GetFileWeight(XSecInfo) / double(GetNEvents()); } std::cout << "[INFO]: Average GENIE XSecWeight = " << fFileWeight << "" << std::endl; + } else if (XSecInfo.has_key("flux")) { + fFileWeight = genietools::GetFileWeightFromReconstructedSplines( + XSecInfo, fGenieNtpl, fInputTree.tree); + } else { + fFileWeight = 1; } } void GENIEInputHandler::GetEntry(ev_index_t idx) const { if (idx >= GetNEvents()) { throw IInputHandler::invalid_entry() << "[ERROR]: Attempted to get entry " << idx << " from an InputHandler with only " << GetNEvents(); } if (fGenieNtpl) { fGenieNtpl->Clear(); } fInputTree.tree->GetEntry(idx); } MinimalEvent const &GENIEInputHandler::GetMinimalEvent(ev_index_t idx) const { GetEntry(idx); genie::GHepRecord *GHep = static_cast(fGenieNtpl->event); if (!GHep) { throw invalid_GENIE_event() << "[ERROR]: GENIE event " << idx << " failed to contain a GHepRecord"; } fReaderEvent.mode = GetEventChannel(*GHep); TObjArrayIter iter(GHep); genie::GHepParticle *p; while ((p = (dynamic_cast((iter).Next())))) { if (!p) { continue; } Particle::Status_t state = GetParticleStatus(*p, fReaderEvent.mode); if (state != Particle::Status_t::kPrimaryInitialState) { continue; } if (!IsNeutralLepton(p->Pdg(), pdgcodes::kMatterAntimatter) && !IsChargedLepton(p->Pdg(), pdgcodes::kMatterAntimatter)) { continue; } fReaderEvent.probe_E = p->E() * 1.E3; fReaderEvent.probe_pdg = p->Pdg(); break; } fReaderEvent.XSecWeight = fFileWeight; if (fWeightCache.size() <= idx) { fWeightCache.push_back(fReaderEvent.XSecWeight); } return fReaderEvent; } FullEvent const &GENIEInputHandler::GetFullEvent(ev_index_t idx) const { (void)GetMinimalEvent(idx); fReaderEvent.ClearParticleStack(); genie::GHepRecord *GHep = static_cast(fGenieNtpl->event); unsigned int npart = GHep->GetEntries(); // Fill Particle Stack genie::GHepParticle *p = 0; TObjArrayIter iter(GHep); // Loop over all particles while ((p = (dynamic_cast((iter).Next())))) { if (!p) { continue; } // Get Status Particle::Status_t state = GetParticleStatus(*p, fReaderEvent.mode); if (!fKeepIntermediates && (state == Particle::Status_t::kIntermediate)) { continue; } if (!fKeepNuclearParticles && IsNuclearPDG(p->Pdg())) { continue; } Particle nuis_part; nuis_part.pdg = p->Pdg(); nuis_part.P4 = TLorentzVector(p->Px(), p->Py(), p->Pz(), p->E()) * 1E3; fReaderEvent.ParticleStack[static_cast(state)].particles.push_back( nuis_part); } return fReaderEvent; } double GENIEInputHandler::GetEventWeight(ev_index_t idx) const { if (idx > fWeightCache.size()) { throw weight_cache_miss() << "[ERROR]: Failed to get cached weight for event index: " << idx; } return fWeightCache[idx]; } double GENIEInputHandler::GetXSecScaleFactor( std::pair const &EnuRange) const { throw input_handler_feature_unimplemented() << "[ERROR]: Flux cuts not yet implemented for GENIE input handler."; } size_t GENIEInputHandler::GetNEvents() const { return fInputTree.tree->GetEntries(); } GeneratorManager::Generator_id_t GENIEInputHandler::GetGeneratorId() const { return GeneratorManager::Get().EnsureGeneratorRegistered("GENIE"); } DECLARE_PLUGIN(IInputHandler, GENIEInputHandler); diff --git a/src/generator/utility/GENIEUtility.cxx b/src/generator/utility/GENIEUtility.cxx index 6011309..741d044 100644 --- a/src/generator/utility/GENIEUtility.cxx +++ b/src/generator/utility/GENIEUtility.cxx @@ -1,479 +1,588 @@ // 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 . *******************************************************************************/ #include "generator/utility/GENIEUtility.hxx" #include "generator/utility/GENIESplineReader.hxx" #include "utility/HistogramUtility.hxx" #include "utility/InteractionChannelUtility.hxx" #include "utility/PDGCodeUtility.hxx" #include "utility/ROOTUtility.hxx" #include "utility/StringUtility.hxx" #ifdef GENIE_V3_INTERFACE #include "Framework/GHEP/GHepParticle.h" #include "Framework/GHEP/GHepRecord.h" #include "Framework/GHEP/GHepUtils.h" #include "Framework/ParticleData/PDGCodes.h" #else #include "GHEP/GHepParticle.h" #include "GHEP/GHepRecord.h" #include "GHEP/GHepUtils.h" #include "PDG/PDGCodes.h" #endif #include "fhiclcpp/ParameterSet.h" #include "string_parsers/from_string.hxx" #include "TGraph.h" namespace nuis { using namespace event; namespace genietools { struct NFSParticleCount { size_t NProton; size_t NNeutron; size_t NPip; size_t NPi0; size_t NPim; size_t NOther; }; NFSParticleCount CountPreFSIParticles(genie::GHepRecord const &ev) { // This code in this method is adapted from the GENIE source code found in // GHep/GHepUtils.cxx This method therefore carries the GENIE copyright // licence as copied below: // /// Copyright (c) 2003-2017, GENIE Neutrino MC Generator Collaboration /// For the full text of the license visit http://copyright.genie-mc.org /// or see $GENIE/LICENSE // genie::Target const &tgt = ev.Summary()->InitState().Tgt(); if (!tgt.HitNucIsSet()) { throw invalid_GENIE_event() << "[ERROR]: Failed to get hit nucleon kinematics as it was not " "included in this GHep event. This is a fatal error."; } genie::GHepParticle *FSLep = ev.FinalStatePrimaryLepton(); genie::GHepParticle *ISLep = ev.Probe(); if (!FSLep || !ISLep) { throw invalid_GENIE_event() << "[ERROR]: Failed to find IS and FS lepton in event: " << ev.Summary()->AsString(); } size_t NPi0 = 0, NPip = 0, NPim = 0, NProton = 0, NNeutron = 0, NOther = 0; bool nuclear_target = tgt.IsNucleus(); TIter event_iter(&ev); genie::GHepParticle *p = 0; while ((p = dynamic_cast(event_iter.Next()))) { genie::GHepStatus_t ghep_ist = (genie::GHepStatus_t)p->Status(); int ghep_pdgc = p->Pdg(); int ghep_fm = p->FirstMother(); int ghep_fmpdgc = (ghep_fm == -1) ? 0 : ev.Particle(ghep_fm)->Pdg(); // For nuclear targets use hadrons marked as 'hadron in the nucleus' // which are the ones passed in the intranuclear rescattering // For free nucleon targets use particles marked as 'final state' // but make an exception for decayed pi0's,eta's (count them and not their // daughters) bool decayed = (ghep_ist == genie::kIStDecayedState && (ghep_pdgc == genie::kPdgPi0 || ghep_pdgc == genie::kPdgEta)); bool parent_included = (ghep_fmpdgc == genie::kPdgPi0 || ghep_fmpdgc == genie::kPdgEta); bool count_it = (nuclear_target && ghep_ist == genie::kIStHadronInTheNucleus) || (!nuclear_target && decayed) || (!nuclear_target && ghep_ist == genie::kIStStableFinalState && !parent_included); if (!count_it) { continue; } if (ghep_pdgc == genie::kPdgPiP) { NPip++; } else if (ghep_pdgc == genie::kPdgPiM) { NPim++; } else if (ghep_pdgc == genie::kPdgPi0) { NPi0++; } else if (ghep_pdgc == genie::kPdgProton) { NProton++; } else if (ghep_pdgc == genie::kPdgNeutron) { NNeutron++; } else if (!utility::IsNeutralLepton( ghep_pdgc, utility::pdgcodes::kMatterAntimatter) && !utility::IsChargedLepton( ghep_pdgc, utility::pdgcodes::kMatterAntimatter)) { NOther++; } } return NFSParticleCount{NProton, NNeutron, NPip, NPi0, NPim, NOther}; } Channel_t GetEventChannel(genie::GHepRecord const &gev) { // Electron Scattering if (gev.Summary()->ProcInfo().IsEM()) { if (gev.Summary()->InitState().ProbePdg() == utility::pdgcodes::kElectron) { if (gev.Summary()->ProcInfo().IsQuasiElastic()) { NFSParticleCount fsparts = CountPreFSIParticles(gev); if (fsparts.NProton) { return Channel_t::kNCELP; } else { return Channel_t::kNCELN; } } else if (gev.Summary()->ProcInfo().IsMEC()) { return Channel_t::kNC2p2h; } else if (gev.Summary()->ProcInfo().IsResonant()) { NFSParticleCount fsparts = CountPreFSIParticles(gev); if (fsparts.NOther || ((fsparts.NPip + fsparts.NPi0 + fsparts.NPim) > 1)) { return Channel_t::kNCTransitionMPi; } else if (fsparts.NPip) { return Channel_t::kNCSPP_NPip; } else if (fsparts.NPi0) { return fsparts.NProton ? Channel_t::kNCSPP_PPi0 : Channel_t::kNCSPP_NPi0; } else if (fsparts.NPim) { return Channel_t::kNCSPP_PPim; } return Channel_t::kNCTransitionMPi; } else if (gev.Summary()->ProcInfo().IsDeepInelastic()) { return Channel_t::kNCDIS; } else { std::cout << "Unknown GENIE Electron Scattering Mode!" << std::endl << "ScatteringTypeId = " << gev.Summary()->ProcInfo().ScatteringTypeId() << " " << "InteractionTypeId = " << gev.Summary()->ProcInfo().InteractionTypeId() << std::endl << genie::ScatteringType::AsString( gev.Summary()->ProcInfo().ScatteringTypeId()) << " " << genie::InteractionType::AsString( gev.Summary()->ProcInfo().InteractionTypeId()) << " " << gev.Summary()->ProcInfo().IsMEC() << std::endl; return Channel_t::kUndefined; } } // Weak CC } else if (gev.Summary()->ProcInfo().IsWeakCC()) { // CC MEC if (gev.Summary()->ProcInfo().IsMEC()) { if (utility::IsNeutralLepton(gev.Summary()->InitState().ProbePdg(), utility::pdgcodes::kMatter)) { return Channel_t::kCC2p2h; } else if (utility::IsNeutralLepton(gev.Summary()->InitState().ProbePdg(), utility::pdgcodes::kAntimatter)) { return Channel_t::kCC2p2h_nub; } // CC OTHER } else { return FromNEUTCode(genie::utils::ghep::NeutReactionCode(&gev)); } // Weak NC } else if (gev.Summary()->ProcInfo().IsWeakNC()) { // NC MEC if (gev.Summary()->ProcInfo().IsMEC()) { if (utility::IsNeutralLepton(gev.Summary()->InitState().ProbePdg(), utility::pdgcodes::kMatter)) { return Channel_t::kNC2p2h; } else if (utility::IsNeutralLepton(gev.Summary()->InitState().ProbePdg(), utility::pdgcodes::kAntimatter)) { return Channel_t::kNC2p2h_nub; } // NC OTHER } else { return FromNEUTCode(genie::utils::ghep::NeutReactionCode(&gev)); } } return Channel_t::kUndefined; } Particle::Status_t GetParticleStatus(genie::GHepParticle const &p, Channel_t chan) { /* kIStUndefined = -1, kIStInitialState = 0, / generator-level initial state / kIStStableFinalState = 1, / generator-level final state: particles to be tracked by detector-level MC / kIStIntermediateState = 2, kIStDecayedState = 3, kIStCorrelatedNucleon = 10, kIStNucleonTarget = 11, kIStDISPreFragmHadronicState = 12, kIStPreDecayResonantState = 13, kIStHadronInTheNucleus = 14, / hadrons inside the nucleus: marked for hadron transport modules to act on / kIStFinalStateNuclearRemnant = 15, / low energy nuclear fragments entering the record collectively as a 'hadronic blob' pseudo-particle / kIStNucleonClusterTarget = 16, // for composite nucleons before phase space decay */ Particle::Status_t state = Particle::Status_t::kUnknown; switch (p.Status()) { case genie::kIStNucleonTarget: case genie::kIStInitialState: case genie::kIStCorrelatedNucleon: case genie::kIStNucleonClusterTarget: { state = Particle::Status_t::kPrimaryInitialState; break; } case genie::kIStStableFinalState: { state = Particle::Status_t::kNuclearLeaving; break; } case genie::kIStHadronInTheNucleus: { state = utility::Is2p2h(chan) ? Particle::Status_t::kPrimaryInitialState : Particle::Status_t::kIntermediate; break; } case genie::kIStPreDecayResonantState: case genie::kIStDISPreFragmHadronicState: case genie::kIStIntermediateState: { state = Particle::Status_t::kIntermediate; break; } case genie::kIStFinalStateNuclearRemnant: case genie::kIStUndefined: case genie::kIStDecayedState: default: { state = Particle::Status_t::kUnknown; } } if (utility::IsNuclearPDG(p.Pdg())) { if (state == Particle::Status_t::kPrimaryInitialState) { state = Particle::Status_t::kPrimaryInitialState; } else if (state == Particle::Status_t::kNuclearLeaving) { state = Particle::Status_t::kPrimaryFinalState; } } return state; } struct TargetSplineBlob { double MolecularWeight; size_t NNucleons; std::string search_term; TGraph TotSpline; }; NEW_NUIS_EXCEPT(Invalid_target_specifier); double GetFileWeight(fhicl::ParameterSet const &xsecinfo, std::set const &spline_list) { double mec_scale = xsecinfo.get("mec_scale", 1); double EMin = xsecinfo.get("EMin", 0); double EMax = xsecinfo.get("EMax", 10); size_t NKnots = xsecinfo.get("NKnots", 100); std::vector TargetSplines; + std::set FoundSplines; int probe = xsecinfo.get("nu_pdg"); for (std::string const &target : xsecinfo.get>("target")) { std::vector splits = nuis::utility::split(target, ";"); if (splits.size() != 1 && splits.size() != 2) { throw Invalid_target_specifier() << "Expected to find 100ZZZAAA0;, e.g. " "1000060120;1, " "but found: \"" << target << "\""; } std::stringstream ss; ss << ".*nu:" << probe << ";tgt:" << splits[0] << ".*"; size_t nuc_pdg = fhicl::string_parsers::str2T(splits[0]); double molwght = splits.size() == 2 ? fhicl::string_parsers::str2T(splits[1]) : 1; TargetSplines.push_back( {molwght, (nuc_pdg % 1000) / 10, ss.str(), TGraph()}); std::cout << "[INFO]:\tGetting splines that match: \"" << ss.str() << "\" with A = " << TargetSplines.back().NNucleons - << " and multiplied by " - << molwght / double(TargetSplines.back().NNucleons) << std::endl; + << " and multiplied by " << molwght << std::endl; } size_t NTotNucleons = 0; std::vector spline_patterns; for (auto const &spb : TargetSplines) { spline_patterns.push_back(spb.search_term); NTotNucleons += spb.MolecularWeight * spb.NNucleons; } double WeightToPerNucleon = 1.0 / double(NTotNucleons); // Read the spline file std::cout << "[INFO]: Reading input file: " << xsecinfo.get("spline_file") << ", this may take a while..." << std::endl; GENIESplineGetter spg(spline_patterns); spg.ParseFile(xsecinfo.get("spline_file").c_str()); std::vector> all_splines = spg.GetTGraphs(); double step = (EMax - EMin) / double(NKnots); // Sum each targets-worth of interaction splines together, weighting for // molecular fraction and possibly fixing the MEC spline oddity for (size_t ts_it = 0; ts_it < TargetSplines.size(); ++ts_it) { TargetSplines[ts_it].TotSpline = TGraph(NKnots); for (size_t p_it = 0; p_it < NKnots; ++p_it) { TargetSplines[ts_it].TotSpline.SetPoint(p_it, EMin + p_it * step, 0); } - double target_weight = (TargetSplines[ts_it].MolecularWeight / - double(TargetSplines[ts_it].NNucleons)); for (size_t c_it = 0; c_it < all_splines[ts_it].size(); ++c_it) { - double weight = target_weight; + double weight = TargetSplines[ts_it].MolecularWeight; if (std::string(all_splines[ts_it][c_it].GetName()).find("MEC") != std::string::npos) { if (mec_scale != 1) { std::cout << "[INFO]: Weighting MEC splines by " << mec_scale << std::endl; weight *= mec_scale; } } + FoundSplines.insert(all_splines[ts_it][c_it].GetName()); for (size_t p_it = 0; p_it < NKnots; ++p_it) { double E, XSec; TargetSplines[ts_it].TotSpline.GetPoint(p_it, E, XSec); - // Copied from GENIE/Convents/Units.h + // Copied from GENIE/Conventions/Units.h static const double gigaelectronvolt = 1.; static const double GeV = gigaelectronvolt; static const double meter = 5.07e+15 / GeV; static const double centimeter = 0.01 * meter; static const double centimeter2 = centimeter * centimeter; XSec += all_splines[ts_it][c_it].Eval(E) * weight / centimeter2; TargetSplines[ts_it].TotSpline.SetPoint(p_it, E, XSec); } } } + for (auto const &sn : spline_list) { + std::cout << "Want: " << sn << std::endl; + if (FoundSplines.count(sn)) { + std::cout << "\tFound!" << std::endl; + } else { + std::cout << "Didn't find." << std::endl; + } + } + // Sum all the correctly weighted per-target splines TGraph MasterSpline(NKnots); + TGraph MasterSpline_pernuc(NKnots); for (size_t p_it = 0; p_it < NKnots; ++p_it) { MasterSpline.SetPoint(p_it, EMin + p_it * step, 0); } for (size_t c_it = 0; c_it < TargetSplines.size(); ++c_it) { for (size_t p_it = 0; p_it < NKnots; ++p_it) { double E, XSec; MasterSpline.GetPoint(p_it, E, XSec); - XSec += TargetSplines[c_it].TotSpline.Eval(E) * WeightToPerNucleon; - + XSec += TargetSplines[c_it].TotSpline.Eval(E); MasterSpline.SetPoint(p_it, E, XSec); } } + for (size_t p_it = 0; p_it < NKnots; ++p_it) { + double E, XSec; + MasterSpline.GetPoint(p_it, E, XSec); + MasterSpline_pernuc.SetPoint(p_it, E, XSec * WeightToPerNucleon); + } + nuis::utility::TFile_ptr outfile(nullptr, [](TFile *) {}); if (xsecinfo.has_key("spline_output_file")) { outfile = nuis::utility::CheckOpenTFile( xsecinfo.get("spline_output_file"), "RECREATE"); for (auto &sp : TargetSplines) { std::string spname = nuis::utility::SanitizeROOTObjectName( sp.search_term.substr(2, sp.search_term.length() - 4)); std::cout << "[INFO]: Dumping spline: " << spname << std::endl; nuis::utility::WriteToTFile(outfile, &sp.TotSpline, spname.c_str()); } nuis::utility::WriteToTFile(outfile, &MasterSpline, "TotalXSec"); + nuis::utility::WriteToTFile(outfile, &MasterSpline_pernuc, + "TotalXSec_pernuc"); } fhicl::ParameterSet fluxps = xsecinfo.get("flux"); // If MonoE if (fluxps.has_key("energy_GeV")) { double monoE = fluxps.get("energy_GeV"); - return MasterSpline.Eval(monoE); + return MasterSpline_pernuc.Eval(monoE); } else { // Use the master spline to get flux*xsec (which the events are thrown // according to) and return the flux-averaged total xsec which is used for // weighting std::unique_ptr Flux = nuis::utility::GetHistogram(fluxps.get("histogram")); if (outfile) { std::unique_ptr Fluxcp = nuis::utility::Clone(Flux); nuis::utility::WriteToTFile(outfile, Fluxcp.get(), "Flux"); } bool per_width = fluxps.get("is_divided_by_bin_width", true); double FluxIntegral = Flux->Integral(per_width ? "width" : ""); for (int bi_it = 0; bi_it < Flux->GetXaxis()->GetNbins(); ++bi_it) { double bc = Flux->GetBinContent(bi_it + 1); double low_e = Flux->GetXaxis()->GetBinLowEdge(bi_it + 1); double up_e = Flux->GetXaxis()->GetBinUpEdge(bi_it + 1); size_t NIntSteps = 100; double step = (up_e - low_e) / double(NIntSteps); double avg_xsec = 0; for (size_t i = 0; i < NIntSteps; ++i) { double e = low_e + (double(i) + 0.5) * step; - avg_xsec += MasterSpline.Eval(e); + if ((e < EMin) || (e > EMax)) { + continue; + } + avg_xsec += MasterSpline_pernuc.Eval(e); } avg_xsec /= double(NIntSteps); Flux->SetBinContent(bi_it + 1, bc * avg_xsec); } if (outfile) { std::unique_ptr Fluxcp = nuis::utility::Clone(Flux); nuis::utility::WriteToTFile(outfile, Fluxcp.get(), "EvRate"); } double EvRateIntegral = Flux->Integral(per_width ? "width" : ""); return EvRateIntegral / FluxIntegral; } } +struct TargetXSecHistBlob { + double MolecularWeight; + size_t NNucleons; + std::string search_term; + TH1 *TotHist; +}; + +std::set GENIEInputHandler::GetFileWeightFromReconstructedSplines( + fhicl::ParameterSet const &xsecinfo, TTree *GENIETree, + genie::NtpMCEventRecord *GenieNtpl) { + + + + double EMin = xsecinfo.get("EMin", 0); + double EMax = xsecinfo.get("EMax", 10); + std::vector TargetXSecHists; + for (std::string const &target : + xsecinfo.get>("target")) { + std::vector splits = nuis::utility::split(target, ";"); + if (splits.size() != 1 && splits.size() != 2) { + throw Invalid_target_specifier() + << "Expected to find 100ZZZAAA0;, e.g. " + "1000060120;1, " + "but found: \"" + << target << "\""; + } + std::stringstream ss; + ss << ".*nu:" << probe << ";tgt:" << splits[0] << ".*"; + + size_t nuc_pdg = fhicl::string_parsers::str2T(splits[0]); + + double molwght = splits.size() == 2 + ? fhicl::string_parsers::str2T(splits[1]) + : 1; + + TargetXSecHists.push_back( + {molwght, (nuc_pdg % 1000) / 10, ss.str(), TGraph()}); + + std::cout << "[INFO]:\tGetting splines that match: \"" << ss.str() + << "\" with A = " << TargetXSecHists.back().NNucleons + << " and multiplied by " << molwght << std::endl; + } + + size_t NTotNucleons = 0; + std::vector spline_patterns; + for (auto const &spb : TargetXSecHists) { + spline_patterns.push_back(spb.search_term); + NTotNucleons += spb.MolecularWeight * spb.NNucleons; + } + double WeightToPerNucleon = 1.0 / double(NTotNucleons); + + std::map splines; + size_t NEvents = GetNEvents(); + size_t ShoutEvery = NEvents / 100; + std::cout << "[INFO]: Rebuilding GENIE splines." << std::endl; + std::cout << "[INFO]: Read " << 0 << "/" << NEvents << " GENIE events." + << std::flush; + for (size_t ev_it = 0; ev_it < NEvents; ++ev_it) { + if (fGenieNtpl) { + fGenieNtpl->Clear(); + } + fInputTree.tree->GetEntry(ev_it); + + if (ShoutEvery && !(ev_it % ShoutEvery)) { + std::cout << "\r[INFO]: Read " << ev_it << "/" << NEvents + << " GENIE events." << std::flush; + } + + genie::GHepRecord *GHep = + static_cast(fGenieNtpl->event); + if (!GHep) { + throw invalid_GENIE_event() << "[ERROR]: GENIE event " << ev_it + << " failed to contain a GHepRecord"; + } + std::string const &spline_name = GHep->Summary()->AsString(); + if (!splinenames.count(spline_name)) { + splinenames.emplace(spline_name, TGraph()); + } + double E_GeV = GHep->Probe()->E(); + double XSec = (genie_record.XSec() / (1E-38 * genie::units::cm2)); + splinenames[spline_name].SetPoint(E_GeV, XSec); + } + std::cout << "[INFO]: Having read GENIE events, found " << splinenames.size() + << " splines: " << std::endl; + for (auto const &sn : splinenames) { + std::cout << "\t" << sn << std::endl; + } + return splinenames; +} + } // namespace genietools } // namespace nuis diff --git a/src/generator/utility/NEUTUtility.cxx b/src/generator/utility/NEUTUtility.cxx index 0b0e04a..d1cf057 100644 --- a/src/generator/utility/NEUTUtility.cxx +++ b/src/generator/utility/NEUTUtility.cxx @@ -1,95 +1,97 @@ #include "generator/utility/NEUTUtility.hxx" #include "generator/input/NEUTInputHandler.hxx" #include "exception/exception.hxx" #include "utility/InteractionChannelUtility.hxx" #include "utility/PDGCodeUtility.hxx" #include "neutpart.h" #include "neutvect.h" using namespace nuis::event; using namespace nuis::utility; namespace nuis { namespace neuttools { NEW_NUIS_EXCEPT(unexpected_NEUT_particle_state); +// #define DEBUG_NEUT_UTILITY + Particle::Status_t GetNeutParticleStatus(NeutPart const &part, Channel_t mode) { #ifdef DEBUG_NEUT_UTILITY std::cout << "[DEBUG]: Mode: " << mode << ", Part: { Status: " << part.fStatus << ", IsAlive: " << part.fIsAlive << ", PDG: " << part.fPID << " }." << std::endl; #endif // Remove Pauli blocked events, probably just single pion events if (part.fStatus == 5) { return Particle::Status_t::kBlocked; // fStatus == -1 means initial state } else if (part.fIsAlive == false && part.fStatus == -1) { return Particle::Status_t::kPrimaryInitialState; // NEUT has a bit of a strange convention for fIsAlive and fStatus // combinations // for NC and neutrino particle isAlive true/false and status 2 means // final state particle // for other particles in NC status 2 means it's an FSI particle // for CC it means it was an FSI particle } else if (part.fStatus == 2) { // NC case is a little strange... The outgoing neutrino might be alive or // not alive. Remaining particles with status 2 are FSI particles that // reinteracted if (IsNC(mode) && IsNeutralLepton(part.fPID)) { return Particle::Status_t::kNuclearLeaving; // The usual CC case } else if (part.fIsAlive == true) { // return Particle::Status_t::kIntermediate; throw unexpected_NEUT_particle_state() << "[ERROR] Found unexpected NEUT particle in neutvect stack: Mode: " << mode << " (IsNC: " << IsNC(mode) << "), Part: { Status: " << part.fStatus << ", IsAlive: " << part.fIsAlive << ", PDG: " << part.fPID << ", IsNeutralLepton: " << IsNeutralLepton(part.fPID) << " }."; } } else if ((part.fIsAlive == true) && (part.fStatus == 2) && IsNeutralLepton(part.fPID)) { return Particle::Status_t::kNuclearLeaving; } else if ((part.fIsAlive == true) && (part.fStatus == 0)) { return Particle::Status_t::kNuclearLeaving; } else if (!part.fIsAlive && ((part.fStatus == 1) || (part.fStatus == 3) || (part.fStatus == 4) || (part.fStatus == 7) || (part.fStatus == 8))) { return Particle::Status_t::kIntermediate; // There's one hyper weird case where fStatus = -3. This apparently // corresponds to a nucleon being ejected via pion FSI when there is "data // available" } else if (!part.fIsAlive && (part.fStatus == -3)) { return Particle::Status_t::kUnknown; // NC neutrino outgoing } else if (!part.fIsAlive && part.fStatus == 0 && IsNeutralLepton(part.fPID)) { return Particle::Status_t::kNuclearLeaving; // Warn if we still find alive particles without classifying them } else if (part.fIsAlive == true) { std::cout << "[WARN]: Undefined NEUT state " << " Alive: " << part.fIsAlive << " Status: " << part.fStatus << " PDG: " << part.fPID << std::endl; return Particle::Status_t::kUnknown; } // Warn if we find dead particles that we haven't classified std::cout << "[WARN]: Undefined NEUT state " << " Alive: " << part.fIsAlive << " Status: " << part.fStatus << " PDG: " << part.fPID << std::endl; return Particle::Status_t::kUnknown; } } // namespace neuttools } // namespace nuis diff --git a/src/samples/nuA/Nuclear/MINERvA/CC0Pi/MINERvA_CC0Pi_CH_xsec_ptpl_nu.cxx b/src/samples/nuA/Nuclear/MINERvA/CC0Pi/MINERvA_CC0Pi_CH_xsec_ptpl_nu.cxx new file mode 100644 index 0000000..e69de29 diff --git a/src/samples/nuA/Nuclear/MINERvA/CCInclusive/CMakeLists.txt b/src/samples/nuA/Nuclear/MINERvA/CCInclusive/CMakeLists.txt new file mode 100644 index 0000000..01227b9 --- /dev/null +++ b/src/samples/nuA/Nuclear/MINERvA/CCInclusive/CMakeLists.txt @@ -0,0 +1,17 @@ +SET(SAMPLES + MINERvA_CCInclusive_CH_xsec_EAvailq3_nu) + +foreach(S ${SAMPLES}) + LIST(APPEND INuADataComparisons ${S}) + LIST(APPEND SAMPLES_src ${S}.cxx) +endforeach() + +add_library(MINERvADataComparisons SHARED ${SAMPLES_src}) +target_link_libraries(MINERvADataComparisons + nuis_event nuis_input nuis_config nuis_persistency) + +install(TARGETS MINERvADataComparisons DESTINATION plugins) + +install(FILES MINERvA.CCInclusive.sample.config.fcl DESTINATION fcl) + +SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) diff --git a/src/samples/nuA/Nuclear/MINERvA/CCInclusive/MINERvA.CCInclusive.sample.config.fcl b/src/samples/nuA/Nuclear/MINERvA/CCInclusive/MINERvA.CCInclusive.sample.config.fcl new file mode 100644 index 0000000..e69de29 diff --git a/src/samples/nuA/Nuclear/MINERvA/CCInclusive/MINERvA_CCInclusive_CH_xsec_EAvailq3_nu.cxx b/src/samples/nuA/Nuclear/MINERvA/CCInclusive/MINERvA_CCInclusive_CH_xsec_EAvailq3_nu.cxx new file mode 100644 index 0000000..e69de29 diff --git a/src/samples/nuA/Nuclear/MINERvA/MINERvA.sample.config.fcl b/src/samples/nuA/Nuclear/MINERvA/MINERvA.sample.config.fcl index 09200c4..afeeb92 100644 --- a/src/samples/nuA/Nuclear/MINERvA/MINERvA.sample.config.fcl +++ b/src/samples/nuA/Nuclear/MINERvA/MINERvA.sample.config.fcl @@ -1 +1,2 @@ #include "MINERvA.CC0Pi.sample.config.fcl" +# include "MINERvA.CCInclusive.sample.config.fcl" diff --git a/src/samples/nuA/Nuclear/MicroBooNE/CCInclusive/CMakeLists.txt b/src/samples/nuA/Nuclear/MicroBooNE/CCInclusive/CMakeLists.txt new file mode 100644 index 0000000..dd82b30 --- /dev/null +++ b/src/samples/nuA/Nuclear/MicroBooNE/CCInclusive/CMakeLists.txt @@ -0,0 +1,17 @@ +SET(SAMPLES +MicroBooNE_CCInclusive_Ar40_xsec_reco_pmu_thetamu_nu) + +foreach(S ${SAMPLES}) + LIST(APPEND INuADataComparisons ${S}) + LIST(APPEND SAMPLES_src ${S}.cxx) +endforeach() + +add_library(MicroBooNEDataComparisons SHARED ${SAMPLES_src}) +target_link_libraries(MicroBooNEDataComparisons + nuis_event nuis_input nuis_config nuis_persistency) + +install(TARGETS MicroBooNEDataComparisons DESTINATION plugins) + +install(FILES MicroBooNE.CCInclusive.sample.config.fcl DESTINATION fcl) + +SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) diff --git a/src/samples/nuA/Nuclear/MicroBooNE/CCInclusive/MicroBooNE.CCInclusive.sample.config.fcl b/src/samples/nuA/Nuclear/MicroBooNE/CCInclusive/MicroBooNE.CCInclusive.sample.config.fcl new file mode 100644 index 0000000..e69de29 diff --git a/src/samples/nuA/Nuclear/MicroBooNE/CCInclusive/MicroBooNE_CCInclusive_Ar40_xsec_reco_pmu_thetamu_nu.cxx b/src/samples/nuA/Nuclear/MicroBooNE/CCInclusive/MicroBooNE_CCInclusive_Ar40_xsec_reco_pmu_thetamu_nu.cxx new file mode 100644 index 0000000..e69de29 diff --git a/src/samples/nuA/Nuclear/MicroBooNE/CMakeLists.txt b/src/samples/nuA/Nuclear/MicroBooNE/CMakeLists.txt new file mode 100644 index 0000000..01ffc5a --- /dev/null +++ b/src/samples/nuA/Nuclear/MicroBooNE/CMakeLists.txt @@ -0,0 +1,7 @@ +add_subdirectory(CCInclusive) + +LIST(APPEND INuADataComparisons_FHiCL MicroBooNE.sample.config.fcl) +install(FILES MicroBooNE.sample.config.fcl DESTINATION fcl) + +SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) +SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE) diff --git a/src/samples/nuA/Nuclear/MicroBooNE/MicroBooNE.sample.config.fcl b/src/samples/nuA/Nuclear/MicroBooNE/MicroBooNE.sample.config.fcl new file mode 100644 index 0000000..675ccaf --- /dev/null +++ b/src/samples/nuA/Nuclear/MicroBooNE/MicroBooNE.sample.config.fcl @@ -0,0 +1 @@ +#include "MiniBooNE.CC0Pi.sample.config.fcl" diff --git a/src/samples/nuA/Nuclear/MiniBooNE/CC0Pi/CMakeLists.txt b/src/samples/nuA/Nuclear/MiniBooNE/CC0Pi/CMakeLists.txt new file mode 100644 index 0000000..d2e924d --- /dev/null +++ b/src/samples/nuA/Nuclear/MiniBooNE/CC0Pi/CMakeLists.txt @@ -0,0 +1,17 @@ +SET(SAMPLES + MiniBooNE_CCQE_CH2_xsec_nu) + +foreach(S ${SAMPLES}) + LIST(APPEND INuADataComparisons ${S}) + LIST(APPEND SAMPLES_src ${S}.cxx) +endforeach() + +add_library(MiniBooNEDataComparisons SHARED ${SAMPLES_src}) +target_link_libraries(MiniBooNEDataComparisons + nuis_event nuis_input nuis_config nuis_persistency) + +install(TARGETS MiniBooNEDataComparisons DESTINATION plugins) + +install(FILES MiniBooNE.CC0Pi.sample.config.fcl DESTINATION fcl) + +SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) diff --git a/src/samples/nuA/Nuclear/MiniBooNE/CC0Pi/MiniBooNE.CC0Pi.sample.config.fcl b/src/samples/nuA/Nuclear/MiniBooNE/CC0Pi/MiniBooNE.CC0Pi.sample.config.fcl new file mode 100644 index 0000000..e69de29 diff --git a/src/samples/nuA/Nuclear/MiniBooNE/CC0Pi/MiniBooNE_CCQE_CH2_xsec_nu.cxx b/src/samples/nuA/Nuclear/MiniBooNE/CC0Pi/MiniBooNE_CCQE_CH2_xsec_nu.cxx new file mode 100644 index 0000000..e69de29 diff --git a/src/samples/nuA/Nuclear/MiniBooNE/CMakeLists.txt b/src/samples/nuA/Nuclear/MiniBooNE/CMakeLists.txt new file mode 100644 index 0000000..f226c0e --- /dev/null +++ b/src/samples/nuA/Nuclear/MiniBooNE/CMakeLists.txt @@ -0,0 +1,7 @@ +add_subdirectory(CC0Pi) + +LIST(APPEND INuADataComparisons_FHiCL MiniBooNE.sample.config.fcl) +install(FILES MiniBooNE.sample.config.fcl DESTINATION fcl) + +SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) +SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE) diff --git a/src/samples/nuA/Nuclear/MiniBooNE/MiniBooNE.sample.config.fcl b/src/samples/nuA/Nuclear/MiniBooNE/MiniBooNE.sample.config.fcl new file mode 100644 index 0000000..675ccaf --- /dev/null +++ b/src/samples/nuA/Nuclear/MiniBooNE/MiniBooNE.sample.config.fcl @@ -0,0 +1 @@ +#include "MiniBooNE.CC0Pi.sample.config.fcl" diff --git a/src/variation/engines/T2KReWeight.cxx b/src/variation/engines/T2KReWeight.cxx index 71eda4e..e9d6d77 100644 --- a/src/variation/engines/T2KReWeight.cxx +++ b/src/variation/engines/T2KReWeight.cxx @@ -1,102 +1,103 @@ #include "variation/engines/T2KReWeight.hxx" #include "fhiclcpp/ParameterSet.h" #include "event/MinimalEvent.hxx" // T2K Engine includes #include "T2KNeutReWeight.h" #include "T2KNIWGReWeight.h" #include "T2KReWeight.h" #include using namespace nuis; using namespace nuis::params; using namespace nuis::event; void T2KReWeightEngine::Initialize(fhicl::ParameterSet const &ps) { fT2KRW = std::make_unique(); TDirectory *dir = gDirectory; fT2KRW->AdoptWghtEngine("NEUT", new t2krew::T2KNeutReWeight); fT2KRW->AdoptWghtEngine("NIWG", new t2krew::T2KNIWGReWeight); dir->cd(); for (fhicl::ParameterSet param_ps : ps.get>("parameters")) { param_ps.put("type", GetName()); std::string const ¶m_name = param_ps.get("name"); T2KSystParam tsp; tsp.t2ksyst = t2krew::T2KSyst::FromString(param_name); if (tsp.t2ksyst == t2krew::kSystNull) { throw invalid_T2K_syst_name() << "[ERROR]: T2KReWeight failed to parse " << std::quoted(param_name) << " as a T2K dial."; } fT2KRW->Systematics().Include(tsp.t2ksyst); tsp.pid = ParameterManager::Get().EnsureParameterRegistered(param_ps); fT2KSysts.push_back(tsp); } Reconfigure(); } void T2KReWeightEngine::Reconfigure() { for (T2KSystParam const &tsp : fT2KSysts) { double val = ParameterManager::Get().GetParameterValue(tsp.pid); fT2KRW->Systematics().SetTwkDial(tsp.t2ksyst, val); + std::cout << "Reconf: " << t2krew::T2KSyst::AsString(tsp.t2ksyst) << " -> " << val << std::endl; } fT2KRW->Reconfigure(); } double T2KReWeightEngine::GetEventWeight(nuis::event::MinimalEvent const &ev) { NeutVect const *nv = static_cast(ev.fGenEvent); if (!nv) { return 1.0; } NeutVect *nv_nc = const_cast(nv); return fT2KRW->CalcWeight(nv_nc); } std::string T2KReWeightEngine::GetName() { return "T2KReWeightEngine"; } std::string T2KReWeightEngine::GetDocumentation() { return ""; } fhicl::ParameterSet T2KReWeightEngine::GetExampleConfiguration() { fhicl::ParameterSet ps; ps.put("name", GetName()); fhicl::ParameterSet dial_maqe; dial_maqe.put("name", "NXSec_MaCCQE"); dial_maqe.put("start", 0); dial_maqe.put("min", -3); dial_maqe.put("max", 3); dial_maqe.put("step", 0.1); fhicl::ParameterSet dial_mares; dial_mares.put("name", "NXSec_MARES"); dial_mares.put("start", 0); dial_mares.put("min", -3); dial_mares.put("max", 3); dial_mares.put("step", 0.1); ps.put>( "parameters", std::vector{{dial_maqe, dial_mares}}); return ps; } GeneratorManager::Generator_id_t T2KReWeightEngine::GetGeneratorId() { return GeneratorManager::Get().EnsureGeneratorRegistered("NEUT"); } DECLARE_PLUGIN(IWeightProvider, T2KReWeightEngine);