No OneTemporary
Actions

Size

45 KB

Subscribers

None

View Options

	diff --git a/src/event/Particle.hxx b/src/event/Particle.hxx
	index 04fcfe8..f71d2e9 100644
	--- a/src/event/Particle.hxx
	+++ b/src/event/Particle.hxx
	@@ -1,78 +1,89 @@
	// 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 <http://www.gnu.org/licenses/>.
	*******************************************************************************/

	-#ifndef EVENT_PARTICLE_HXX_SEEN
	-#define EVENT_PARTICLE_HXX_SEEN
	+#pragma once

	#include "event/types.hxx"

	#include "TLorentzVector.h"

	namespace nuis {
	namespace event {
	class Particle {
	public:
	NEW_NUIS_EXCEPT(invalid_particle);

	#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;
	TLorentzVector P4;

	bool operator!() const { return (pdg == std::numeric_limits<PDG_t>::max()); }

	double E() const { return P4.E(); }
	double KE() const { return P4.E() - P4.M(); }
	double P() const { return P4.Vect().Mag(); }
	TVector3 P3() const { return P4.Vect(); }
	+ TVector3 Dir() const { return P3().Unit(); }
	double M() const { return P4.M(); }
	double Theta() const { return P4.Vect().Theta(); }
	double CosTheta() const { return P4.Vect().CosTheta(); }
	};
	} // namespace event
	} // namespace nuis

	#define X(A, B) \
	case nuis::event::Particle::Status_t::A: { \
	return os << #A; \
	}

	inline std::ostream &operator<<(std::ostream &os,
	nuis::event::Particle::Status_t te) {
	switch (te) { STATUS_LIST }
	return os;
	}
	#undef X
	#undef STATUS_LIST

	-#endif
	+inline bool operator==(nuis::event::Particle const &l,
	+ nuis::event::Particle const &r) {
	+ if (l.pdg != r.pdg) {
	+ return false;
	+ }
	+ for (size_t i = 0; i < 4; ++i) {
	+ if (l.P4[i] != r.P4[i]) {
	+ return false;
	+ }
	+ }
	+ return true;
	+}
	diff --git a/src/samples/MCTools/CMakeLists.txt b/src/samples/MCTools/CMakeLists.txt
	index 236fdfa..6d1455b 100644
	--- a/src/samples/MCTools/CMakeLists.txt
	+++ b/src/samples/MCTools/CMakeLists.txt
	@@ -1,16 +1,16 @@
	if(USE_NUWRO)
	LIST(APPEND MC_TOOL_IMPL NuisToNuWro.cxx)
	endif(USE_NUWRO)

	-LIST(APPEND MC_TOOL_IMPL VerboseEventSummary.cxx)
	+LIST(APPEND MC_TOOL_IMPL VerboseEventSummary.cxx NuisToCAFTruth)

	add_library(MCTools SHARED ${MC_TOOL_IMPL})
	target_link_libraries(MCTools ${SAMPLES_MCTOOLS_LIBS})

	if(USE_NUWRO)
	target_compile_options(MCTools PUBLIC ${NUWRO_CXX_FLAGS})
	include_directories(${NUWRO_INCLUDE_DIRS})
	target_link_libraries(MCTools -Wl,--no-as-needed ${NUWRO_IMPORTED_TARGETS})
	endif()

	install(TARGETS MCTools DESTINATION plugins)
	diff --git a/src/samples/MCTools/NuisToCAFTruth.cxx b/src/samples/MCTools/NuisToCAFTruth.cxx
	new file mode 100644
	index 0000000..5ba6935
	--- /dev/null
	+++ b/src/samples/MCTools/NuisToCAFTruth.cxx
	@@ -0,0 +1,253 @@
	+// 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 <http://www.gnu.org/licenses/>.
	+ *******************************************************************************/
	+
	+//********************************************************************
	+
	+#include "samples/SimpleMCStudy.hxx"
	+
	+#include "utility/EventTopologyUtility.hxx"
	+#include "utility/FullEventUtility.hxx"
	+#include "utility/KinematicUtility.hxx"
	+#include "utility/PDGCodeUtility.hxx"
	+#include "utility/ROOTUtility.hxx"
	+
	+#include "persistency/ROOTOutput.hxx"
	+
	+using namespace nuis::event;
	+using namespace nuis::utility;
	+
	+class NuisToCAFTruth : public SimpleMCStudy {
	+
	+ TreeFile outputTree;
	+
	+ // int mode;
	+
	+ int isFD;
	+ int isFHC;
	+ int isCC;
	+
	+ int nuPDG;
	+ int nuPDGunosc;
	+ int LepPDG;
	+
	+ double NuMomX;
	+ double NuMomY;
	+ double NuMomZ;
	+ double Ev;
	+ double LepMomX;
	+ double LepMomY;
	+ double LepMomZ;
	+ double LepE;
	+ double LepNuAngle;
	+
	+ double Q2;
	+ double W;
	+ double X;
	+ double Y;
	+
	+ int nP;
	+ int nN;
	+ int nipip;
	+ int nipim;
	+ int nipi0;
	+ int nikp;
	+ int nikm;
	+ int nik0;
	+ int niem;
	+ int niother;
	+ int nNucleus;
	+ int nUNKNOWN;
	+
	+ double eP;
	+ double eN;
	+ double ePip;
	+ double ePim;
	+ double ePi0;
	+ double eOther;
	+
	+ double xsecWeight;
	+
	+ bool isOscSwap;
	+
	+public:
	+ NuisToCAFTruth() { ReadGlobalConfigDefaults(); }
	+
	+ void Initialize(fhicl::ParameterSet const &instance_sample_configuration) {
	+
	+ outputTree = AddNewTTreeToFile(nuis::persistency::GetOutputFile().get(),
	+ "cafTruthTree");
	+
	+ isFD = instance_sample_configuration.get<bool>("is_fd", false);
	+ isFHC = instance_sample_configuration.get<bool>("is_numode", true);
	+ isOscSwap = instance_sample_configuration.get<bool>("is_nueswap", false);
	+
	+ // Perform any per-sample configuration in the base class
	+ SimpleMCStudy::Initialize(instance_sample_configuration);
	+
	+ // outputTree->Branch("mode", &mode, "mode/I");
	+ outputTree->Branch("isFD", &isFD, "isFD/I");
	+ outputTree->Branch("isFHC", &isFHC, "isFHC/I");
	+ outputTree->Branch("isCC", &isCC, "isCC/I");
	+ outputTree->Branch("nuPDG", &nuPDG, "nuPDG/I");
	+ outputTree->Branch("nuPDGunosc", &nuPDGunosc, "nuPDGunosc/I");
	+ outputTree->Branch("LepPDG", &LepPDG, "LepPDG/I");
	+ outputTree->Branch("NuMomX", &NuMomX, "NuMomX/D");
	+ outputTree->Branch("NuMomY", &NuMomY, "NuMomY/D");
	+ outputTree->Branch("NuMomZ", &NuMomZ, "NuMomZ/D");
	+ outputTree->Branch("Ev", &Ev, "Ev/D");
	+ outputTree->Branch("LepMomX", &LepMomX, "LepMomX/D");
	+ outputTree->Branch("LepMomY", &LepMomY, "LepMomY/D");
	+ outputTree->Branch("LepMomZ", &LepMomZ, "LepMomZ/D");
	+ outputTree->Branch("LepE", &LepE, "LepE/D");
	+ outputTree->Branch("LepNuAngle", &LepNuAngle, "LepNuAngle/D");
	+ outputTree->Branch("Q2", &Q2, "Q2/D");
	+ outputTree->Branch("W", &W, "W/D");
	+ outputTree->Branch("X", &X, "X/D");
	+ outputTree->Branch("Y", &Y, "Y/D");
	+ outputTree->Branch("nP", &nP, "nP/I");
	+ outputTree->Branch("nN", &nN, "nN/I");
	+ outputTree->Branch("nipip", &nipip, "nipip/I");
	+ outputTree->Branch("nipim", &nipim, "nipim/I");
	+ outputTree->Branch("nipi0", &nipi0, "nipi0/I");
	+ outputTree->Branch("nikp", &nikp, "nikp/I");
	+ outputTree->Branch("nikm", &nikm, "nikm/I");
	+ outputTree->Branch("nik0", &nik0, "nik0/I");
	+ outputTree->Branch("niem", &niem, "niem/I");
	+ outputTree->Branch("niother", &niother, "niother/I");
	+ outputTree->Branch("nNucleus", &nNucleus, "nNucleus/I");
	+ outputTree->Branch("nUNKNOWN", &nUNKNOWN, "nUNKNOWN/I");
	+ outputTree->Branch("eP", &eP, "eP/D");
	+ outputTree->Branch("eN", &eN, "eN/D");
	+ outputTree->Branch("ePip", &ePip, "ePip/D");
	+ outputTree->Branch("ePim", &ePim, "ePim/D");
	+ outputTree->Branch("ePi0", &ePi0, "ePi0/D");
	+ outputTree->Branch("eOther", &eOther, "eOther/D");
	+ outputTree->Branch("xsecWeight", &xsecWeight, "xsecWeight/D");
	+
	+ //! Here you do whatever you want to do on a per-event basis.
	+ //! This method will be run for every event given by the input handler,
	+ //! the weight allows differential xsecs to be plotted easily and also
	+ //! contains any requested reweightable parameter variations.
	+ ProcessEventFunction = [&](nuis::event::FullEvent const &ev,
	+ double weight) -> void {
	+ xsecWeight = weight;
	+
	+ nP = 0;
	+ nN = 0;
	+ nipip = 0;
	+ nipim = 0;
	+ nipi0 = 0;
	+ nikp = 0;
	+ nikm = 0;
	+ nik0 = 0;
	+ niem = 0;
	+ niother = 0;
	+ nNucleus = 0;
	+ nUNKNOWN = 0;
	+ eP = 0;
	+ eN = 0;
	+ ePip = 0;
	+ ePim = 0;
	+ ePi0 = 0;
	+ eOther = 0;
	+
	+ Particle ISLep = GetHMISNeutralLepton(ev);
	+ Particle FSLep = GetHMFSLepton(ev);
	+
	+ isCC = IsChargedLepton(FSLep.pdg);
	+
	+ nuPDG = ISLep.pdg;
	+ nuPDGunosc = isOscSwap
	+ ? (ISLep.pdg < 0 ? pdgcodes::kNuMuBar : pdgcodes::kNuMu)
	+ : ISLep.pdg;
	+ LepPDG = FSLep.pdg;
	+
	+ NuMomX = ISLep.P4.X();
	+ NuMomY = ISLep.P4.Y();
	+ NuMomZ = ISLep.P4.Z();
	+ Ev = ISLep.P4.E();
	+
	+ LepMomX = FSLep.P4.X();
	+ LepMomY = FSLep.P4.Y();
	+ LepMomZ = FSLep.P4.Z();
	+ LepE = FSLep.P4.E();
	+
	+ LepNuAngle = FSLep.Dir().Dot(ISLep.Dir());
	+
	+ TLorentzVector FourMomTransfer = GetEnergyMomentumTransfer(ev);
	+
	+ Q2 = -FourMomTransfer.Mag2();
	+ W = sqrt(-Q2 +
	+ 2 * pdgmasses::kNucleonAverageMass_MeV * FourMomTransfer.E() +
	+ pdgmasses::kNucleonAverageMass_MeV *
	+ pdgmasses::kNucleonAverageMass_MeV);
	+ X = Q2 / (2 * pdgmasses::kNucleonAverageMass_MeV * FourMomTransfer.E());
	+ Y = 1 - LepE / Ev;
	+
	+ for (Particle const &p : GetNuclearLeavingParticles(ev)) {
	+ if (p == FSLep) {
	+ continue;
	+ }
	+
	+ if (IsProton(p.pdg)) {
	+ nP++;
	+ eP += p.KE();
	+ } else if (IsNeutron(p.pdg)) {
	+ nN++;
	+ eN += p.KE();
	+ } else if (IsPositivePion(p.pdg)) {
	+ nipip++;
	+ ePip += p.KE();
	+ } else if (IsNegativePion(p.pdg)) {
	+ nipim++;
	+ ePim += p.KE();
	+ } else if (IsNeutralPion(p.pdg)) {
	+ nipi0++;
	+ ePi0 += p.KE();
	+ } else if (IsPositiveKaon(p.pdg)) {
	+ nikp++;
	+ eOther += p.KE();
	+ } else if (IsNegativeKaon(p.pdg)) {
	+ nikm++;
	+ eOther += p.KE();
	+ } else if (IsNeutralKaon(p.pdg)) {
	+ nik0++;
	+ eOther += p.KE();
	+ } else if (IsGamma(p.pdg)) {
	+ niem++;
	+ eOther += p.KE();
	+ } else if (IsNuclearPDG(p.pdg)) {
	+ nNucleus++;
	+ } else {
	+ niother++;
	+ eOther += p.KE();
	+ }
	+ }
	+
	+ outputTree->Fill();
	+ };
	+ }
	+
	+ std::string Name() { return "NuisToCAFTruth"; }
	+
	+ void Write() {}
	+};
	+
	+//! These declarations allow your class to be loaded dynamically by NUISANCE
	+DECLARE_PLUGIN(IEventProcessor, NuisToCAFTruth);
	diff --git a/src/utility/FullEventUtility.cxx b/src/utility/FullEventUtility.cxx
	index 1406c6b..7046b2e 100644
	--- a/src/utility/FullEventUtility.cxx
	+++ b/src/utility/FullEventUtility.cxx
	@@ -1,169 +1,172 @@
	#include "utility/FullEventUtility.hxx"

	#include "event/FullEvent.hxx"

	#include "utility/PDGCodeUtility.hxx"

	using namespace nuis::event;

	namespace nuis {
	namespace utility {

	event::Particle GetHMParticle(std::vector<event::Particle> const &parts) {
	if (parts.size()) {
	return event::Particle();
	}

	return *std::max_element(
	parts.begin(), parts.end(),
	[](event::Particle const &l, event::Particle const &r) {
	return l.P() < r.P();
	});
	}

	std::vector<Particle> GetParticles(FullEvent const &ev,
	std::vector<PDG_t> const &pdgs,
	Particle::Status_t status,
	bool include_matching_pdg) {
	std::vector<Particle> selected_particles;
	for (auto const &parts : ev.ParticleStack) {
	if (parts.status != status) {
	continue;
	}
	for (Particle const &part : parts.particles) {
	bool matched_pdg = false;
	for (PDG_t pdg : pdgs) {
	matched_pdg = matched_pdg \|\| (part.pdg == pdg);
	}
	bool keep = ((include_matching_pdg && matched_pdg) \|\|
	(!include_matching_pdg && !matched_pdg));
	if (!keep) {
	continue;
	}
	selected_particles.push_back(part);
	}
	}
	return selected_particles;
	}

	std::vector<Particle> const &GetISParticles(FullEvent const &ev) {
	return ev
	.ParticleStack[static_cast<size_t>(
	Particle::Status_t::kPrimaryInitialState)]
	.particles;
	}
	std::vector<Particle> const &GetPrimaryFSParticles(FullEvent const &ev) {
	return ev
	.ParticleStack[static_cast<size_t>(
	Particle::Status_t::kPrimaryFinalState)]
	.particles;
	}
	std::vector<Particle> const &GetNuclearLeavingParticles(FullEvent const &ev) {
	return ev
	.ParticleStack[static_cast<size_t>(Particle::Status_t::kNuclearLeaving)]
	.particles;
	}

	Particle GetHMParticle(FullEvent const &ev, std::vector<PDG_t> const &pdgs,
	Particle::Status_t status, bool include_matching_pdg) {
	Particle HMParticle;
	for (auto const &parts : ev.ParticleStack) {
	if (parts.status != status) {
	continue;
	}
	for (Particle const &part : parts.particles) {
	bool matched_pdg = false;
	for (PDG_t pdg : pdgs) {
	matched_pdg = matched_pdg \|\| (part.pdg == pdg);
	}
	bool keep = ((include_matching_pdg && matched_pdg) \|\|
	(!include_matching_pdg && !matched_pdg));
	if (!keep) {
	continue;
	}
	if (part.P() > HMParticle.P()) {
	HMParticle = part;
	}
	}
	}
	return HMParticle;
	}

	event::Particle GetHMISParticle(event::FullEvent const &ev,
	std::vector<event::PDG_t> const &pdgs) {
	return GetHMParticle(ev, pdgs, Particle::Status_t::kPrimaryInitialState);
	}

	event::Particle GetHMFSParticle(event::FullEvent const &ev,
	std::vector<event::PDG_t> const &pdgs) {
	return GetHMParticle(ev, pdgs);
	}

	std::vector<Particle> GetFSChargedLeptons(FullEvent const &ev) {
	return GetParticles(ev, pdgcodes::ChargedLeptons);
	}
	std::vector<Particle> GetFSNeutralLeptons(FullEvent const &ev) {
	return GetParticles(ev, pdgcodes::NeutralLeptons);
	}
	std::vector<Particle> GetISNeutralLeptons(FullEvent const &ev) {
	return GetParticles(ev, pdgcodes::NeutralLeptons,
	Particle::Status_t::kPrimaryInitialState);
	}
	std::vector<Particle> GetFSChargedPions(FullEvent const &ev) {
	return GetParticles(ev, pdgcodes::ChargedPions);
	}
	std::vector<Particle> GetFSNeutralPions(FullEvent const &ev) {
	return GetParticles(ev, pdgcodes::NeutralPions);
	}
	std::vector<Particle> GetFSPions(FullEvent const &ev) {
	return GetParticles(ev, pdgcodes::Pions);
	}
	std::vector<Particle> GetFSProtons(FullEvent const &ev) {
	return GetParticles(ev, pdgcodes::Protons);
	}
	std::vector<Particle> GetFSNeutrons(FullEvent const &ev) {
	return GetParticles(ev, pdgcodes::Neutron);
	}
	std::vector<Particle> GetFSNucleons(FullEvent const &ev) {
	return GetParticles(ev, pdgcodes::Nucleons);
	}
	std::vector<Particle> GetFSOthers(FullEvent const &ev) {
	return GetParticles(ev, pdgcodes::CommonParticles,
	Particle::Status_t::kNuclearLeaving, false);
	}

	Particle GetHMFSChargedLepton(FullEvent const &ev) {
	return GetHMParticle(ev, pdgcodes::ChargedLeptons);
	}
	Particle GetHMFSNeutralLepton(FullEvent const &ev) {
	return GetHMParticle(ev, pdgcodes::NeutralLeptons);
	}
	+Particle GetHMFSLepton(event::FullEvent const &ev) {
	+ return GetHMParticle(ev, pdgcodes::Leptons);
	+}
	Particle GetHMISNeutralLepton(FullEvent const &ev) {
	return GetHMParticle(ev, pdgcodes::NeutralLeptons,
	Particle::Status_t::kPrimaryInitialState);
	}
	Particle GetHMFSChargedPion(FullEvent const &ev) {
	return GetHMParticle(ev, pdgcodes::ChargedPions);
	}
	Particle GetHMFSNeutralPion(FullEvent const &ev) {
	return GetHMParticle(ev, pdgcodes::NeutralPions);
	}
	Particle GetHMFSPion(FullEvent const &ev) {
	return GetHMParticle(ev, pdgcodes::Pions);
	}
	Particle GetHMFSProton(FullEvent const &ev) {
	return GetHMParticle(ev, pdgcodes::Protons);
	}
	Particle GetHMFSNeutron(FullEvent const &ev) {
	return GetHMParticle(ev, pdgcodes::Neutron);
	}
	Particle GetHMFSNucleon(FullEvent const &ev) {
	return GetHMParticle(ev, pdgcodes::Nucleons);
	}
	Particle GetHMFSOther(FullEvent const &ev) {
	return GetHMParticle(ev, pdgcodes::CommonParticles,
	Particle::Status_t::kNuclearLeaving, false);
	}

	} // namespace utility
	} // namespace nuis
	diff --git a/src/utility/FullEventUtility.hxx b/src/utility/FullEventUtility.hxx
	index aa25bea..22c7ac9 100644
	--- a/src/utility/FullEventUtility.hxx
	+++ b/src/utility/FullEventUtility.hxx
	@@ -1,84 +1,85 @@
	// 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 <http://www.gnu.org/licenses/>.
	*******************************************************************************/

	#pragma once

	#include "event/types.hxx"
	#include "event/Particle.hxx"

	#include <vector>

	namespace nuis {
	namespace event {
	class FullEvent;
	} // namespace event
	} // namespace nuis

	namespace nuis {
	namespace utility {

	event::Particle GetHMParticle(std::vector<event::Particle>);

	std::vector<event::Particle>
	GetParticles(event::FullEvent const &, std::vector<event::PDG_t> const &,
	event::Particle::Status_t status =
	event::Particle::Status_t::kNuclearLeaving,
	bool include_matching_pdg = true);

	std::vector<event::Particle> const &GetISParticles(event::FullEvent const &);
	std::vector<event::Particle> const &
	GetPrimaryFSParticles(event::FullEvent const &);
	std::vector<event::Particle> const &
	GetNuclearLeavingParticles(event::FullEvent const &);

	event::Particle GetHMParticle(event::FullEvent const &,
	std::vector<event::PDG_t> const &,
	event::Particle::Status_t status =
	event::Particle::Status_t::kNuclearLeaving,
	bool include_matching_pdg = true);

	event::Particle GetHMISParticle(event::FullEvent const &,
	std::vector<event::PDG_t> const &);
	event::Particle GetHMFSParticle(event::FullEvent const &,
	std::vector<event::PDG_t> const &);

	std::vector<event::Particle> GetFSChargedLeptons(event::FullEvent const &);
	std::vector<event::Particle> GetFSNeutralLeptons(event::FullEvent const &);
	std::vector<event::Particle> GetISNeutralLeptons(event::FullEvent const &);
	std::vector<event::Particle> GetFSChargedPions(event::FullEvent const &);
	std::vector<event::Particle> GetFSNeutralPions(event::FullEvent const &);
	std::vector<event::Particle> GetFSPions(event::FullEvent const &);
	std::vector<event::Particle> GetFSProtons(event::FullEvent const &);
	std::vector<event::Particle> GetFSNeutrons(event::FullEvent const &);
	std::vector<event::Particle> GetFSNucleons(event::FullEvent const &);
	std::vector<event::Particle> GetFSOthers(event::FullEvent const &);

	event::Particle GetHMFSChargedLepton(event::FullEvent const &);
	event::Particle GetHMFSNeutralLepton(event::FullEvent const &);
	+event::Particle GetHMFSLepton(event::FullEvent const &);
	event::Particle GetHMISNeutralLepton(event::FullEvent const &);
	event::Particle GetHMFSChargedPion(event::FullEvent const &);
	event::Particle GetHMFSNeutralPion(event::FullEvent const &);
	event::Particle GetHMFSPion(event::FullEvent const &);
	event::Particle GetHMFSProton(event::FullEvent const &);
	event::Particle GetHMFSNeutron(event::FullEvent const &);
	event::Particle GetHMFSNucleon(event::FullEvent const &);
	event::Particle GetHMFSOther(event::FullEvent const &);

	} // namespace utility
	} // namespace nuis
	diff --git a/src/utility/KinematicUtility.cxx b/src/utility/KinematicUtility.cxx
	index cfdbd69..be2a305 100644
	--- a/src/utility/KinematicUtility.cxx
	+++ b/src/utility/KinematicUtility.cxx
	@@ -1,228 +1,228 @@
	#include "utility/KinematicUtility.hxx"

	#include "event/FullEvent.hxx"
	#include "event/Particle.hxx"

	#include "utility/FullEventUtility.hxx"
	#include "utility/PDGCodeUtility.hxx"

	using namespace nuis::event;

	namespace nuis {
	namespace utility {

	double GetNeutrinoEQERec(FullEvent const &fev, double SeparationEnergy_MeV) {

	Particle const &neutrino = GetHMISNeutralLepton(fev);
	if (!neutrino) {
	throw Particle::invalid_particle()
	<< "[ERROR]: In GetNeutrinoEQERec, expected to be able to get IS "
	"neutral lepton, but found none:"
	<< "\n"
	<< fev.to_string();
	}

	Particle const &charged_lepton = GetHMFSChargedLepton(fev);
	if (!charged_lepton) {
	throw Particle::invalid_particle()
	<< "[ERROR]: In GetNeutrinoEQERec, expected to be able to get FS "
	"charged lepton, but found none:"
	<< "\n"
	<< fev.to_string();
	}

	double mass_ISN =
	(IsMatter(neutrino.pdg)
	? (pdgmasses::kNeutronMass_MeV - SeparationEnergy_MeV)
	: (pdgmasses::kProtonMass_MeV - SeparationEnergy_MeV)) /
	1000.0;

	double mass_FSN = (IsMatter(neutrino.pdg) ? pdgmasses::kProtonMass_MeV
	: pdgmasses::kNeutronMass_MeV) /
	1000.0;

	double el_GeV = charged_lepton.E() / 1000.0;
	double pl_GeV = charged_lepton.P() / 1000.0; // momentum of lepton
	double ml_GeV = charged_lepton.M() / 1000.0; // lepton mass

	double Theta_nu_mu = neutrino.P3().Angle(charged_lepton.P3());

	return (2.0 * mass_ISN * el_GeV - ml_GeV * ml_GeV + mass_FSN * mass_FSN -
	mass_ISN * mass_ISN) /
	(2.0 * (mass_ISN - el_GeV + pl_GeV * cos(Theta_nu_mu)));
	}

	double GetNeutrinoQ2QERec(FullEvent const &fev, double SeparationEnergy_MeV) {

	Particle const &neutrino = GetHMISNeutralLepton(fev);
	if (!neutrino) {
	throw Particle::invalid_particle()
	<< "[ERROR]: In GetNeutrinoQ2QERec, expected to be able to get IS "
	"neutral lepton, but found none:"
	<< "\n"
	<< fev.to_string();
	}

	Particle const &charged_lepton = GetHMFSChargedLepton(fev);
	if (!charged_lepton) {
	throw Particle::invalid_particle()
	<< "[ERROR]: In GetNeutrinoQ2QERec, expected to be able to get FS "
	"charged lepton, but found none:"
	<< "\n"
	<< fev.to_string();
	}

	double el_GeV = charged_lepton.E() / 1000.0;
	double pl_GeV = charged_lepton.P() / 1000.0; // momentum of lepton
	double ml_GeV = charged_lepton.M() / 1000.0; // lepton mass

	double Theta_nu_mu = neutrino.P3().Angle(charged_lepton.P3());

	return -ml_GeV * ml_GeV + 2.0 * GetNeutrinoEQERec(fev, SeparationEnergy_MeV) *
	(el_GeV - pl_GeV * cos(Theta_nu_mu));
	}

	TVector3 GetVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal) {
	TVector3 pnUnit = planarNormal.Unit();
	double inpProjectPN = inp.Dot(pnUnit);

	TVector3 InPlanarInput = inp - (inpProjectPN * pnUnit);
	return InPlanarInput;
	}

	TVector3 GetUnitVectorInTPlane(const TVector3 &inp,
	const TVector3 &planarNormal) {
	return GetVectorInTPlane(inp, planarNormal).Unit();
	}

	double GetDeltaPhiT(TVector3 const &V_lepton, TVector3 const &V_other,
	TVector3 const &Normal, bool PiMinus) {
	TVector3 V_lepton_T = GetUnitVectorInTPlane(V_lepton, Normal);

	TVector3 V_other_T = GetUnitVectorInTPlane(V_other, Normal);

	return PiMinus ? acos(V_lepton_T.Dot(V_other_T))
	: (M_PI - acos(V_lepton_T.Dot(V_other_T)));
	}

	TVector3 GetDeltaPT(TVector3 const &V_lepton, TVector3 const &V_other,
	TVector3 const &Normal) {
	TVector3 V_lepton_T = GetVectorInTPlane(V_lepton, Normal);

	TVector3 V_other_T = GetVectorInTPlane(V_other, Normal);

	return V_lepton_T + V_other_T;
	}

	double GetDeltaAlphaT(TVector3 const &V_lepton, TVector3 const &V_other,
	TVector3 const &Normal, bool PiMinus) {
	TVector3 DeltaPT = GetDeltaPT(V_lepton, V_other, Normal);

	return GetDeltaPhiT(V_lepton, DeltaPT, Normal, PiMinus);
	}

	TVector3 GetDeltaPT_CC0PiN(event::FullEvent const &fev, event::PDG_t islep_pdg,
	event::PDG_t fslep_pdg, event::PDG_t fsnuc_pdg) {
	if (fslep_pdg == nuis::utility::pdgcodes::kDefault) { // Assume neutrino CC
	fslep_pdg = islep_pdg > 0 ? islep_pdg - 1 : islep_pdg + 1;
	}
	if (fsnuc_pdg == nuis::utility::pdgcodes::kDefault) {
	fsnuc_pdg = islep_pdg > 0 ? nuis::utility::pdgcodes::kProton
	: nuis::utility::pdgcodes::kNeutron;
	}

	Particle ISSLep = GetHMISParticle(fev, {islep_pdg});
	if (!ISSLep) {
	return {-std::numeric_limits<double>::max(),
	-std::numeric_limits<double>::max(),
	-std::numeric_limits<double>::max()};
	}
	Particle FSLep = GetHMFSParticle(fev, {fslep_pdg});
	if (!FSLep) {
	return {-std::numeric_limits<double>::max(),
	-std::numeric_limits<double>::max(),
	-std::numeric_limits<double>::max()};
	}
	Particle FSNuc = GetHMFSParticle(fev, {fsnuc_pdg});
	if (!FSNuc) {
	return {-std::numeric_limits<double>::max(),
	-std::numeric_limits<double>::max(),
	-std::numeric_limits<double>::max()};
	}

	return GetDeltaPT(FSLep.P3(), FSNuc.P3(), ISSLep.P3());
	}
	double GetDeltaPhiT_CC0PiN(event::FullEvent const &fev, event::PDG_t islep_pdg,
	event::PDG_t fslep_pdg, event::PDG_t fsnuc_pdg) {
	if (fslep_pdg == nuis::utility::pdgcodes::kDefault) { // Assume neutrino CC
	fslep_pdg = islep_pdg > 0 ? islep_pdg - 1 : islep_pdg + 1;
	}
	if (fsnuc_pdg == nuis::utility::pdgcodes::kDefault) {
	fsnuc_pdg = islep_pdg > 0 ? nuis::utility::pdgcodes::kProton
	: nuis::utility::pdgcodes::kNeutron;
	}

	Particle ISSLep = GetHMISParticle(fev, {islep_pdg});
	if (!ISSLep) {
	return -std::numeric_limits<double>::max();
	}
	Particle FSLep = GetHMFSParticle(fev, {fslep_pdg});
	if (!FSLep) {
	return -std::numeric_limits<double>::max();
	}
	Particle FSNuc = GetHMFSParticle(fev, {fsnuc_pdg});
	if (!FSNuc) {
	return -std::numeric_limits<double>::max();
	}

	return GetDeltaPhiT(FSLep.P3(), FSNuc.P3(), ISSLep.P3());
	}
	double GetDeltaAlphaT_CC0PiN(event::FullEvent const &fev,
	event::PDG_t islep_pdg, event::PDG_t fslep_pdg,
	event::PDG_t fsnuc_pdg) {
	if (fslep_pdg == nuis::utility::pdgcodes::kDefault) { // Assume neutrino CC
	fslep_pdg = islep_pdg > 0 ? islep_pdg - 1 : islep_pdg + 1;
	}
	if (fsnuc_pdg == nuis::utility::pdgcodes::kDefault) {
	fsnuc_pdg = islep_pdg > 0 ? nuis::utility::pdgcodes::kProton
	: nuis::utility::pdgcodes::kNeutron;
	}

	Particle ISSLep = GetHMISParticle(fev, {islep_pdg});
	if (!ISSLep) {
	return -std::numeric_limits<double>::max();
	}
	Particle FSLep = GetHMFSParticle(fev, {fslep_pdg});
	if (!FSLep) {
	return -std::numeric_limits<double>::max();
	}
	Particle FSNuc = GetHMFSParticle(fev, {fsnuc_pdg});
	if (!FSNuc) {
	return -std::numeric_limits<double>::max();
	}

	return GetDeltaAlphaT(FSLep.P3(), FSNuc.P3(), ISSLep.P3());
	}

	TLorentzVector GetEnergyMomentumTransfer(event::FullEvent const &fev) {
	Particle const &neutrino = GetHMISNeutralLepton(fev);
	if (!neutrino) {
	throw Particle::invalid_particle()
	<< "[ERROR]: In GetEnergyMomentumTransfer, expected to be able to get "
	"IS neutral lepton, but found none: \n"
	<< fev.to_string();
	}

	- Particle const &charged_lepton = GetHMFSChargedLepton(fev);
	- if (!charged_lepton) {
	+ Particle const &fs_lepton = GetHMFSLepton(fev);
	+ if (!fs_lepton) {
	throw Particle::invalid_particle()
	<< "[ERROR]: In GetEnergyMomentumTransfer, expected to be able to get "
	- "FS charged lepton, but found none: \n"
	+ "FS lepton, but found none: \n"
	<< fev.to_string();
	}

	- return (neutrino.P4 - charged_lepton.P4);
	+ return (neutrino.P4 - fs_lepton.P4);
	}

	} // namespace utility
	} // namespace nuis
	diff --git a/src/utility/PDGCodeUtility.cxx b/src/utility/PDGCodeUtility.cxx
	index f777df9..d2d7e4b 100644
	--- a/src/utility/PDGCodeUtility.cxx
	+++ b/src/utility/PDGCodeUtility.cxx
	@@ -1,118 +1,128 @@
	#include "utility/PDGCodeUtility.hxx"

	#include <algorithm>
	#include <map>

	using namespace nuis::event;
	using namespace nuis::utility::pdgcodes;
	using namespace nuis::utility::pdgmasses;

	namespace nuis {
	namespace utility {

	static std::map<PDG_t, double> const PDGMasses{
	{kNuMu, kNuMuMass_MeV}, {kNuMuBar, kNuMuBarMass_MeV},
	{kMu, kMuMass_MeV}, {kMuPlus, kMuPlusMass_MeV},
	{kNue, kNueMass_MeV}, {kNueBar, kNueBarMass_MeV},
	{kPiPlus, kPiPlusMass_MeV}, {kPiMinus, kPiMinusMass_MeV},
	{kPi0, kPi0Mass_MeV}, {kProton, kNeutronMass_MeV},
	{kNeutron, kProtonMass_MeV}};

	double GetPDGMass(PDG_t pdg) {
	if (PDGMasses.find(pdg) == PDGMasses.end()) {
	throw unhandled_pdg_code()
	<< "[ERROR]: Unknown mass for particle with PDG code: " << pdg
	<< ", please add it to src/utility/PDGCodeUtility.cxx:PDGMasses";
	}
	return PDGMasses.at(pdg);
	}

	NEW_NUIS_EXCEPT(invalid_MatterType);

	bool IsInPDGList(PDG_t pdg, std::vector<PDG_t> const &MatterList,
	std::vector<PDG_t> const &AntiMatterList, MatterType type) {
	switch (type) {
	case kMatter: {
	return std::count(MatterList.begin(), MatterList.end(), pdg);
	}
	case kMatterAntimatter: {
	return std::count(MatterList.begin(), MatterList.end(), pdg) \|\|
	std::count(AntiMatterList.begin(), AntiMatterList.end(), pdg);
	}
	case kAntimatter: {
	return std::count(AntiMatterList.begin(), AntiMatterList.end(), pdg);
	}
	default: { throw invalid_MatterType(); }
	}
	}
	bool IsNeutralLepton(PDG_t pdg, MatterType type) {
	return IsInPDGList(pdg, NeutralLeptons_matter, NeutralLeptons_antimatter,
	type);
	}
	bool IsChargedLepton(PDG_t pdg, MatterType type) {
	return IsInPDGList(pdg, ChargedLeptons_matter, ChargedLeptons_antimatter,
	type);
	}
	bool IsProton(PDG_t pdg, MatterType type) {
	return IsInPDGList(pdg, Proton_matter, Proton_antimatter, type);
	}
	bool IsNeutron(PDG_t pdg) { return pdg == Neutron[0]; }
	bool IsChargedPion(PDG_t pdg) {
	return std::count(ChargedPions.begin(), ChargedPions.end(), pdg);
	}
	-bool IsNeutralPion(PDG_t pdg) {
	- return std::count(NeutralPions.begin(), NeutralPions.end(), pdg);
	-}
	+bool IsNeutralPion(PDG_t pdg) { return (pdg == kPi0); }
	+bool IsPositivePion(PDG_t pdg) { return (pdg == kPiPlus); }
	+bool IsNegativePion(PDG_t pdg) { return (pdg == kPiMinus); }
	bool IsPion(PDG_t pdg) { return std::count(Pions.begin(), Pions.end(), pdg); }
	+bool IsChargedKaon(PDG_t pdg) {
	+ return std::count(ChargedKaons.begin(), ChargedKaons.end(), pdg);
	+}
	+bool IsNeutralKaon(PDG_t pdg) {
	+ return std::count(NeutralKaons.begin(), NeutralKaons.end(), pdg);
	+}
	+bool IsPositiveKaon(PDG_t pdg) { return (pdg == kKPlus); }
	+bool IsNegativeKaon(PDG_t pdg) { return (pdg == kKMinus); }
	+bool IsKaon(PDG_t pdg) { return std::count(Kaons.begin(), Kaons.end(), pdg); }
	+bool IsGamma(PDG_t pdg) { return (pdg == kGamma); }
	bool IsOther(PDG_t pdg) {
	return !std::count(CommonParticles.begin(), CommonParticles.end(), pdg);
	}

	bool IsMatter(PDG_t pdg) {
	// Special cases
	switch (pdg) {
	case kPiPlus: {
	return true;
	}
	case kPiMinus: {
	return true;
	}
	case kPi0: {
	return true;
	}
	case kNeutron: {
	return true;
	}
	}
	return (pdg > 0);
	}
	bool IsAntiMatter(PDG_t pdg) {
	// Special cases
	switch (pdg) {
	case kPiPlus: {
	return true;
	}
	case kPiMinus: {
	return true;
	}
	case kPi0: {
	return true;
	}
	case kNeutron: {
	return true;
	}
	}
	return (pdg < 0);
	}

	bool IsNuclearPDG(event::PDG_t pdg) { return (pdg > 1000000000); }

	PDG_t MakeNuclearPDG(size_t A, size_t Z) {
	return 1000 * Z + 10 * A + 1000000000;
	}

	size_t GetA(PDG_t pdg) { return ((pdg / 10) % 1000); }
	size_t GetZ(PDG_t pdg) { return ((pdg / 1000) % 1000); }

	} // namespace utility
	} // namespace nuis
	diff --git a/src/utility/PDGCodeUtility.hxx b/src/utility/PDGCodeUtility.hxx
	index 3cf290c..6cf5c7f 100644
	--- a/src/utility/PDGCodeUtility.hxx
	+++ b/src/utility/PDGCodeUtility.hxx
	@@ -1,133 +1,157 @@
	// 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 <http://www.gnu.org/licenses/>.
	*******************************************************************************/

	#pragma once

	#include "event/types.hxx"

	#include "exception/exception.hxx"

	#include <vector>

	namespace nuis {
	namespace utility {
	namespace pdgcodes {

	NEW_NUIS_EXCEPT(unhandled_pdg_code);

	enum MatterType { kMatter = 1, kMatterAntimatter = 0, kAntimatter = -1 };

	static event::PDG_t const kDefault = 0;

	static event::PDG_t const kNuMu = 14;
	static event::PDG_t const kNuMuBar = -14;

	static event::PDG_t const kMu = 13;
	static event::PDG_t const kMuPlus = -13;

	static event::PDG_t const kElectron = 11;
	static event::PDG_t const kPositron = -11;

	static event::PDG_t const kNue = 12;
	static event::PDG_t const kNueBar = 12;

	static event::PDG_t const kPiPlus = 211;
	static event::PDG_t const kPiMinus = -211;
	static event::PDG_t const kPi0 = 111;

	+static event::PDG_t const kK0Long = 130;
	+static event::PDG_t const kK0Short = 130;
	+static event::PDG_t const kK0 = 311;
	+static event::PDG_t const kKPlus = 321;
	+static event::PDG_t const kKMinus = -321;
	+
	static event::PDG_t const kProton = 2212;
	static event::PDG_t const kNeutron = 2112;

	+static event::PDG_t const kGamma = 22;
	+
	static std::vector<event::PDG_t> const ChargedLeptons{kElectron, kMu, 15,
	kPositron, kMuPlus, -15};
	static std::vector<event::PDG_t> const ChargedLeptons_matter{kElectron, kMu,
	15};
	static std::vector<event::PDG_t> const ChargedLeptons_antimatter{kPositron,
	kMuPlus, -15};

	static std::vector<event::PDG_t> const NeutralLeptons{kNue, kNuMu, 16,
	kNueBar, kNuMuBar, -16};
	static std::vector<event::PDG_t> const NeutralLeptons_matter{kNue, kNuMu, 16};
	static std::vector<event::PDG_t> const NeutralLeptons_antimatter{kNueBar,
	kNuMuBar, -16};
	+static std::vector<event::PDG_t> const Leptons{
	+ kElectron, kMu, 15, kPositron, kMuPlus, -15,
	+ kNue, kNuMu, 16, kNueBar, kNuMuBar, -16};

	static std::vector<event::PDG_t> const ChargedPions{kPiPlus, kPiMinus};
	static std::vector<event::PDG_t> const NeutralPions{kPi0};
	static std::vector<event::PDG_t> const Pions{kPiPlus, kPiMinus, kPi0};
	+static std::vector<event::PDG_t> const ChargedKaons{kKPlus, kKMinus};
	+static std::vector<event::PDG_t> const NeutralKaons{kK0, kK0Long, kK0Short};
	+static std::vector<event::PDG_t> const Kaons{kKPlus, kKMinus, kK0, kK0Long,
	+ kK0Short};
	static std::vector<event::PDG_t> const Protons{kProton, -kProton};
	static std::vector<event::PDG_t> const Proton_matter{kProton};
	static std::vector<event::PDG_t> const Proton_antimatter{-kProton};
	static std::vector<event::PDG_t> const Neutron{kNeutron};
	static std::vector<event::PDG_t> const Nucleons{kProton, kNeutron, -kProton};
	static std::vector<event::PDG_t> const Nucleons_matter{kProton, kNeutron};
	static std::vector<event::PDG_t> const Nucleons_antimatter{-kProton, kNeutron};

	static std::vector<event::PDG_t> const CommonParticles{
	kElectron, kMu, 15, kPositron, kMuPlus, -15, kNue, kNuMu, 16,
	kNueBar, kNuMuBar, -16, kPiPlus, kPiMinus, kPi0, kProton, kNeutron};
	} // namespace pdgcodes

	namespace pdgmasses {
	static double const kNuMuMass_MeV = 0;
	static double const kNuMuBarMass_MeV = 0;

	static double const kMuMass_MeV = 105.65;
	static double const kMuPlusMass_MeV = 105.65;

	static double const kNueMass_MeV = 0;
	static double const kNueBarMass_MeV = 0;

	static double const kPiPlusMass_MeV = 139.57;
	static double const kPiMinusMass_MeV = 139.57;
	static double const kPi0Mass_MeV = 134.97;

	static double const kNeutronMass_MeV = 939.56;
	static double const kProtonMass_MeV = 938.27;
	+static double const kNucleonAverageMass_MeV = 938.27;
	} // namespace pdgmasses

	double GetPDGMass(event::PDG_t);

	bool IsInPDGList(event::PDG_t pdg, std::vector<event::PDG_t> const &MatterList,
	std::vector<event::PDG_t> const &AntiMatterList,
	pdgcodes::MatterType type = pdgcodes::kMatterAntimatter);
	bool IsNeutralLepton(event::PDG_t pdg,
	pdgcodes::MatterType type = pdgcodes::kMatterAntimatter);
	bool IsChargedLepton(event::PDG_t pdg,
	pdgcodes::MatterType type = pdgcodes::kMatterAntimatter);
	bool IsProton(event::PDG_t pdg,
	pdgcodes::MatterType type = pdgcodes::kMatterAntimatter);

	bool IsNeutron(event::PDG_t pdg);
	bool IsChargedPion(event::PDG_t pdg);
	+bool IsPositivePion(event::PDG_t pdg);
	+bool IsNegativePion(event::PDG_t pdg);
	bool IsNeutralPion(event::PDG_t pdg);
	bool IsPion(event::PDG_t pdg);
	+bool IsChargedKaon(event::PDG_t pdg);
	+bool IsPositiveKaon(event::PDG_t pdg);
	+bool IsNegativeKaon(event::PDG_t pdg);
	+bool IsNeutralKaon(event::PDG_t pdg);
	+bool IsKaon(event::PDG_t pdg);
	+bool IsGamma(event::PDG_t pdg);
	bool IsOther(event::PDG_t pdg);

	bool IsMatter(event::PDG_t pdg);
	bool IsAntiMatter(event::PDG_t pdg);

	bool IsNuclearPDG(event::PDG_t pdg);

	event::PDG_t MakeNuclearPDG(size_t A, size_t Z);

	size_t GetA(event::PDG_t);
	size_t GetZ(event::PDG_t);

	} // namespace utility
	} // namespace nuis
	diff --git a/src/utility/ROOTUtility.hxx b/src/utility/ROOTUtility.hxx
	index bc62623..60551e8 100644
	--- a/src/utility/ROOTUtility.hxx
	+++ b/src/utility/ROOTUtility.hxx
	@@ -1,142 +1,155 @@
	// 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 <http://www.gnu.org/licenses/>.
	*******************************************************************************/

	#ifndef UTILITY_ROOTUTILITY_HXX_SEEN
	#define UTILITY_ROOTUTILITY_HXX_SEEN

	#include "exception/exception.hxx"

	#include "TFile.h"
	#include "TTree.h"

	#include <iomanip>
	#include <iostream>
	#include <memory>
	#include <string>

	namespace nuis {
	namespace utility {

	NEW_NUIS_EXCEPT(failed_to_open_TFile);
	NEW_NUIS_EXCEPT(failed_to_get_TTree);

	inline void CloseTFile(TFile *f = nullptr) {
	if (f) {
	std::cout << "[INFO]: Shutting TFile: " << f->GetName() << ", " << f
	<< std::endl;
	f->Close();
	}
	delete f;
	}

	using TFile_ptr = std::unique_ptr<TFile, decltype(&CloseTFile)>;
	inline TFile_ptr make_unique_TFile(TFile *f) {
	return TFile_ptr(f, &CloseTFile);
	}

	inline TFile_ptr CheckOpenTFile(std::string const &fname,
	char const *opts = "") {

	TDirectory *ogDir = gDirectory;

	TFile *inpF = new TFile(fname.c_str(), opts);
	if (!inpF \|\| !inpF->IsOpen()) {
	throw failed_to_open_TFile()
	<< "[ERROR]: Couldn't open input file: " << std::quoted(fname);
	}

	if (ogDir) {
	ogDir->cd();
	}

	return make_unique_TFile(inpF);
	}

	struct TreeFile {
	TFile_ptr file;
	TTree *tree;
	bool file_owned;
	TreeFile()
	: file(make_unique_TFile(nullptr)), tree(nullptr), file_owned(false) {}
	TreeFile(TreeFile const &other) = delete;
	TreeFile &operator=(TreeFile &&other) {
	file = std::move(other.file);
	tree = other.tree;
	file_owned = other.file_owned;

	other.file = nullptr;
	other.tree = nullptr;
	other.file_owned = false;

	return *this;
	}

	TreeFile(TreeFile &&other)
	: file(std::move(other.file)), tree(other.tree),
	file_owned(other.file_owned) {
	other.file = nullptr;
	other.tree = nullptr;
	other.file_owned = false;
	}
	+
	+ TTree *operator->() { return tree; }
	+
	~TreeFile() {
	if (!file_owned) {
	file.release();
	}
	}
	};

	inline TreeFile CheckGetTTree(TFile *file, std::string const &tname) {
	TreeFile tf;
	tf.file = make_unique_TFile(file);
	tf.tree = dynamic_cast<TTree *>(tf.file->Get(tname.c_str()));
	tf.file_owned = false;

	if (!tf.tree) {
	throw failed_to_get_TTree()
	<< "[ERROR]: Failed to get TTree named: " << std::quoted(tname)
	<< " from TFile: " << std::quoted(file->GetName());
	}
	std::cout << "[INFO]: Opened TFile: " << file->GetName() << ", " << file
	<< std::endl;
	return tf;
	}

	inline TreeFile CheckGetTTree(std::string const &fname,
	std::string const &tname, char const *opts = "") {
	TreeFile tf = CheckGetTTree(CheckOpenTFile(fname, opts).release(), tname);
	tf.file_owned = true;
	return tf;
	}

	inline TreeFile MakeNewTTree(std::string const &fname, std::string const &tname,
	char const *opts = "") {
	TreeFile tf;
	tf.file = CheckOpenTFile(fname, opts);
	tf.tree = new TTree(tname.c_str(), "");
	tf.tree->SetDirectory(tf.file.get());
	tf.file_owned = true;
	return tf;
	}

	+inline TreeFile AddNewTTreeToFile(TFile *f, std::string const &tname) {
	+ TreeFile tf;
	+ // This is pretty dumb...
	+ tf.file = make_unique_TFile(f);
	+ tf.tree = new TTree(tname.c_str(), "");
	+ tf.tree->SetDirectory(tf.file.get());
	+ tf.file_owned = false;
	+ return tf;
	+}
	+
	inline std::string SanitizeROOTObjectName(std::string name) { return name; }

	} // namespace utility
	} // namespace nuis

	#endif

File Metadata

Mime Type: text/x-diff
Expires: Wed, May 14, 10:58 AM (20 h, 50 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 5084395
Default Alt Text: (45 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions