diff --git a/src/MCStudies/OscStudies.md b/src/MCStudies/OscStudies.md
index 6dc1941..9cda387 100644
--- a/src/MCStudies/OscStudies.md
+++ b/src/MCStudies/OscStudies.md
@@ -1,150 +1,150 @@
# Oscillation studies in NUISANCE
This simple module allows the construction of simple NOvA-style unfold/propagate/refold/compare style neutrino oscillation studies in NUSIANCE. It is implemented as a standard NUSIANCE sample, although with signficantly more complicated markup than usual.
## Full example
```xml
-
+
-
+
-
+
```
## Details
### NUISANCE Parameters
Explained in detail elsewhere, but briefly repeated here for completeness.
A reweightable parameter is specified by an element like:
```xml
```
where the attributes are hopefully self explanatory. The state should be `FREE`
or `FIXED` depending on whether in a minimisation or chi2 scan the parameter
should be allowed to vary. The valid names for oscillation parameters currently
implemented through Prob3++ are:
* `sinsq_theta12`
* `sinsq_theta13`
* `sinsq_theta23`
* `dm12`
* `dm23`
* `dcp`
In general, all parameters should be specified if any are, as the default values
may not be sensible.
The baseline is calculated from the assumption of a beam passing through the
earth to be detected on the surface.
### Analysis structure
The oscillation analysis inputs are near detector observations, far detector
observations and neutrino energy migration matricies that encode how detector
effects cause interactions of some true neutrino energy to be reconstructed
with some other neutrino energy. As oscillations happen as a function of
neutrino energy, being able to corectly identify the energy spectrum of
neutrinos before and after oscillation is the most key aspect of an oscillation
analysis.
The obseved near detector reconstructed energy spectrum is unfolded through the
near detector migration matrix to give the data-motivated, estimated true energy
spectrum of interactions at the near detector. This is then oscillated according
to the current PMNS parameters and corrected for far/near beam divergence and
detector mass. The true energy spectrum at the far detector is then smeared with
the far detector migration matrix, which may or may not be the same as the near
detector one, to give the predicted far detector reconstructed energy spectrum.
A test statistic is evaluated between the predicted and observed far detector
spectrum for the current set of oscillation, and potentially nuisance,
parameters.
In general the unfolding can be used to estimate the true energy spectrum of
the observed interactions, the true energy spectrum of the interactions that
occured (including efficiency), or the true energy spectrum of the neutrino
flux (also including cross-section). As the far over near correction can
currently only be a single number to account for beam divergence, it is advised
that the supplied migration matrices should account for the migration from true
interactions to true energy spectra (i.e. include the efficiency
correction/application at the near and far detectors). Currently, the nue and
numu cross sections are assumed identical.
**N.B.:** The near and far detector samples are assigned a neutrino PDG code, *i.e.* they must represent a CCInclusive sample of a single neutrino flavor. There is no support for interaction model-fitting style (T2K) analyses.
### The Sample element
```xml
...
```
Every oscillation analysis needs at least one near and one far detector sample.
To keep the structure of NUISANCE sample elements, the first near
detector sample is specified slightly differently to subsequent samples.
The sample element specifies the near detector mass through the attributes `DetectorVolume` and `DetectorDensity`, the separation is purely for your bookkeeping, they must multiply to give a mass in the same units as the far detector masses that will be specified. This is purely used for near-to-far extrapolation, the input histograms are expected to be well-normalised event rate predictions in the relevant near and far detector configurations.
The `ScalePOT` attribute is applied to all input histograms as a re-normalisation factor, this allows event rate predictions to be recalculated by a simple scale factor at runtime.
The `SetErrorsFromRate` attribute specifies whether the input histogram errors should be respected or they should be rescaled to the poisson error on the fully normalised event rate (included the `ScalePOT` attributed) in each bin.
**N.B.:** Only a single near detector can be specified, it would be fairly trivial to add multiple near detector setups, but that is currently not implemented. However, multiple near detector CCInc samples may be declared.
### ND sample descriptors
```xml
```
The input histograms for each ND sample should be given in the order: [observed rate in ERec], [ERec:ETrue migration matrix]. These can come from a single ROOT file, as shown above, or from two separate ROOT files, specified like: `ObsInput="FILE1.root[ObsRateERecHistName],FILE2.root[ERecETrueMigrationMatrixTH2DName]"`. The first ND sample, which is specified in the sample element, as opposed to in a child element named `NDObs`, must be prefixed with `HISTO:` to satisfy NUISANCE's requirement for input types.
The `NuPDG` attribute specifies the neutrinp PDG code which is observed sample corresponds to, most likely one of `-12,12,-14,14`, depending on the analysis setup.
The `TruncateStart` and `TruncateUpTo` attributes specify the number of singular values that may be set to 0 in the SVD decoposition to ensure that the unfolded spectrum at the near detector is >= 0 within the union of specified `FitRegion`s in the far detector samples (see the next section). If there are still negative values after the maximum allowed regularisation has been attempted, the analysis will throw an exception.
### FD sample descriptors
```xml
-
+
```
Each far detector sample has the same inputs as each near detector sample, an observed histogram (correctly normalised) and a migration matrix. For the far detector samples, the migration matrix is not inverted, but used to forward-fold the true spectrum prediction at the far detector.
-The two elements `FitRegion_min` and `FitRegion_max` are used to specify the range of the observation histogram that is included in the Chi2 calculation. These ranges are also used to inform the regularisation of the SVD inversion used ont he near detector samples.
+The two elements `FitRegion_Min` and `FitRegion_Max` are used to specify the range of the observation histogram that is included in the Chi2 calculation. These ranges are also used to inform the regularisation of the SVD inversion used ont he near detector samples.
Each far detector sample receives contributions from every near detector sample, these may well be 0, but the extrapolation is performed independently for each far detector sample. Therefore, each far detector sample element must set up its detector mass, even if they are all the same. In the example shown above the far detector is set up to have a 40 kilotonne fiducial mass.
The `OscillateToPDG` attribtue functions similarly to the `NuPDG` attribute for the near detector samples and specifies what species of neutrino the sample characterises, muon neutrino, electron anti-neutrino, etc... This is used to calculate the correction oscillation probability from each of the near detector samples. *e.g.* For a electron neutrino far detector sample, the muon neutrino to electron neutrino appearance probability is applied during the propagation to a muon neutrino near detector sample, and the electron neutrino survival probability is applied to an intrinsic electron neutrino near detector sample.
To inform the correct propagation of each near detector sample, a child `FDNDRatio` tag must be supplied. This specifies the beam divergence factor for each near detector sample --- linked via its `NuPDG` --- this is applied as a simple scale factor during near-to-far propagation. It would be possible to allow this a functional form to account for the relative difference in shape of the unoscillated flux at the near and far detector, but this is not yet implemented. Any near detector samples which do not have a corresponding `FDNDRatio` tag are given a divergence factor of `0` and effectively do not contribute to the far detector prediction. *e.g.* In the example at the top of this note, the far detector muon neutrino disapeerance sample recieves no contribution from the intrinsic electron neutrino sample, this is trivial to enable, but makes a completely insignificant difference to the far detector prediction.
## Input generation details
The input histograms are simplest to prepare through the `nuissmear` application, which by default will output both that are required for each input.
The observed ERec distribution for each sample should be self-explanatory. The response matrix should be supplied as a TH2D, with the true energy on the X axis and the reconstructed energy on the Y axis. The generation shape should be divided out of this. There is a small choice to be made here, either the efficiency correction can be included as part of the unfolding, this would involve dividing every bin by the total number of interaction generated in that bin of true neutrino energy -- events that are not reconstructed and thus do not make it into the histogram are accounted for by the sum of each Y slice of bins not totalling unity. Alternatively, you could re-normalised each Y slice so that the sum of bins totals to unity and this would mean that the unfolding only accounted for the migration and not the detector/selection efficiency. This implicitly would include the assumption into any susequent analysis that the detectors are functionally identical. It would also be possible to apply the efficiency ratio during the propagation, but this can only be a flat scale calculated from the ratio of detector masses and the beam divergence factor, it is not currently possible.
diff --git a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx
index 7b6df60..8b11a3f 100644
--- a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx
+++ b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx
@@ -1,823 +1,913 @@
// 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];
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) {
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/Smear_SVDUnfold_Propagation_Osc.h b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.h
index 874b081..752888f 100644
--- a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.h
+++ b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.h
@@ -1,118 +1,119 @@
// 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 .
*******************************************************************************/
#ifndef Smear_SVDUnfold_Propagation_Osc_H_SEEN
#define Smear_SVDUnfold_Propagation_Osc_H_SEEN
#include "Measurement1D.h"
/// Class for building simple oscillation analyses.
///
/// An example tag which runs a FHC beam and performs a numudisp and nue app
/// measurement might look like:
///
///
///
///
///
///
class Smear_SVDUnfold_Propagation_Osc : public Measurement1D {
+public:
struct FDSample {
TH1D *FDDataHist;
std::vector FD_Propagated_Spectrum_Hist_NDSamples;
std::vector NDFD_Corrected_Spectrum_Hist_NDSamples;
TH1D *FD_Propagated_Spectrum_Hist;
TH1D *NDFD_Corrected_Spectrum_Hist;
TH1D *FD_Smeared_Spectrum_Hist;
TH2D *SpectrumToFDSmearingMatrix_TH2;
TMatrixD SpectrumToFDSmearingMatrix;
Int_t OscillateToPDG;
double FitRegion_Min;
double FitRegion_Max;
std::map FDNDRatios;
double FDNDMassRatio;
std::pair DetectorInfo;
};
struct NDSample {
TH1D *NDDataHist;
Int_t TruncateStart;
Int_t TruncateUpTo;
TH1D *ND_Unfolded_Spectrum_Hist;
TH2D *NDToSpectrumSmearingMatrix;
TMatrixD NDToSpectrumResponseMatrix;
Int_t NuPDG;
double XSecToEvRateScale;
};
public:
Smear_SVDUnfold_Propagation_Osc(nuiskey samplekey);
virtual ~Smear_SVDUnfold_Propagation_Osc(){};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
void UnfoldToNDETrueSpectrum(size_t);
void ConvertEventRates(void);
void Write(std::string drawOpt);
void ReadExtraConfig(nuiskey &samplekey);
void SetupChi2Hists();
void UpdateChi2Hists();
// ----- Near detector things
void AddNDInputs(nuiskey &samplekey);
void SetupNDInputs();
std::vector NDSamples;
std::pair NDetectorInfo;
double FitRegion_Min;
double FitRegion_Max;
// ----- Far detector things
Int_t NFDAnalysisBins;
Int_t GetFDSampleNAnalysisBins(size_t fds_it);
void AddFDTarget(nuiskey &nk);
void PropagateFDSample(size_t fds_it);
void FinaliseFDSamples();
std::vector FDSamples;
double ScalePOT;
// ----- Other config
bool UseRateErrors;
};
#endif