diff --git a/cmake/InputHandlerSetup.cmake b/cmake/InputHandlerSetup.cmake
index 88ca096..e69e227 100644
--- a/cmake/InputHandlerSetup.cmake
+++ b/cmake/InputHandlerSetup.cmake
@@ -1,48 +1,56 @@
# 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 .
################################################################################
################################# NuWro ######################################
if(USE_NuWro)
include(${CMAKE_SOURCE_DIR}/cmake/NuWroSetup.cmake)
- cmessage(STATUS "Using NuWro Reweight engine.")
+ cmessage(STATUS "Using NuWro.")
set(USE_NuWro TRUE CACHE BOOL "Whether to enable NuWro support. " FORCE)
endif()
+################################# NEUT #######################################
+if(USE_NEUT)
+ include(${CMAKE_SOURCE_DIR}/cmake/NEUTSetup.cmake)
+ cmessage(STATUS "Using NEUT.")
+ set(USE_NEUT TRUE CACHE BOOL "Whether to enable NEUT support. " FORCE)
+endif()
+
+
if(NEED_ROOTEVEGEN)
cmessage(STATUS "Require ROOT eve generation libraries")
LIST(REVERSE ROOT_LIBS)
LIST(APPEND ROOT_LIBS
Gui
Ged
Geom
TreePlayer
EG
Eve)
LIST(REVERSE ROOT_LIBS)
endif()
if(NEED_ROOTPYTHIA6)
cmessage(STATUS "Require ROOT Pythia6 libraries")
LIST(APPEND ROOT_LIBS
EGPythia6
Pythia6)
endif()
LIST(APPEND EXTRA_LIBS ${ROOT_LIBS})
diff --git a/cmake/NEUTSetup.cmake b/cmake/NEUTSetup.cmake
new file mode 100644
index 0000000..619d09f
--- /dev/null
+++ b/cmake/NEUTSetup.cmake
@@ -0,0 +1,114 @@
+# 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(NOT IS_NEUT_54)
+ 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 EXTRA_CXX_FLAGS -DNEUT_ENABLED )
+
+LIST(APPEND EXTRA_CXX_FLAGS
+ -I${NEUT_ROOT}/include
+ -I${NEUT_ROOT}/src/neutclass
+ -I${NEUT_ROOT}/src/reweight)
+
+LIST(APPEND EXTRA_LINK_DIRS
+ ${NEUT_LIB_DIR}
+ ${CERN}/${CERN_LEVEL}/lib
+ ${NEUT_ROOT}/src/reweight)
+
+if(IS_NEUT_54)
+ LIST(APPEND EXTRA_LIBS
+ NReWeight
+ 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
+ jetset74
+ pdflib804
+ mathlib
+ packlib
+ pawlib)
+else()
+ LIST(APPEND EXTRA_LIBS
+ NReWeight
+ neutcore
+ nuccorrspl
+ nuceff
+ partnuck
+ skmcsvc
+ tauola
+ jetset74
+ pdflib804
+ mathlib
+ packlib
+ pawlib)
+endif()
+
+set(NEUT_ROOT_LIBS)
+
+LIST(APPEND NEUT_ROOT_LIBS
+ ${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 EXTRA_CXX_FLAGS -DNEUT_NUCFSI_ENABLED)
+
+ LIST(APPEND NEUT_ROOT_LIBS
+ ${NEUT_CLASS}/neutnucfsistep.so
+ ${NEUT_CLASS}/neutnucfsivert.so
+ )
+endif()
+
+if(NOT IS_NEUT_54)
+ LIST(APPEND NEUT_ROOT_LIBS
+ ${NEUT_CLASS}/neutrootTreeSingleton.so)
+endif()
+
+LIST(APPEND NEUT_ROOT_LIBS
+ ${NEUT_CLASS}/neutvtx.so
+ ${NEUT_CLASS}/neutfsipart.so
+ ${NEUT_CLASS}/neutpart.so
+ ${NEUT_CLASS}/neutvect.so
+ )
+
+foreach(OBJ ${NEUT_ROOT_LIBS})
+ LIST(APPEND EXTRA_SHAREDOBJS ${OBJ})
+endforeach()
diff --git a/cmake/NuWroSetup.cmake b/cmake/NuWroSetup.cmake
index 9440310..332c32e 100644
--- a/cmake/NuWroSetup.cmake
+++ b/cmake/NuWroSetup.cmake
@@ -1,48 +1,48 @@
# 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 EXTRA_CXX_FLAGS -D__NUWRO_ENABLED__)
+LIST(APPEND EXTRA_CXX_FLAGS -DNUWRO_ENABLED)
LIST(APPEND EXTRA_CXX_FLAGS -I${NUWRO}/src)
if(NOT EXISTS ${NUWRO}/bin/event1.so)
if(EXISTS ${NUWRO}/build/${CMAKE_SYSTEM_NAME}/lib)
if(NUWRO_INC STREQUAL "")
cmessage(FATAL_ERROR "Variable NUWRO_INC is not defined. "
"This must be set to point to an installed NuWro instance.")
endif()
LIST(APPEND EXTRA_LINK_DIRS ${NUWRO}/build/${CMAKE_SYSTEM_NAME}/lib)
LIST(APPEND EXTRA_LIBS event)
else()
cmessage(FATAL_ERROR "Expected to find the NuWro event library in: ${NUWRO}/bin/event1.so, or if using NuWro with reweight support: ${NUWRO}/build/${CMAKE_SYSTEM_NAME}/lib/libevent.a. Is NuWro built?")
endif()
else()
LIST(APPEND EXTRA_SHAREDOBJS ${NUWRO}/bin/event1.so)
endif()
set(NEED_PYTHIA6 TRUE)
set(NEED_ROOTPYTHIA6 TRUE)
diff --git a/cmake/cacheVariables.cmake b/cmake/cacheVariables.cmake
index 64ea088..11f3ae4 100644
--- a/cmake/cacheVariables.cmake
+++ b/cmake/cacheVariables.cmake
@@ -1,89 +1,99 @@
# 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 .
################################################################################
function(CheckAndSetDefaultEnv VARNAME DEFAULT CACHETYPE DOCSTRING ENVNAME)
#cmessage(DEBUG "Trying to assign variable ${VARNAME} into the cache.")
if(NOT DEFINED ${VARNAME})
if(DEFINED ENV{${ENVNAME}} AND NOT $ENV{${ENVNAME}} STREQUAL "")
set(${VARNAME} $ENV{${ENVNAME}} CACHE ${CACHETYPE} ${DOCSTRING})
cmessage(DEBUG " Read ${VARNAME} from ENVVAR ${ENVNAME} as $ENV{${ENVNAME}}.")
else()
set(${VARNAME} ${DEFAULT} CACHE ${CACHETYPE} ${DOCSTRING})
endif()
else()
set(${VARNAME} ${${VARNAME}} CACHE ${CACHETYPE} ${DOCSTRING})
unset(${VARNAME})
endif()
cmessage(CACHE "--Set cache variable: \"${VARNAME}\" to \"${${VARNAME}}\", in cache ${CACHETYPE}.")
endfunction()
function(CheckAndSetDefaultCache VARNAME DEFAULT CACHETYPE DOCSTRING)
# cmessage(DEBUG "Trying to assign variable ${VARNAME} into the cache.")
if(NOT DEFINED ${VARNAME})
set(${VARNAME} ${DEFAULT} CACHE ${CACHETYPE} ${DOCSTRING})
else()
set(${VARNAME} ${${VARNAME}} CACHE ${CACHETYPE} ${DOCSTRING})
unset(${VARNAME})
endif()
cmessage(CACHE "--Set cache variable: \"${VARNAME}\" to \"${${VARNAME}}\", in cache ${CACHETYPE}.")
endfunction()
function(CheckAndSetDefault VARNAME DEFAULT)
# cmessage(DEBUG "Trying to assign variable ${VARNAME}.")
if(NOT DEFINED ${VARNAME})
set(${VARNAME} ${DEFAULT} PARENT_SCOPE)
set(${VARNAME} ${DEFAULT})
endif()
cmessage(CACHE "--Set variable: \"${VARNAME}\" to \"${${VARNAME}}\".")
endfunction()
CheckAndSetDefaultCache(VERBOSE TRUE BOOL "Whether to configure loudly.")
set (CMAKE_SKIP_BUILD_RPATH TRUE)
#Changes default install path to be a subdirectory of the build dir.
#Can set build dir at configure time with -DCMAKE_INSTALL_PREFIX=/install/path
if(CMAKE_INSTALL_PREFIX STREQUAL "" OR CMAKE_INSTALL_PREFIX STREQUAL
"/usr/local")
set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}")
elseif(NOT DEFINED CMAKE_INSTALL_PREFIX)
set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}")
endif()
if(CMAKE_BUILD_TYPE STREQUAL "")
set(CMAKE_BUILD_TYPE DEBUG)
elseif(NOT DEFINED CMAKE_BUILD_TYPE)
set(CMAKE_BUILD_TYPE DEBUG)
endif()
+# Misc
CheckAndSetDefaultCache(EXTRA_SETUP_SCRIPT "" PATH "The path to an extra script to inject into the NUISANCE setup script. <>")
CheckAndSetDefaultCache(USE_MINIMIZER TRUE INTERNAL "Whether we are using the ROOT minimization libraries. ")
CheckAndSetDefaultCache(USE_ROOT6 FALSE INTERNAL "Whether we are using the ROOT 6. ")
+# NuWro
CheckAndSetDefaultCache(USE_NuWro FALSE BOOL "Whether to enable NuWro support. ")
CheckAndSetDefaultEnv(NUWRO "" PATH "Path to NuWro source tree root directory. Overrides environment variable \$NUWRO <>" NUWRO)
CheckAndSetDefaultEnv(NUWRO_INC "" PATH "Path to NuWro installed includes directory, needs to contain \"params_all.h\". Overrides environment variable \$NUWRO_INC <>" NUWRO_INC)
+# NEUT
+CheckAndSetDefaultCache(USE_NEUT FALSE BOOL "Whether to enable NEUT (reweight) support. Requires external libraries. ")
+CheckAndSetDefaultCache(IS_NEUT_54 FALSE BOOL "Whether to enabled NEUT is version 5.4 or greater. ")
+CheckAndSetDefaultEnv(NEUT_ROOT "" PATH "Path to NEUT source tree root directory. Overrides environment variable \$NEUT_ROOT <>" NEUT_ROOT)
+CheckAndSetDefaultEnv(CERN "" PATH "Path to CERNLIB source tree root directory that NEUT was built against. Overrides environment variable \$CERN <>" CERN)
+CheckAndSetDefaultEnv(CERN_LEVEL "" STRING "CERNLIB Library version. Overrides environment variable \$CERN_LEVEL <>" CERN_LEVEL)
+
+# Pythia
CheckAndSetDefaultEnv(PYTHIA6 "" PATH "Path to directory containing libPythia6.so. Overrides environment variable \$PYTHIA6 <>" PYTHIA6)
CheckAndSetDefault(NEED_PYTHIA6 FALSE)
CheckAndSetDefault(NEED_ROOTEVEGEN FALSE)
CheckAndSetDefault(NEED_ROOTPYTHIA6 FALSE)
diff --git a/src/core/FullEvent.cxx b/src/core/FullEvent.cxx
index 83ebcf4..f0fdad0 100644
--- a/src/core/FullEvent.cxx
+++ b/src/core/FullEvent.cxx
@@ -1,13 +1,24 @@
#include "core/FullEvent.hxx"
namespace nuis {
namespace core {
FullEvent::FullEvent() : MinimalEvent() {
- ParticleStack.resize(static_cast(Particle::Status_t::kNParticleStatus));
+
+ for (size_t status_it = 0;
+ status_it < static_cast(Particle::Status_t::kNParticleStatus);
+ ++status_it) {
+ ParticleStack.push_back({static_cast(status_it), {}});
+ }
}
FullEvent::FullEvent(FullEvent &&other)
: MinimalEvent(std::move(other)),
ParticleStack(std::move(other.ParticleStack)) {}
+void FullEvent::ClearParticleStack() {
+ for (auto &status_stack : ParticleStack) {
+ status_stack.particles.clear();
+ }
+}
+
} // namespace core
} // namespace nuis
diff --git a/src/core/FullEvent.hxx b/src/core/FullEvent.hxx
index 5835c10..d351550 100644
--- a/src/core/FullEvent.hxx
+++ b/src/core/FullEvent.hxx
@@ -1,46 +1,48 @@
// 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 .
*******************************************************************************/
#ifndef CORE_FULLEVENT_HXX_SEEN
#define CORE_FULLEVENT_HXX_SEEN
#include "core/MinimalEvent.hxx"
#include "core/Particle.hxx"
namespace nuis {
namespace core {
///\brief The full, internal event format.
class FullEvent : public MinimalEvent {
public:
struct StatusParticles {
Particle::Status_t status;
std::vector particles;
};
FullEvent();
FullEvent(FullEvent const&) = delete;
FullEvent(FullEvent&&);
std::vector ParticleStack;
+
+ void ClearParticleStack();
};
} // namespace core
} // namespace nuis
#endif
diff --git a/src/core/MinimalEvent.cxx b/src/core/MinimalEvent.cxx
index 1ab1fd9..31a5332 100644
--- a/src/core/MinimalEvent.cxx
+++ b/src/core/MinimalEvent.cxx
@@ -1,21 +1,22 @@
#include "core/MinimalEvent.hxx"
namespace nuis {
namespace core {
MinimalEvent::MinimalEvent()
- : mode(0), probe_E(0), probe_pdg(0), XSecWeight(1), RWWeight(1) {
+ : mode(Channel_t::kUndefined), probe_E(0), probe_pdg(0), XSecWeight(1),
+ RWWeight(1) {
#ifdef __NUWRO_ENABLED__
fNuWroEvent = nullptr;
#endif
}
MinimalEvent::MinimalEvent(MinimalEvent &&other)
: mode(other.mode), probe_E(other.probe_E), probe_pdg(other.probe_pdg),
XSecWeight(other.XSecWeight), RWWeight(other.RWWeight) {
#ifdef __NUWRO_ENABLED__
fNuWroEvent = other.fNuWroEvent;
other.fNuWroEvent = nullptr;
#endif
}
} // namespace core
} // namespace nuis
diff --git a/src/core/MinimalEvent.hxx b/src/core/MinimalEvent.hxx
index 4c83e11..92e0dc2 100644
--- a/src/core/MinimalEvent.hxx
+++ b/src/core/MinimalEvent.hxx
@@ -1,64 +1,72 @@
// 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 .
*******************************************************************************/
#ifndef CORE_MINIMALEVENT_HXX_SEEN
#define CORE_MINIMALEVENT_HXX_SEEN
-#ifdef __NUWRO_ENABLED__
+#ifdef NUWRO_ENABLED
#include "event1.h"
#endif
+#ifdef NEUT_ENABLED
+#include "neutpart.h"
+#include "neutvect.h"
+#endif
+
#include "core/types.hxx"
namespace nuis {
namespace core {
///\brief The minimal event information needed to perform reweights.
///
/// Most often, event selections cannot be applied using this reduced format.
class MinimalEvent {
public:
MinimalEvent();
MinimalEvent(MinimalEvent const &) = delete;
MinimalEvent(MinimalEvent &&);
/// True interaction mode
Channel_t mode;
/// True probe energy
double probe_E;
/// True probe particle code
PDG_t probe_pdg;
+ PDG_t target_pdg;
+
/// Event-weight that can be used to scale to a cross-section prediction
double XSecWeight;
/// Event weight incurred from current reweight engine state.
double RWWeight;
-#ifdef __NUWRO_ENABLED__
+#ifdef NUWRO_ENABLED
///\brief Pointer to Nuwro event
///
/// This will usually be tied to a TTree and so we are not responsible for
/// deleting it
event *fNuWroEvent;
#endif
- // Use this to check if
- bool operator!() { return mode == 0; }
+#ifdef NEUT_ENABLED
+ NeutVect *fNeutVect;
+#endif
};
} // namespace core
} // namespace nuis
#endif
diff --git a/src/core/Particle.hxx b/src/core/Particle.hxx
index 0bb047e..88a9528 100644
--- a/src/core/Particle.hxx
+++ b/src/core/Particle.hxx
@@ -1,48 +1,67 @@
// 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 .
*******************************************************************************/
#ifndef CORE_PARTICLE_HXX_SEEN
#define CORE_PARTICLE_HXX_SEEN
#include "core/types.hxx"
#include "TLorentzVector.h"
namespace nuis {
namespace core {
class Particle {
public:
- enum class Status_t {
- kNuclearLeaving = 0,
- kPrimaryInitialState,
- kPrimaryFinalState,
- kIntermediate,
- kUnknown,
- kNParticleStatus
- };
+#define STATUS_LIST \
+ X(kNuclearLeaving, 0) \
+ X(kPrimaryInitialState, 1) \
+ X(kPrimaryFinalState, 2) \
+ X(kIntermediate, 3) \
+ X(kUnknown, 4) \
+ X(kBlocked, 5) \
+ X(kNParticleStatus, 6)
+
+#define X(A, B) A = B,
+ enum class Status_t { STATUS_LIST };
+#undef X
+
Particle();
Particle(Particle const &);
PDG_t pdg;
Status_t status;
TLorentzVector P4;
};
} // namespace core
} // namespace nuis
+
+#define X(A, B) \
+ case nuis::core::Particle::Status_t::A: { \
+ return os << #A; \
+ }
+
+inline std::ostream &operator<<(std::ostream &os,
+ nuis::core::Particle::Status_t te) {
+ switch (te) { STATUS_LIST }
+ return os;
+}
+#undef X
+#undef STATUS_LIST
+
#endif
diff --git a/src/core/types.hxx b/src/core/types.hxx
index b4b2cb6..3891185 100644
--- a/src/core/types.hxx
+++ b/src/core/types.hxx
@@ -1,28 +1,104 @@
// 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 .
*******************************************************************************/
#ifndef CORE_TYPES_HXX_SEEN
#define CORE_TYPES_HXX_SEEN
+#include "exception/exception.hxx"
+
namespace nuis {
namespace core {
-typedef long Channel_t;
+
+#define NUIS_INTERACTION_CHANNEL_LIST \
+ X(kCCQE, 1) \
+ X(kCC2p2h, 2) \
+ X(kCCSPP_PPip, 11) \
+ X(kCCSPP_PPi0, 12) \
+ X(kCCSPP_NPip, 13) \
+ X(kCCCohPi, 16) \
+ X(kCCResGamma, 17) \
+ X(kCCTransitionMPi, 21) \
+ X(kCCResEta0, 22) \
+ X(kCCResK, 23) \
+ X(kCCDIS, 26) \
+ \
+ X(kNCSPP_NPi0, 31) \
+ X(kNCSPP_PPi0, 32) \
+ X(kNCSPP_PPim, 33) \
+ X(kNCSPP_NPip, 34) \
+ X(kNCCohPi, 36) \
+ X(kNCResNGamma, 38) \
+ X(kNCResPGamma, 39) \
+ X(kNCTransitionMPi, 41) \
+ X(kNCResNEta0, 42) \
+ X(kNCResPEta0, 43) \
+ X(kNCResK0, 44) \
+ X(kNCResKp, 45) \
+ X(kNCDIS, 46) \
+ X(kNCELP, 51) \
+ X(kNCELN, 52) \
+ \
+ X(kCCQE_nub, -1) \
+ X(kCC2p2h_nub, -2) \
+ X(kCCSPP_NPim_nub, -11) \
+ X(kCCSPP_NPi0_nub, -12) \
+ X(kCCSPP_PPim_nub, -13) \
+ X(kCCCohPi_nub, -16) \
+ X(kCCResGamma_nub, -17) \
+ X(kCCTransitionMPi_nub, -21) \
+ X(kCCResEta0_nub, -22) \
+ X(kCCResK_nub, -23) \
+ X(kCCDIS_nub, -26) \
+ \
+ X(kNCSPP_NPi0_nub, -31) \
+ X(kNCSPP_PPi0_nub, -32) \
+ X(kNCSPP_PPim_nub, -33) \
+ X(kNCSPP_NPip_nub, -34) \
+ X(kNCCohPi_nub, -36) \
+ X(kNCResNGamma_nub, -38) \
+ X(kNCResPGamma_nub, -39) \
+ X(kNCTransitionMPi_nub, -41) \
+ X(kNCResNEta0_nub, -42) \
+ X(kNCResPEta0_nub, -43) \
+ X(kNCResK0_nub, -44) \
+ X(kNCResKp_nub, -45) \
+ X(kNCDIS_nub, -46) \
+ X(kNCELP_nub, -51) \
+ X(kNCELN_nub, -52)\
+ \
+ X(kUndefined, 0)
+
+#define X(A, B) A = B,
+enum class Channel_t { NUIS_INTERACTION_CHANNEL_LIST };
+#undef X
+
typedef long PDG_t;
} // namespace core
} // namespace nuis
+
+#define X(A, B) \
+ case nuis::core::Channel_t::A: { \
+ return os << #A; \
+ }
+
+inline std::ostream &operator<<(std::ostream &os, nuis::core::Channel_t te) {
+ switch (te) { NUIS_INTERACTION_CHANNEL_LIST }
+ return os;
+}
+#undef X
#endif
diff --git a/src/generator/CMakeLists.txt b/src/generator/CMakeLists.txt
index 45c16fe..0b35aad 100644
--- a/src/generator/CMakeLists.txt
+++ b/src/generator/CMakeLists.txt
@@ -1,6 +1,12 @@
if(USE_NuWro)
- add_library(InputHandlers SHARED NuWroInputHandler.cxx)
- target_link_libraries(InputHandlers nuis_core nuis_config)
-
- install(TARGETS InputHandlers DESTINATION plugins)
+ LIST(APPEND INPUT_HANDLERS_IMPL NuWroInputHandler.cxx)
endif(USE_NuWro)
+
+if(USE_NEUT)
+ LIST(APPEND INPUT_HANDLERS_IMPL NEUTInputHandler.cxx)
+endif(USE_NEUT)
+
+add_library(InputHandlers SHARED ${INPUT_HANDLERS_IMPL})
+target_link_libraries(InputHandlers nuis_core nuis_config)
+
+install(TARGETS InputHandlers DESTINATION plugins)
diff --git a/src/generator/NEUTInputHandler.cxx b/src/generator/NEUTInputHandler.cxx
new file mode 100644
index 0000000..70bf5ce
--- /dev/null
+++ b/src/generator/NEUTInputHandler.cxx
@@ -0,0 +1,85 @@
+#include "NEUTInputHandler.hxx"
+
+#include "core/FullEvent.hxx"
+
+#include "utility/ROOTUtility.hxx"
+
+#include "generator/NEUTUtility.hxx"
+
+#include "fhiclcpp/ParameterSet.h"
+
+using namespace nuis::core;
+using namespace nuis::utility;
+
+NEUTInputHandler::NEUTInputHandler() : fInputTree(nullptr) {}
+NEUTInputHandler::NEUTInputHandler(NEUTInputHandler &&other)
+ : fInputTree(std::move(other.fInputTree)),
+ fReaderEvent(std::move(other.fReaderEvent)) {}
+
+void NEUTInputHandler::Initialize(fhicl::ParameterSet const &ps) {
+
+ fInputTree = CheckGetTTree(ps.get("file"), "neuttree");
+
+ fReaderEvent.fNeutVect = nullptr;
+ fInputTree->tree->SetBranchAddress("vectorbranch", &fReaderEvent.fNeutVect);
+}
+
+MinimalEvent const &NEUTInputHandler::GetMinimalEvent(ev_index_t idx) const {
+ if (idx >= GetNEvents()) {
+ throw IInputHandler::invalid_entry()
+ << "[ERROR]: Attempted to get entry " << idx
+ << " from an InputHandler with only " << GetNEvents();
+ }
+ fInputTree->tree->GetEntry(idx);
+
+ fReaderEvent.mode = IntToChannel(fReaderEvent.fNeutVect->Mode);
+
+ size_t NPart = fReaderEvent.fNeutVect->Npart();
+ for (size_t part_it = 0; part_it < NPart; part_it++) {
+ NeutPart const &part = (*fReaderEvent.fNeutVect->PartInfo(part_it));
+ if ((part.fIsAlive == false) && (part.fStatus == -1) &&
+ IsNeutralLepton(part.fPID)) {
+ fReaderEvent.probe_E = part.fP.T();
+ fReaderEvent.probe_pdg = part.fPID;
+ break;
+ }
+ }
+
+ return fReaderEvent;
+}
+
+FullEvent const &NEUTInputHandler::GetFullEvent(ev_index_t idx) const {
+ (void)GetMinimalEvent(idx);
+
+ fReaderEvent.ClearParticleStack();
+
+ if (fReaderEvent.fNeutVect->Ibound) {
+ fReaderEvent.target_pdg = MakeNuclearPDG(fReaderEvent.fNeutVect->TargetA,
+ fReaderEvent.fNeutVect->TargetZ);
+ } else {
+ fReaderEvent.target_pdg = MakeNuclearPDG(1, 1);
+ }
+
+ size_t NPart = fReaderEvent.fNeutVect->Npart();
+ for (size_t part_it = 0; part_it < NPart; part_it++) {
+ NeutPart const &part = (*fReaderEvent.fNeutVect->PartInfo(part_it));
+
+ nuis::core::Particle nuis_part;
+ nuis_part.pdg = part.fPID;
+ nuis_part.status = GetNeutParticleStatus(part, fReaderEvent.mode);
+ nuis_part.P4 = part.fP;
+
+ size_t state_int = static_cast(nuis_part.status);
+
+ fReaderEvent.ParticleStack[state_int].particles.push_back(
+ std::move(nuis_part));
+ }
+
+ return fReaderEvent;
+}
+
+size_t NEUTInputHandler::GetNEvents() const {
+ return fInputTree->tree->GetEntries();
+}
+
+DECLARE_PLUGIN(IInputHandler, NEUTInputHandler);
diff --git a/src/samples/ISample.hxx b/src/generator/NEUTInputHandler.hxx
similarity index 61%
copy from src/samples/ISample.hxx
copy to src/generator/NEUTInputHandler.hxx
index 2333cbd..9426658 100644
--- a/src/samples/ISample.hxx
+++ b/src/generator/NEUTInputHandler.hxx
@@ -1,53 +1,56 @@
// 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 .
*******************************************************************************/
-#ifndef SAMPLES_ISAMPLE_HXX_SEEN
-#define SAMPLES_ISAMPLE_HXX_SEEN
+#ifndef GENERATOR_NEUTINPUTHANDLER_HXX_SEEN
+#define GENERATOR_NEUTINPUTHANDLER_HXX_SEEN
-#include "plugins/traits.hxx"
+#include "core/IInputHandler.hxx"
+#include "core/FullEvent.hxx"
-#include "exception/exception.hxx"
+#include
namespace fhicl {
class ParameterSet;
}
namespace nuis {
namespace core {
-class FullEvent;
class MinimalEvent;
} // namespace core
+namespace utility {
+class TreeFile;
+}
} // namespace nuis
-class ISample {
-public:
- NEW_NUIS_EXCEPT(uninitialized_ISample);
-
- virtual void Initialize(fhicl::ParameterSet const &) = 0;
-
- virtual void ProcessSample() = 0;
+class NEUTInputHandler : public IInputHandler {
+ mutable std::unique_ptr fInputTree;
+ mutable nuis::core::FullEvent fReaderEvent;
- virtual void Write() = 0;
-
- virtual ~ISample(){}
+public:
+ NEUTInputHandler();
+ NEUTInputHandler(NEUTInputHandler const &) = delete;
+ NEUTInputHandler(NEUTInputHandler &&);
+
+ void Initialize(fhicl::ParameterSet const &);
+ nuis::core::MinimalEvent const &GetMinimalEvent(ev_index_t idx) const;
+ nuis::core::FullEvent const &GetFullEvent(ev_index_t idx) const;
+ size_t GetNEvents() const;
};
-DECLARE_PLUGIN_INTERFACE(ISample);
-
#endif
diff --git a/src/generator/NEUTUtility.hxx b/src/generator/NEUTUtility.hxx
new file mode 100644
index 0000000..7d35c45
--- /dev/null
+++ b/src/generator/NEUTUtility.hxx
@@ -0,0 +1,81 @@
+#ifndef GENERATOR_NEUTUTILITY_HXX_SEEN
+#define GENERATOR_NEUTUTILITY_HXX_SEEN
+
+#include "core/Particle.hxx"
+
+#include "utility/ChannelUtility.hxx"
+#include "utility/PDGCodeUtility.hxx"
+
+#include "exception/exception.hxx"
+
+#include "neutpart.h"
+#include "neutvect.h"
+
+NEW_NUIS_EXCEPT(unexpected_NEUT_particle_state);
+
+inline nuis::core::Particle::Status_t
+GetNeutParticleStatus(NeutPart const &part, nuis::core::Channel_t mode) {
+ // Remove Pauli blocked events, probably just single pion events
+ if (part.fStatus == 5) {
+ return nuis::core::Particle::Status_t::kBlocked;
+
+ // fStatus == -1 means initial state
+ } else if (part.fIsAlive == false && part.fStatus == -1) {
+ return nuis::core::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 (nuis::utility::IsNC(mode) &&
+ nuis::utility::IsNeutralLepton(part.fPID)) {
+ return nuis::core::Particle::Status_t::kNuclearLeaving;
+ // The usual CC case
+ } else if (part.fIsAlive == true) {
+ return nuis::core::Particle::Status_t::kIntermediate;
+ }
+
+ } else if ((part.fIsAlive == true) && (part.fStatus == 2) &&
+ nuis::utility::IsNeutralLepton(part.fPID)) {
+ return nuis::core::Particle::Status_t::kNuclearLeaving;
+
+ } else if ((part.fIsAlive == true) && (part.fStatus == 0)) {
+ return nuis::core::Particle::Status_t::kNuclearLeaving;
+
+ } else if (!part.fIsAlive && ((part.fStatus == 1) || (part.fStatus == 3) ||
+ (part.fStatus == 4) || (part.fStatus == 7) ||
+ (part.fStatus == 8))) {
+ return nuis::core::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 nuis::core::Particle::Status_t::kUnknown;
+ // NC neutrino outgoing
+ } else if (!part.fIsAlive && part.fStatus == 0 &&
+ (abs(part.fPID) == 16 || abs(part.fPID) == 14 ||
+ abs(part.fPID) == 12)) {
+ return nuis::core::Particle::Status_t::kNuclearLeaving;
+
+ // Warn if we still find alive particles without classifying them
+ } else if (part.fIsAlive == true) {
+ throw unexpected_NEUT_particle_state()
+ << "[ERROR]: Undefined NEUT state "
+ << " Alive: " << part.fIsAlive << " Status: " << part.fStatus
+ << " PDG: " << part.fPID;
+ }
+ // Warn if we find dead particles that we haven't classified
+ throw unexpected_NEUT_particle_state()
+ << "[ERROR]: Undefined NEUT state "
+ << " Alive: " << part.fIsAlive << " Status: " << part.fStatus
+ << " PDG: " << part.fPID;
+}
+
+#endif
diff --git a/src/samples/ISample.hxx b/src/samples/ISample.hxx
index 2333cbd..23df06a 100644
--- a/src/samples/ISample.hxx
+++ b/src/samples/ISample.hxx
@@ -1,53 +1,65 @@
// 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 .
*******************************************************************************/
#ifndef SAMPLES_ISAMPLE_HXX_SEEN
#define SAMPLES_ISAMPLE_HXX_SEEN
#include "plugins/traits.hxx"
#include "exception/exception.hxx"
+#include
+#include
+
namespace fhicl {
class ParameterSet;
}
namespace nuis {
namespace core {
class FullEvent;
class MinimalEvent;
} // namespace core
} // namespace nuis
class ISample {
public:
NEW_NUIS_EXCEPT(uninitialized_ISample);
+ NEW_NUIS_EXCEPT(unimplemented_ISample_optional_method);
virtual void Initialize(fhicl::ParameterSet const &) = 0;
- virtual void ProcessSample() = 0;
+ virtual void ProcessEvent(nuis::core::FullEvent const &) {
+ throw unimplemented_ISample_optional_method()
+ << "[ERROR]: ISample " << Name() << " does not implement ProcessEvent.";
+ }
+
+ virtual void
+ ProcessSample(size_t nmax = std::numeric_limits::max()) = 0;
virtual void Write() = 0;
- virtual ~ISample(){}
+ virtual std::string Name() = 0;
+
+ virtual ~ISample() {}
};
DECLARE_PLUGIN_INTERFACE(ISample);
#endif
diff --git a/src/samples/MCTools/CMakeLists.txt b/src/samples/MCTools/CMakeLists.txt
index 0f760ae..cff9823 100644
--- a/src/samples/MCTools/CMakeLists.txt
+++ b/src/samples/MCTools/CMakeLists.txt
@@ -1,8 +1,10 @@
if(USE_NuWro)
LIST(APPEND MC_TOOL_IMPL NuisToNuWro.cxx)
endif(USE_NuWro)
+LIST(APPEND MC_TOOL_IMPL VerboseEventSummary.cxx)
+
add_library(MCTools SHARED ${MC_TOOL_IMPL})
target_link_libraries(MCTools nuis_core nuis_config)
install(TARGETS MCTools DESTINATION plugins)
diff --git a/src/samples/MCTools/NuisToNuWro.cxx b/src/samples/MCTools/NuisToNuWro.cxx
index 68ae343..eef38c0 100644
--- a/src/samples/MCTools/NuisToNuWro.cxx
+++ b/src/samples/MCTools/NuisToNuWro.cxx
@@ -1,40 +1,45 @@
#include "samples/ISample.hxx"
-#include "core/InputManager.hxx"
#include "core/FullEvent.hxx"
+#include "core/InputManager.hxx"
-#include
#include
+#include
using namespace nuis::core;
class NuisToNuWro : public ISample {
public:
InputManager::Input_id_t fIH_id;
NuisToNuWro()
: fIH_id(std::numeric_limits::max()) {}
void Initialize(fhicl::ParameterSet const &ps) {
fIH_id = InputManager::Get().EnsureInputLoaded(ps);
}
- void Process(FullEvent const &ps) {
+ void ProcessEvent(FullEvent const &ps) {
std::cout << ps.fNuWroEvent->dyn << std::endl;
}
- void ProcessSample() {
+ void ProcessSample(size_t nmax) {
if (fIH_id == std::numeric_limits::max()) {
throw uninitialized_ISample();
}
IInputHandler const &IH = InputManager::Get().GetInputHandler(fIH_id);
+ size_t n = 0;
for (auto const &fe : IH) {
- Process(fe);
+ if (++n > nmax) {
+ break;
+ }
+ ProcessEvent(fe);
}
}
void Write() {}
+ std::string Name() { return "NuisToNuWro"; }
};
DECLARE_PLUGIN(ISample, NuisToNuWro);
diff --git a/src/samples/MCTools/VerboseEventSummary.cxx b/src/samples/MCTools/VerboseEventSummary.cxx
new file mode 100644
index 0000000..73470ce
--- /dev/null
+++ b/src/samples/MCTools/VerboseEventSummary.cxx
@@ -0,0 +1,58 @@
+#include "samples/ISample.hxx"
+
+#include "core/FullEvent.hxx"
+#include "core/InputManager.hxx"
+
+#include
+#include
+
+using namespace nuis::core;
+
+class VerboseEventSummary : public ISample {
+public:
+ InputManager::Input_id_t fIH_id;
+
+ VerboseEventSummary()
+ : fIH_id(std::numeric_limits::max()) {}
+
+ void Initialize(fhicl::ParameterSet const &ps) {
+ fIH_id = InputManager::Get().EnsureInputLoaded(ps);
+ }
+
+ void ProcessEvent(FullEvent const &ps) {
+ std::cout << "Event: Interaction mode = " << ps.mode
+ << ", probe: { PDG: " << ps.probe_pdg
+ << ", Energy: " << ps.probe_E << "}." << std::endl;
+ for (auto &status_stack : ps.ParticleStack) {
+ std::cout << "\t[" << status_stack.status << "]" << std::endl;
+
+ for (nuis::core::Particle const &part : status_stack.particles) {
+ std::cout << "\t\t{ PDG: " << part.pdg << ", P3: [ " << part.P4[0]
+ << ", " << part.P4[1] << ", " << part.P4[2]
+ << "], E: " << part.P4[3] << ", M: " << part.P4.M()
+ << std::endl
+ << std::endl;
+ }
+ }
+ }
+
+ void ProcessSample(size_t nmax) {
+ if (fIH_id == std::numeric_limits::max()) {
+ throw uninitialized_ISample();
+ }
+
+ IInputHandler const &IH = InputManager::Get().GetInputHandler(fIH_id);
+ size_t n = 0;
+ for (auto const &fe : IH) {
+ if (++n > nmax) {
+ break;
+ }
+ ProcessEvent(fe);
+ }
+ }
+
+ void Write() {}
+ std::string Name() { return "VerboseEventSummary"; }
+};
+
+DECLARE_PLUGIN(ISample, VerboseEventSummary);
diff --git a/src/utility/ChannelUtility.hxx b/src/utility/ChannelUtility.hxx
new file mode 100644
index 0000000..9aaa523
--- /dev/null
+++ b/src/utility/ChannelUtility.hxx
@@ -0,0 +1,58 @@
+#ifndef UTILITY_CHANNELUTILITY_HXX_SEEN
+#define UTILITY_CHANNELUTILITY_HXX_SEEN
+
+#include "core/types.hxx"
+
+#include "exception/exception.hxx"
+
+#include
+
+namespace nuis {
+
+namespace utility {
+NEW_NUIS_EXCEPT(invalid_channel_conversion);
+
+#define X(A, B) \
+ case B: { \
+ return nuis::core::Channel_t::A; \
+ }
+
+inline core::Channel_t IntToChannel(int mode) {
+ switch (mode) {
+ NUIS_INTERACTION_CHANNEL_LIST
+ default: {
+ throw invalid_channel_conversion()
+ << "[ERROR]: Failed to parse mode integer " << mode
+ << " as a nuis::core::Channel_t.";
+ }
+ }
+}
+
+#undef X
+
+#define X(A, B) \
+ case core::Channel_t::A: { \
+ return B; \
+ }
+
+inline int ChannelToInt(core::Channel_t mode) {
+ switch (mode) {
+ NUIS_INTERACTION_CHANNEL_LIST
+ default: {
+ throw invalid_channel_conversion()
+ << "[ERROR]: Attempting to convert "
+ "undefined nuis::core::Channel_t to an "
+ "integer.";
+ }
+ }
+}
+
+inline bool IsNC(core::Channel_t mode) { return abs(ChannelToInt(mode) > 30); }
+inline bool IsCC(core::Channel_t mode) { return !IsNC(mode); }
+inline bool IsNu(core::Channel_t mode) { return ChannelToInt(mode) > 0; }
+inline bool IsNub(core::Channel_t mode) { return !IsNu(mode); }
+
+} // namespace utility
+} // namespace nuis
+
+#endif
diff --git a/src/utility/PDGCodeUtility.hxx b/src/utility/PDGCodeUtility.hxx
index 62cbd5e..f62044c 100644
--- a/src/utility/PDGCodeUtility.hxx
+++ b/src/utility/PDGCodeUtility.hxx
@@ -1,48 +1,116 @@
// 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 .
*******************************************************************************/
#ifndef UTILITY_PDGCODEUTILITY_HXX_SEEN
#define UTILITY_PDGCODEUTILITY_HXX_SEEN
+#include "core/types.hxx"
+
namespace nuis {
namespace utility {
namespace pdgcodes {
-static std::vector const ChargedLeptons{11, 13, 15, -11, -13, -15};
-static std::vector const ChargedLeptons_matter{11, 13, 15};
-static std::vector const ChargedLeptons_antimatter{-11, -13, -15};
-
-static std::vector const NeutralLeptons{12, 14, 16, -12, -14, -16};
-static std::vector const NeutralLeptons_matter{12, 14, 16};
-static std::vector const NeutralLeptons_antimatter{-12, -14, -16};
-
-static std::vector const ChargedPions{211, -211};
-static std::vector const NeutralPions{111};
-static std::vector const Pions{211, -211, 111};
-static std::vector const Protons{2212};
-static std::vector const Neutron{2112};
-static std::vector const Nucleons{2212, 2112};
-
-static std::vector const CommonParticles{11, 13, 15, -11, -13, -15,
- 12, 14, 16, -12, -14, -16,
- 211, -211, 111, 2212, 2112};
+
+enum MatterType { kMatter = 1, kMatterAntimatter = 0, kAntimatter = -1 };
+
+static std::vector const ChargedLeptons{11, 13, 15, -11, -13, -15};
+static std::vector const ChargedLeptons_matter{11, 13, 15};
+static std::vector const ChargedLeptons_antimatter{-11, -13, -15};
+
+static std::vector const NeutralLeptons{12, 14, 16, -12, -14, -16};
+static std::vector const NeutralLeptons_matter{12, 14, 16};
+static std::vector const NeutralLeptons_antimatter{-12, -14, -16};
+
+static std::vector const ChargedPions{211, -211};
+static std::vector const NeutralPions{111};
+static std::vector const Pions{211, -211, 111};
+static std::vector const Protons{2212, -2122};
+static std::vector const Proton_matter{2212};
+static std::vector const Proton_antimatter{-2212};
+static std::vector const Neutron{2112};
+static std::vector const Nucleons{2212, 2112, -2212};
+static std::vector const Nucleons_matter{2212, 2112};
+static std::vector const Nucleons_antimatter{-2212};
+
+static std::vector const CommonParticles{
+ 11, 13, 15, -11, -13, -15, 12, 14, 16,
+ -12, -14, -16, 211, -211, 111, 2212, 2112};
} // namespace pdgcodes
+
+inline bool
+IsInPDGList(core::PDG_t pdg, std::vector const &MatterList,
+ std::vector const &AntiMatterList,
+ pdgcodes::MatterType type = pdgcodes::kMatterAntimatter) {
+ switch (type) {
+ case pdgcodes::kMatter: {
+ return std::count(MatterList.begin(), MatterList.end(), pdg);
+ }
+ case pdgcodes::kMatterAntimatter: {
+ return std::count(MatterList.begin(), MatterList.end(), pdg) ||
+ std::count(AntiMatterList.begin(), AntiMatterList.end(), pdg);
+ }
+ case pdgcodes::kAntimatter: {
+ return std::count(AntiMatterList.begin(), AntiMatterList.end(), pdg);
+ }
+ }
+}
+inline bool
+IsNeutralLepton(core::PDG_t pdg,
+ pdgcodes::MatterType type = pdgcodes::kMatterAntimatter) {
+ return IsInPDGList(pdg, pdgcodes::NeutralLeptons_matter,
+ pdgcodes::NeutralLeptons_antimatter, type);
+}
+inline bool
+IsChargedLepton(core::PDG_t pdg,
+ pdgcodes::MatterType type = pdgcodes::kMatterAntimatter) {
+ return IsInPDGList(pdg, pdgcodes::ChargedLeptons_matter,
+ pdgcodes::ChargedLeptons_antimatter, type);
+}
+inline bool IsProton(core::PDG_t pdg,
+ pdgcodes::MatterType type = pdgcodes::kMatterAntimatter) {
+ return IsInPDGList(pdg, pdgcodes::Proton_matter, pdgcodes::Proton_antimatter,
+ type);
+}
+inline bool IsNeutron(core::PDG_t pdg) { return pdg == pdgcodes::Neutron[0]; }
+inline bool IsChargedPion(core::PDG_t pdg) {
+ return std::count(pdgcodes::ChargedPions.begin(),
+ pdgcodes::ChargedPions.end(), pdg);
+}
+inline bool IsNeutralPion(core::PDG_t pdg) {
+ return std::count(pdgcodes::NeutralPions.begin(),
+ pdgcodes::NeutralPions.end(), pdg);
+}
+inline bool IsPion(core::PDG_t pdg) {
+ return std::count(pdgcodes::Pions.begin(), pdgcodes::Pions.end(), pdg);
+}
+inline bool IsOther(core::PDG_t pdg) {
+ return !std::count(pdgcodes::CommonParticles.begin(),
+ pdgcodes::CommonParticles.end(), pdg);
+}
+
+core::PDG_t MakeNuclearPDG(size_t A, size_t Z) {
+ return 1000 * Z + 10 * A + 1000000000;
+}
+
+size_t GetA(core::PDG_t pdg) { return ((pdg / 10) % 1000); }
+size_t GetZ(core::PDG_t pdg) { return ((pdg / 1000) % 1000); }
+
} // namespace utility
} // namespace nuis
#endif