diff --git a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx
index 014f927..4d56759 100644
--- a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx
+++ b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx
@@ -1,924 +1,924 @@
// 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 .
*******************************************************************************/
#include "Smear_SVDUnfold_Propagation_Osc.h"
#include "SmearceptanceUtils.h"
#include "OscWeightEngine.h"
#include "HistogramInputHandler.h"
void Smear_SVDUnfold_Propagation_Osc::AddNDInputs(nuiskey &samplekey) {
NDSample nds;
// Plot Setup -------------------------------------------------------
// Check that we have the relevant near detector histograms specified.
if (!NDSamples.size()) { // If this is the first ND sample, take from the
// sample input
InputHandlerBase *InputBase = GetInput();
if (InputBase->GetType() != kHISTO) {
THROW(
"Smear_SVDUnfold_Propagation_Osc expects a Histogram input that "
"contains the ND observed spectrum.");
}
HistoInputHandler *HInput = dynamic_cast(InputBase);
if (!HInput) {
THROW(
"Smear_SVDUnfold_Propagation_Osc expects a Histogram input that "
"contains the ND observed spectrum.");
}
if (HInput->NHistograms() != 2) {
THROW(
"Input expected to contain 2 histograms. "
"HISTO:input.root[NDObs_TH1D,NDSmear_TH2D]");
}
nds.NDDataHist = dynamic_cast(HInput->GetHistogram(0));
nds.NDToSpectrumSmearingMatrix =
dynamic_cast(HInput->GetHistogram(1));
if (!nds.NDDataHist) {
THROW("Expected a valid TH1D input for the ND observed spectrum.");
}
if (!nds.NDToSpectrumSmearingMatrix) {
THROW("Expected a valid TH2D input for the ND observed smearing.");
}
} else {
std::vector NDObsInputs =
PlotUtils::GetTH1sFromRootFile(samplekey.GetS("ObsInput"));
if (NDObsInputs.size() < 2) {
THROW(
"Near detector sample must contain the observed ERec spectrum and "
"the "
"ND ETrue/ERec smearing matrix. e.g. "
"ObsInput=\"input.root[NDObs_species,NDSmearing_species]\"");
}
nds.NDDataHist = dynamic_cast(NDObsInputs[0]);
if (!nds.NDDataHist) {
ERROR(FTL,
"First histogram from ObsInput attribute was not a TH1D containing "
"the near detector observed ERec spectrum ("
<< samplekey.GetS("ObsInput") << ").");
THROW(
"Near detector sample must contain the observed ERec spectrum and "
"the "
"ND ETrue/ERec smearing matrix. e.g. "
"ObsInput=\"input.root[FDObs_species,FDSmearing_species]\"");
}
nds.NDToSpectrumSmearingMatrix = dynamic_cast(NDObsInputs[1]);
if (!nds.NDToSpectrumSmearingMatrix) {
ERROR(
FTL,
"Second histogram from ObsInput attribute was not a TH2D containing "
"the near detector ETrue/ERec smearing matrix ("
<< samplekey.GetS("ObsInput") << ").");
THROW(
"Near detector sample must contain the observed ERec spectrum and "
"the "
"ND ETrue/ERec smearing matrix. e.g. "
"ObsInput=\"input.root[FDObs_species,FDSmearing_species]\"");
}
}
nds.NDDataHist->Scale(ScalePOT);
if (UseRateErrors) {
for (Int_t bi_it = 1; bi_it < nds.NDDataHist->GetXaxis()->GetNbins() + 1;
++bi_it) {
nds.NDDataHist->SetBinError(bi_it,
sqrt(nds.NDDataHist->GetBinContent(bi_it)));
}
}
nds.TruncateStart = 0;
if (samplekey.Has("TruncateStart")) {
nds.TruncateStart = samplekey.GetI("TruncateStart");
}
nds.TruncateUpTo = 0;
if (samplekey.Has("TruncateUpTo")) {
nds.TruncateUpTo = samplekey.GetI("TruncateUpTo");
QLOG(SAM, "\tAllowed to truncate unfolding matrix by up to "
<< nds.TruncateUpTo
<< " singular values to limit negative ENu spectra.");
}
nds.NuPDG = 14;
if (samplekey.Has("NuPDG")) {
nds.NuPDG = samplekey.GetI("NuPDG");
}
NDSamples.push_back(nds);
}
void Smear_SVDUnfold_Propagation_Osc::SetupNDInputs() {
for (size_t nd_it = 0; nd_it < NDSamples.size(); ++nd_it) {
NDSample &nds = NDSamples[nd_it];
TMatrixD NDToSpectrumResponseMatrix_l = SmearceptanceUtils::GetMatrix(
SmearceptanceUtils::SVDGetInverse(nds.NDToSpectrumSmearingMatrix));
nds.NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l);
nds.NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l;
if (nds.TruncateStart != 0) {
TMatrixD NDToSpectrumResponseMatrix_l =
SmearceptanceUtils::GetMatrix(SmearceptanceUtils::SVDGetInverse(
nds.NDToSpectrumSmearingMatrix, nds.TruncateStart));
nds.NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l);
nds.NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l;
}
if (nds.TruncateStart >= nds.TruncateUpTo) {
nds.TruncateUpTo = nds.TruncateStart + 1;
}
UnfoldToNDETrueSpectrum(nd_it);
}
}
void Smear_SVDUnfold_Propagation_Osc::ReadExtraConfig(nuiskey &samplekey) {
UseRateErrors = false;
if (samplekey.Has("SetErrorsFromRate")) {
UseRateErrors = samplekey.GetI("SetErrorsFromRate");
}
NDetectorInfo.first = 0xdeadbeef;
if (samplekey.Has("DetectorVolume") && samplekey.Has("DetectorDensity")) {
NDetectorInfo.first = samplekey.GetD("DetectorVolume");
NDetectorInfo.second = samplekey.GetD("DetectorDensity");
double TargetMass_kg = NDetectorInfo.first * NDetectorInfo.second;
QLOG(SAM, "\tND sample detector mass : ");
QLOG(SAM, "\t\tTarget volume : " << NDetectorInfo.first);
QLOG(SAM, "\t\tTarget density : " << NDetectorInfo.second);
QLOG(SAM, "\t\tTarget mass : " << TargetMass_kg << " kg");
}
ScalePOT = 1;
if (samplekey.Has("ScalePOT")) {
ScalePOT = samplekey.GetD("ScalePOT");
}
}
void Smear_SVDUnfold_Propagation_Osc::AddFDTarget(nuiskey &nk) {
FDSample fds;
fds.FitRegion_Min = 0xdeadbeef;
if (nk.Has("FitRegion_Min")) {
fds.FitRegion_Min = nk.GetD("FitRegion_Min");
QLOG(SAM, "FD Sample [" << FDSamples.size() << "] imposes FitRegion E_nu > "
<< fds.FitRegion_Min);
}
if ((FitRegion_Min == 0xdeadbeef) || FitRegion_Min > fds.FitRegion_Min) {
FitRegion_Min = fds.FitRegion_Min;
}
nk.Print();
fds.FitRegion_Max = 0xdeadbeef;
if (nk.Has("FitRegion_Max")) {
fds.FitRegion_Max = nk.GetD("FitRegion_Max");
QLOG(SAM, "FD Sample [" << FDSamples.size() << "] imposes FitRegion E_nu < "
<< fds.FitRegion_Max);
}
if ((FitRegion_Max == 0xdeadbeef) || FitRegion_Max < fds.FitRegion_Max) {
FitRegion_Max = fds.FitRegion_Max;
}
fds.OscillateToPDG = 0;
if (nk.Has("OscillateToPDG")) {
fds.OscillateToPDG = nk.GetD("OscillateToPDG");
}
std::vector FDNDRatioElements = nk.GetListOfChildNodes("FDNDRatio");
for (size_t fdnd_it = 0; fdnd_it < FDNDRatioElements.size(); ++fdnd_it) {
nuiskey &fnr = FDNDRatioElements[fdnd_it];
if (fnr.Has("FromPDG") && fnr.Has("DivergenceFactor")) {
fds.FDNDRatios[fnr.GetI("FromPDG")] = fnr.GetD("DivergenceFactor");
QLOG(SAM, "FDND DivergenceFactor for far detector sample index: "
<< FDSamples.size() << " for PDG: " << fnr.GetI("FromPDG")
<< " -> " << fds.OscillateToPDG << " = "
<< fnr.GetD("DivergenceFactor"));
} else {
THROW(
"Far detector sample contained FDNDRatio element, but couldn't find "
"both FromPDG and Factor attributes.");
}
}
fds.FDNDMassRatio = 1;
if (NDetectorInfo.first != 0xdeadbeef) {
if ((!nk.Has("DetectorVolume")) || (!nk.Has("DetectorDensity"))) {
THROW(
"Near detector sample has specified volume but FD doesn't. This is "
"needed to scale the predicted event rate by the mass ratio.");
}
fds.DetectorInfo.first = nk.GetD("DetectorVolume");
fds.DetectorInfo.second = nk.GetD("DetectorDensity");
double TargetMass_kg = fds.DetectorInfo.first * fds.DetectorInfo.second;
fds.FDNDMassRatio =
TargetMass_kg / (NDetectorInfo.first * NDetectorInfo.second);
QLOG(SAM, "\tFD[" << FDSamples.size() << "] Event rate prediction : ");
QLOG(SAM, "\t\tTarget volume : " << fds.DetectorInfo.first);
QLOG(SAM, "\t\tTarget density : " << fds.DetectorInfo.second);
QLOG(SAM, "\t\tFD/ND mass : " << fds.FDNDMassRatio);
}
if (!nk.Has("ObsInput")) {
THROW("Far detector sample must specify at least ObsInput.");
}
std::vector FDObsInputs =
PlotUtils::GetTH1sFromRootFile(nk.GetS("ObsInput"));
if (FDObsInputs.size() < 2) {
THROW(
"Far detector sample must contain the observed ERec spectrum and the "
"FD ETrue/ERec smearing matrix. "
"ObsInput=\"input.root[FDObs_species,FDSmearing_species]\"");
}
fds.FDDataHist = NULL;
for (size_t hist_it = 0; hist_it < FDObsInputs.size() - 1; ++hist_it) {
if (!dynamic_cast(FDObsInputs[hist_it])) {
ERROR(FTL, "Input spectrum index "
<< hist_it
<< " from ObsInput attribute was not a TH1D containing "
"a far detector observed ERec spectrum ("
<< nk.GetS("ObsInput") << ").");
THROW(
"Far detector sample must contain the observed ERec spectrum and the "
"FD ETrue/ERec smearing matrix. "
"ObsInput=\"input.root[FDObs_species,(FDObs_species2),FDSmearing_"
"species]\"");
}
FDObsInputs[hist_it]->Scale(ScalePOT);
if (!fds.FDDataHist) {
fds.FDDataHist = dynamic_cast(FDObsInputs[hist_it]);
} else {
fds.FDDataHist->Add(dynamic_cast(FDObsInputs[hist_it]));
}
QLOG(SAM, "Added " << (FDObsInputs.size() - 1)
<< " far detector component spectra to form Observed "
"spectra for sample index "
<< FDSamples.size() << ".");
}
fds.SpectrumToFDSmearingMatrix_TH2 = dynamic_cast(FDObsInputs.back());
if (!fds.SpectrumToFDSmearingMatrix_TH2) {
ERROR(FTL,
"last histogram from ObsInput attribute was not a TH2D containing "
"the far detector ETrue/ERec smearing matrix ("
<< nk.GetS("ObsInput") << ").");
THROW(
"Far detector sample must contain the observed ERec spectrum and the "
"FD ETrue/ERec smearing matrix. "
"ObsInput=\"input.root[FDObs_species,FDSmearing_species]\"");
}
TMatrixD SpectrumToFDSmearingMatrix_l =
SmearceptanceUtils::GetMatrix(fds.SpectrumToFDSmearingMatrix_TH2);
fds.SpectrumToFDSmearingMatrix.ResizeTo(SpectrumToFDSmearingMatrix_l);
fds.SpectrumToFDSmearingMatrix = SpectrumToFDSmearingMatrix_l;
FDSamples.push_back(fds);
}
void Smear_SVDUnfold_Propagation_Osc::FinaliseFDSamples() {
std::stringstream ss("");
for (size_t fds_it = 0; fds_it < FDSamples.size(); ++fds_it) {
FDSample &fds = FDSamples[fds_it];
// Set up FD histograms.
// ==============================
for (size_t nd_it = 0; nd_it < NDSamples.size(); ++nd_it) {
NDSample &nds = NDSamples[nd_it];
TH1D *sampleHist =
static_cast(nds.ND_Unfolded_Spectrum_Hist->Clone());
sampleHist->Reset();
ss.str("");
ss << "FD_Propagated_Spectrum_Hist_" << fds_it << "_NDSample_" << nd_it;
sampleHist->SetName(ss.str().c_str());
fds.FD_Propagated_Spectrum_Hist_NDSamples.push_back(sampleHist);
}
fds.FD_Propagated_Spectrum_Hist = static_cast(
fds.FD_Propagated_Spectrum_Hist_NDSamples[0]->Clone());
fds.FD_Propagated_Spectrum_Hist->Reset();
ss.str("");
ss << "FD_Propagated_Spectrum_Hist_" << fds_it;
fds.FD_Propagated_Spectrum_Hist->SetName(ss.str().c_str());
for (size_t nd_it = 0; nd_it < NDSamples.size(); ++nd_it) {
- NDSample &nds = NDSamples[nd_it];
+ //NDSample &nds = NDSamples[nd_it];
TH1D *sampleHist =
static_cast(fds.FD_Propagated_Spectrum_Hist->Clone());
sampleHist->Reset();
ss.str("");
ss << "NDFD_Corrected_Spectrum_Hist_" << fds_it << "_NDSample_" << nd_it;
sampleHist->SetName(ss.str().c_str());
fds.NDFD_Corrected_Spectrum_Hist_NDSamples.push_back(sampleHist);
}
fds.NDFD_Corrected_Spectrum_Hist =
static_cast(fds.FD_Propagated_Spectrum_Hist->Clone());
fds.NDFD_Corrected_Spectrum_Hist->Reset();
ss.str("");
ss << "NDFD_Corrected_Spectrum_Hist_" << fds_it;
fds.NDFD_Corrected_Spectrum_Hist->SetName(ss.str().c_str());
fds.FD_Smeared_Spectrum_Hist =
new TH1D(*fds.SpectrumToFDSmearingMatrix_TH2->ProjectionX());
fds.FD_Smeared_Spectrum_Hist->SetDirectory(NULL);
fds.FD_Smeared_Spectrum_Hist->Reset();
ss.str("");
ss << "FD_Smeared_Spectrum_Hist_" << fds_it;
fds.FD_Smeared_Spectrum_Hist->SetName(ss.str().c_str());
for (size_t nd_it = 0; nd_it < NDSamples.size(); ++nd_it) {
NDSample &nds = NDSamples[nd_it];
if (!fds.FDNDRatios.count(nds.NuPDG)) {
ERROR(WRN, "Have an ND sample that provides PDG:"
<< nds.NuPDG
<< " neutrinos but far detector sample index " << fds_it
<< " doesn't have an NDFD ratio for this sample. "
"Setting to 0 (contribution from sample ignored.)");
fds.FDNDRatios[nds.NuPDG] = 0;
}
}
}
}
Int_t Smear_SVDUnfold_Propagation_Osc::GetFDSampleNAnalysisBins(size_t fds_it) {
if (fds_it >= FDSamples.size()) {
THROW("Requested FD sample index " << fds_it << " but only initialised "
<< FDSamples.size());
}
FDSample &fds = FDSamples[fds_it];
Int_t NAnalysisBins = 0;
for (Int_t bi_it = 1; bi_it < fds.FDDataHist->GetXaxis()->GetNbins() + 1;
++bi_it) {
if ((fds.FitRegion_Min != 0xdeadbeef) &&
(fds.FDDataHist->GetXaxis()->GetBinUpEdge(bi_it) <=
fds.FitRegion_Min)) {
continue;
}
if ((fds.FitRegion_Max != 0xdeadbeef) &&
(fds.FDDataHist->GetXaxis()->GetBinLowEdge(bi_it) >
fds.FitRegion_Max)) {
continue;
}
NAnalysisBins++;
}
return NAnalysisBins;
}
void Smear_SVDUnfold_Propagation_Osc::SetupChi2Hists() {
fDataHist =
new TH1D("SmearSVDUnfold", "", NFDAnalysisBins, 0, NFDAnalysisBins);
fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str());
Int_t CurrAnalysisBin = 1;
for (size_t fds_it = 0; fds_it < FDSamples.size(); ++fds_it) {
FDSample &fds = FDSamples[fds_it];
for (Int_t bi_it = 1; bi_it < fds.FDDataHist->GetXaxis()->GetNbins() + 1;
++bi_it) {
if ((fds.FitRegion_Min != 0xdeadbeef) &&
(fds.FDDataHist->GetXaxis()->GetBinUpEdge(bi_it) <=
fds.FitRegion_Min)) {
continue;
}
if ((fds.FitRegion_Max != 0xdeadbeef) &&
(fds.FDDataHist->GetXaxis()->GetBinLowEdge(bi_it) >
fds.FitRegion_Max)) {
continue;
}
fDataHist->SetBinContent(CurrAnalysisBin,
fds.FDDataHist->GetBinContent(bi_it));
if (UseRateErrors) {
fDataHist->SetBinError(CurrAnalysisBin,
sqrt(fds.FDDataHist->GetBinContent(bi_it)));
} else {
fDataHist->SetBinError(CurrAnalysisBin,
fds.FDDataHist->GetBinError(bi_it));
}
CurrAnalysisBin++;
}
}
fMCHist = static_cast(fDataHist->Clone());
fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(),
(fSettings.GetFullTitles()).c_str());
fMCHist->Reset();
}
void Smear_SVDUnfold_Propagation_Osc::UpdateChi2Hists() {
Int_t CurrAnalysisBin = 1;
for (size_t fds_it = 0; fds_it < FDSamples.size(); ++fds_it) {
FDSample &fds = FDSamples[fds_it];
for (Int_t bi_it = 1;
bi_it < fds.FD_Smeared_Spectrum_Hist->GetXaxis()->GetNbins() + 1;
++bi_it) {
if ((fds.FitRegion_Min != 0xdeadbeef) &&
(fds.FD_Smeared_Spectrum_Hist->GetXaxis()->GetBinUpEdge(bi_it) <=
fds.FitRegion_Min)) {
continue;
}
if ((fds.FitRegion_Max != 0xdeadbeef) &&
(fds.FD_Smeared_Spectrum_Hist->GetXaxis()->GetBinLowEdge(bi_it) >
fds.FitRegion_Max)) {
continue;
}
fMCHist->SetBinContent(
CurrAnalysisBin, fds.FD_Smeared_Spectrum_Hist->GetBinContent(bi_it));
fMCHist->SetBinError(CurrAnalysisBin,
fds.FD_Smeared_Spectrum_Hist->GetBinError(bi_it));
CurrAnalysisBin++;
}
}
}
//********************************************************************
Smear_SVDUnfold_Propagation_Osc::Smear_SVDUnfold_Propagation_Osc(
nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip =
"Simple measurement class for doing fake data oscillation studies.\n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("Osc Studies");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("XXX");
fSettings.SetYTitle("Number of events");
fSettings.SetEnuRange(0.0, 50);
fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG");
fSettings.DefineAllowedTargets("*");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") / (fNEvents + 0.);
fMCHist = NULL;
fDataHist = NULL;
ReadExtraConfig(samplekey);
AddNDInputs(samplekey);
std::vector NDTargets = samplekey.GetListOfChildNodes("NDObs");
for (size_t nd_it = 0; nd_it < NDTargets.size(); ++nd_it) {
AddNDInputs(NDTargets[nd_it]);
}
std::vector FDTargets = samplekey.GetListOfChildNodes("FDObs");
NFDAnalysisBins = 0;
if (!FDTargets.size()) { // If no child elements, assume that everything is
// contained on the sample element
AddFDTarget(samplekey);
} else {
for (size_t fd_it = 0; fd_it < FDTargets.size(); ++fd_it) {
AddFDTarget(FDTargets[fd_it]);
NFDAnalysisBins += GetFDSampleNAnalysisBins(fd_it);
}
}
if ((FitRegion_Min != 0xdeadbeef) || (FitRegion_Max != 0xdeadbeef)) {
QLOG(SAM, "When unfolding, interested region limited to:");
if (FitRegion_Min != 0xdeadbeef) {
QLOG(SAM, "\tE_nu > " << FitRegion_Min);
}
if (FitRegion_Max != 0xdeadbeef) {
QLOG(SAM, "\tE_nu < " << FitRegion_Max);
}
}
SetupNDInputs();
FinaliseFDSamples();
QLOG(SAM, "Set up " << FDSamples.size()
<< " samples for oscillation analysis with "
<< NFDAnalysisBins << " analysis bins.");
SetupChi2Hists();
SetCovarFromDiagonal();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void Smear_SVDUnfold_Propagation_Osc::FillEventVariables(FitEvent *event){};
bool Smear_SVDUnfold_Propagation_Osc::isSignal(FitEvent *event) {
return false;
}
void DumpUnfoldDebugInfo(Smear_SVDUnfold_Propagation_Osc::NDSample &nds,
size_t nd_it, size_t truncations) {
TDirectory *ogDir = gDirectory;
std::stringstream ss;
ss.str("");
ss << "DEBUG_FailedInvert_NDSample_" << nd_it;
Config::Get().out->mkdir(ss.str().c_str());
Config::Get().out->cd(ss.str().c_str());
ss.str("");
ss << "ND_Smearing_Matrix_" << nd_it;
nds.NDToSpectrumSmearingMatrix->Write(ss.str().c_str(), TObject::kOverwrite);
ss.str("");
ss << "ND_Inverse_Smearing_Matrix_" << nd_it;
SmearceptanceUtils::SVDGetInverse(nds.NDToSpectrumSmearingMatrix, 0)
->Write(ss.str().c_str(), TObject::kOverwrite);
ss.str("");
ss << "ND_Inverse_Smearing_Matrix_" << nd_it << "_Trunc_" << truncations;
SmearceptanceUtils::SVDGetInverse(nds.NDToSpectrumSmearingMatrix, truncations)
->Write(ss.str().c_str(), TObject::kOverwrite);
ss.str("");
ss << "ND_Obs_" << nd_it;
nds.NDDataHist->Write(ss.str().c_str(), TObject::kOverwrite);
TMatrixD NDToSpectrumResponseMatrix_notrunc = SmearceptanceUtils::GetMatrix(
SmearceptanceUtils::SVDGetInverse(nds.NDToSpectrumSmearingMatrix, 0));
nds.NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_notrunc);
nds.NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_notrunc;
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
nds.NDDataHist, nds.ND_Unfolded_Spectrum_Hist,
nds.NDToSpectrumResponseMatrix, 1000, false);
ss.str("");
ss << "ND_Unfolded_" << nd_it;
nds.ND_Unfolded_Spectrum_Hist->Write(ss.str().c_str(), TObject::kOverwrite);
TMatrixD NDToSpectrumResponseMatrix_trunc =
SmearceptanceUtils::GetMatrix(SmearceptanceUtils::SVDGetInverse(
nds.NDToSpectrumSmearingMatrix, truncations));
nds.NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_trunc);
nds.NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_trunc;
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
nds.NDDataHist, nds.ND_Unfolded_Spectrum_Hist,
nds.NDToSpectrumResponseMatrix, 1000, false);
ss.str("");
ss << "ND_Unfolded_" << nd_it << "_Trunc_" << truncations;
nds.ND_Unfolded_Spectrum_Hist->Write(ss.str().c_str(), TObject::kOverwrite);
if (ogDir) {
ogDir->cd();
} else {
Config::Get().out->cd();
}
}
void Smear_SVDUnfold_Propagation_Osc::UnfoldToNDETrueSpectrum(size_t nd_it) {
if (nd_it >= NDSamples.size()) {
THROW("Attempting to unfold ND sample index "
<< nd_it << " but only have " << NDSamples.size() << " ND samples.");
}
NDSample &nds = NDSamples[nd_it];
nds.ND_Unfolded_Spectrum_Hist =
new TH1D(*nds.NDToSpectrumSmearingMatrix->ProjectionY());
nds.ND_Unfolded_Spectrum_Hist->SetDirectory(NULL);
nds.ND_Unfolded_Spectrum_Hist->Clear();
std::stringstream ss("");
ss << "ND_Unfolded_Spectrum_Hist_" << nd_it;
nds.ND_Unfolded_Spectrum_Hist->SetName(ss.str().c_str());
bool HasNegValue = false;
int truncations = nds.TruncateStart;
do {
if (truncations >= nds.TruncateUpTo) {
DumpUnfoldDebugInfo(nds, nd_it, truncations);
THROW("Unfolded enu spectrum had negative values even after "
<< truncations << " SVD singular value truncations.");
}
// Unfold ND ERec -> Enu spectrum
// ------------------------------
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
nds.NDDataHist, nds.ND_Unfolded_Spectrum_Hist,
nds.NDToSpectrumResponseMatrix, 1000, false);
HasNegValue = false;
for (Int_t bi_it = 1;
bi_it < nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetNbins() + 1;
++bi_it) {
if ((FitRegion_Min != 0xdeadbeef) &&
(nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinUpEdge(bi_it) <=
FitRegion_Min)) {
continue;
}
if ((FitRegion_Max != 0xdeadbeef) &&
(nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinLowEdge(bi_it) >
FitRegion_Max)) {
continue;
}
if (nds.ND_Unfolded_Spectrum_Hist->GetBinContent(bi_it) < 0) {
HasNegValue = true;
QLOG(SAM,
"After "
<< truncations << " truncations, bin " << (bi_it - 1) << " ["
<< nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinLowEdge(
bi_it)
<< " -- "
<< nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinUpEdge(
bi_it)
<< " ] has value: "
<< nds.ND_Unfolded_Spectrum_Hist->GetBinContent(bi_it));
break;
}
}
if (HasNegValue) {
TMatrixD NDToSpectrumResponseMatrix_l =
SmearceptanceUtils::GetMatrix(SmearceptanceUtils::SVDGetInverse(
nds.NDToSpectrumSmearingMatrix, truncations));
nds.NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l);
nds.NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l;
}
truncations++;
} while (HasNegValue);
}
void Smear_SVDUnfold_Propagation_Osc::PropagateFDSample(size_t fds_it) {
if (fds_it >= FDSamples.size()) {
THROW("Requested FD sample index " << fds_it << " but only initialised "
<< FDSamples.size());
}
FDSample &fds = FDSamples[fds_it];
// Apply Oscillations
// ------------------------------
FitWeight *fw = FitBase::GetRW();
OscWeightEngine *oscWE =
dynamic_cast(fw->GetRWEngine(kOSCILLATION));
if (!oscWE) {
THROW(
"Couldn't load oscillation weight engine for sample: "
"Smear_SVDUnfold_Propagation_Osc.");
}
for (Int_t bi_it = 1;
bi_it < fds.NDFD_Corrected_Spectrum_Hist->GetXaxis()->GetNbins() + 1;
++bi_it) {
double content_osc = 0;
double error_osc = 0;
double content_fdnd = 0;
double error_fdnd = 0;
// Oscillate each ND sample to FD neutrino
for (size_t nd_it = 0; nd_it < NDSamples.size(); ++nd_it) {
NDSample &nds = NDSamples[nd_it];
double oscWeight = oscWE->CalcWeight(
nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinCenter(bi_it),
nds.NuPDG, fds.OscillateToPDG);
double sample_content_osc =
nds.ND_Unfolded_Spectrum_Hist->GetBinContent(bi_it) * oscWeight;
double sample_error_osc =
nds.ND_Unfolded_Spectrum_Hist->GetBinError(bi_it) * oscWeight;
fds.FD_Propagated_Spectrum_Hist_NDSamples[nd_it]->SetBinContent(
bi_it, sample_content_osc);
fds.FD_Propagated_Spectrum_Hist_NDSamples[nd_it]->SetBinError(
bi_it, sample_error_osc);
double sample_content_fdnd =
sample_content_osc * fds.FDNDRatios[nds.NuPDG] * fds.FDNDMassRatio;
double sample_error_fdnd =
sample_error_osc * fds.FDNDRatios[nds.NuPDG] * fds.FDNDMassRatio;
fds.NDFD_Corrected_Spectrum_Hist_NDSamples[nd_it]->SetBinContent(
bi_it, sample_content_fdnd);
fds.NDFD_Corrected_Spectrum_Hist_NDSamples[nd_it]->SetBinError(
bi_it, sample_error_fdnd);
content_osc += sample_content_osc;
error_osc += sample_error_osc * sample_error_osc;
content_fdnd += sample_content_fdnd;
error_fdnd += sample_error_fdnd * sample_error_fdnd;
}
fds.FD_Propagated_Spectrum_Hist->SetBinContent(bi_it, content_osc);
fds.FD_Propagated_Spectrum_Hist->SetBinError(bi_it, sqrt(error_osc));
fds.NDFD_Corrected_Spectrum_Hist->SetBinContent(bi_it, content_fdnd);
fds.NDFD_Corrected_Spectrum_Hist->SetBinError(bi_it, sqrt(error_fdnd));
}
// Forward fold Spectrum -> ERec FD
// ------------------------------
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
fds.NDFD_Corrected_Spectrum_Hist, fds.FD_Smeared_Spectrum_Hist,
fds.SpectrumToFDSmearingMatrix, 1000, true);
}
void Smear_SVDUnfold_Propagation_Osc::ConvertEventRates(void) {
/// Get topology weights
/// for each ND sample
/// /// Build NDSample::NDToSpectrumSmearingMatrix from sum of topology smearing histos scaled by topology weights
///
/// SetupNDInputs();
///
/// for each FD sample
/// /// Build FDSample::SpectrumToFDSmearingMatrix_TH2 as above.
/// /// Convert to TMatrixD, FDSample::SpectrumToFDSmearingMatrix
for (size_t fds_it = 0; fds_it < FDSamples.size(); ++fds_it) {
PropagateFDSample(fds_it);
}
UpdateChi2Hists();
}
void Smear_SVDUnfold_Propagation_Osc::Write(std::string drawOpt) {
TDirectory *ogDir = gDirectory;
ConvertEventRates();
FitWeight *fw = FitBase::GetRW();
OscWeightEngine *oscWE =
dynamic_cast(fw->GetRWEngine(kOSCILLATION));
if (!oscWE) {
THROW(
"Couldn't load oscillation weight engine for sample: "
"Smear_SVDUnfold_Propagation_Osc.");
}
oscWE->Print();
// Write ND samples
//----------------
for (size_t nd_it = 0; nd_it < NDSamples.size(); ++nd_it) {
NDSample &nds = NDSamples[nd_it];
std::stringstream ss("");
ss << "Smear_SVDUnfold_Propagation_Osc_NDSample_" << nd_it;
if (ogDir) {
ogDir->mkdir(ss.str().c_str());
ogDir->cd(ss.str().c_str());
} else {
Config::Get().out->mkdir(ss.str().c_str());
Config::Get().out->cd(ss.str().c_str());
}
nds.NDToSpectrumSmearingMatrix->Write("SmearingMatrix_ND",
TObject::kOverwrite);
nds.NDDataHist->Write("Obs_ND", TObject::kOverwrite);
nds.ND_Unfolded_Spectrum_Hist->Write(
nds.ND_Unfolded_Spectrum_Hist->GetName(), TObject::kOverwrite);
if (ogDir) {
ogDir->cd();
} else {
Config::Get().out->cd();
}
}
// Write FD samples
//----------------
for (size_t fds_it = 0; fds_it < FDSamples.size(); ++fds_it) {
FDSample &fds = FDSamples[fds_it];
std::stringstream ss("");
ss << "Smear_SVDUnfold_Propagation_Osc_FDSample_" << fds_it;
if (ogDir) {
ogDir->mkdir(ss.str().c_str());
ogDir->cd(ss.str().c_str());
} else {
Config::Get().out->mkdir(ss.str().c_str());
Config::Get().out->cd(ss.str().c_str());
}
// Calc oscillation probability
// ------------------------------
FitWeight *fw = FitBase::GetRW();
OscWeightEngine *oscWE =
dynamic_cast(fw->GetRWEngine(kOSCILLATION));
if (!oscWE) {
THROW(
"Couldn't load oscillation weight engine for sample: "
"Smear_SVDUnfold_Propagation_Osc.");
}
for (size_t nd_it = 0; nd_it < NDSamples.size(); ++nd_it) {
NDSample &nds = NDSamples[nd_it];
ss.str("");
ss << "NDSample_" << nd_it << "_contribution";
fds.FD_Propagated_Spectrum_Hist_NDSamples[nd_it]->Write(
ss.str().c_str(), TObject::kOverwrite);
ss.str("");
ss << "NDSample_" << nd_it << "_FDNDCorrected_contribution";
fds.NDFD_Corrected_Spectrum_Hist_NDSamples[nd_it]->Write(
ss.str().c_str(), TObject::kOverwrite);
ss.str("");
ss << "NDSample_" << nd_it << "_OscillationProb";
TGraph POsc;
POsc.Set(1E4 - 1);
double min = nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinLowEdge(1);
double step =
(nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinUpEdge(
nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetNbins()) -
nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinLowEdge(1)) /
double(1E4);
for (size_t i = 1; i < 1E4; ++i) {
double enu = min + i * step;
double ow = oscWE->CalcWeight(enu, nds.NuPDG, fds.OscillateToPDG);
if (ow != ow) {
std::cout << "Bad osc weight for ENu: " << enu << std::endl;
}
POsc.SetPoint(i - 1, enu, ow);
}
POsc.Write(ss.str().c_str(), TObject::kOverwrite);
}
fds.FDDataHist->Write("Obs_FD", TObject::kOverwrite);
fds.FD_Propagated_Spectrum_Hist->Write(
fds.FD_Propagated_Spectrum_Hist->GetName(), TObject::kOverwrite);
fds.NDFD_Corrected_Spectrum_Hist->Write(
fds.NDFD_Corrected_Spectrum_Hist->GetName(), TObject::kOverwrite);
fds.SpectrumToFDSmearingMatrix_TH2->Write("SmearingMatrix_FD",
TObject::kOverwrite);
fds.FD_Smeared_Spectrum_Hist->Write(fds.FD_Smeared_Spectrum_Hist->GetName(),
TObject::kOverwrite);
if (ogDir) {
ogDir->cd();
} else {
Config::Get().out->cd();
}
}
fMCHist->Write("Pred_FD", TObject::kOverwrite);
fDataHist->Write("Obs_FD", TObject::kOverwrite);
Measurement1D::Write(drawOpt);
}
diff --git a/src/MCStudies/Smearceptance_Tester.cxx b/src/MCStudies/Smearceptance_Tester.cxx
index e19e54b..2904770 100644
--- a/src/MCStudies/Smearceptance_Tester.cxx
+++ b/src/MCStudies/Smearceptance_Tester.cxx
@@ -1,881 +1,881 @@
// 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 .
*******************************************************************************/
#include "Smearceptance_Tester.h"
#include "SmearceptanceUtils.h"
#include "Smearcepterton.h"
//#define DEBUG_SMEARTESTER
//********************************************************************
/// @brief Class to perform smearceptance MC Studies on a custom measurement
Smearceptance_Tester::Smearceptance_Tester(nuiskey samplekey) {
//********************************************************************
samplekey.Print();
// Sample overview ---------------------------------------------------
std::string descrip =
"Simple measurement class for producing an event summary tree of smeared "
"events.\n";
if (Config::HasPar("NPOT")) {
samplekey.SetS("NPOT", Config::GetParS("NPOT"));
}
if (Config::HasPar("FluxIntegralOverride")) {
samplekey.SetS("FluxIntegralOverride",
Config::GetParS("FluxIntegralOverride"));
}
if (Config::HasPar("TargetVolume")) {
samplekey.SetS("TargetVolume", Config::GetParS("TargetVolume"));
}
if (Config::HasPar("TargetMaterialDensity")) {
samplekey.SetS("TargetMaterialDensity",
Config::GetParS("TargetMaterialDensity"));
}
OutputSummaryTree = true;
if (Config::HasPar("smear.OutputSummaryTree")) {
OutputSummaryTree = Config::GetParI("smear.OutputSummaryTree");
}
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("Smearceptance Studies");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("XXX");
fSettings.SetYTitle("Number of events");
fSettings.SetEnuRange(0.0, 1E5);
fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG");
fSettings.DefineAllowedTargets("*");
fSettings.DefineAllowedSpecies("*");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor =
(GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) /
TotalIntegratedFlux();
// Measurement Details
std::vector splitName = GeneralUtils::ParseToStr(fName, "_");
- size_t firstUS = fName.find_first_of("_");
+ //size_t firstUS = fName.find_first_of("_");
std::string smearceptorName = samplekey.GetS("smearceptor");
QLOG(SAM, "Using smearceptor: " << smearceptorName
<< " (parsed from: " << fName << ").");
fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1);
SetupDefaultHist();
fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist);
covar = StatUtils::GetInvert(fFullCovar);
eventVariables = NULL;
QLOG(SAM, "Smearceptance Flux Scaling Factor = " << fScaleFactor);
if (fScaleFactor <= 0.0) {
ERROR(WRN, "SCALE FACTOR TOO LOW ");
sleep(20);
}
// Setup our TTrees
AddEventVariablesToTree();
smearceptor = &Smearcepterton::Get().GetSmearcepter(smearceptorName);
Int_t RecNBins = 20, TrueNBins = 20;
double RecBinL = 0xdeadbeef, TrueBinL = 0, RecBinH = 10, TrueBinH = 10;
if (Config::HasPar("smear.reconstructed.binning")) {
std::vector args = GeneralUtils::ParseToStr(
Config::GetParS("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::HasPar("smear.true.binning")) {
std::vector args =
GeneralUtils::ParseToStr(Config::GetParS("smear.true.binning"), ",");
TrueNBins = GeneralUtils::StrToInt(args[0]);
TrueBinL = GeneralUtils::StrToDbl(args[1]);
TrueBinH = GeneralUtils::StrToDbl(args[2]);
}
SVDTruncation = 0;
if (Config::HasPar("smear.true.binning")) {
SVDTruncation = Config::GetParI("smear.SVD.truncation");
QLOG(SAM, "Applying SVD truncation of: " << SVDTruncation)
}
ETrueDistrib = NULL;
ETrueDistrib_noweight = NULL;
ERecDistrib = NULL;
RecoSmear = NULL;
if (RecBinL != 0xdeadbeef) {
QLOG(SAM, "Using binning True: " << TrueNBins << ", [" << TrueBinL << " -- "
<< TrueBinH << "], Rec: " << RecNBins
<< ", [" << RecBinL << " -- " << RecBinH
<< "]");
ETrueDistrib = new TH1D("ELep_rate", ";True E_{#nu};Count", TrueNBins,
TrueBinL, TrueBinH);
ETrueDistrib_noweight =
new TH1D("ELep_rate_noweight", ";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", ";True E_{#nu};Recon. E_{#nu}", RecNBins,
RecBinL, RecBinH, TrueNBins, TrueBinL, TrueBinH);
RecoSmear->Sumw2();
}
// Final setup ---------------------------------------------------
FinaliseMeasurement();
}
void Smearceptance_Tester::AddEventVariablesToTree() {
if (OutputSummaryTree) {
// Setup the TTree to save everything
if (!eventVariables) {
Config::Get().out->cd();
eventVariables =
new TTree((fName + "_VARS").c_str(), (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_clep_true", &HMFS_clep_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_pi0_true", &HMFS_pi0_true);
eventVariables->Branch("HMFS_cK_true", &HMFS_cK_true);
eventVariables->Branch("HMFS_K0_true", &HMFS_K0_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_cK_true", &KEFSHad_cK_true,
"KEFSHad_cK_true/F");
eventVariables->Branch("KEFSHad_K0_true", &KEFSHad_K0_true,
"KEFSHad_K0_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("NcK_true", &NcK_true, "NcK_true/I");
eventVariables->Branch("NK0_true", &NK0_true, "NK0_true/I");
eventVariables->Branch("HMFS_clep_rec", &HMFS_clep_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_pi0_rec", &HMFS_pi0_rec);
eventVariables->Branch("HMFS_cK_rec", &HMFS_cK_rec);
eventVariables->Branch("HMFS_K0_rec", &HMFS_K0_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_cK_rec", &KEFSHad_cK_rec,
"KEFSHad_cK_rec/F");
eventVariables->Branch("KEFSHad_K0_rec", &KEFSHad_K0_rec,
"KEFSHad_K0_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_cK", &EFSVis_cK, "EFSVis_cK/F");
eventVariables->Branch("EFSVis_K0", &EFSVis_K0, "EFSVis_K0/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("NcK_seen", &NcK_seen, "NcK_seen/I");
eventVariables->Branch("NK0_seen", &NK0_seen, "NK0_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("NcK_contributed", &NcK_contributed,
"NcK_contributed/I");
eventVariables->Branch("NK0_contributed", &NK0_contributed,
"NK0_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");
xsecScaling = fScaleFactor;
eventVariables->Branch("xsecScaling", &xsecScaling, "xsecScaling/F");
eventVariables->Branch("flagCCINC_true", &flagCCINC_true,
"flagCCINC_true/O");
eventVariables->Branch("flagCC0K_true", &flagCC0K_true, "flagCC0K_true/O");
eventVariables->Branch("flagCC0Pi_true", &flagCC0Pi_true,
"flagCC0Pi_true/O");
eventVariables->Branch("flagCC1Pi_true", &flagCC1Pi_true,
"flagCC1Pi_true/O");
eventVariables->Branch("flagCCINC_rec", &flagCCINC_rec, "flagCCINC_rec/O");
eventVariables->Branch("flagCC0K_rec", &flagCC0K_rec, "flagCC0K_rec/O");
eventVariables->Branch("flagCC0Pi_rec", &flagCC0Pi_rec, "flagCC0Pi_rec/O");
eventVariables->Branch("flagCC1Pi_rec", &flagCC1Pi_rec, "flagCC1Pi_rec/O");
}
PredEvtRateWeight = 1;
if (fEvtRateScaleFactor != 0xdeadbeef) {
if (OutputSummaryTree) {
eventVariables->Branch("PredEvtRateWeight", &PredEvtRateWeight,
"PredEvtRateWeight/F");
}
PredEvtRateWeight = fScaleFactor * fEvtRateScaleFactor;
}
}
template
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
int CountNNotPdgsSeen(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
for (size_t p_it = 0; p_it < ri.RecObjClass.size(); ++p_it) {
if (!std::count(pdgs, pdgs + N, ri.RecObjClass[p_it])) {
sum++;
}
}
return sum;
}
template
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
int CountNNotPdgsContributed(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
for (size_t p_it = 0; p_it < ri.TrueContribPDGs.size(); ++p_it) {
if (!std::count(pdgs, pdgs + N, ri.TrueContribPDGs[p_it])) {
sum++;
}
}
return sum;
}
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.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
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
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
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
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) {
//********************************************************************
static int const cpipPDG[] = {211};
static int const cpimPDG[] = {-211};
static int const pi0PDG[] = {111};
static int const cKPDG[] = {321, -321};
static int const K0PDG[] = {311, 310, 130};
static int const ProtonPDG[] = {2212};
static int const NeutronPDG[] = {2112};
static int const GammaPDG[] = {22};
static int const CLeptonPDGs[] = {11, 13, 15, -11, -13, -15};
static int const ExplicitPDGs[] = {211, -211, 111, 321, -321, 311, 310, 130,
2212, 2112, 22, 11, 13, 15, 12, 14,
16, -11, -13, -15, -12, -14, -16};
RecoInfo *ri = smearceptor->Smearcept(event);
//** START Pions
HMFS_clep_true = TLorentzVector(0, 0, 0, 0);
HMFS_clep_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsCLep = event->GetHMFSParticle(CLeptonPDGs);
if (fsCLep) {
HMFS_clep_true = fsCLep->P4();
HMFS_clep_rec = GetHMFSRecParticles(*ri, fsCLep->PDG());
}
//** END Charged leptons
//** START Pions
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, fsPip->PDG());
}
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, fsPim->PDG());
}
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_pi0_true = TLorentzVector(0, 0, 0, 0);
HMFS_pi0_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsPi0 = event->GetHMFSPiZero();
if (fsPi0) {
HMFS_pi0_true = fsPi0->P4();
HMFS_pi0_rec = GetHMFSRecParticles(*ri, fsPi0->PDG());
}
//** END Pions
//** START Kaons
HMFS_cK_true = TLorentzVector(0, 0, 0, 0);
HMFS_cK_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fscK = event->GetHMFSParticle(cKPDG);
if (fscK) {
HMFS_cK_true = fscK->P4();
HMFS_cK_rec = GetHMFSRecParticles(*ri, fscK->PDG());
}
HMFS_K0_true = TLorentzVector(0, 0, 0, 0);
HMFS_K0_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsK0 = event->GetHMFSParticle(K0PDG);
if (fsK0) {
HMFS_K0_true = fsK0->P4();
HMFS_K0_rec = GetHMFSRecParticles(*ri, fsK0->PDG());
}
//** END Kaons
//** START Nucleons
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, fsP->PDG());
}
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_cK_true = FitUtils::SumTE_PartVect(event->GetAllFSParticle(cKPDG));
KEFSHad_K0_true = FitUtils::SumTE_PartVect(event->GetAllFSParticle(K0PDG));
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 +
KEFSHad_cK_true + KEFSHad_K0_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();
NcK_true = event->GetAllFSParticle(cKPDG).size();
NK0_true = event->GetAllFSParticle(K0PDG).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_cK_rec =
SumKE_RecoInfo(*ri, cKPDG, PhysConst::mass_cK * PhysConst::mass_MeV);
KEFSHad_K0_rec =
SumKE_RecoInfo(*ri, K0PDG, PhysConst::mass_K0 * 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;
TLorentzVector FSLepMom_rec(0, 0, 0, 0);
if (event->GetHMFSAnyLeptons()) {
FSLepMom_rec = GetHMFSRecParticles(*ri, event->GetHMFSAnyLeptons()->PDG());
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_cK = SumVisE_RecoInfo(*ri, cKPDG);
EFSVis_K0 = SumVisE_RecoInfo(*ri, K0PDG);
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_cpi + EFSVis_pi0 + EFSVis_p + EFSVis_n + EFSVis_gamma +
EFSVis_cK + EFSVis_K0;
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);
NcK_seen = CountNPdgsSeen(*ri, cKPDG);
NK0_seen = CountNPdgsSeen(*ri, K0PDG);
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);
NcK_contributed = CountNPdgsContributed(*ri, cKPDG);
NK0_contributed = CountNPdgsContributed(*ri, K0PDG);
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;
flagCCINC_true = PDGFSLep_true & 1;
flagCC0K_true = (NcK_true + NK0_true) == 0;
flagCC0Pi_true =
flagCCINC_true && flagCC0K_true && ((Ncpi_true + Npi0_true) == 0);
flagCC1Pi_true =
flagCCINC_true && flagCC0K_true && ((Ncpi_true + Npi0_true) == 1);
flagCCINC_rec = FSCLep_seen && flagCCINC_true;
flagCC0K_rec = (NcK_seen + NK0_seen) == 0;
flagCC0Pi_rec =
flagCCINC_rec && flagCC0K_rec && ((Ncpi_seen + Npi0_seen) == 0);
flagCC1Pi_rec =
flagCCINC_rec && flagCC0K_rec && ((Ncpi_seen + Npi0_seen) == 1);
if (OutputSummaryTree) {
// Fill the eventVariables Tree
eventVariables->Fill();
}
if (RecoSmear) {
RecoSmear->Fill(EISLep_true / 1000.0,
flagCCINC_rec ? EISLep_LepHadVis_rec / 1000.0 : -1, Weight);
ETrueDistrib_noweight->Fill(EISLep_true / 1000.0,
flagCCINC_true ? Weight : 0);
ETrueDistrib->Fill(EISLep_true / 1000.0,
flagCCINC_true ? Weight * PredEvtRateWeight : 0);
ERecDistrib->Fill(EISLep_LepHadVis_rec / 1000.0,
flagCCINC_rec ? Weight * PredEvtRateWeight : 0);
}
};
//********************************************************************
void Smearceptance_Tester::Write(std::string drawOpt) {
//********************************************************************
if (OutputSummaryTree) {
// First save the TTree
eventVariables->Write();
}
// Save Flux and Event Histograms too
GetInput()->GetFluxHistogram()->Write();
GetInput()->GetEventHistogram()->Write();
if (!RecoSmear) {
return;
}
TH2D *SmearMatrix_ev =
static_cast(RecoSmear->Clone("ELepHadVis_Smear_ev"));
for (Int_t trueAxis_it = 1;
trueAxis_it < RecoSmear->GetXaxis()->GetNbins() + 1; ++trueAxis_it) {
double NEISLep = ETrueDistrib_noweight->GetBinContent(trueAxis_it);
for (Int_t recoAxis_it = 1;
recoAxis_it < RecoSmear->GetYaxis()->GetNbins() + 1; ++recoAxis_it) {
if (NEISLep > std::numeric_limits::epsilon()) {
SmearMatrix_ev->SetBinContent(
trueAxis_it, recoAxis_it,
SmearMatrix_ev->GetBinContent(trueAxis_it, recoAxis_it) / NEISLep);
}
}
}
ETrueDistrib_noweight->Write();
ETrueDistrib->Write();
ERecDistrib->Write();
RecoSmear->Write();
SmearMatrix_ev->Write();
TH2D *ResponseMatrix_ev =
SmearceptanceUtils::SVDGetInverse(SmearMatrix_ev, SVDTruncation);
ResponseMatrix_ev->SetName("ResponseMatrix_ev");
ResponseMatrix_ev->Write();
#ifdef DEBUG_SMEARTESTER
TMatrixD SmearMatrix_ev_md = SmearceptanceUtils::GetMatrix(SmearMatrix_ev);
TH1D *SmearedEvt = static_cast(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
#ifdef __PROB3PP_ENABLED__
FitWeight *fw = FitBase::GetRW();
if (fw->HasRWEngine(kOSCILLATION)) {
OscWeightEngine *oscWE =
dynamic_cast(fw->GetRWEngine(kOSCILLATION));
TGraph POsc;
POsc.Set(1E4 - 1);
double min = ETrueDistrib->GetXaxis()->GetBinLowEdge(1);
double step = (ETrueDistrib->GetXaxis()->GetBinUpEdge(
ETrueDistrib->GetXaxis()->GetNbins()) -
ETrueDistrib->GetXaxis()->GetBinLowEdge(1)) /
double(1E4);
for (size_t i = 1; i < 1E4; ++i) {
double enu = min + i * step;
double ow = oscWE->CalcWeight(enu, 14);
if (ow != ow) {
std::cout << "Bad osc weight for ENu: " << enu << std::endl;
}
POsc.SetPoint(i - 1, enu, ow);
}
POsc.Write("POsc", TObject::kOverwrite);
}
#endif
TMatrixD ResponseMatrix_evt_md =
SmearceptanceUtils::GetMatrix(ResponseMatrix_ev);
TH1D *Unfolded_enu_obs = static_cast(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
fCurrentNorm = norm;
return;
}
//********************************************************************
void Smearceptance_Tester::FillHistograms() {
//********************************************************************
// No Histograms need filling........
return;
}
//********************************************************************
void Smearceptance_Tester::ResetAll() {
//********************************************************************
if (OutputSummaryTree) {
eventVariables->Reset();
}
return;
}
//********************************************************************
float Smearceptance_Tester::GetChi2() {
//********************************************************************
// No Likelihood to test, purely MC
return 0.0;
}