No OneTemporary
Actions

Size

110 KB

Subscribers

None

View Options

	diff --git a/src/MCStudies/Smearceptance_Tester.cxx b/src/MCStudies/Smearceptance_Tester.cxx
	index 468278e..6110443 100644
	--- a/src/MCStudies/Smearceptance_Tester.cxx
	+++ b/src/MCStudies/Smearceptance_Tester.cxx
	@@ -1,646 +1,690 @@
	// 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 <http://www.gnu.org/licenses/>.
	*******************************************************************************/

	#include "Smearceptance_Tester.h"

	#include "SmearceptanceUtils.h"

	#include "Smearcepterton.h"

	//********************************************************************
	/// @brief Class to perform smearceptance MC Studies on a custom measurement
	Smearceptance_Tester::Smearceptance_Tester(std::string name,
	std::string inputfile, FitWeight *rw,
	std::string type,
	std::string fakeDataFile) {
	//********************************************************************

	// Measurement Details
	std::vector<std::string> splitName = GeneralUtils::ParseToStr(name, "_");
	size_t firstUS = name.find_first_of("_");

	std::string smearceptorName = name.substr(firstUS + 1);

	fName = name.substr(0, firstUS);
	eventVariables = NULL;

	// Define our energy range for flux calcs
	EnuMin = 0.;
	EnuMax = 100.; // Arbritrarily high energy limit

	// Set default fitter flags
	fIsDiag = true;
	fIsShape = false;
	fIsRawEvents = false;

	// This function will sort out the input files automatically and parse all the
	// inputs,flags,etc.
	// There may be complex cases where you have to do this by hand, but usually
	// this will do.
	Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile);

	eventVariables = NULL;

	// Setup fDataHist as a placeholder
	this->fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1);
	this->SetupDefaultHist();
	fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist);
	covar = StatUtils::GetInvert(fFullCovar);

	this->fScaleFactor =
	(GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) /
	this->TotalIntegratedFlux();

	LOG(SAM) << "Smearceptance Flux Scaling Factor = " << fScaleFactor
	<< std::endl;

	if (fScaleFactor <= 0.0) {
	ERR(WRN) << "SCALE FACTOR TOO LOW " << std::endl;
	sleep(20);
	}

	// Setup our TTrees
	this->AddEventVariablesToTree();

	smearceptor = &Smearcepterton::Get().GetSmearcepter(smearceptorName);

	Int_t RecNBins = 20, TrueNBins = 20;
	double RecBinL = 0, TrueBinL = 0, RecBinH = 10, TrueBinH = 10;

	if (Config::Get().GetConfigNode("smear.reconstructed.binning")) {
	std::vector<std::string> args = GeneralUtils::ParseToStr(
	Config::Get().ConfS("smear.reconstructed.binning"), ",");
	RecNBins = GeneralUtils::StrToInt(args[0]);
	RecBinL = GeneralUtils::StrToDbl(args[1]);
	RecBinH = GeneralUtils::StrToDbl(args[2]);
	TrueNBins = RecNBins;
	TrueBinL = RecBinL;
	TrueBinH = RecBinH;
	}

	if (Config::Get().GetConfigNode("smear.true.binning")) {
	std::vector<std::string> args = GeneralUtils::ParseToStr(
	Config::Get().ConfS("smear.true.binning"), ",");
	TrueNBins = GeneralUtils::StrToInt(args[0]);
	TrueBinL = GeneralUtils::StrToDbl(args[1]);
	TrueBinH = GeneralUtils::StrToDbl(args[2]);
	}
	SVDTruncation = 0;
	if (Config::Get().GetConfigNode("smear.true.binning")) {
	SVDTruncation = Config::Get().ConfI("smear.SVD.truncation");
	QLOG(SAM, "Applying SVD truncation of: " << SVDTruncation)
	}

	QLOG(SAM, "Using binning True: " << TrueNBins << ", [" << TrueBinL << " -- "
	<< TrueBinH << "], Rec: " << RecNBins
	<< ", [" << RecBinL << " -- " << RecBinH
	<< "]");

	ETrueDistrib = new TH1D("ELep_rate", ";True E_{#nu};Count", TrueNBins,
	TrueBinL, TrueBinH);
	ERecDistrib = new TH1D("ELepRec_rate", ";Rec E_{#nu};Count", RecNBins,
	RecBinL, RecBinH);
	ETrueDistrib->Sumw2();
	ERecDistrib->Sumw2();

	RecoSmear =
	new TH2D("ELepHadVis_Recon", ";Recon. E_{#nu};True E_{#nu}", RecNBins,
	RecBinL, RecBinH, TrueNBins, TrueBinL, TrueBinH);
	}

	void Smearceptance_Tester::AddEventVariablesToTree() {
	// Setup the TTree to save everything
	if (!eventVariables) {
	FitPar::Config().out->cd();
	eventVariables = new TTree((this->fName + "_VARS").c_str(),
	(this->fName + "_VARS").c_str());
	}

	LOG(SAM) << "Adding Event Variables" << std::endl;

	eventVariables->Branch("Omega_true", &Omega_true, "Omega_true/F");
	eventVariables->Branch("Q2_true", &Q2_true, "Q2_true/F");
	eventVariables->Branch("Mode_true", &Mode_true, "Mode_true/I");

	eventVariables->Branch("EISLep_true", &EISLep_true, "EISLep_true/F");

	+ eventVariables->Branch("HMFS_mu_true", &HMFS_mu_true);
	+ eventVariables->Branch("HMFS_pip_true", &HMFS_pip_true);
	+ eventVariables->Branch("HMFS_pim_true", &HMFS_pim_true);
	+ eventVariables->Branch("HMFS_cpi_true", &HMFS_cpi_true);
	+ eventVariables->Branch("HMFS_p_true", &HMFS_p_true);
	+
	eventVariables->Branch("KEFSHad_cpip_true", &KEFSHad_cpip_true,
	"KEFSHad_cpip_true/F");
	eventVariables->Branch("KEFSHad_cpim_true", &KEFSHad_cpim_true,
	"KEFSHad_cpim_true/F");
	eventVariables->Branch("KEFSHad_cpi_true", &KEFSHad_cpi_true,
	"KEFSHad_cpi_true/F");
	eventVariables->Branch("TEFSHad_pi0_true", &TEFSHad_pi0_true,
	"TEFSHad_pi0_true/F");
	eventVariables->Branch("KEFSHad_p_true", &KEFSHad_p_true, "KEFSHad_p_true/F");
	eventVariables->Branch("KEFSHad_n_true", &KEFSHad_n_true, "KEFSHad_n_true/F");

	eventVariables->Branch("EFSHad_true", &EFSHad_true, "EFSHad_true/F");
	eventVariables->Branch("EFSChargedEMHad_true", &EFSChargedEMHad_true,
	"EFSChargedEMHad_true/F");

	eventVariables->Branch("EFSLep_true", &EFSLep_true, "EFSLep_true/F");
	eventVariables->Branch("EFSgamma_true", &EFSgamma_true, "EFSgamma_true/F");

	eventVariables->Branch("PDGISLep_true", &PDGISLep_true, "PDGISLep_true/I");
	eventVariables->Branch("PDGFSLep_true", &PDGFSLep_true, "PDGFSLep_true/I");

	eventVariables->Branch("Nprotons_true", &Nprotons_true, "Nprotons_true/I");
	eventVariables->Branch("Nneutrons_true", &Nneutrons_true, "Nneutrons_true/I");
	eventVariables->Branch("Ncpiplus_true", &Ncpiplus_true, "Ncpiplus_true/I");
	eventVariables->Branch("Ncpiminus_true", &Ncpiminus_true, "Ncpiminus_true/I");
	eventVariables->Branch("Ncpi_true", &Ncpi_true, "Ncpi_true/I");
	eventVariables->Branch("Npi0_true", &Npi0_true, "Npi0_true/I");

	+ eventVariables->Branch("HMFS_mu_rec", &HMFS_mu_rec);
	+ eventVariables->Branch("HMFS_pip_rec", &HMFS_pip_rec);
	+ eventVariables->Branch("HMFS_pim_rec", &HMFS_pim_rec);
	+ eventVariables->Branch("HMFS_cpi_rec", &HMFS_cpi_rec);
	+ eventVariables->Branch("HMFS_p_rec", &HMFS_p_rec);
	+
	eventVariables->Branch("KEFSHad_cpip_rec", &KEFSHad_cpip_rec,
	"KEFSHad_cpip_rec/F");
	eventVariables->Branch("KEFSHad_cpim_rec", &KEFSHad_cpim_rec,
	"KEFSHad_cpim_rec/F");
	eventVariables->Branch("KEFSHad_cpi_rec", &KEFSHad_cpi_rec,
	"KEFSHad_cpi_rec/F");
	eventVariables->Branch("TEFSHad_pi0_rec", &TEFSHad_pi0_rec,
	"TEFSHad_pi0_rec/F");
	eventVariables->Branch("KEFSHad_p_rec", &KEFSHad_p_rec, "KEFSHad_p_rec/F");
	eventVariables->Branch("KEFSHad_n_rec", &KEFSHad_n_rec, "KEFSHad_n_rec/F");

	eventVariables->Branch("EFSHad_rec", &EFSHad_rec, "EFSHad_rec/F");
	eventVariables->Branch("EFSLep_rec", &EFSLep_rec, "EFSLep_rec/F");

	eventVariables->Branch("EFSVis_cpip", &EFSVis_cpip, "EFSVis_cpip/F");
	eventVariables->Branch("EFSVis_cpim", &EFSVis_cpim, "EFSVis_cpim/F");
	eventVariables->Branch("EFSVis_cpi", &EFSVis_cpi, "EFSVis_cpi/F");
	eventVariables->Branch("EFSVis_pi0", &EFSVis_pi0, "EFSVis_pi0/F");
	eventVariables->Branch("EFSVis_p", &EFSVis_p, "EFSVis_p/F");
	eventVariables->Branch("EFSVis_n", &EFSVis_n, "EFSVis_n/F");
	eventVariables->Branch("EFSVis_gamma", &EFSVis_gamma, "EFSVis_gamma/F");
	eventVariables->Branch("EFSVis_other", &EFSVis_other, "EFSVis_other/F");
	eventVariables->Branch("EFSVis", &EFSVis, "EFSVis/F");

	eventVariables->Branch("FSCLep_seen", &FSCLep_seen, "FSCLep_seen/I");
	eventVariables->Branch("Nprotons_seen", &Nprotons_seen, "Nprotons_seen/I");
	eventVariables->Branch("Nneutrons_seen", &Nneutrons_seen, "Nneutrons_seen/I");
	eventVariables->Branch("Ncpip_seen", &Ncpip_seen, "Ncpip_seen/I");
	eventVariables->Branch("Ncpim_seen", &Ncpim_seen, "Ncpim_seen/I");
	eventVariables->Branch("Ncpi_seen", &Ncpi_seen, "Ncpi_seen/I");
	eventVariables->Branch("Npi0_seen", &Npi0_seen, "Npi0_seen/I");
	eventVariables->Branch("Nothers_seen", &Nothers_seen, "Nothers_seen/I");

	eventVariables->Branch("EISLep_QE_rec", &EISLep_QE_rec, "EISLep_QE_rec/F");
	eventVariables->Branch("EISLep_LepHad_rec", &EISLep_LepHad_rec,
	"EISLep_LepHad_rec/F");
	eventVariables->Branch("EISLep_LepHadVis_rec", &EISLep_LepHadVis_rec,
	"EISLep_LepHadVis_rec/F");

	eventVariables->Branch("Nprotons_contributed", &Nprotons_contributed,
	"Nprotons_contributed/I");
	eventVariables->Branch("Nneutrons_contributed", &Nneutrons_contributed,
	"Nneutrons_contributed/I");
	eventVariables->Branch("Ncpip_contributed", &Ncpip_contributed,
	"Ncpip_contributed/I");
	eventVariables->Branch("Ncpim_contributed", &Ncpim_contributed,
	"Ncpim_contributed/I");
	eventVariables->Branch("Ncpi_contributed", &Ncpi_contributed,
	"Ncpi_contributed/I");
	eventVariables->Branch("Npi0_contributed", &Npi0_contributed,
	"Npi0_contributed/I");
	eventVariables->Branch("Ngamma_contributed", &Ngamma_contributed,
	"Ngamma_contributed/I");
	eventVariables->Branch("Nothers_contibuted", &Nothers_contibuted,
	"Nothers_contibuted/I");

	eventVariables->Branch("Weight", &Weight, "Weight/F");
	eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F");
	eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F");
	eventVariables->Branch("FluxWeight", &FluxWeight, "FluxWeight/F");
	eventVariables->Branch("EffWeight", &EffWeight, "EffWeight/F");

	eventVariables->Branch("xsecScaling", &xsecScaling, "xsecScaling/F");

	eventVariables->Branch("flagCCINC_true", &flagCCINC_true, "flagCCINC_true/O");
	eventVariables->Branch("flagCC0Pi_true", &flagCC0Pi_true, "flagCC0Pi_true/O");

	eventVariables->Branch("flagCCINC_rec", &flagCCINC_rec, "flagCCINC_rec/O");
	eventVariables->Branch("flagCC0Pi_rec", &flagCC0Pi_rec, "flagCC0Pi_rec/O");
	-
	-#ifdef DEBUG_SMEARTESTER
	- eventVariables->Branch("FSMuon_True", &FSMuon_True);
	- eventVariables->Branch("FSMuon_Smeared", &FSMuon_Smeared);
	-#endif
	}

	template <size_t N>
	-int CountNPdgsSeen(RecoInfo ri, int (&pdgs)[N]) {
	+int CountNPdgsSeen(RecoInfo ri, int const (&pdgs)[N]) {
	int sum = 0;
	for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
	sum +=
	std::count(ri.RecObjClass.begin(), ri.RecObjClass.end(), pdgs[pdg_it]);
	}
	return sum;
	}

	template <size_t N>
	-int CountNNotPdgsSeen(RecoInfo ri, int (&pdgs)[N]) {
	+int CountNNotPdgsSeen(RecoInfo ri, int const (&pdgs)[N]) {
	int sum = 0;
	for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
	sum +=
	(std::count(ri.RecObjClass.begin(), ri.RecObjClass.end(), pdgs[pdg_it])
	? 0
	: 1);
	}
	return sum;
	}

	template <size_t N>
	-int CountNPdgsContributed(RecoInfo ri, int (&pdgs)[N]) {
	+int CountNPdgsContributed(RecoInfo ri, int const (&pdgs)[N]) {
	int sum = 0;
	for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
	sum += std::count(ri.TrueContribPDGs.begin(), ri.TrueContribPDGs.end(),
	pdgs[pdg_it]);
	}
	return sum;
	}

	template <size_t N>
	-int CountNNotPdgsContributed(RecoInfo ri, int (&pdgs)[N]) {
	+int CountNNotPdgsContributed(RecoInfo ri, int const (&pdgs)[N]) {
	int sum = 0;
	for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
	sum += (std::count(ri.TrueContribPDGs.begin(), ri.TrueContribPDGs.end(),
	pdgs[pdg_it])
	? 0
	: 1);
	}
	return sum;
	}

	-TVector3 GetHMFSRecParticles(RecoInfo ri, int pdg) {
	- TVector3 mom(0, 0, 0);
	+TLorentzVector GetHMFSRecParticles(RecoInfo ri, int pdg) {
	+ TLorentzVector mom(0, 0, 0, 0);
	for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) {
	if ((ri.RecObjClass[p_it] == pdg) &&
	(mom.Mag() < ri.RecObjMom[p_it].Mag())) {
	- mom = ri.RecObjMom[p_it];
	+ mom.SetXYZM(ri.RecObjMom[p_it].X(), ri.RecObjMom[p_it].Y(),
	+ ri.RecObjMom[p_it].Z(),
	+ PhysConst::GetMass(ri.RecObjClass[p_it]) * 1.0E3);
	}
	}
	return mom;
	}

	template <size_t N>
	-double SumKE_RecoInfo(RecoInfo ri, int (&pdgs)[N], double mass) {
	+double SumKE_RecoInfo(RecoInfo ri, int const (&pdgs)[N], double mass) {
	double sum = 0;
	for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) {
	if (!std::count(pdgs, pdgs + N,
	ri.RecObjClass[p_it])) { // If we don't care about this
	// particle type.
	continue;
	}
	sum += sqrt(ri.RecObjMom[p_it].Mag2() + mass * mass) - mass;
	}

	return sum;
	}

	template <size_t N>
	-double SumTE_RecoInfo(RecoInfo ri, int (&pdgs)[N], double mass) {
	+double SumTE_RecoInfo(RecoInfo ri, int const (&pdgs)[N], double mass) {
	double sum = 0;
	for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) {
	if (!std::count(pdgs, pdgs + N,
	ri.RecObjClass[p_it])) { // If we don't care about this
	// particle type.
	continue;
	}
	sum += sqrt(ri.RecObjMom[p_it].Mag2() + mass * mass);
	}

	return sum;
	}

	template <size_t N>
	-double SumVisE_RecoInfo(RecoInfo ri, int (&pdgs)[N]) {
	+double SumVisE_RecoInfo(RecoInfo ri, int const (&pdgs)[N]) {
	double sum = 0;

	for (size_t p_it = 0; p_it < ri.RecVisibleEnergy.size(); ++p_it) {
	if (!std::count(pdgs, pdgs + N,
	ri.TrueContribPDGs[p_it])) { // If we don't care about this
	// particle type.
	continue;
	}
	sum += ri.RecVisibleEnergy[p_it];
	}

	return sum;
	}

	template <size_t N>
	-double SumVisE_RecoInfo_NotPdgs(RecoInfo ri, int (&pdgs)[N]) {
	+double SumVisE_RecoInfo_NotPdgs(RecoInfo ri, int const (&pdgs)[N]) {
	double sum = 0;

	for (size_t p_it = 0; p_it < ri.RecVisibleEnergy.size(); ++p_it) {
	if (std::count(pdgs, pdgs + N,
	ri.TrueContribPDGs[p_it])) { // If we know about this
	// particle type.
	continue;
	}
	sum += ri.RecVisibleEnergy[p_it];
	}

	return sum;
	}

	//********************************************************************
	void Smearceptance_Tester::FillEventVariables(FitEvent *event) {
	//********************************************************************

	- int cpipPDG[] = {211};
	- int cpimPDG[] = {-211};
	- int pi0PDG[] = {111};
	- int ProtonPDG[] = {2212};
	- int NeutronPDG[] = {2112};
	- int GammaPDG[] = {22};
	- int CLeptonPDGs[] = {11, 13, 15};
	- int ExplicitPDGs[] = {211, -211, 11, 2212, 2112, 22, 11, 13, 15, 12, 14, 16};
	+ static int const cpipPDG[] = {211};
	+ static int const cpimPDG[] = {-211};
	+ static int const pi0PDG[] = {111};
	+ static int const ProtonPDG[] = {2212};
	+ static int const NeutronPDG[] = {2112};
	+ static int const GammaPDG[] = {22};
	+ static int const CLeptonPDGs[] = {11, 13, 15};
	+ static int const ExplicitPDGs[] = {211, -211, 11, 2212, 2112, 22,
	+ 11, 13, 15, 12, 14, 16};

	RecoInfo *ri = smearceptor->Smearcept(event);

	-#ifdef DEBUG_SMEARTESTER
	-
	- FSMuon_True = TVector3(0, 0, 0);
	- FSMuon_Smeared = TVector3(0, 0, 0);
	+ HMFS_mu_true = TLorentzVector(0, 0, 0, 0);
	+ HMFS_mu_rec = TLorentzVector(0, 0, 0, 0);
	FitParticle *fsMu = event->GetHMFSMuon();
	if (fsMu) {
	- FSMuon_True = fsMu->P3();
	- FSMuon_Smeared = GetHMFSRecParticles(*ri, 13);
	+ HMFS_mu_true = fsMu->P4();
	+ HMFS_mu_rec = GetHMFSRecParticles(*ri, 13);
	+ }
	+
	+ HMFS_pip_true = TLorentzVector(0, 0, 0, 0);
	+ HMFS_pip_rec = TLorentzVector(0, 0, 0, 0);
	+ FitParticle *fsPip = event->GetHMFSPiPlus();
	+ if (fsPip) {
	+ HMFS_pip_true = fsPip->P4();
	+ HMFS_pip_rec = GetHMFSRecParticles(*ri, 211);
	+ }
	+
	+ HMFS_pim_true = TLorentzVector(0, 0, 0, 0);
	+ HMFS_pim_rec = TLorentzVector(0, 0, 0, 0);
	+ FitParticle *fsPim = event->GetHMFSPiMinus();
	+ if (fsPim) {
	+ HMFS_pim_true = fsPim->P4();
	+ HMFS_pim_rec = GetHMFSRecParticles(*ri, -211);
	+ }
	+
	+ HMFS_cpi_true = TLorentzVector(0, 0, 0, 0);
	+ HMFS_cpi_rec = TLorentzVector(0, 0, 0, 0);
	+ if (fsPip \|\| fsPim) {
	+ if (!fsPip) {
	+ HMFS_cpi_true = HMFS_pim_true;
	+ HMFS_cpi_rec = HMFS_pim_rec;
	+ } else if (!fsPim) {
	+ HMFS_cpi_true = HMFS_pip_true;
	+ HMFS_cpi_rec = HMFS_pip_rec;
	+ } else {
	+ HMFS_cpi_true =
	+ (fsPip->p2() > fsPim->p2()) ? HMFS_pip_true : HMFS_pim_true;
	+ HMFS_cpi_rec = (fsPip->p2() > fsPim->p2()) ? HMFS_pip_rec : HMFS_pim_rec;
	+ }
	+ }
	+
	+ HMFS_p_true = TLorentzVector(0, 0, 0, 0);
	+ HMFS_p_rec = TLorentzVector(0, 0, 0, 0);
	+ FitParticle *fsP = event->GetHMFSProton();
	+ if (fsP) {
	+ HMFS_p_true = fsP->P4();
	+ HMFS_p_rec = GetHMFSRecParticles(*ri, 2212);
	}
	-#endif

	TLorentzVector FourMomentumTransfer =
	(event->GetHMISAnyLeptons()->P4() - event->GetHMFSAnyLeptons()->P4());

	Omega_true = FourMomentumTransfer.E();
	Q2_true = -1 * FourMomentumTransfer.Mag2();
	Mode_true = event->Mode;

	EISLep_true = event->GetHMISAnyLeptons()->E();
	KEFSHad_cpip_true = FitUtils::SumTE_PartVect(event->GetAllFSPiPlus());
	KEFSHad_cpim_true = FitUtils::SumTE_PartVect(event->GetAllFSPiMinus());
	KEFSHad_cpi_true = KEFSHad_cpip_true + KEFSHad_cpim_true;
	TEFSHad_pi0_true = FitUtils::SumTE_PartVect(event->GetAllFSPiZero());
	KEFSHad_p_true = FitUtils::SumKE_PartVect(event->GetAllFSProton());
	KEFSHad_n_true = FitUtils::SumKE_PartVect(event->GetAllFSNeutron());
	EFSHad_true =
	KEFSHad_cpi_true + TEFSHad_pi0_true + KEFSHad_p_true + KEFSHad_n_true;
	EFSChargedEMHad_true = KEFSHad_cpi_true + TEFSHad_pi0_true + KEFSHad_p_true;

	EFSLep_true = event->GetHMFSAnyLeptons()->E();
	EFSgamma_true = FitUtils::SumTE_PartVect(event->GetAllFSPhoton());

	PDGISLep_true = event->GetHMISAnyLeptons()->PDG();
	PDGFSLep_true = event->GetHMFSAnyLeptons()->PDG();

	Nprotons_true = event->GetAllFSProton().size();
	Nneutrons_true = event->GetAllFSNeutron().size();
	Ncpiplus_true = event->GetAllFSPiPlus().size();
	Ncpiminus_true = event->GetAllFSPiMinus().size();
	Ncpi_true = Ncpiplus_true + Ncpiminus_true;
	Npi0_true = event->GetAllFSPiZero().size();

	KEFSHad_cpip_rec =
	SumKE_RecoInfo(ri, cpipPDG, PhysConst::mass_cpi PhysConst::mass_MeV);
	KEFSHad_cpim_rec =
	SumKE_RecoInfo(ri, cpimPDG, PhysConst::mass_cpi PhysConst::mass_MeV);
	KEFSHad_cpi_rec = KEFSHad_cpip_rec + KEFSHad_cpim_rec;

	TEFSHad_pi0_rec =
	SumTE_RecoInfo(ri, pi0PDG, PhysConst::mass_pi0 PhysConst::mass_MeV);
	KEFSHad_p_rec = SumKE_RecoInfo(*ri, ProtonPDG,
	PhysConst::mass_proton * PhysConst::mass_MeV);
	KEFSHad_n_rec = SumKE_RecoInfo(*ri, NeutronPDG,
	PhysConst::mass_neutron * PhysConst::mass_MeV);
	EFSHad_rec =
	KEFSHad_cpi_rec + TEFSHad_pi0_rec + KEFSHad_p_rec + KEFSHad_n_rec;

	- TVector3 FSLepMom_rec(0, 0, 0);
	+ TLorentzVector FSLepMom_rec(0, 0, 0, 0);
	if (event->GetHMFSAnyLeptons()) {
	- double massLep = event->GetHMFSAnyLeptons()->M();
	FSLepMom_rec = GetHMFSRecParticles(*ri, event->GetHMFSAnyLeptons()->PDG());
	- EFSLep_rec = (FSLepMom_rec.Mag() > 1E-5)
	- ? sqrt(FSLepMom_rec * FSLepMom_rec + massLep * massLep)
	- : 0;
	+ EFSLep_rec = FSLepMom_rec.E();
	} else {
	EFSLep_rec = 0;
	}

	EFSVis_cpip = SumVisE_RecoInfo(*ri, cpipPDG);
	EFSVis_cpim = SumVisE_RecoInfo(*ri, cpimPDG);
	EFSVis_cpi = EFSVis_cpip + EFSVis_cpim;
	EFSVis_pi0 = SumVisE_RecoInfo(*ri, pi0PDG);
	EFSVis_p = SumVisE_RecoInfo(*ri, ProtonPDG);
	EFSVis_n = SumVisE_RecoInfo(*ri, NeutronPDG);
	EFSVis_gamma = SumVisE_RecoInfo(*ri, GammaPDG);
	EFSVis_other = SumVisE_RecoInfo_NotPdgs(*ri, ExplicitPDGs);
	EFSVis = EFSVis_cpip + EFSVis_cpim + EFSVis_cpi + EFSVis_pi0 + EFSVis_p +
	EFSVis_n + EFSVis_gamma;

	FSCLep_seen = CountNPdgsSeen(*ri, CLeptonPDGs);
	Nprotons_seen = CountNPdgsSeen(*ri, ProtonPDG);
	Nneutrons_seen = CountNPdgsSeen(*ri, NeutronPDG);
	Ncpip_seen = CountNPdgsSeen(*ri, cpipPDG);
	Ncpim_seen = CountNPdgsSeen(*ri, cpimPDG);
	Ncpi_seen = Ncpip_seen + Ncpim_seen;
	Npi0_seen = CountNPdgsSeen(*ri, pi0PDG);
	Nothers_seen = CountNNotPdgsSeen(*ri, ExplicitPDGs);

	if (FSCLep_seen && (FSLepMom_rec.Mag() > 1E-8)) {
	EISLep_QE_rec =
	FitUtils::EnuQErec(FSLepMom_rec.Mag() / 1000.0, FSLepMom_rec.CosTheta(),
	34, PDGFSLep_true > 0) *
	1000.0;
	} else {
	EISLep_QE_rec = 0;
	}
	EISLep_LepHad_rec = EFSLep_rec + EFSHad_rec;
	EISLep_LepHadVis_rec = EFSLep_rec + EFSHad_rec + EFSVis;

	Nprotons_contributed = CountNPdgsContributed(*ri, ProtonPDG);
	Nneutrons_contributed = CountNPdgsContributed(*ri, NeutronPDG);
	Ncpip_contributed = CountNPdgsContributed(*ri, cpipPDG);
	Ncpim_contributed = CountNPdgsContributed(*ri, cpimPDG);
	Ncpi_contributed = Ncpip_contributed + Ncpim_contributed;
	Npi0_contributed = CountNPdgsContributed(*ri, pi0PDG);
	Ngamma_contributed = CountNPdgsContributed(*ri, GammaPDG);
	Nothers_contibuted = CountNNotPdgsContributed(*ri, ExplicitPDGs);

	Weight = event->RWWeight * event->InputWeight;
	RWWeight = event->RWWeight;
	InputWeight = event->InputWeight;
	FluxWeight = GetFluxHistogram()->GetBinContent(
	GetFluxHistogram()->FindBin(EISLep_true)) /
	GetFluxHistogram()->Integral();
	EffWeight = ri->Weight;

	xsecScaling = fScaleFactor;

	flagCCINC_true = PDGFSLep_true & 1;
	flagCC0Pi_true = (Ncpi_true + Npi0_true) == 0;

	flagCCINC_rec = FSCLep_seen && PDGFSLep_true & 1;
	flagCC0Pi_rec = ((Ncpi_seen + Npi0_seen) == 0) && flagCCINC_rec;

	// Fill the eventVariables Tree
	eventVariables->Fill();

	RecoSmear->Fill(flagCCINC_rec ? EISLep_LepHadVis_rec / 1000.0 : -1,
	EISLep_true / 1000.0, Weight);
	ETrueDistrib->Fill(EISLep_true / 1000.0, flagCCINC_true ? Weight : 0);

	ERecDistrib->Fill(EISLep_LepHadVis_rec / 1000.0, flagCCINC_rec ? Weight : 0);

	return;
	};

	//********************************************************************
	void Smearceptance_Tester::Write(std::string drawOpt) {
	//********************************************************************

	// First save the TTree
	eventVariables->Write();

	// Save Flux and Event Histograms too
	GetInput()->GetFluxHistogram()->Write();
	GetInput()->GetEventHistogram()->Write();

	TH2D *SmearMatrix_ev =
	static_cast<TH2D *>(RecoSmear->Clone("ELepHadVis_Smear_ev"));

	for (Int_t trueAxis_it = 1;
	trueAxis_it < RecoSmear->GetYaxis()->GetNbins() + 1; ++trueAxis_it) {
	double NEISLep = ETrueDistrib->GetBinContent(trueAxis_it);

	for (Int_t recoAxis_it = 1;
	recoAxis_it < RecoSmear->GetXaxis()->GetNbins() + 1; ++recoAxis_it) {
	if (NEISLep > std::numeric_limits<double>::epsilon()) {
	SmearMatrix_ev->SetBinContent(
	recoAxis_it, trueAxis_it,
	SmearMatrix_ev->GetBinContent(recoAxis_it, trueAxis_it) / NEISLep);
	}
	}
	}

	ETrueDistrib->Write();
	ERecDistrib->Write();

	RecoSmear->Write();

	SmearMatrix_ev->Write();

	TH2D *ResponseMatrix_ev =
	SmearceptanceUtils::SVDGetInverse(SmearMatrix_ev, SVDTruncation);
	ResponseMatrix_ev = SmearceptanceUtils::SwapXYTH2D(ResponseMatrix_ev);
	ResponseMatrix_ev->SetName("ResponseMatrix_ev");
	ResponseMatrix_ev->Write();

	#ifdef DEBUG_SMEARTESTER

	TMatrixD SmearMatrix_ev_md = SmearceptanceUtils::GetMatrix(
	SmearceptanceUtils::SwapXYTH2D(SmearMatrix_ev));

	TH1D SmearedEvt = static_cast<TH1D >(ERecDistrib->Clone());
	SmearedEvt->SetNameTitle("SmearedEvt", ";Rec E_{#nu}; count");

	SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
	ETrueDistrib, SmearedEvt, SmearMatrix_ev_md, 5000, false);

	SmearedEvt->Write();

	SmearedEvt->Scale(1, "width");
	SmearedEvt->SetName("SmearedEvt_bw");
	SmearedEvt->Write();

	#endif

	TMatrixD ResponseMatrix_evt_md =
	SmearceptanceUtils::GetMatrix(ResponseMatrix_ev);

	TH1D Unfolded_enu_obs = static_cast<TH1D >(ETrueDistrib->Clone());
	Unfolded_enu_obs->SetNameTitle("UnfoldedENu_evt", ";True E_{#nu};count");

	SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
	ERecDistrib, Unfolded_enu_obs, ResponseMatrix_evt_md, 5000, false);

	Unfolded_enu_obs->Write();

	Unfolded_enu_obs->Scale(1, "width");
	Unfolded_enu_obs->SetName("UnfoldedENu_evt_bw");
	Unfolded_enu_obs->Write();

	ETrueDistrib->Scale(1, "width");
	ETrueDistrib->SetName("ELep_rate_bw");
	ETrueDistrib->Write();

	ERecDistrib->Scale(1, "width");
	ERecDistrib->SetName("ELepRec_rate_bw");
	ERecDistrib->Write();
	}

	// -------------------------------------------------------------------
	// Purely MC Plot
	// Following functions are just overrides to handle this
	// -------------------------------------------------------------------
	//********************************************************************
	/// Everything is classed as signal...
	bool Smearceptance_Tester::isSignal(FitEvent *event) {
	//********************************************************************
	(void)event;
	return true;
	};

	//********************************************************************
	void Smearceptance_Tester::ScaleEvents() {
	//********************************************************************
	// Saving everything to a TTree so no scaling required
	return;
	}

	//********************************************************************
	void Smearceptance_Tester::ApplyNormScale(float norm) {
	//********************************************************************

	// Saving everything to a TTree so no scaling required
	this->fCurrentNorm = norm;
	return;
	}

	//********************************************************************
	void Smearceptance_Tester::FillHistograms() {
	//********************************************************************
	// No Histograms need filling........
	return;
	}

	//********************************************************************
	void Smearceptance_Tester::ResetAll() {
	//********************************************************************
	eventVariables->Reset();
	return;
	}

	//********************************************************************
	float Smearceptance_Tester::GetChi2() {
	//********************************************************************
	// No Likelihood to test, purely MC
	return 0.0;
	}
	diff --git a/src/MCStudies/Smearceptance_Tester.h b/src/MCStudies/Smearceptance_Tester.h
	index ef5a87f..f0abbaa 100644
	--- a/src/MCStudies/Smearceptance_Tester.h
	+++ b/src/MCStudies/Smearceptance_Tester.h
	@@ -1,171 +1,176 @@
	// 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 <http://www.gnu.org/licenses/>.
	*******************************************************************************/

	#ifndef Smearceptance_Tester_H_SEEN
	#define Smearceptance_Tester_H_SEEN

	#include "Measurement1D.h"

	#include "ISmearcepter.h"

	#ifdef __PROB3PP_ENABLED__
	#include "OscWeightEngine.h"
	#endif

	-#define DEBUG_SMEARTESTER
	-
	//********************************************************************
	class Smearceptance_Tester : public Measurement1D {
	//********************************************************************

	public:
	Smearceptance_Tester(std::string name, std::string inputfile, FitWeight *rw,
	std::string type, std::string fakeDataFile);
	virtual ~Smearceptance_Tester(){};

	//! Grab info from event
	void FillEventVariables(FitEvent *event);

	//! Fill Custom Histograms
	void FillHistograms();

	//! ResetAll
	void ResetAll();

	//! Scale
	void ScaleEvents();

	//! Norm
	void ApplyNormScale(float norm);

	//! Define this samples signal
	bool isSignal(FitEvent *nvect);

	//! Write Files
	void Write(std::string drawOpt);

	//! Get Chi2
	float GetChi2();

	void AddEventVariablesToTree();

	private:
	ISmearcepter *smearceptor;

	TTree *eventVariables;

	float Omega_true;
	float Q2_true;
	int Mode_true;

	float EISLep_true;

	+ TLorentzVector HMFS_mu_true;
	+ TLorentzVector HMFS_pip_true;
	+ TLorentzVector HMFS_pim_true;
	+ TLorentzVector HMFS_cpi_true;
	+ TLorentzVector HMFS_p_true;
	+
	float KEFSHad_cpip_true;
	float KEFSHad_cpim_true;
	float KEFSHad_cpi_true;
	float TEFSHad_pi0_true;
	float KEFSHad_p_true;
	float KEFSHad_n_true;

	float EFSHad_true;
	float EFSChargedEMHad_true;

	float EFSLep_true;
	float EFSgamma_true;

	int PDGISLep_true;
	int PDGFSLep_true;

	int Nprotons_true;
	int Nneutrons_true;
	int Ncpiplus_true;
	int Ncpiminus_true;
	int Ncpi_true;
	int Npi0_true;

	+ TLorentzVector HMFS_mu_rec;
	+ TLorentzVector HMFS_pip_rec;
	+ TLorentzVector HMFS_pim_rec;
	+ TLorentzVector HMFS_cpi_rec;
	+ TLorentzVector HMFS_p_rec;
	+
	float KEFSHad_cpip_rec;
	float KEFSHad_cpim_rec;
	float KEFSHad_cpi_rec;
	float TEFSHad_pi0_rec;
	float KEFSHad_p_rec;
	float KEFSHad_n_rec;

	float EFSHad_rec;
	float EFSLep_rec;

	float EFSVis_cpip;
	float EFSVis_cpim;
	float EFSVis_cpi;
	float EFSVis_pi0;
	float EFSVis_p;
	float EFSVis_n;
	float EFSVis_gamma;
	float EFSVis_other;
	float EFSVis;

	int FSCLep_seen;
	int Nprotons_seen;
	int Nneutrons_seen;
	int Ncpip_seen;
	int Ncpim_seen;
	int Ncpi_seen;
	int Npi0_seen;
	int Nothers_seen;

	float EISLep_QE_rec;
	float EISLep_LepHad_rec;
	float EISLep_LepHadVis_rec;

	int Nprotons_contributed;
	int Nneutrons_contributed;
	int Ncpip_contributed;
	int Ncpim_contributed;
	int Ncpi_contributed;
	int Npi0_contributed;
	int Ngamma_contributed;
	int Nothers_contibuted;

	float Weight;
	float RWWeight;
	float InputWeight;
	float FluxWeight;
	float EffWeight;

	float xsecScaling;

	bool flagCCINC_true;
	bool flagCC0Pi_true;

	bool flagCCINC_rec;
	bool flagCC0Pi_rec;

	-#ifdef DEBUG_SMEARTESTER
	- TVector3 FSMuon_True;
	- TVector3 FSMuon_Smeared;
	-#endif

	int SVDTruncation;

	TH2D *RecoSmear;
	TH1D *ETrueDistrib;
	TH1D *ERecDistrib;
	-
	};

	#endif
	diff --git a/src/Smearceptance/EfficiencyApplicator.cxx b/src/Smearceptance/EfficiencyApplicator.cxx
	index 5c43594..b046cb5 100644
	--- a/src/Smearceptance/EfficiencyApplicator.cxx
	+++ b/src/Smearceptance/EfficiencyApplicator.cxx
	@@ -1,325 +1,388 @@
	// 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 <http://www.gnu.org/licenses/>.
	*******************************************************************************/

	#include "EfficiencyApplicator.h"

	#include "TEfficiency.h"
	#include "TH2.h"
	#include "TH3.h"

	//#define DEBUG_EFFAPP

	EfficiencyApplicator::DependVar GetVarType(std::string const &axisvar) {
	if (axisvar == "kMomentum") {
	return EfficiencyApplicator::kMomentum;
	} else if (axisvar == "kKE") {
	return EfficiencyApplicator::kKE;
	} else if (axisvar == "kTheta") {
	return EfficiencyApplicator::kTheta;
	} else if (axisvar == "kCosTheta") {
	return EfficiencyApplicator::kCosTheta;
	} else if (axisvar == "kPhi") {
	return EfficiencyApplicator::kPhi;
	}
	return EfficiencyApplicator::kNoAxis;
	}

	TH1 GetEffHist(TFile inputFile, std::string const &HistName) {
	TH1 hist = dynamic_cast<TH1 >(inputFile->Get(HistName.c_str()));
	if (hist) {
	return hist;
	}

	TEfficiency *effHist =
	dynamic_cast<TEfficiency *>(inputFile->Get(HistName.c_str()));

	if (effHist) {
	TH1D numer = dynamic_cast<TH1D >(effHist->GetCopyPassedHisto());
	TH1D denom = dynamic_cast<TH1D >(effHist->GetCopyTotalHisto());

	if (numer && denom) {
	numer->Divide(denom);

	denom->SetDirectory(NULL);
	delete denom;

	// Gonna be a memory leak, but I'll get over it
	numer->SetDirectory(NULL);
	return numer;
	}
	ERROR(FTL, "TEfficiency internal histograms were not TH1Ds.");
	}

	THROW("Couldn't get appropriate efficiency object named "
	<< HistName << " from input file " << inputFile->GetName());
	}

	/// Reads particle efficiency nodes
	///
	/// Nodes look like:
	/// <EfficiencyApplicator Name="D00N_ND_LAr">
	/// <EfficiencyCurve PDG="211,-211" InputFile="effs.root"
	/// HistName="cpion_eff_mom_ctheta" NDims="2" XAxis="kMomentum"
	/// YAxis="kCosTheta" Interpolate="false" />
	/// <!-- Can also contain ThresholdAccepter elements as below-->
	/// <RecoThreshold PDG="2212" RecoThresholdAbsCosTheta_Min="0" />
	/// <VisThreshold PDG="2212" VisThresholdKE_MeV="10" Contrib="K" />
	/// </EfficiencyApplicator>
	void EfficiencyApplicator::SpecifcSetup(nuiskey &nk) {
	rand.~TRandom3();
	new (&rand) TRandom3();

	std::vector<nuiskey> effDescriptors =
	nk.GetListOfChildNodes("EfficiencyCurve");

	for (size_t t_it = 0; t_it < effDescriptors.size(); ++t_it) {
	std::string inputFileName = effDescriptors[t_it].GetS("InputFile");
	std::string HistName = effDescriptors[t_it].GetS("HistName");

	TFile inputFile(inputFileName.c_str());
	if (!inputFile.IsOpen()) {
	THROW("Couldn't open specified input root file: " << inputFileName);
	}

	TH1 *inpHist = GetEffHist(&inputFile, HistName);
	if (!inpHist) {
	THROW("Couldn't get TH1D named: " << HistName << " from input root file: "
	<< inputFileName);
	}

	int NDims = effDescriptors[t_it].GetI("NDims");

	if (NDims < 1 \|\| NDims > 3) {
	THROW("Read NDims attribute as: " << NDims << ", efficiency curve can "
	"currently have between 1 "
	"and 3 dimensions.");
	}

	EfficiencyApplicator::DependVar XVar =
	GetVarType(effDescriptors[t_it].GetS("XAxis"));
	double XAxisScale = effDescriptors[t_it].Has("XAxisScaleToInternal")
	? effDescriptors[t_it].GetD("XAxisScaleToInternal")
	: 1;
	EfficiencyApplicator::DependVar YVar =
	NDims > 1 ? GetVarType(effDescriptors[t_it].GetS("YAxis"))
	: EfficiencyApplicator::kNoAxis;
	double YAxisScale = effDescriptors[t_it].Has("YAxisScaleToInternal")
	? effDescriptors[t_it].GetD("YAxisScaleToInternal")
	: 1;
	EfficiencyApplicator::DependVar ZVar =
	NDims > 2 ? GetVarType(effDescriptors[t_it].GetS("ZAxis"))
	: EfficiencyApplicator::kNoAxis;
	double ZAxisScale = effDescriptors[t_it].Has("ZAxisScaleToInternal")
	? effDescriptors[t_it].GetD("ZAxisScaleToInternal")
	: 1;

	bool Interpolate = effDescriptors[t_it].Has("Interpolate") &&
	effDescriptors[t_it].GetI("Interpolate");

	// bool ApplyAsWeight = effDescriptors[t_it].Has("ApplyAsWeight") &&
	// effDescriptors[t_it].GetI("ApplyAsWeight");

	std::string pdgs_s = effDescriptors[t_it].GetS("PDG");
	std::vector<int> pdgs_i = GeneralUtils::ParseToInt(pdgs_s, ",");
	for (size_t pdg_it = 0; pdg_it < pdgs_i.size(); ++pdg_it) {
	if (Efficiencies.count(pdgs_i[pdg_it])) {
	ERROR(WRN, "Smearceptor " << ElementName << ":" << InstanceName
	<< " already has a efficiency for PDG: "
	<< pdgs_i[pdg_it]);
	}

	EffMap em;
	em.EffCurve = static_cast<TH1 *>(inpHist->Clone());
	em.EffCurve->SetDirectory(NULL);
	em.Interpolate = Interpolate;
	// em.ApplyAsWeight = ApplyAsWeight;
	em.NDims = NDims;
	em.DependVars[0] = XVar;
	em.DependVars[1] = YVar;
	em.DependVars[2] = ZVar;
	em.AxisScales[0] = XAxisScale;
	em.AxisScales[1] = YAxisScale;
	em.AxisScales[2] = ZAxisScale;

	Efficiencies[pdgs_i[pdg_it]] = em;

	QLOG(SAM,
	"Added reconstruction efficiency curve for PDG: " << pdgs_i[pdg_it]);
	}
	}

	SlaveTA.Setup(nk);
	}

	RecoInfo EfficiencyApplicator::Smearcept(FitEvent fe) {
	RecoInfo *ri = new RecoInfo();

	for (size_t p_it = 0; p_it < fe->NParticles(); ++p_it) {
	FitParticle *fp = fe->GetParticle(p_it);
	#ifdef DEBUG_EFFAPP
	std::cout << std::endl;
	std::cout << "[" << p_it << "]: " << fp->PDG() << ", " << fp->Status()
	<< ", " << fp->E() << " -- KE:" << fp->KE()
	<< " Mom: " << fp->P3().Mag() << std::flush;
	#endif

	if (fp->Status() != kFinalState) {
	#ifdef DEBUG_EFFAPP
	std::cout << " -- Not final state." << std::flush;
	#endif
	continue;
	}

	if (!Efficiencies.count(fp->PDG())) {
	SlaveTA.SmearceptOneParticle(ri, fp
	#ifdef DEBUG_THRESACCEPT
	,
	p_it
	#endif
	);
	continue;
	}

	EffMap &em = Efficiencies[fp->PDG()];

	double kineProps[3];
	for (Int_t dim_it = 0; dim_it < em.NDims; ++dim_it) {
	switch (em.DependVars[dim_it]) {
	case kMomentum: {
	kineProps[dim_it] = fp->P3().Mag();
	break;
	}
	case kKE: {
	kineProps[dim_it] = fp->KE();
	break;
	}
	case kTheta: {
	kineProps[dim_it] = fp->P3().Theta();
	break;
	}
	case kCosTheta: {
	kineProps[dim_it] = fp->P3().CosTheta();
	break;
	}
	case kPhi: {
	kineProps[dim_it] = fp->P3().Phi();
	break;
	}
	default: { THROW("Trying to find particle value for a kNoAxis."); }
	}
	kineProps[dim_it] /= em.AxisScales[dim_it];
	}

	double effProb = 0;

	switch (em.NDims) {
	case 1: {
	TH1 *hist = em.EffCurve;
	if (em.Interpolate &&
	(hist->GetXaxis()->GetBinCenter(1) < kineProps[0]) &&
	(hist->GetXaxis()->GetBinCenter(hist->GetXaxis()->GetNbins()) >
	kineProps[0])) {
	effProb = hist->Interpolate(kineProps[0]);
	} else {
	Int_t xbin = hist->GetXaxis()->FindFixBin(kineProps[0]);
	+
	+ if (!xbin \|\| ((hist->GetXaxis()->GetNbins() + 1) == xbin)) {
	+ ERROR(WRN, "Tried to apply effiency but XBin: "
	+ << xbin << " is outside range (/"
	+ << hist->GetXaxis()->GetNbins() << "): Prop "
	+ << kineProps[0] << ", ["
	+ << hist->GetXaxis()->GetBinLowEdge(1) << " -- "
	+ << hist->GetXaxis()->GetBinUpEdge(
	+ hist->GetXaxis()->GetNbins()));
	+ }
	+
	effProb = hist->GetBinContent(xbin);
	}
	break;
	}
	case 2: {
	TH2 hist = static_cast<TH2 >(em.EffCurve);

	if (em.Interpolate &&
	(hist->GetXaxis()->GetBinCenter(1) < kineProps[0]) &&
	(hist->GetXaxis()->GetBinCenter(hist->GetXaxis()->GetNbins()) >
	kineProps[0]) &&
	(hist->GetYaxis()->GetBinCenter(1) < kineProps[1]) &&
	(hist->GetYaxis()->GetBinCenter(hist->GetYaxis()->GetNbins()) >
	kineProps[1])) {
	effProb = hist->Interpolate(kineProps[0], kineProps[1]);
	} else {
	Int_t xbin = hist->GetXaxis()->FindFixBin(kineProps[0]);
	Int_t ybin = hist->GetYaxis()->FindFixBin(kineProps[1]);
	+
	+ if (!xbin \|\| ((hist->GetXaxis()->GetNbins() + 1) == xbin)) {
	+ ERROR(WRN, "Tried to apply effiency but XBin: "
	+ << xbin << " is outside range (/"
	+ << hist->GetXaxis()->GetNbins() << "): Prop "
	+ << kineProps[0] << ", ["
	+ << hist->GetXaxis()->GetBinLowEdge(1) << " -- "
	+ << hist->GetXaxis()->GetBinUpEdge(
	+ hist->GetXaxis()->GetNbins()));
	+ }
	+
	+ if (!ybin \|\| ((hist->GetYaxis()->GetNbins() + 1) == ybin)) {
	+ ERROR(WRN, "Tried to apply effiency but XBin: "
	+ << ybin << " is outside range (/"
	+ << hist->GetYaxis()->GetNbins() << "): Prop "
	+ << kineProps[0] << ", ["
	+ << hist->GetYaxis()->GetBinLowEdge(1) << " -- "
	+ << hist->GetYaxis()->GetBinUpEdge(
	+ hist->GetYaxis()->GetNbins()));
	+ }
	+
	effProb = hist->GetBinContent(xbin, ybin);

	#ifdef DEBUG_EFFAPP
	std::cout << "\t\t: XProp: " << kineProps[0]
	<< ", YProp: " << kineProps[1] << " x/y bins: " << xbin
	<< "/" << ybin << ". Prop ? " << effProb << std::endl;
	#endif
	}
	break;
	}
	case 3: {
	TH3 hist = static_cast<TH3 >(em.EffCurve);

	if (em.Interpolate &&
	(hist->GetXaxis()->GetBinCenter(1) < kineProps[0]) &&
	(hist->GetXaxis()->GetBinCenter(hist->GetXaxis()->GetNbins()) >
	kineProps[0]) &&
	(hist->GetYaxis()->GetBinCenter(1) < kineProps[1]) &&
	(hist->GetYaxis()->GetBinCenter(hist->GetYaxis()->GetNbins()) >
	kineProps[2]) &&
	(hist->GetZaxis()->GetBinCenter(hist->GetZaxis()->GetNbins()) >
	kineProps[2])) {
	effProb = hist->Interpolate(kineProps[0], kineProps[1], kineProps[2]);
	} else {
	Int_t xbin = hist->GetXaxis()->FindFixBin(kineProps[0]);
	Int_t ybin = hist->GetYaxis()->FindFixBin(kineProps[1]);
	Int_t zbin = hist->GetZaxis()->FindFixBin(kineProps[2]);
	+
	+ if (!xbin \|\| ((hist->GetXaxis()->GetNbins() + 1) == xbin)) {
	+ ERROR(WRN, "Tried to apply effiency but XBin: "
	+ << xbin << " is outside range (/"
	+ << hist->GetXaxis()->GetNbins() << "): Prop "
	+ << kineProps[0] << ", ["
	+ << hist->GetXaxis()->GetBinLowEdge(1) << " -- "
	+ << hist->GetXaxis()->GetBinUpEdge(
	+ hist->GetXaxis()->GetNbins()));
	+ }
	+
	+ if (!ybin \|\| ((hist->GetYaxis()->GetNbins() + 1) == ybin)) {
	+ ERROR(WRN, "Tried to apply effiency but XBin: "
	+ << ybin << " is outside range (/"
	+ << hist->GetYaxis()->GetNbins() << "): Prop "
	+ << kineProps[0] << ", ["
	+ << hist->GetYaxis()->GetBinLowEdge(1) << " -- "
	+ << hist->GetYaxis()->GetBinUpEdge(
	+ hist->GetYaxis()->GetNbins()));
	+ }
	+
	+ if (!zbin \|\| ((hist->GetZaxis()->GetNbins() + 1) == zbin)) {
	+ ERROR(WRN, "Tried to apply effiency but ZBin: "
	+ << zbin << " is outside range (/"
	+ << hist->GetZaxis()->GetNbins() << "): Prop "
	+ << kineProps[0] << ", ["
	+ << hist->GetZaxis()->GetBinLowEdge(1) << " -- "
	+ << hist->GetZaxis()->GetBinUpEdge(
	+ hist->GetZaxis()->GetNbins()));
	+ }
	+
	effProb = hist->GetBinContent(xbin, ybin, zbin);
	}
	break;
	}
	}

	bool accepted = (rand.Uniform() < effProb);

	if (accepted) {
	#ifdef DEBUG_EFFAPP
	std::cout << " -- Reconstructed with probability: " << effProb
	<< std::flush;
	#endif
	ri->RecObjMom.push_back(fp->P3());
	ri->RecObjClass.push_back(fp->PDG());

	continue;
	}

	#ifdef DEBUG_EFFAPP
	std::cout << " -- Rejected with probability: " << effProb << std::flush;
	#endif

	SlaveTA.SmearceptOneParticle(ri, fp
	#ifdef DEBUG_THRESACCEPT
	,
	p_it
	#endif
	);
	}
	#ifdef DEBUG_EFFAPP
	std::cout << std::endl;
	#endif

	#ifdef DEBUG_EFFAPP
	std::cout << "Reconstructed " << ri->RecObjMom.size() << " particles. "
	<< std::endl;
	#endif
	return ri;
	}
	diff --git a/src/Smearceptance/GaussianSmearer.cxx b/src/Smearceptance/GaussianSmearer.cxx
	index 61fecf3..5f57177 100644
	--- a/src/Smearceptance/GaussianSmearer.cxx
	+++ b/src/Smearceptance/GaussianSmearer.cxx
	@@ -1,477 +1,529 @@
	// 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 <http://www.gnu.org/licenses/>.
	*******************************************************************************/

	#include "GaussianSmearer.h"

	-// #define DEBUG_THRESACCEPT
	-
	namespace {
	GaussianSmearer::GSmearType GetVarType(std::string const &type) {
	if (type == "Absolute") {
	return GaussianSmearer::kAbsolute;
	} else if (type == "Fractional") {
	return GaussianSmearer::kFractional;
	} else if (type == "Function") {
	return GaussianSmearer::kFunction;
	}
	return GaussianSmearer::kNoType;
	}

	GaussianSmearer::DependVar GetKineType(std::string const &axisvar) {
	if (axisvar == "Momentum") {
	return GaussianSmearer::kMomentum;
	} else if (axisvar == "KE") {
	return GaussianSmearer::kKE;
	} else if (axisvar == "TEVis") {
	return GaussianSmearer::kTEVis;
	} else if (axisvar == "KEVis") {
	return GaussianSmearer::kKEVis;
	} else if (axisvar == "CosTheta") {
	return GaussianSmearer::kCosTheta;
	} else if (axisvar == "Theta") {
	return GaussianSmearer::kTheta;
	}
	THROW("Failed to parse smear type from \"" << axisvar << "\"");
	}

	std::string GetKineTypeName(GaussianSmearer::DependVar dv) {
	switch (dv) {
	case GaussianSmearer::kMomentum: {
	return "Momentum";
	}
	case GaussianSmearer::kKE: {
	return "KE";
	}
	case GaussianSmearer::kTEVis: {
	return "TEVis";
	}
	case GaussianSmearer::kKEVis: {
	return "KEVis";
	}
	case GaussianSmearer::kCosTheta: {
	return "CosTheta";
	}
	case GaussianSmearer::kTheta: {
	return "Theta";
	}
	default: { THROW("NO VAR!"); }
	}
	}
	}

	/// Nodes look like:
	/// Function attribute is given to a TF1, where any "V"s are replaced with the
	-/// selected kinematic property on an event-by-event basis. e.g Function="V +
	-/// gaus(0.2*V),..." should give the same result as Type="kFractional"
	-/// Width="0.2".
	+/// selected kinematic property on an event-by-event basis. e.g Function="{V} +
	+/// gaus({P1}),..." with P1="0.2", should give the same result as
	+/// Type="kFractional" and Width="0.2".
	/// <GaussianSmearer Name="D00N_ND_LAr">
	-/// <Smear PDG="211" Type="[kAbsolute,kFractional,kFunction]"
	-/// Width="2" Kinematics="[KE\|Momentum\|VisKE\|VisTE]" Function="V +
	-/// gaus(1/V),<lowlim>,<highlim>" />
	+/// <Smear PDG="211" Type="[Absolute\|Fractional\|Function]"
	+/// Kinematics="[KE\|Momentum\|VisKE\|VisTE\|CosTheta\|Theta]" (Width="2")
	+/// (Function="V +
	+/// gaus(1/{V}),<lowlim>,<highlim>") (AllowNeg="0") />
	/// </GaussianSmearer>
	void GaussianSmearer::SpecifcSetup(nuiskey &nk) {
	rand.~TRandom3();
	new (&rand) TRandom3();

	std::vector<nuiskey> smearDescriptors = nk.GetListOfChildNodes("Smear");

	for (size_t t_it = 0; t_it < smearDescriptors.size(); ++t_it) {
	std::string pdgs_s = smearDescriptors[t_it].GetS("PDG");
	std::vector<int> pdgs_i = GeneralUtils::ParseToInt(pdgs_s, ",");

	double Width = smearDescriptors[t_it].Has("Width")
	? smearDescriptors[t_it].GetD("Width")
	: 0xdeadbeef;

	GaussianSmearer::GSmearType Type =
	GetVarType(smearDescriptors[t_it].GetS("Type"));
	GaussianSmearer::DependVar Kinematics =
	GetKineType(smearDescriptors[t_it].GetS("Kinematics"));
	bool IsVisSmear = (Kinematics == GaussianSmearer::kKEVis) \|\|
	(Kinematics == GaussianSmearer::kTEVis);

	TF1 *sf = NULL;

	if (Type == GaussianSmearer::kFunction) {
	std::string funcDescriptor = smearDescriptors[t_it].Has("Function")
	? smearDescriptors[t_it].GetS("Function")
	: "";
	if (funcDescriptor.size()) {
	std::vector<std::string> funcP =
	- GeneralUtils::ParseToStr(funcDescriptor, ",");
	+ GeneralUtils::ParseToStr(funcDescriptor, "$");
	if (funcP.size() != 3) {
	THROW(
	"Expected Function attribute to contain 3 comma separated "
	- "entries. e.g. Function=\"1/x,<low lim>,<high lim>\". ");
	+ "entries. e.g. Function=\"1/{V}$<low lim>$<high lim>\". ");
	+ }
	+ bool FoundParam;
	+ int pCtr = 1;
	+ std::map<std::string, std::string> PVals;
	+ do {
	+ std::stringstream pv_str("");
	+ pv_str << "P" << pCtr++;
	+ if (smearDescriptors[t_it].Has(pv_str.str())) {
	+ PVals.insert(
	+ std::make_pair(std::string("{") + pv_str.str() + "}",
	+ smearDescriptors[t_it].GetS(pv_str.str())));
	+ FoundParam = true;
	+ } else {
	+ FoundParam = false;
	+ }
	+ } while (FoundParam);
	+
	+ for (std::map<std::string, std::string>::iterator v_it = PVals.begin();
	+ v_it != PVals.end(); ++v_it) {
	+ funcP[0] =
	+ GeneralUtils::ReplaceAll(funcP[0], v_it->first, v_it->second);
	}
	- funcP[0] = GeneralUtils::ReplaceAll(funcP[0], "V", "[0]");
	- QLOG(FIT, "Added smearing func: " << funcP[0]);
	+
	+ funcP[0] = GeneralUtils::ReplaceAll(funcP[0], "{V}", "[0]");
	+ QLOG(FIT, "Added smearing func: "
	+ << funcP[0] << ", [" << GeneralUtils::StrToDbl(funcP[1])
	+ << " -- " << GeneralUtils::StrToDbl(funcP[2]) << "].");
	sf = new TF1("smear_dummy", funcP[0].c_str(),
	GeneralUtils::StrToDbl(funcP[1]),
	GeneralUtils::StrToDbl(funcP[2]));
	} else {
	THROW(
	"Expected Function attribute with 3 comma separated "
	"entries. e.g. Function=\"1/x,<low lim>,<high lim>\". ");
	}
	}

	for (size_t pdg_it = 0; pdg_it < pdgs_i.size(); ++pdg_it) {
	if (IsVisSmear && VisGausSmears.count(pdgs_i[pdg_it])) {
	ERROR(WRN,
	"Smearceptor "
	<< ElementName << ":" << InstanceName
	<< " already has a Visible Energy smearing function for PDG: "
	<< pdgs_i[pdg_it]);
	}

	GSmear gs;

	gs.type = Type;
	gs.smearVar = Kinematics;
	gs.width = Width;
	gs.func = sf ? static_cast<TF1 *>(sf->Clone()) : NULL;
	if (sf) {
	std::stringstream ss("");
	ss << "GausSmear"
	<< "_PDG" << pdgs_i[pdg_it];
	gs.func->SetName(ss.str().c_str());
	}

	if (IsVisSmear) {
	VisGausSmears[pdgs_i[pdg_it]] = gs;
	} else {
	TrackedGausSmears[pdgs_i[pdg_it]].push_back(gs);
	}

	QLOG(SAM, "Added gaussian "
	<< GetKineTypeName(gs.smearVar)
	<< " smearing function for PDG: " << pdgs_i[pdg_it]);
	}
	delete sf;
	}
	}

	-RecoInfo GaussianSmearer::Smearcept(FitEvent fe) {
	- RecoInfo *ri = new RecoInfo();
	-
	- for (size_t p_it = 0; p_it < fe->NParticles(); ++p_it) {
	- FitParticle *fp = fe->GetParticle(p_it);
	-#ifdef DEBUG_THRESACCEPT
	- std::cout << std::endl;
	- std::cout << "[" << p_it << "]: " << fp->PDG() << ", " << fp->Status()
	- << ", " << fp->E() << " -- KE:" << fp->KE()
	- << " Mom: " << fp->P3().Mag() << std::flush;
	+void GaussianSmearer::SmearceptOneParticle(RecoInfo ri, FitParticle fp
	+#ifdef DEBUG_GAUSSSMEAR
	+ ,
	+ size_t p_it
	#endif
	-
	- if (fp->Status() != kFinalState) {
	-#ifdef DEBUG_THRESACCEPT
	- std::cout << " -- Not final state." << std::flush;
	-#endif
	- continue;
	- }
	-
	- if ((TrackedGausSmears.count(fp->PDG()) + VisGausSmears.count(fp->PDG())) ==
	- 0) {
	-#ifdef DEBUG_THRESACCEPT
	- std::cout << " -- Undetectable." << std::flush;
	-#endif
	- continue;
	- }
	-
	- if (TrackedGausSmears.count(fp->PDG())) {
	- TVector3 ThreeMom = fp->P3();
	- for (size_t sm_it = 0; sm_it < TrackedGausSmears[fp->PDG()].size();
	- ++sm_it) {
	- GSmear &sm = TrackedGausSmears[fp->PDG()][sm_it];
	-
	- double kineProp = 0;
	-
	- switch (sm.smearVar) {
	- case GaussianSmearer::kMomentum: {
	- kineProp = fp->P3().Mag();
	- break;
	- }
	- case GaussianSmearer::kKE: {
	- kineProp = fp->KE();
	- break;
	- }
	- case GaussianSmearer::kCosTheta: {
	- kineProp = fp->P3().CosTheta();
	- break;
	- }
	- case GaussianSmearer::kTheta: {
	- kineProp = fp->P3().Theta();
	- break;
	- }
	- default: { THROW("Trying to find particle value for a kNoVar."); }
	- }
	-
	- double Smeared;
	- bool ok = false;
	- while (!ok) {
	- if (sm.type == GaussianSmearer::kFunction) {
	- sm.func->SetParameter(0, kineProp);
	- Smeared = sm.func->GetRandom();
	- } else {
	- double sThrow =
	- rand.Gaus(0, sm.width * ((sm.type == GaussianSmearer::kAbsolute)
	- ? 1
	- : kineProp));
	- Smeared = kineProp + sThrow;
	- }
	- switch (
	- sm.smearVar) { // Different kinematics have different truncation.
	- case GaussianSmearer::kMomentum:
	- case GaussianSmearer::kKE: {
	- ok = (Smeared > 0);
	- break;
	- }
	- case GaussianSmearer::kCosTheta: {
	- ok = ((Smeared >= -1) && (Smeared <= 1));
	- break;
	- }
	- case GaussianSmearer::kTheta: {
	- ok = true;
	- break;
	- }
	-
	- default: { THROW("SHOULDN'T BE HERE."); }
	- }
	- }
	-
	- switch (sm.smearVar) {
	- case GaussianSmearer::kMomentum: {
	- ThreeMom = (ThreeMom.Unit() * Smeared);
	- break;
	- }
	- case GaussianSmearer::kKE: {
	- double mass = fp->P4().M();
	- double TE = mass + Smeared;
	- double magP = sqrt(TE * TE - mass * mass);
	- ThreeMom = (ThreeMom.Unit() * magP);
	- break;
	- }
	- case GaussianSmearer::kCosTheta: {
	- ThreeMom.SetTheta(acos(Smeared));
	- break;
	- }
	- case GaussianSmearer::kTheta: {
	- ThreeMom.SetTheta(Smeared);
	- break;
	- }
	- default: {}
	- }
	- }
	-#ifdef DEBUG_THRESACCEPT
	- std::cout << " -- momentum reconstructed as Mom: "
	- << ri->RecObjMom.back().Mag()
	- << ", CT: " << ri->RecObjMom.back().CosTheta() << " from "
	- << ThreeMom.Mag() << ", " << fp->P3().CosTheta() << " true."
	- << std::endl;
	+ ) {
	+#ifdef DEBUG_GAUSSSMEAR
	+ std::cout << std::endl;
	+ std::cout << "[" << p_it << "]: " << fp->PDG() << ", " << fp->Status() << ", "
	+ << fp->E() << " -- KE:" << fp->KE() << " Mom: " << fp->P3().Mag()
	+ << std::flush;
	#endif
	- ri->RecObjMom.push_back(ThreeMom);
	- ri->RecObjClass.push_back(fp->PDG());
	- } else { // Smear to EVis

	- GSmear &sm = VisGausSmears[fp->PDG()];
	-
	- double kineProp = 0;
	-
	- switch (sm.smearVar) {
	- case GaussianSmearer::kKEVis: {
	- kineProp = fp->KE();
	- break;
	- }
	- case GaussianSmearer::kTEVis: {
	- kineProp = fp->E();
	- break;
	- }
	- default: { THROW("Trying to find particle value for a kNoVar."); }
	- }
	-
	- double Smeared;
	- if (sm.type == GaussianSmearer::kFunction) {
	- sm.func->SetParameter(0, kineProp);
	- Smeared = sm.func->GetRandom();
	- } else {
	- double sThrow = rand.Gaus(
	- 0, sm.width *
	- ((sm.type == GaussianSmearer::kAbsolute) ? 1.0 : kineProp));
	- Smeared = kineProp + sThrow;
	- }
	- Smeared = (Smeared < 0) ? 0 : Smeared;
	-#ifdef DEBUG_THRESACCEPT
	- std::cout << " -- Saw " << Smeared << " visible energy from " << kineProp
	- << " available. [PDG: " << fp->PDG() << "]" << std::endl;
	+ if (fp->Status() != kFinalState) {
	+#ifdef DEBUG_GAUSSSMEAR
	+ std::cout << " -- Not final state." << std::flush;
	#endif
	+ return;
	+ }

	- ri->RecVisibleEnergy.push_back(Smeared);
	- ri->TrueContribPDGs.push_back(fp->PDG());
	- }
	-#ifdef DEBUG_THRESACCEPT
	- std::cout << std::endl;
	+ if ((TrackedGausSmears.count(fp->PDG()) + VisGausSmears.count(fp->PDG())) ==
	+ 0) {
	+#ifdef DEBUG_GAUSSSMEAR
	+ std::cout << " -- Undetectable." << std::flush;
	#endif
	+ return;
	}

	- return ri;
	-}
	-
	-void GaussianSmearer::SmearRecoInfo(RecoInfo *ri) {
	- // Smear tracked particles
	- for (size_t p_it = 0; p_it < ri->RecObjMom.size(); ++p_it) {
	- if (!TrackedGausSmears.count(ri->RecObjClass[p_it])) {
	- continue;
	- }
	- TVector3 ThreeMom = ri->RecObjMom[p_it];
	- TVector3 OriginalKP = ThreeMom;
	- for (size_t sm_it = 0;
	- sm_it < TrackedGausSmears[ri->RecObjClass[p_it]].size(); ++sm_it) {
	- GSmear &sm = TrackedGausSmears[ri->RecObjClass[p_it]][sm_it];
	+ if (TrackedGausSmears.count(fp->PDG())) {
	+ TVector3 ThreeMom = fp->P3();
	+ for (size_t sm_it = 0; sm_it < TrackedGausSmears[fp->PDG()].size();
	+ ++sm_it) {
	+ GSmear &sm = TrackedGausSmears[fp->PDG()][sm_it];

	double kineProp = 0;

	switch (sm.smearVar) {
	case GaussianSmearer::kMomentum: {
	- kineProp = ri->RecObjMom[p_it].Mag();
	+ kineProp = fp->P3().Mag();
	break;
	}
	case GaussianSmearer::kKE: {
	- double mass = PhysConst::GetMass(ri->RecObjClass[p_it]) * 1.0E3;
	- kineProp = sqrt(ri->RecObjMom[p_it].Mag2() + mass * mass) - mass;
	+ kineProp = fp->KE();
	break;
	}
	case GaussianSmearer::kCosTheta: {
	- kineProp = OriginalKP.CosTheta();
	+ kineProp = fp->P3().CosTheta();
	break;
	}
	case GaussianSmearer::kTheta: {
	- kineProp = OriginalKP.Theta();
	+ kineProp = fp->P3().Theta();
	break;
	}
	default: { THROW("Trying to find particle value for a kNoVar."); }
	}

	double Smeared;
	+ size_t attempt = 0;
	bool ok = false;
	- int attempt = 0;
	while (!ok) {
	if (sm.type == GaussianSmearer::kFunction) {
	sm.func->SetParameter(0, kineProp);
	Smeared = sm.func->GetRandom();
	} else {
	- double sThrow =
	- rand.Gaus(0, sm.width * ((sm.type == GaussianSmearer::kAbsolute)
	- ? 1.0
	- : kineProp));
	+ double sThrow = rand.Gaus(
	+ 0, sm.width *
	+ ((sm.type == GaussianSmearer::kAbsolute) ? 1 : kineProp));
	Smeared = kineProp + sThrow;
	-#ifdef DEBUG_THRESACCEPT
	- std::cout << "*** [" << attempt << "] Gaus(0,"
	- << (sm.width * (sm.type == GaussianSmearer::kAbsolute)
	- ? 1
	- : kineProp)
	- << "[" << sm.width << "]) = " << sThrow << ": " << kineProp
	- << " -> " << Smeared << std::endl;
	-#endif
	}
	switch (
	sm.smearVar) { // Different kinematics have different truncation.
	case GaussianSmearer::kMomentum:
	case GaussianSmearer::kKE: {
	ok = (Smeared > 0);
	break;
	}
	case GaussianSmearer::kCosTheta: {
	ok = ((Smeared >= -1) && (Smeared <= 1));
	break;
	}
	case GaussianSmearer::kTheta: {
	ok = true;
	break;
	}
	+
	default: { THROW("SHOULDN'T BE HERE."); }
	}
	attempt++;
	+ if (attempt > 1000) {
	+ THROW("Didn't get a good smeared value after " << attempt
	+ << " attempts.");
	+ }
	}

	switch (sm.smearVar) {
	case GaussianSmearer::kMomentum: {
	ThreeMom = (ThreeMom.Unit() * Smeared);
	break;
	}
	case GaussianSmearer::kKE: {
	- double mass = PhysConst::GetMass(ri->RecObjClass[p_it]) * 1.0E3;
	+ double mass = fp->P4().M();
	double TE = mass + Smeared;
	double magP = sqrt(TE * TE - mass * mass);
	ThreeMom = (ThreeMom.Unit() * magP);
	break;
	}
	case GaussianSmearer::kCosTheta: {
	ThreeMom.SetTheta(acos(Smeared));
	break;
	}
	case GaussianSmearer::kTheta: {
	ThreeMom.SetTheta(Smeared);
	break;
	}
	default: {}
	}
	}
	- ri->RecObjMom[p_it] = ThreeMom;
	-
	-#ifdef DEBUG_THRESACCEPT
	+#ifdef DEBUG_GAUSSSMEAR
	std::cout << " -- momentum reconstructed as Mom: "
	- << ri->RecObjMom[p_it].Mag()
	- << ", CT: " << ri->RecObjMom[p_it].CosTheta() << " from "
	- << OriginalKP.Mag() << ", " << OriginalKP.CosTheta() << " true."
	+ << ri->RecObjMom.back().Mag()
	+ << ", CT: " << ri->RecObjMom.back().CosTheta() << " from "
	+ << ThreeMom.Mag() << ", " << fp->P3().CosTheta() << " true."
	<< std::endl;
	#endif
	- }
	+ ri->RecObjMom.push_back(ThreeMom);
	+ ri->RecObjClass.push_back(fp->PDG());
	+ } else { // Smear to EVis

	- for (size_t ve_it = 0; ve_it < ri->RecVisibleEnergy.size(); ++ve_it) {
	- if (!VisGausSmears.count(ri->TrueContribPDGs[ve_it])) {
	- continue;
	+ GSmear &sm = VisGausSmears[fp->PDG()];
	+
	+ double kineProp = 0;
	+
	+ switch (sm.smearVar) {
	+ case GaussianSmearer::kKEVis: {
	+ kineProp = fp->KE();
	+ break;
	+ }
	+ case GaussianSmearer::kTEVis: {
	+ kineProp = fp->E();
	+ break;
	+ }
	+ default: { THROW("Trying to find particle value for a kNoVar."); }
	}
	- GSmear &sm = VisGausSmears[ri->TrueContribPDGs[ve_it]];
	- double kineProp = ri->RecVisibleEnergy[ve_it];

	double Smeared;
	if (sm.type == GaussianSmearer::kFunction) {
	sm.func->SetParameter(0, kineProp);
	Smeared = sm.func->GetRandom();
	} else {
	double sThrow = rand.Gaus(
	0, sm.width *
	((sm.type == GaussianSmearer::kAbsolute) ? 1.0 : kineProp));
	Smeared = kineProp + sThrow;
	}
	Smeared = (Smeared < 0) ? 0 : Smeared;
	-#ifdef DEBUG_THRESACCEPT
	+#ifdef DEBUG_GAUSSSMEAR
	std::cout << " -- Saw " << Smeared << " visible energy from " << kineProp
	- << " available. [PDG: " << ri->TrueContribPDGs[ve_it] << "]"
	- << std::endl;
	+ << " available. [PDG: " << fp->PDG() << "]" << std::endl;
	+#endif
	+
	+ ri->RecVisibleEnergy.push_back(Smeared);
	+ ri->TrueContribPDGs.push_back(fp->PDG());
	+ }
	+#ifdef DEBUG_GAUSSSMEAR
	+ std::cout << std::endl;
	+#endif
	+}
	+
	+RecoInfo GaussianSmearer::Smearcept(FitEvent fe) {
	+ RecoInfo *ri = new RecoInfo();
	+
	+ for (size_t p_it = 0; p_it < fe->NParticles(); ++p_it) {
	+ FitParticle *fp = fe->GetParticle(p_it);
	+ SmearceptOneParticle(ri, fp
	+#ifdef DEBUG_GAUSSSMEAR
	+ ,
	+ p_it
	+#endif
	+ );
	+ }
	+
	+ return ri;
	+}
	+
	+void GaussianSmearer::SmearceptOneParticle(TVector3 &RecObjMom,
	+ int RecObjClass) {
	+ if (!TrackedGausSmears.count(RecObjClass)) {
	+ return;
	+ }
	+ TVector3 ThreeMom = RecObjMom;
	+ TVector3 OriginalKP = ThreeMom;
	+ for (size_t sm_it = 0; sm_it < TrackedGausSmears[RecObjClass].size();
	+ ++sm_it) {
	+ GSmear &sm = TrackedGausSmears[RecObjClass][sm_it];
	+
	+ double kineProp = 0;
	+
	+ switch (sm.smearVar) {
	+ case GaussianSmearer::kMomentum: {
	+ kineProp = RecObjMom.Mag();
	+ break;
	+ }
	+ case GaussianSmearer::kKE: {
	+ double mass = PhysConst::GetMass(RecObjClass) * 1.0E3;
	+ kineProp = sqrt(RecObjMom.Mag2() + mass * mass) - mass;
	+ break;
	+ }
	+ case GaussianSmearer::kCosTheta: {
	+ kineProp = OriginalKP.CosTheta();
	+ break;
	+ }
	+ case GaussianSmearer::kTheta: {
	+ kineProp = OriginalKP.Theta();
	+ break;
	+ }
	+ default: { THROW("Trying to find particle value for a kNoVar."); }
	+ }
	+
	+ double Smeared;
	+ bool ok = false;
	+ int attempt = 0;
	+ while (!ok) {
	+ if (sm.type == GaussianSmearer::kFunction) {
	+ sm.func->SetParameter(0, kineProp);
	+ Smeared = sm.func->GetRandom();
	+ } else {
	+ double sThrow = rand.Gaus(
	+ 0, sm.width *
	+ ((sm.type == GaussianSmearer::kAbsolute) ? 1.0 : kineProp));
	+ Smeared = kineProp + sThrow;
	+#ifdef DEBUG_GAUSSSMEAR
	+ std::cout << "*** [" << attempt << "] Gaus(0,"
	+ << (sm.width * (sm.type == GaussianSmearer::kAbsolute)
	+ ? 1
	+ : kineProp)
	+ << "[" << sm.width << "]) = " << sThrow << ": " << kineProp
	+ << " -> " << Smeared << std::endl;
	#endif
	+ }
	+ switch (sm.smearVar) { // Different kinematics have different truncation.
	+ case GaussianSmearer::kMomentum:
	+ case GaussianSmearer::kKE: {
	+ ok = (Smeared > 0);
	+ break;
	+ }
	+ case GaussianSmearer::kCosTheta: {
	+ ok = ((Smeared >= -1) && (Smeared <= 1));
	+ break;
	+ }
	+ case GaussianSmearer::kTheta: {
	+ ok = true;
	+ break;
	+ }
	+ default: { THROW("SHOULDN'T BE HERE."); }
	+ }
	+ attempt++;
	+ if (attempt > 1000) {
	+ THROW("Didn't get a good smeared value after " << attempt
	+ << " attempts.");
	+ }
	+ }

	- ri->RecVisibleEnergy[ve_it] = Smeared;
	+ switch (sm.smearVar) {
	+ case GaussianSmearer::kMomentum: {
	+ ThreeMom = (ThreeMom.Unit() * Smeared);
	+ break;
	+ }
	+ case GaussianSmearer::kKE: {
	+ double mass = PhysConst::GetMass(RecObjClass) * 1.0E3;
	+ double TE = mass + Smeared;
	+ double magP = sqrt(TE * TE - mass * mass);
	+ ThreeMom = (ThreeMom.Unit() * magP);
	+ break;
	+ }
	+ case GaussianSmearer::kCosTheta: {
	+ ThreeMom.SetTheta(acos(Smeared));
	+ break;
	+ }
	+ case GaussianSmearer::kTheta: {
	+ ThreeMom.SetTheta(Smeared);
	+ break;
	+ }
	+ default: {}
	+ }
	+ }
	+ RecObjMom = ThreeMom;
	+
	+#ifdef DEBUG_GAUSSSMEAR
	+ std::cout << " -- momentum reconstructed as Mom: " << RecObjMom.Mag()
	+ << ", CT: " << RecObjMom.CosTheta() << " from " << OriginalKP.Mag()
	+ << ", " << OriginalKP.CosTheta() << " true." << std::endl;
	+#endif
	+}
	+
	+void GaussianSmearer::SmearceptOneParticle(double &RecVisibleEnergy,
	+ int TrueContribPDG) {
	+ if (!VisGausSmears.count(TrueContribPDG)) {
	+ return;
	+ }
	+ GSmear &sm = VisGausSmears[TrueContribPDG];
	+ double kineProp = RecVisibleEnergy;
	+
	+ double Smeared;
	+ if (sm.type == GaussianSmearer::kFunction) {
	+ sm.func->SetParameter(0, kineProp);
	+ Smeared = sm.func->GetRandom();
	+ } else {
	+ double sThrow = rand.Gaus(
	+ 0,
	+ sm.width * ((sm.type == GaussianSmearer::kAbsolute) ? 1.0 : kineProp));
	+ Smeared = kineProp + sThrow;
	+ }
	+ Smeared = (Smeared < 0) ? 0 : Smeared;
	+#ifdef DEBUG_GAUSSSMEAR
	+ std::cout << " -- Saw " << Smeared << " visible energy from " << kineProp
	+ << " available. [PDG: " << TrueContribPDG << "]" << std::endl;
	+#endif
	+
	+ RecVisibleEnergy = Smeared;
	+}
	+
	+void GaussianSmearer::SmearRecoInfo(RecoInfo *ri) {
	+ // Smear tracked particles
	+ for (size_t p_it = 0; p_it < ri->RecObjMom.size(); ++p_it) {
	+ SmearceptOneParticle(ri->RecObjMom[p_it], ri->RecObjClass[p_it]);
	+ }
	+
	+ for (size_t ve_it = 0; ve_it < ri->RecVisibleEnergy.size(); ++ve_it) {
	+ SmearceptOneParticle(ri->RecVisibleEnergy[ve_it],
	+ ri->TrueContribPDGs[ve_it]);
	}

	-#ifdef DEBUG_THRESACCEPT
	+#ifdef DEBUG_GAUSSSMEAR
	std::cout << std::endl;
	#endif
	}
	diff --git a/src/Smearceptance/GaussianSmearer.h b/src/Smearceptance/GaussianSmearer.h
	index 4b9a62a..8add0cf 100644
	--- a/src/Smearceptance/GaussianSmearer.h
	+++ b/src/Smearceptance/GaussianSmearer.h
	@@ -1,55 +1,67 @@
	// 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 <http://www.gnu.org/licenses/>.
	*******************************************************************************/

	#ifndef GAUSSIANSMEARER_HXX_SEEN
	#define GAUSSIANSMEARER_HXX_SEEN

	#include "ISmearcepter.h"

	#include <map>

	+// #define DEBUG_GAUSSSMEAR
	+
	class GaussianSmearer : public ISmearcepter {
	public:
	enum GSmearType { kAbsolute, kFractional, kFunction, kNoType };
	enum DependVar { kMomentum, kKE, kKEVis, kTEVis, kCosTheta, kTheta, kNoVar };

	private:
	struct GSmear {
	GSmearType type;
	DependVar smearVar;
	double width;
	TF1 *func;
	- ~GSmear() { delete func; }
	};

	std::map<int, std::vector<GSmear> > TrackedGausSmears;
	std::map<int, GSmear> VisGausSmears;

	TRandom3 rand;

	void SpecifcSetup(nuiskey &);

	public:
	RecoInfo Smearcept(FitEvent );
	+
	+ void SmearceptOneParticle(RecoInfo ri, FitParticle fp
	+#ifdef DEBUG_GAUSSSMEAR
	+ ,
	+ size_t p_it
	+#endif
	+ );
	/// Helper method for using this class as a component in a more complex
	/// smearer
	void SmearRecoInfo(RecoInfo *);
	+
	+ void SmearceptOneParticle(TVector3 &RecObjMom, int RecObjClass);
	+
	+ void SmearceptOneParticle(double &RecVisibleEnergy, int TrueContribPDGs);
	};

	#endif
	diff --git a/src/Smearceptance/TrackedMomentumMatrixSmearer.cxx b/src/Smearceptance/TrackedMomentumMatrixSmearer.cxx
	index bb5251a..25f1f5f 100644
	--- a/src/Smearceptance/TrackedMomentumMatrixSmearer.cxx
	+++ b/src/Smearceptance/TrackedMomentumMatrixSmearer.cxx
	@@ -1,289 +1,431 @@
	// 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 <http://www.gnu.org/licenses/>.
	*******************************************************************************/

	#include "TrackedMomentumMatrixSmearer.h"

	-// #define DEBUG_EFFAPP
	namespace {
	TrackedMomentumMatrixSmearer::DependVar GetVarType(std::string const &axisvar) {
	if (axisvar == "Momentum") {
	return TrackedMomentumMatrixSmearer::kMomentum;
	} else if (axisvar == "KE") {
	return TrackedMomentumMatrixSmearer::kKE;
	+ } else if (axisvar == "TE") {
	+ return TrackedMomentumMatrixSmearer::kTE;
	}
	return TrackedMomentumMatrixSmearer::kNoVar;
	}
	}

	TH1D const *TrackedMomentumMatrixSmearer::SmearMap::GetRecoSlice(double val) {
	- // std::vector<std::pair<std::pair<double, double>, TH1D *> >
	-
	if ((val < RecoSlices.front().first.first) \|\|
	(val > RecoSlices.back().first.second)) {
	ERROR(WRN,
	"Kinematic property: " << val << ", not within smearable range: ["
	<< RecoSlices.front().first.first << " -- "
	<< RecoSlices.back().first.second << "].");
	return NULL;
	}

	int L = 0, U = RecoSlices.size();
	+#ifdef DEBUG_MATSMEAR
	+ std::cout << "\nBin search for " << val
	+ << " within: " << RecoSlices.front().first.first << " -- "
	+ << RecoSlices.back().first.second << std::endl;
	+#endif
	while (true) {
	+#ifdef DEBUG_MATSMEAR
	+ std::cout << "\t\t U:" << U << ", L: " << L << " (/" << RecoSlices.size()
	+ << ")" << std::endl;
	+#endif
	if (U == L) {
	+#ifdef DEBUG_MATSMEAR
	+ std::cout << "Found: " << L << std::endl;
	+#endif
	return RecoSlices[L].second;
	}
	int R = (U - L);
	int m = L + (R / 2);

	- if (val < RecoSlices[m].first.first) {
	+#ifdef DEBUG_MATSMEAR
	+ std::cout << "Bin: " << m << "[ " << RecoSlices[m].first.first << " -- "
	+ << RecoSlices[m].first.second << "] (Search: " << val << ")"
	+ << std::endl;
	+#endif
	+
	+ if (val <= RecoSlices[m].first.first) {
	U = m - 1;
	continue;
	}
	- if (val > RecoSlices[m].first.first) {
	+ if (val > RecoSlices[m].first.second) {
	L = m + 1;
	continue;
	}
	if ((val > RecoSlices[m].first.first) &&
	- (val < RecoSlices[m].first.first)) {
	+ (val <= RecoSlices[m].first.second)) {
	+#ifdef DEBUG_MATSMEAR
	+ std::cout << "Found: " << m << std::endl;
	+#endif
	return RecoSlices[m].second;
	}
	THROW("Binary smearing search failed. Check logic.");
	}
	}

	TH1D GetMapSlice(TH2D mp, int SliceBin, bool AlongX) {
	int NBins = (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetNbins();
	int NOtherBins = (AlongX ? mp->GetYaxis() : mp->GetXaxis())->GetNbins();
	if (SliceBin >= NOtherBins) {
	THROW("Asked for slice " << SliceBin << " but the " << (AlongX ? 'Y' : 'X')
	<< " axis only has " << NOtherBins);
	}
	- std::stringstream ss("");

	+ if ((AlongX ? mp->GetXaxis() : mp->GetYaxis())->IsVariableBinSize() &&
	+ ((NBins + 1) !=
	+ (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetXbins()->GetSize())) {
	+ THROW(
	+ "Attemping to take binning slice of variable binning, but NBins+1 != "
	+ "NBinEdges: "
	+ << (NBins + 1) << " != "
	+ << (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetXbins()->GetSize());
	+ }
	+
	+ std::stringstream ss("");
	ss << mp->GetName() << (AlongX ? 'Y' : 'X') << "Slice_" << SliceBin;
	- TH1D *Ret = new TH1D(
	- ss.str().c_str(), mp->GetTitle(), NBins,
	- (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetXbins()->GetArray());
	+ std::stringstream st("");
	+ st << mp->GetTitle() << ";"
	+ << (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetTitle() << ";"
	+ << "Count";
	+ TH1D *Ret;
	+ if ((AlongX ? mp->GetXaxis() : mp->GetYaxis())->IsVariableBinSize()) {
	+ Ret = new TH1D(
	+ ss.str().c_str(), st.str().c_str(), NBins,
	+ (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetXbins()->GetArray());
	+ } else {
	+ Ret = new TH1D(ss.str().c_str(), st.str().c_str(), NBins,
	+ (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetXmin(),
	+ (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetXmax());
	+ }
	+
	for (int bi_it = 0; bi_it < NBins + 2; ++bi_it) {
	int X = AlongX ? bi_it : SliceBin + 1;
	int Y = AlongX ? SliceBin + 1 : bi_it;
	int GBin = mp->GetBin(X, Y);
	Ret->SetBinContent(bi_it, mp->GetBinContent(GBin));
	Ret->SetBinError(bi_it, mp->GetBinError(GBin));
	}
	+
	+#ifdef DEBUG_MATSMEAR
	+ std::cout << "Took slice: " << SliceBin << " spanning ["
	+ << Ret->GetXaxis()->GetBinLowEdge(1) << " -- "
	+ << Ret->GetXaxis()->GetBinUpEdge(Ret->GetXaxis()->GetNbins())
	+ << "] (Orignal span ["
	+ << (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetBinLowEdge(1)
	+ << " -- "
	+ << (AlongX ? mp->GetXaxis() : mp->GetYaxis())
	+ ->GetBinUpEdge(
	+ (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetNbins())
	+ << "]) " << (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetXmin()
	+ << " -- " << (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetXmax()
	+ << " IsVBin: "
	+ << (AlongX ? mp->GetXaxis() : mp->GetYaxis())->IsVariableBinSize()
	+ << std::endl;
	+#endif
	+
	Ret->SetDirectory(NULL);
	return Ret;
	}

	void TrackedMomentumMatrixSmearer::SmearMap::SetSlicesFromMap(TH2D *map,
	bool TruthIsY) {
	int NSlices = (TruthIsY ? map->GetYaxis() : map->GetXaxis())->GetNbins();

	for (Int_t TrueSlice_it = 0; TrueSlice_it < NSlices; ++TrueSlice_it) {
	std::pair<double, double> BinEdges;
	BinEdges.first = (TruthIsY ? map->GetYaxis() : map->GetXaxis())
	->GetBinLowEdge(TrueSlice_it + 1);
	BinEdges.second = (TruthIsY ? map->GetYaxis() : map->GetXaxis())
	->GetBinUpEdge(TrueSlice_it + 1);

	TH1D *slice = GetMapSlice(map, TrueSlice_it, TruthIsY);

	RecoSlices.push_back(std::make_pair(BinEdges, slice));
	}
	QLOG(FIT, "\tAdded " << RecoSlices.size() << " reco slices.");
	}

	/// Reads particle efficiency nodes
	///
	/// Nodes look like:
	/// <TrackedMomentumMatrixSmearer Name="D00N_ND_LAr">
	-/// <SmearMatrix PDG="211,-211" InputFile="effs.root"
	-/// HistName="cpion_eff_mom_ctheta" Kinematics="[KE\|Momentum]"
	-/// GeVToMapUnits="1E3" YIsTrue="true"/>
	+/// <SmearMatrix PDG="211,-211" InputFile="smear.root"
	+/// HistName="cpion_mom_smear" Kinematics="[KE\|Momentum]"
	+/// MatrixToInternal="1E3" YIsTrue="true"/>
	/// </TrackedMomentumMatrixSmearer>
	void TrackedMomentumMatrixSmearer::SpecifcSetup(nuiskey &nk) {
	std::vector<nuiskey> effDescriptors = nk.GetListOfChildNodes("SmearMatrix");

	for (size_t t_it = 0; t_it < effDescriptors.size(); ++t_it) {
	std::string inputFileName = effDescriptors[t_it].GetS("InputFile");
	std::string HistName = effDescriptors[t_it].GetS("HistName");
	bool YIsTrue = effDescriptors[t_it].Has("YIsTrue")
	? effDescriptors[t_it].GetI("YIsTrue")
	: true;

	- double UnitsScale = effDescriptors[t_it].Has("GeVToMapUnits")
	- ? effDescriptors[t_it].GetD("GeVToMapUnits")
	+ double UnitsScale = effDescriptors[t_it].Has("MatrixToInternal")
	+ ? effDescriptors[t_it].GetD("MatrixToInternal")
	: 1;

	TFile inputFile(inputFileName.c_str());
	if (!inputFile.IsOpen()) {
	THROW("Couldn't open specified input root file: " << inputFileName);
	}

	TH2D inpHist = dynamic_cast<TH2D >(inputFile.Get(HistName.c_str()));
	if (!inpHist) {
	THROW("Couldn't get TH2D named: " << HistName << " from input root file: "
	<< inputFileName);
	}

	TrackedMomentumMatrixSmearer::DependVar var =
	GetVarType(effDescriptors[t_it].GetS("Kinematics"));

	std::string pdgs_s = effDescriptors[t_it].GetS("PDG");
	std::vector<int> pdgs_i = GeneralUtils::ParseToInt(pdgs_s, ",");
	for (size_t pdg_it = 0; pdg_it < pdgs_i.size(); ++pdg_it) {
	if (ParticleMappings.count(pdgs_i[pdg_it])) {
	ERROR(WRN, "Smearceptor " << ElementName << ":" << InstanceName
	- << " already has a efficiency for PDG: "
	+ << " already has a smearing for PDG: "
	<< pdgs_i[pdg_it]);
	}

	SmearMap sm;
	sm.SetSlicesFromMap(inpHist, YIsTrue);
	sm.SmearVar = var;
	sm.UnitsScale = UnitsScale;

	ParticleMappings[pdgs_i[pdg_it]] = sm;

	QLOG(SAM, "Added smearing map for PDG: " << pdgs_i[pdg_it]);
	}
	}
	+ SlaveGS.Setup(nk);
	}

	RecoInfo TrackedMomentumMatrixSmearer::Smearcept(FitEvent fe) {
	RecoInfo *ri = new RecoInfo();

	for (size_t p_it = 0; p_it < fe->NParticles(); ++p_it) {
	FitParticle *fp = fe->GetParticle(p_it);
	-#ifdef DEBUG_EFFAPP
	+#ifdef DEBUG_MATSMEAR
	std::cout << std::endl;
	std::cout << "[" << p_it << "]: " << fp->PDG() << ", " << fp->Status()
	<< ", " << fp->E() << " -- KE:" << fp->KE()
	<< " Mom: " << fp->P3().Mag() << std::flush;
	#endif

	if (fp->Status() != kFinalState) {
	-#ifdef DEBUG_EFFAPP
	+#ifdef DEBUG_MATSMEAR
	std::cout << " -- Not final state." << std::flush;
	#endif
	continue;
	}

	if (!ParticleMappings.count(fp->PDG())) {
	-#ifdef DEBUG_EFFAPP
	- std::cout << " -- Undetectable." << std::flush;
	+ SlaveGS.SmearceptOneParticle(ri, fp
	+#ifdef DEBUG_GAUSSSMEAR
	+ ,
	+ p_it
	#endif
	+ );
	continue;
	}

	SmearMap &sm = ParticleMappings[fp->PDG()];
	double kineProp = 0;

	switch (sm.SmearVar) {
	case kMomentum: {
	kineProp = fp->P3().Mag();
	break;
	}
	case kKE: {
	kineProp = fp->KE();
	break;
	}
	+ case kTE: {
	+ kineProp = fp->E();
	+ break;
	+ }
	default: { THROW("Trying to find particle value for a kNoAxis."); }
	}

	- TH1 const recoDistrib = sm.GetRecoSlice(kineProp sm.UnitsScale);
	+ TH1 const *recoDistrib = sm.GetRecoSlice(kineProp / sm.UnitsScale);
	+
	+#ifdef DEBUG_MATSMEAR
	+ std::cout << " -- Got slice spanning ["
	+ << recoDistrib->GetXaxis()->GetBinLowEdge(1) << " -- "
	+ << recoDistrib->GetXaxis()->GetBinUpEdge(
	+ recoDistrib->GetXaxis()->GetNbins())
	+ << "]" << std::endl;
	+#endif

	if (!recoDistrib) {
	-#ifdef DEBUG_EFFAPP
	+#ifdef DEBUG_MATSMEAR
	std::cout << " -- outside smearable range." << std::flush;
	#endif
	continue;
	}

	- double Smeared = recoDistrib->GetRandom() / sm.UnitsScale;
	+ if (recoDistrib->Integral() == 0) {
	+ ERROR(WRN, "True slice has no reconstructed events. Not smearing.")
	+ continue;
	+ }
	+
	+ double Smeared = recoDistrib->GetRandom() * sm.UnitsScale;
	+#ifdef DEBUG_MATSMEAR
	+ std::cout << "GotRandom: " << Smeared << ", MPV: "
	+ << recoDistrib->GetXaxis()->GetBinCenter(
	+ recoDistrib->GetMaximumBin()) *
	+ sm.UnitsScale
	+ << std::endl;
	+#endif

	switch (sm.SmearVar) {
	case kMomentum: {
	ri->RecObjMom.push_back(fp->P3().Unit() * Smeared);
	+
	+#ifdef DEBUG_MATSMEAR
	+ std::cout << " -- Smeared: " << fp->p() << " -> " << Smeared << "."
	+ << std::flush;
	+#endif
	+
	break;
	}
	case kKE: {
	double mass = fp->P4().M();
	double TE = mass + Smeared;
	double magP = sqrt(TE * TE - mass * mass);
	+
	ri->RecObjMom.push_back(fp->P3().Unit() * magP);
	+
	+#ifdef DEBUG_MATSMEAR
	+ std::cout << " -- Smeared: " << fp->KE() << " (mass: " << mass
	+ << ") -> " << Smeared
	+ << ". Smear Mom: " << ri->RecObjMom.back().Mag() << "."
	+ << std::flush;
	+#endif
	+ break;
	+ }
	+ case kTE: {
	+ double mass = fp->P4().M();
	+ double TE = Smeared;
	+ double magP = sqrt(TE * TE - mass * mass);
	+
	+ ri->RecObjMom.push_back(fp->P3().Unit() * magP);
	+
	+#ifdef DEBUG_MATSMEAR
	+ std::cout << " -- Smeared: " << fp->E() << " (mass: " << mass
	+ << ") -> " << Smeared
	+ << ". Smear Mom: " << ri->RecObjMom.back().Mag() << "."
	+ << std::flush;
	+#endif
	break;
	}
	default: {}
	}
	-#ifdef DEBUG_EFFAPP
	- std::cout << " -- momentum reconstructed as " << ri->RecObjMom.back() << "."
	- << std::flush;
	+#ifdef DEBUG_MATSMEAR
	+ std::cout << " -- momentum reconstructed as " << ri->RecObjMom.back().Mag()
	+ << "." << std::endl;
	#endif
	- ri->RecObjClass.push_back(fp->PDG());
	+ if (ri->RecObjMom.back().Mag() != ri->RecObjMom.back().Mag()) {
	+ ERROR(WRN, "Invalid particle built.");
	+ ri->RecObjMom.pop_back();
	+
	+#include "TCanvas.h"
	+ TCanvas *Test = new TCanvas("c1", "");
	+ static_cast<TH1 *>(recoDistrib->Clone())->Draw();
	+ Test->SaveAs("Fail.png");
	+ delete Test;
	+ THROW("ARGH");
	+ } else {
	+ ri->RecObjClass.push_back(fp->PDG());
	+ }
	}
	-#ifdef DEBUG_EFFAPP
	+#ifdef DEBUG_MATSMEAR
	std::cout << std::endl;
	#endif
	return ri;
	}

	void TrackedMomentumMatrixSmearer::SmearRecoInfo(RecoInfo *ri) {
	for (size_t p_it = 0; p_it < ri->RecObjMom.size(); ++p_it) {
	if (!ParticleMappings.count(ri->RecObjClass[p_it])) {
	+ SlaveGS.SmearceptOneParticle(ri->RecObjMom[p_it], ri->RecObjClass[p_it]);
	continue;
	}
	SmearMap &sm = ParticleMappings[ri->RecObjClass[p_it]];
	double kineProp = 0;

	switch (sm.SmearVar) {
	case kMomentum: {
	kineProp = ri->RecObjMom[p_it].Mag();
	break;
	}
	case kKE: {
	double mass = PhysConst::GetMass(ri->RecObjClass[p_it]) * 1E3;
	kineProp = sqrt(ri->RecObjMom[p_it].Mag2() + mass * mass) - mass;
	break;
	}
	+ case kTE: {
	+ double mass = PhysConst::GetMass(ri->RecObjClass[p_it]) * 1E3;
	+ kineProp = sqrt(ri->RecObjMom[p_it].Mag2() + mass * mass);
	+ break;
	+ }
	default: { THROW("Trying to find particle value for a kNoAxis."); }
	}
	- TH1 const recoDistrib = sm.GetRecoSlice(kineProp sm.UnitsScale);
	+ TH1 const *recoDistrib = sm.GetRecoSlice(kineProp / sm.UnitsScale);
	if (!recoDistrib) {
	continue;
	}

	- double Smeared = recoDistrib->GetRandom() / sm.UnitsScale;
	+ double Smeared = recoDistrib->GetRandom() * sm.UnitsScale;

	switch (sm.SmearVar) {
	case kMomentum: {
	ri->RecObjMom[p_it] = ri->RecObjMom[p_it].Unit() * Smeared;
	break;
	}
	case kKE: {
	double mass = PhysConst::GetMass(ri->RecObjClass[p_it]) * 1E3;
	double TE = mass + Smeared;
	double magP = sqrt(TE * TE - mass * mass);
	ri->RecObjMom[p_it] = ri->RecObjMom[p_it].Unit() * magP;
	break;
	}
	+ case kTE: {
	+ double mass = PhysConst::GetMass(ri->RecObjClass[p_it]) * 1E3;
	+ double TE = Smeared;
	+ double magP = sqrt(TE * TE - mass * mass);
	+ ri->RecObjMom[p_it] = ri->RecObjMom[p_it].Unit() * magP;
	+ break;
	+ }
	default: {}
	}
	}
	}
	diff --git a/src/Smearceptance/TrackedMomentumMatrixSmearer.h b/src/Smearceptance/TrackedMomentumMatrixSmearer.h
	index 0dfc9a3..c84698c 100644
	--- a/src/Smearceptance/TrackedMomentumMatrixSmearer.h
	+++ b/src/Smearceptance/TrackedMomentumMatrixSmearer.h
	@@ -1,66 +1,72 @@
	// 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 <http://www.gnu.org/licenses/>.
	*******************************************************************************/

	#ifndef TRACKEDMOMENTUMMATRIXSMEARER_HXX_SEEN
	#define TRACKEDMOMENTUMMATRIXSMEARER_HXX_SEEN

	#include "ISmearcepter.h"

	+#include "GaussianSmearer.h"
	+
	#include "TRandom3.h"
	#include "TH2D.h"
	#include "TH1D.h"

	#include <vector>
	#include <utility>

	+// #define DEBUG_MATSMEAR
	+
	class TrackedMomentumMatrixSmearer : public ISmearcepter {
	public:
	- enum DependVar { kMomentum, kKE, kNoVar };
	+ enum DependVar { kMomentum, kKE, kTE, kNoVar };

	private:
	class SmearMap {
	/// Input True -> Reco mapping.
	std::vector<std::pair<std::pair<double, double>, TH1D *> > RecoSlices;

	public:
	TH1D const *GetRecoSlice(double val);
	void SetSlicesFromMap(TH2D *, bool TruthIsY);
	/// Particle variable to smear: Momentum/KE
	///
	/// In the future should be able to smear multi-kinematic property
	/// distributions, but requires the definition of axis mappings that I
	/// cannot be bothered with at this second.
	TrackedMomentumMatrixSmearer::DependVar SmearVar;

	double UnitsScale;
	};
	std::map<int, SmearMap> ParticleMappings;

	+ GaussianSmearer SlaveGS;
	+
	void SpecifcSetup(nuiskey &);

	public:
	/// Will reject any particle that is not known about.
	RecoInfo Smearcept(FitEvent );
	/// Helper method for using this class as a component in a more complex
	/// smearer
	void SmearRecoInfo(RecoInfo *);
	~TrackedMomentumMatrixSmearer();
	};

	#endif
	diff --git a/src/Smearceptance/smearceptance.md b/src/Smearceptance/smearceptance.md
	index a9875fa..d96cb7f 100644
	--- a/src/Smearceptance/smearceptance.md
	+++ b/src/Smearceptance/smearceptance.md
	@@ -1,187 +1,262 @@
	# NUISANCE Smearceptance

	1. Introduction
	2. Quickstart: Applying thresholds
	3. Included smearcepters
	4. Writing your own

	## 1. Introduction
	-A generalised 'Fast MC' interface for NUISANCE. Smearers, accepters, and
	-smearcepters can be written and configured to mock up the effects of using a
	-'realistic' detection techniques on any input event format that NUISANCE can
	-read. A simple example of this would be a threshold accepter, which lets the
	-user define charged particle tracking thresholds for different types of final
	-state particles. These thresholds would be applied so that any downstream
	+A generalised 'Fast MC' interface for NUISANCE. Smearers, accepters, and
	+smearcepters can be written and configured to mock up the effects of using a
	+'realistic' detection techniques on any input event format that NUISANCE can
	+read. A simple example of this would be a threshold accepter, which lets the
	+user define charged particle tracking thresholds for different types of final
	+state particles. These thresholds would be applied so that any downstream
	analysis can only see particles above threshold.

	-The eventual aim of such a module would be to expand the use of NUISANCE as a
	-data-release platform so that analysers could add their
	-reconstruction-level/uncorrected results, along with the information to
	-forward-fold, or 'smearcept' theory inputs for comparison and model tuning. It
	-is not currently clear that this is feasible, but it is an interesting
	+The eventual aim of such a module would be to expand the use of NUISANCE as a
	+data-release platform so that analysers could add their
	+reconstruction-level/uncorrected results, along with the information to
	+forward-fold, or 'smearcept' theory inputs for comparison and model tuning. It
	+is not currently clear that this is feasible, but it is an interesting
	possibility.

	-The current use cases include generating fake data for different types of
	-detector technology, running simple mock experiments as a first gauge of
	-potential sensitivities, and simple model overlays for uncorrected data from
	+The current use cases include generating fake data for different types of
	+detector technology, running simple mock experiments as a first gauge of
	+potential sensitivities, and simple model overlays for uncorrected data from
	both nu-A and e-A scattering experiments.

	-The interface is designed to be abstract, and very simple. With any real data
	-and analyses, the devil is often in the detail and it is likely that tailored
	-Smearcepters would have to be written for each analysis. However, a number of
	-examples and simple smearcepters already exist and are useful for applying
	-thresholds, multi-dimensional efficiency curves, and simple detector smearing
	+The interface is designed to be abstract, and very simple. With any real data
	+and analyses, the devil is often in the detail and it is likely that tailored
	+Smearcepters would have to be written for each analysis. However, a number of
	+examples and simple smearcepters already exist and are useful for applying
	+thresholds, multi-dimensional efficiency curves, and simple detector smearing
	effects.

	## 2. Quickstart: Applying thresholds

	Steps to acceptance:
	* Write a smearcepter card file, `WaterCherenkovThresholds.xml`:
	```xml
	<nuisance>
	<smearcepters>
	<ThresholdAccepter Name="WCThresh">
	<RecoThreshold PDG="13" RecoThresholdMomentum_MeV="160" />
	<RecoThreshold PDG="2212" RecoThresholdMomentum_MeV="1400" />
	</ThresholdAccepter>
	</smearcepters>
	</nuisance>
	```
	-* Produce a smeared event summary tree:
	+* Produce a smeared event summary tree:
	`$ nuissmear -i InputVector.root -c WaterCherenkovThresholds.xml -t WCThresh -o WaterCherenkov_summary.root`
	-* Draw plots from the event tree:
	+* Draw plots from the event tree:
	`[root] FlatTree_VARS->Draw("EISLep_true:EISLep_LepHad_rec >> (100,0,2000,100,0,2000)","Weight*(flagCCINC_rec==1)","COLZ")`

	-The branchs in the event summary tree are described in detail the class
	+The branchs in the event summary tree are described in detail the class
	documentation of `src/MCStudies/Smearceptance_Tester.cxx`

	## 3. Included smearcepters

	### General

	-Smearcepters are defined in a NUISANCE xml configuration file within a tag of
	-the `<nuisance>` root element named `<smearcepters>`. Any number can be defined,
	-and each should be distinctly named with a `Name="<name>"` attribute.
	+Smearcepters are defined in a NUISANCE xml configuration file within a tag of
	+the `<nuisance>` root element named `<smearcepters>`. Any number can be defined,
	+and each should be distinctly named with a `Name="<name>"` attribute.

	-An example XML tag is shown below that highlights some of the notation used in
	+An example XML tag is shown below that highlights some of the notation used in
	this section.
	```xml
	<tag A[B\|C]="<name>" (D="<value>") E="[F\|G]" />
	```
	-* `AB` and `AC` are both valid attribute names, at least one of which is a
	+* `AB` and `AC` are both valid attribute names, at least one of which is a
	required attribute.
	* `D` is an optional attribute.
	* `<name>` is an example attribute value.
	* The attribute `E` must take either the value `F` or `G`.

	### ThresholdAccepter

	Applies tracking and visible energy thresholds to an input vector.

	#### Full example
	```xml
	<nuisance>
	<smearcepters>
	<ThresholdAccepter Name="ForwardMuonHighHadronThresholdDetector">
	<RecoThreshold PDG="13" RecoThresholdMomentum_MeV="300" />
	<RecoThreshold PDG="13" RecoThresholdCosTheta_Min="0.5" />
	<RecoThreshold PDG="2212" RecoThresholdMomentum_MeV="500" />
	<RecoThreshold PDG="2212" RecoThresholdAbsCosTheta_Min="0.5" />
	<RecoThreshold PDG="221,-211" RecoThresholdKE_MeV="120" />
	-
	+
	<VisThreshold PDG="2212" VisThresholdKE_MeV="10" Contrib="K" />
	<VisThreshold PDG="111" VisThresholdKE_MeV="0" Contrib="T" />
	<VisThreshold PDG="221,-211" VisThresholdKE_MeV="30" Contrib="K" Fraction="0.9" />
	</ThresholdAccepter>
	</smearcepters>
	</nuisance>
	```

	#### Details
	-A tracking threshold can be placed on a particles production zenith angle by
	+A tracking threshold can be placed on a particles production zenith angle by
	the element:
	* `<RecoThreshold PDG="<211,-211,2212>" RecoThresholdCosTheta_[Min,Max]="<value>" />`
	* `<RecoThreshold PDG="<X,Y,Z>" RecoThresholdAbsCosTheta_[Min,Max]="<value>" />`

	-Tracking and visible energy threshold values can be applied by the following
	+Tracking and visible energy threshold values can be applied by the following
	elements:
	* Particle momentum: Use `<[Reco\|Vis]Threshold PDG="<X,Y,Z>" [Reco\|Vis]ThresholdMomentum_MeV="<value>" (Contrib="[T\|K]" Fraction="<value>") />`
	-* Particle kinetic energy: Use `<[Reco\|Vis]Threshold PDG="<X,Y,Z>" [Reco\|Vis]ThresholdKE_MeV="<value>" (Contrib="[T\|K]" Fraction="<value>") />`
	+* Particle kinetic energy: Use `<[Reco\|Vis]Threshold PDG="<X,Y,Z>" [Reco\|Vis]ThresholdKE_MeV="<value>" (Contrib="[T\|K]" Fraction="<value>") />`

	-When defining visible energy thresholds, the additional attributes
	-`Contrib="[<T>,K]"` and `Fraction="<value>"` can be specified to detail whether
	-the energy deposit contribution comes from the particle total or kinetic energy
	-only and what fraction of that deposit is visible. The deposited energy defaults
	+When defining visible energy thresholds, the additional attributes
	+`Contrib="[<T>,K]"` and `Fraction="<value>"` can be specified to detail whether
	+the energy deposit contribution comes from the particle total or kinetic energy
	+only and what fraction of that deposit is visible. The deposited energy defaults
	to the particle total energy and the fraction defaults to 1.

	-If a particle kinematics do not exceed a momentum or kinetic energy threshold
	-for tracking they are tested against the visible energy deposit threshold. If a
	-particle is rejected due to an angular threshold, it is not allowed to deposit
	+If a particle kinematics do not exceed a momentum or kinetic energy threshold
	+for tracking they are tested against the visible energy deposit threshold. If a
	+particle is rejected due to an angular threshold, it is not allowed to deposit
	energy.

	-Multiple thresholds can target the same input PDG, however it is undefined
	-behavior to use more than one energy/momentum tracking threshold, or use both
	+Multiple thresholds can target the same input PDG, however it is undefined
	+behavior to use more than one energy/momentum tracking threshold, or use both
	`Abs` and non-`Abs` zenith angle cuts.

	-N.B. If a PDG code is encountered for which there are no defined thresholds,
	+N.B. If a PDG code is encountered for which there are no defined thresholds,
	it is not accepted.

	### GaussianSmearer

	Documentation to follow. Sorry.

	### EfficiencyApplicator

	-Applies tracking efficiencies from one, two, or three dimensional efficiencies
	+Applies tracking efficiencies from one, two, or three dimensional efficiencies
	passed in as input.

	#### Full example
	```xml
	<EfficiencyApplicator Name="SimpleLArLike">
	<EfficiencyCurve PDG="13" InputFile="uBooNEEffs.root" HistName="Tracking/Muon/p_costheta" NDims="2" XAxis="kMomentum" YAxis="kCosTheta" Interpolate="false" XAxisScaleToMeV="1E3" />
	<EfficiencyCurve PDG="211,-211" InputFile="uBooNEEffs.root" HistName="Tracking/Pion/p_costheta" NDims="2" XAxis="kMomentum" YAxis="kCosTheta" Interpolate="false" XAxisScaleToMeV="1E3" />
	<EfficiencyCurve PDG="2212" InputFile="uBooNEEffs.root" HistName="Tracking/Proton/p_costheta" NDims="2" XAxis="kMomentum" YAxis="kCosTheta" Interpolate="false" XAxisScaleToMeV="1E3" />
	-
	+
	<VisThreshold PDG="2212" VisThresholdKE_MeV="10" Contrib="K" />
	<VisThreshold PDG="211,-211" VisThresholdKE_MeV="30" Contrib="K" />
	<VisThreshold PDG="111" VisThresholdKE_MeV="0" Contrib="T" />
	</EfficiencyApplicator>
	```
	#### Details

	Each `<EfficiencyCurve>` tag requires at least:
	* `PDG="<211,-211>"`: The particle IDs to apply this efficiency to.
	-* `InputFile="<file.root>"`: The location of the ROOT file containing the
	+* `InputFile="<file.root>"`: The location of the ROOT file containing the
	efficiency histogram.
	-* `HistName="<histname>"`: The name (may include directories) of the `TH1` or
	+* `HistName="<histname>"`: The name (may include directories) of the `TH1` or
	`TEfficiency` that describes the efficiency.
	* `NDims="[1\|2\|3]"`: The number of kinematic axes for the efficiency curve.
	-* `XAxis="[kMomentum\|kKE\|kCosTheta\|kTheta\|kPhi]"`: The kinematic variable
	+* `XAxis="[kMomentum\|kKE\|kCosTheta\|kTheta\|kPhi]"`: The kinematic variable
	spanned by the first dimension of the efficiency histogram.

	optional attributes:
	-* `YAxis="[kMomentum\|kKE\|kCosTheta\|kTheta\|kPhi]"`: The kinematic variable
	+* `YAxis="[kMomentum\|kKE\|kCosTheta\|kTheta\|kPhi]"`: The kinematic variable
	spanned by the second dimension of the efficiency histogram.
	-* `ZAxis="[kMomentum\|kKE\|kCosTheta\|kTheta\|kPhi]"`: The kinematic variable
	+* `ZAxis="[kMomentum\|kKE\|kCosTheta\|kTheta\|kPhi]"`: The kinematic variable
	spanned by the third dimension of the efficiency histogram.
	-* `Interpolate="[true\|false]"`: Whether to use TH1:Interpolate or to take the
	+* `Interpolate="[true\|false]"`: Whether to use TH1:Interpolate or to take the
	bin-averaged efficiency.
	-* `[X\|Y\|Z]AxisScaleToInternal="<value>"`: A scale factor to translate the axes
	-to the internal NUISANCE units (MeV and radians). e.g. Use
	-`XAxisScaleToInternal="1E3"` if the input histogram uses kinetic energy in units
	+* `[X\|Y\|Z]AxisScaleToInternal="<value>"`: A scale factor to translate the axes
	+to the internal NUISANCE units (MeV and radians). e.g. Use
	+`XAxisScaleToInternal="1E3"` if the input histogram uses kinetic energy in units
	of GeV.

	-An `EfficiencyApplicator` Can also contain `<RecoThreshold>` and
	-`<VisThreshold>` tags, as described in `ThresholdAccepter`. Particles that are
	-rejected due to inefficiency get passed to a configured `ThresholdAccepter`
	-instance. N.B. it is unlikely that you want to have an `<EfficiencyCurve>` tag
	-and a `<RecoThreshold>` tag for the same PDG. However, it is possible that
	-particles that failed to be tracked, left some visible energy.
	+An `EfficiencyApplicator` Can also contain `<RecoThreshold>` and
	+`<VisThreshold>` tags, as described in `ThresholdAccepter`. Particles that are
	+rejected due to inefficiency get passed to a configured `ThresholdAccepter`
	+instance. N.B. it is unlikely that you want to have an `<EfficiencyCurve>` tag
	+and a `<RecoThreshold>` tag for the same PDG. However, it is possible that
	+particles that failed to be tracked, left some visible energy.
	+
	+### GaussianSmearer
	+
	+Applies simple smearing to tracked particle kinematics or visible energy deposits.
	+
	+#### Full example
	+```xml
	+<GaussianSmearer Name="D00N_ND_LAr">
	+ <Smear PDG="211" Type="Fractional" Kinematics="KE" Width="0.1" />
	+ <Smear PDG="13" Type="Absolute" Kinematics="Momentum" Width="0.025" />
	+ <Smear PDG="2212" Type="Function" Kinematics="Momentum" Function="{P1}*TMath::Landau(x - {V})$0$1E4" P1="0.2" />
	+</GaussianSmearer>
	+```
	+#### Details
	+
	+Each `<Smear>` tag requires at least:
	+* `PDG="<211,-211>"`: The particle IDs to apply this smearing to.
	+* `Type="[Fractional\|Absolute\|Function]"`: The type of smearer to use. `Fractional` and `Absolute` use a Gaussian smearer. For `Fractional` the Gaussian width is set to `<WidthAttribute>*ParticleKinematicPropertyValue` and for Absolute it is simple `<WidthAttribute>`.
	+* `Kinematics="[KE\|Momentum\|TEVis\|KEVis\|Theta\|CosTheta]"`: The particle kinematics to smear. For the majority of practical uses, where an accepter is also in use, `TEVis` and `KEVis` produce the same results. As only a single property smearer can be used per PDG, the use of `Theta` and `CosTheta` are of very limited use.
	+
	+Additional attributes that depend on the value of the `Type` attribute are as follows:
	+* `Width="<value>"`: The width of the Gaussian used for `Fractional` or `Absolute` type smearers.
	+* `Function="<TMath::Gaus(x - {V},0,0.2*{V})>"`: The function to throw the smeared values from for `Function` type smearers. The example function here should give the same result as `<Smear Type="Fractional" Width="0.2" />`.
	+* `P[1\|2\|3\|...]="<value>"`: Extra parameters to replace in the `Function` attribute value. e.g.: `<Smear Function="{V} + TMath::Gaus(x - {V},0,{P1}*{V})" P1="0.2" />`. Mostly just for clarity. An element can contain any number of numbered parameters, but they must be numbered in increasing order.

	### TrackedMomentumMatrixSmearer

	-Documentation to follow. Sorry.
	+Applies momentum or kinetic energy smearing to tracked particles.
	+
	+#### Full example
	+```xml
	+<TrackedMomentumMatrixSmearer Name="D00N_ND_LAr">
	+ <SmearMatrix PDG="211,-211" InputFile="smear.root" HistName="cpion_mom_smear" Kinematics="Momentum" MatrixToInternal="1E3" YIsTrue="true"/>
	+ <Smear PDG="13" Type="Absolute" Kinematics="KEVis" Width="0.025" />
	+</TrackedMomentumMatrixSmearer>
	+```
	+#### Details
	+
	+Each `<Smear>` tag requires at least:
	+* `PDG="<211,-211>"`: The particle IDs to apply this smearing to.
	+* `InputFile="<ROOT file location>"`: The root file that contains the input smearing matrix.
	+* `HistName="<TH1D name>"`: The name of the histogram within in the input root file. Can be located within a subdirectory of the TFile.
	+* `Kinematics="[KE\|TE\|Momentum]"`: The kinematics to smear.
	+* `MatrixToInternal="<value>"`: The scale factor used to scale the units used in the smearing matrix to MeV.
	+* YIsTrue="[1\|0]": Whether the Y or X axis of the input histogram denotes the true kinematic property axis.
	+
	+The `<TrackedMomentummatrixSmearer>` element can also contain any number of `<Smear>` element, which are used for PDG codes that do not have a `<SmearMatrix>` element. Their behavior is as in `<GaussianSmearer>`. As this smearer can only smear tracked particles, using `<Smear>` elements to add visible energy smearing is the expected use case.
	+
	+### VisECoalescer
	+
	+### EnergyShuffler

	### MetaSimpleSmearcepter

	+#### Full example
	+
	+```xml
	+<MetaSimpleSmearcepter Name="D00N_ND_LAr">
	+ <EfficiencyApplicator>
	+ <EfficiencyCurve PDG="13" InputFile="uBooNEEffs.root" HistName="Tracking/Muon/p_costheta" NDims="2" XAxis="kMomentum" YAxis="kCosTheta" Interpolate="false" XAxisScaleToMeV="1E3" />
	+ <EfficiencyCurve PDG="211,-211" InputFile="uBooNEEffs.root" HistName="Tracking/Pion/p_costheta" NDims="2" XAxis="kMomentum" YAxis="kCosTheta" Interpolate="false" XAxisScaleToMeV="1E3" />
	+ <EfficiencyCurve PDG="2212" InputFile="uBooNEEffs.root" HistName="Tracking/Proton/p_costheta" NDims="2" XAxis="kMomentum" YAxis="kCosTheta" Interpolate="false" XAxisScaleToMeV="1E3" />
	+
	+ <VisThreshold PDG="2212" VisThresholdKE_MeV="10" Contrib="K" />
	+ <VisThreshold PDG="211,-211" VisThresholdKE_MeV="30" Contrib="K" />
	+ <VisThreshold PDG="111" VisThresholdKE_MeV="0" Contrib="T" />
	+ </EfficiencyApplicator>
	+ <VisECoaslescer/> <!-- Obscure information about VisEContributor multiplicity >
	+ <EnergyShuffler> <!-- Apply event-by-event energy shuffling outside the FSI/interaction model >
	+ <Shuffle From="2212" To="2112" Amount="0.1" />
	+ <Shuffle From="211,-211" Amount="0.05" />
	+ <Shuffle From="111" Amount="0.25" />
	+ </EnergyShuffler>
	+ <TrackedMomentumMatrixSmearer>
	+ <SmearMatrix PDG="211,-211" InputFile="smear.root" HistName="cpion_mom_smear" Kinematics="Momentum" MatrixToInternal="1E3" YIsTrue="true"/>
	+ <SmearMatrix PDG="13" InputFile="smear.root" HistName="cpion_mom_smear" Kinematics="Momentum" MatrixToInternal="1E3" YIsTrue="true"/>
	+ <Smear PDG="2212" Type="Absolute" Kinematics="KEVis" Width="0.025" />
	+ <Smear PDG="111" Type="Fractional" Kinematics="KEVis" Width="0.1" />
	+ </TrackedMomentumMatrixSmearer>
	+</MetaSimpleSmearcepter>
	+```
	+
	Documentation to follow. Sorry.

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 3:53 PM (1 d, 20 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4023357
Default Alt Text: (110 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions