diff --git a/.gitignore b/.gitignore index 145d03de..09737d8e 100644 --- a/.gitignore +++ b/.gitignore @@ -24,4 +24,5 @@ gsl/ cmake-build-debug/ cmake-build/ *.sif -cmake-build-forprof/ \ No newline at end of file +cmake-build-forprof/ +cmake-build-valgrind/ \ No newline at end of file diff --git a/examples/MCMC.cpp b/examples/MCMC.cpp index 9079607a..379509ce 100644 --- a/examples/MCMC.cpp +++ b/examples/MCMC.cpp @@ -24,6 +24,8 @@ in these sorts of circumstances, which is valuable for us to know. #include #include #include +#include +#include // global consts constexpr double N_SIGNAL = 100.; // Truth signal rate in dataset constexpr double N_BACK = 300.; // Truth background rate in dataset @@ -31,7 +33,8 @@ constexpr double SIGMA_BACK = 50.; // Constraint uncertainty on background norma constexpr double SIGMA_ESCALE = 0.01; constexpr double SIGMA_ESMEAR = 0.05; constexpr double SIGMA_RSCALE = 0.01; -constexpr size_t N_STEPS = 10000; +constexpr double SIGMA_RSMEAR = 0.02; +constexpr size_t N_STEPS = 100; const std::string outfilename_root = "mcmc_example_output.root"; std::vector load_mc(const AxisCollection& ax, const std::vector& observables) { @@ -118,6 +121,7 @@ void add_systematics(BinnedNLLH& lh_function) { * - Energy scale * - Energy smear * - Radial scale + * - Radial smear */ // 1: Add energy scale Scale* escale = new Scale("energy_scale"); @@ -131,21 +135,24 @@ void add_systematics(BinnedNLLH& lh_function) { lh_function.AddSystematic(escale); lh_function.SetConstraint("energy_scale_factor", 1.0, SIGMA_ESCALE); - // 2: Add energy smear (a bit complicated!) + // 2: Add energy smear (less complicated than it used to be!) // A Gaussian kernel... - VaryingCDF* smearer = new VaryingCDF("energy_smear"); - Gaussian* gaus = new Gaussian(0, 0.01, "e_gaus"); // temp sigma value - gaus->RenameParameter("means_0", "mean"); - gaus->RenameParameter("stddevs_0", "sigma"); - smearer->SetKernel(gaus); - // ...With a width scaling by the square root of the energy... - SquareRootScale* smear_sigma_func = new SquareRootScale("e_smear_sigma_func"); - smear_sigma_func->RenameParameter("grad", "energy_smear_param"); - smear_sigma_func->SetGradient(0.03); // temp gradient value - smearer->SetDependance("sigma", smear_sigma_func); - // ...which smears the PDFs along the energy axis. - Convolution* esmear = new Convolution("esmear"); - esmear->SetConditionalPDF(smearer); + // VaryingCDF* smearer = new VaryingCDF("energy_smear"); + // Gaussian* gaus = new Gaussian(0, 0.01, "e_gaus"); // temp sigma value + // gaus->RenameParameter("means_0", "mean"); + // gaus->RenameParameter("stddevs_0", "sigma"); + // smearer->SetKernel(gaus); + // // ...With a width scaling by the square root of the energy... + // SquareRootScale* smear_sigma_func = new SquareRootScale("e_smear_sigma_func"); + // smear_sigma_func->RenameParameter("grad", "energy_smear_param"); + // smear_sigma_func->SetGradient(0.03); // temp gradient value + // smearer->SetDependance("sigma", smear_sigma_func); + // // ...which smears the PDFs along the energy axis. + // Convolution* esmear = new Convolution("esmear"); + // esmear->SetConditionalPDF(smearer); + + GaussianSqrtConvolution* esmear = new GaussianSqrtConvolution("esmear"); + esmear->RenameSigma("energy_smear_param"); esmear->SetAxes(lh_function.GetDataDist().GetAxes()); esmear->SetDistributionObs(ObsSet(std::vector({"energy", "r3"}))); esmear->SetTransformationObs(ObsSet("energy")); @@ -154,7 +161,7 @@ void add_systematics(BinnedNLLH& lh_function) { lh_function.AddSystematic(esmear); lh_function.SetConstraint("energy_smear_param", 0., SIGMA_ESMEAR); - // 3: Add radial smear + // 3: Add radial scale Scale* rscale = new Scale("radial_scale"); rscale->RenameParameter("scaleFactor", "radial_scale_factor"); rscale->SetScaleFactor(1.0); @@ -165,6 +172,17 @@ void add_systematics(BinnedNLLH& lh_function) { lh_function.AddSystematic(rscale); lh_function.SetConstraint("radial_scale_factor", 1.0, SIGMA_RSCALE); + + // 4: Add radial smear + GaussianConvolution* rsmear = new GaussianConvolution("radial_smear"); + rsmear->RenameSigma("radial_smear_sigma"); + rsmear->SetAxes(lh_function.GetDataDist().GetAxes()); + rsmear->SetDistributionObs(ObsSet(std::vector({"energy", "r3"}))); + rsmear->SetTransformationObs(ObsSet("r3")); + rsmear->Construct(); + + lh_function.AddSystematic(rsmear); + lh_function.SetConstraint("radial_smear_sigma", 0., SIGMA_RSMEAR); } std::pair calc_scale_bounds(BinnedNLLH& lh_function, @@ -222,6 +240,11 @@ void setup_metaparams(ParameterDict& minima, ParameterDict& maxima, maxima["radial_scale_factor"] = rscale_bounds.second; initial_vals["radial_scale_factor"] = 1.0; initial_err["radial_scale_factor"] = SIGMA_RSCALE; + // radial_smear_sigma + minima["radial_smear_sigma"] = 0.; + maxima["radial_smear_sigma"] = 5.*SIGMA_RSMEAR; + initial_vals["radial_smear_sigma"] = 0.00001; + initial_err["radial_smear_sigma"] = SIGMA_RSMEAR; } TTree* run_mcmc(BinnedNLLH& lh_function, const ParameterDict& minima, const ParameterDict& maxima, diff --git a/src/function/Gaussian.cpp b/src/function/Gaussian.cpp index 907e3d16..282de54b 100644 --- a/src/function/Gaussian.cpp +++ b/src/function/Gaussian.cpp @@ -36,7 +36,7 @@ Gaussian::Gaussian(const std::vector &mean_, Initialise(mean_, stdDev_, name_); } -Gaussian::Gaussian(const Gaussian ©_) : fFitter(this, copy_.GetMeanNames(), copy_.GetStDevNames()) +Gaussian::Gaussian(const Gaussian ©_) : fFitter(this, copy_.GetMeanNames(), copy_.GetStDevNames(), copy_.GetHideMeanParameters()) { fMeans = copy_.fMeans; fStdDevs = copy_.fStdDevs; @@ -53,7 +53,7 @@ Gaussian::operator=(const Gaussian ©_) fCdfCutOff = copy_.fCdfCutOff; fNDims = copy_.fNDims; fName = std::string(copy_.fName + "_copy"); - fFitter = GaussianFitter(this, copy_.GetMeanNames(), copy_.GetStDevNames()); + fFitter = GaussianFitter(this, copy_.GetMeanNames(), copy_.GetStDevNames(), copy_.GetHideMeanParameters()); return *this; } @@ -171,6 +171,16 @@ Gaussian::GetNDims() const return fNDims; } +void Gaussian::HideMeanParameters() +{ + fFitter.HideMeanParameters(); +} + +bool Gaussian::GetHideMeanParameters() const +{ + return fFitter.GetHideMeanParameters(); +} + ///////////////// // Probability // ///////////////// diff --git a/src/function/Gaussian.h b/src/function/Gaussian.h index 3a30096b..2449899f 100644 --- a/src/function/Gaussian.h +++ b/src/function/Gaussian.h @@ -42,6 +42,9 @@ class Gaussian : public PDF void SetCdfCutOff(double); size_t GetNDims() const; + void HideMeanParameters(); + bool GetHideMeanParameters() const; + // Make this object fittable void SetParameter(const std::string &name_, double value); double GetParameter(const std::string &name_) const; diff --git a/src/function/GaussianFitter.cpp b/src/function/GaussianFitter.cpp index 2b96f137..04269d89 100644 --- a/src/function/GaussianFitter.cpp +++ b/src/function/GaussianFitter.cpp @@ -10,7 +10,7 @@ using ContainerTools::ToString; -GaussianFitter::GaussianFitter(Gaussian *gaus, const size_t &nDims_) +GaussianFitter::GaussianFitter(Gaussian *gaus, const size_t &nDims_) : fHideMeans(false) { fOrignalFunc = gaus; std::stringstream ss; @@ -25,7 +25,7 @@ GaussianFitter::GaussianFitter(Gaussian *gaus, const size_t &nDims_) } } -GaussianFitter::GaussianFitter(Gaussian *gaus, const std::vector &meanNames_, const std::vector &stdDevNames_) +GaussianFitter::GaussianFitter(Gaussian *gaus, const std::vector &meanNames_, const std::vector &stdDevNames_, bool hideMeans_) : fHideMeans(hideMeans_) { if (meanNames_.size() != stdDevNames_.size()) throw OXSXException(Formatter() << "GaussianFitter:: #meanName != #stdDevNames"); @@ -43,11 +43,14 @@ void GaussianFitter::RenameParameter(const std::string &old_, const std::string if (find(fMeansNames.begin(), fMeansNames.end(), old_) == fMeansNames.end() && find(fStdDevsNames.begin(), fStdDevsNames.end(), old_) == fStdDevsNames.end()) throw NotFoundError(Formatter() << "GaussianFitter:: When attempting to renaming the parameter " << old_ << ", it wasn't found. Available names: " << ToString(GetParameterNames())); - it = find(fMeansNames.begin(), fMeansNames.end(), old_); - while (it != fMeansNames.end()) + if (!fHideMeans) { - *it = new_; - it = find(it++, fMeansNames.end(), old_); + it = find(fMeansNames.begin(), fMeansNames.end(), old_); + while (it != fMeansNames.end()) + { + *it = new_; + it = find(it++, fMeansNames.end(), old_); + } } it = find(fStdDevsNames.begin(), fStdDevsNames.end(), old_); @@ -73,11 +76,14 @@ GaussianFitter::GetStdDevNames() const void GaussianFitter::SetParameter(const std::string &name_, double value_) { std::vector::iterator it; - it = find(fMeansNames.begin(), fMeansNames.end(), name_); - while (it != fMeansNames.end()) + if (!fHideMeans) { - fOrignalFunc->SetMean(it - fMeansNames.begin(), value_); - it = find(++it, fMeansNames.end(), name_); + it = find(fMeansNames.begin(), fMeansNames.end(), name_); + while (it != fMeansNames.end()) + { + fOrignalFunc->SetMean(it - fMeansNames.begin(), value_); + it = find(++it, fMeansNames.end(), name_); + } } it = find(fStdDevsNames.begin(), fStdDevsNames.end(), name_); while (it != fStdDevsNames.end()) @@ -96,7 +102,7 @@ GaussianFitter::GetParameter(const std::string &name_) const std::vector::const_iterator it; it = find(fMeansNames.begin(), fMeansNames.end(), name_); - if (it == fMeansNames.end()) + if (fHideMeans || it == fMeansNames.end()) { it = find(fStdDevsNames.begin(), fStdDevsNames.end(), name_); if (it == fStdDevsNames.end()) @@ -122,12 +128,12 @@ GaussianFitter::GetParameters() const std::vector values; values.reserve(means.size() + stddevs.size()); // preallocate memory - values.insert(values.end(), means.begin(), means.end()); + if (!fHideMeans) { values.insert(values.end(), means.begin(), means.end()); } values.insert(values.end(), stddevs.begin(), stddevs.end()); std::vector names; names.reserve(fMeansNames.size() + fStdDevsNames.size()); // preallocate memory - names.insert(names.end(), fMeansNames.begin(), fMeansNames.end()); + if (!fHideMeans) { names.insert(names.end(), fMeansNames.begin(), fMeansNames.end()); } names.insert(names.end(), fStdDevsNames.begin(), fStdDevsNames.end()); return ContainerTools::CreateMap(names, values); @@ -136,7 +142,7 @@ GaussianFitter::GetParameters() const size_t GaussianFitter::GetParameterCount() const { - return fMeansNames.size() + fStdDevsNames.size(); + return fHideMeans ? fStdDevsNames.size() : fMeansNames.size() + fStdDevsNames.size(); } std::set @@ -145,7 +151,7 @@ GaussianFitter::GetParameterNames() const std::set names; for (size_t i = 0; i < fMeansNames.size(); ++i) { - names.insert(fMeansNames.at(i)); + if (!fHideMeans) { names.insert(fMeansNames.at(i)); } names.insert(fStdDevsNames.at(i)); } return names; diff --git a/src/function/GaussianFitter.h b/src/function/GaussianFitter.h index 9b6378f3..05589b5c 100644 --- a/src/function/GaussianFitter.h +++ b/src/function/GaussianFitter.h @@ -12,7 +12,7 @@ class GaussianFitter { public: GaussianFitter(Gaussian *gaus, const size_t &nDims); - GaussianFitter(Gaussian *gaus, const std::vector &, const std::vector &); + GaussianFitter(Gaussian *gaus, const std::vector &, const std::vector &, bool hideMeans_); void SetParameter(const std::string &name_, double value); double GetParameter(const std::string &name_) const; @@ -27,9 +27,16 @@ class GaussianFitter std::vector GetMeanNames() const; std::vector GetStdDevNames() const; + // Option to hide the means from being observed as a parameter + // (Useful when we want to use the Gaussian as a kernel, and only + // want the fit to be able to modify the sigmas and not the means) + void HideMeanParameters() { fHideMeans = true; } + bool GetHideMeanParameters() const { return fHideMeans; } + private: Gaussian *fOrignalFunc; std::vector fMeansNames; std::vector fStdDevsNames; + bool fHideMeans; }; #endif diff --git a/src/systematic/CMakeLists.txt b/src/systematic/CMakeLists.txt index 04af0ee8..190ba47c 100644 --- a/src/systematic/CMakeLists.txt +++ b/src/systematic/CMakeLists.txt @@ -1,5 +1,7 @@ set(systematic_src systematic/Convolution.cpp + systematic/GaussianConvolution.cpp + systematic/GaussianSqrtConvolution.cpp systematic/Renorm.cpp systematic/Scale.cpp systematic/ScaleFunction.cpp @@ -11,6 +13,8 @@ set(systematic_src set(systematic_headers Convolution.h + GaussianConvolution.h + GaussianSqrtConvolution.h Renorm.h Scale.h ScaleFunction.h diff --git a/src/systematic/Convolution.cpp b/src/systematic/Convolution.cpp index b22318f3..920e0d95 100644 --- a/src/systematic/Convolution.cpp +++ b/src/systematic/Convolution.cpp @@ -40,6 +40,7 @@ void Convolution::ConstructSubmatrix(std::vector &column column_indices.reserve(fSubMapAxes.GetNBins() * fSubMapAxes.GetNBins()); row_indices.reserve(fSubMapAxes.GetNBins() * fSubMapAxes.GetNBins()); vals.reserve(fSubMapAxes.GetNBins() * fSubMapAxes.GetNBins()); + // Loop over all entries of the sub-matrix to determine their values for (long long unsigned int origBin = 0; origBin < fSubMapAxes.GetNBins(); origBin++) { @@ -103,7 +104,7 @@ void Convolution::MakeFullMatrix(const std::vector &colu void Convolution::Construct() { /* - * Method that constructs the response matrix associated with this convlution systematic. + * Method that constructs the response matrix associated with this convolution systematic. * Attempts to be clever and speedy about building this, because this method can get computationally- * expensive, fast, when dealing with BinnedED objects with numerous dimensions and many bins! * diff --git a/src/systematic/Convolution.h b/src/systematic/Convolution.h index 5cc7a2d8..4b09e16a 100644 --- a/src/systematic/Convolution.h +++ b/src/systematic/Convolution.h @@ -39,7 +39,7 @@ class Convolution : public Systematic std::string GetName() const; void SetName(const std::string &); -private: +protected: ConditionalPDF *fDist; // kernel used in convolution std::string fName; // name of this object // Members needed for efficient calculations: diff --git a/src/systematic/GaussianConvolution.cpp b/src/systematic/GaussianConvolution.cpp new file mode 100644 index 00000000..8483626b --- /dev/null +++ b/src/systematic/GaussianConvolution.cpp @@ -0,0 +1,155 @@ +#include +#include +#include +#include + +GaussianConvolution::GaussianConvolution(const std::string &name_) : +Convolution(name_), fSigmaName("stddevs_0") +{ + // Set the kernel to be a Gaussian. + // SetFunction() clones this Gaussian object, packaging it within + // a JumpPDF object (a kind of ConditionalPDF object). + Gaussian gauss_tmp(1, ""); + gauss_tmp.HideMeanParameters(); // we don't want kernel's mean to be a fit parameter + SetFunction(&gauss_tmp); +} + +GaussianConvolution::GaussianConvolution(const std::string &name_, double cutoff_) : +Convolution(name_), fSigmaName("stddevs_0") +{ + // Set the kernel to be a Gaussian. + // SetFunction() clones this Gaussian object, packaging it within + // a JumpPDF object (a kind of ConditionalPDF object). + Gaussian gauss_tmp(1, ""); + gauss_tmp.SetCdfCutOff(cutoff_); // set a custom CDF cutoff + gauss_tmp.HideMeanParameters(); // we don't want kernel's mean to be a fit parameter + SetFunction(&gauss_tmp); +} + +void GaussianConvolution::ConstructSubmatrix(std::vector &column_indices, std::vector &row_indices, + std::vector &vals) const +{ + /* + * Construct the sub-matrix associated with the convolution in the transformed axes. + * Returns this information, by reference, in the form of column & row indices for non-zero entries of this submatrix, + * and their associated values. This form is what Armadillo prefers for sparse matrices! + */ + // variables storing the axes bin information + std::vector binCentres(fSubMapAxes.GetNDimensions()); + std::vector lowEdges(fSubMapAxes.GetNDimensions()); + std::vector highEdges(fSubMapAxes.GetNDimensions()); + // Pre-allocate memory for the sub-matrix data + // likely to need only a fraction of this memory! + column_indices.reserve(fSubMapAxes.GetNBins() * fSubMapAxes.GetNBins()); + row_indices.reserve(fSubMapAxes.GetNBins() * fSubMapAxes.GetNBins()); + vals.reserve(fSubMapAxes.GetNBins() * fSubMapAxes.GetNBins()); + + // Some tricks we are pulling here to make things go even faster: + // - In most use cases, the binning along the smearing axis is ~equal, + // so the amount of smearing will become translation-invariant! + // Cache calculations of the smearing Integral() (which is expensive) + // in a map object based on the integration endpoints (the output bin edges) + // After the first row, we expect ~all integrals to have already been calculated! + // - For typical smearing kernels, the contribution necessarily goes monotonically down + // for bins further away from the original bin. In other words, we expect this submatrix + // to be quite diagonal-heavy, with no random spikes far away. + // Also, beyond some distance we expect the contribution to go to zero. + // Instead of blithely scanning as usual over the full 2D grid of (origBin, destBin), + // we now start ~along the diagonal and work outwards - first forwards, then backwards. + // Whenever a zero is reached, we immediately exit the loop! + std::map, double> integral_cache; + // Loop over all entries of the sub-matrix to determine their values + for (long long unsigned int origBin = 0; origBin < fSubMapAxes.GetNBins(); origBin++) + { + // get the centre of the bin. Need to offset by this for a convolution + fSubMapAxes.GetBinCentres(origBin, binCentres); + + // loop over the bins it can be smeared into, going forwards + // from the bin after the diagonal first + for (long long unsigned int destBin = origBin+1; destBin < fSubMapAxes.GetNBins(); destBin++) + { + fSubMapAxes.GetBinLowEdges(destBin, lowEdges); + fSubMapAxes.GetBinHighEdges(destBin, highEdges); + + // Calculate destination bin edges along smearing axis, relative + // to origBin's centre + const double xlo = lowEdges.at(0) - binCentres.at(0); + const double xhi = lowEdges.at(0) - binCentres.at(0); + const std::pair edges {xlo, xhi}; + // Has an integral already been calculated for this pair of (relative) bin edges? + const auto it = std::find_if(integral_cache.begin(), integral_cache.end(), [edges](const std::pair, double>& p){ return p.first == edges; }); + double integral = 0.; + if (it == integral_cache.end()) + { + // Nope, need to calculate the integral (and add it to the cache) + integral = fDist->Integral(lowEdges, highEdges, binCentres); + integral_cache[edges] = integral; + } else + { + // Yep, use that result! + integral = it->second; + } + + // Only bother adding to matrix if non-zero! + // No point looking at any further destination bins in this loop + // as we know they'll be zero too + if (integral == 0.) + { + break; + } + column_indices.push_back(origBin); + row_indices.push_back(destBin); + vals.push_back(integral); + } + // Now do the same inner loop, but going backwards, starting from + // the diagonal - this ensures we don't accidentally try negative bins initially + for (long long unsigned int destBin = origBin; destBin < fSubMapAxes.GetNBins(); destBin--) + { + fSubMapAxes.GetBinLowEdges(destBin, lowEdges); + fSubMapAxes.GetBinHighEdges(destBin, highEdges); + + const double xlo = lowEdges.at(0) - binCentres.at(0); + const double xhi = lowEdges.at(0) - binCentres.at(0); + const std::pair edges {xlo, xhi}; + const auto it = std::find_if(integral_cache.begin(), integral_cache.end(), [edges](const std::pair, double>& p){ return p.first == edges; }); + double integral = 0.; + if (it == integral_cache.end()) + { + integral = fDist->Integral(lowEdges, highEdges, binCentres); + integral_cache[edges] = integral; + } else + { + integral = it->second; + } + + if (integral == 0.) + { + break; + } + column_indices.push_back(origBin); + row_indices.push_back(destBin); + vals.push_back(integral); + } + } +} + +double GaussianConvolution::GetSigma() const +{ + return fDist->GetParameter(fSigmaName); +} + +void GaussianConvolution::SetSigma(double sigma_) +{ + fDist->SetParameter(fSigmaName, sigma_); +} + +std::string GaussianConvolution::GetSigmaName() const +{ + return fSigmaName; +} + +void GaussianConvolution::RenameSigma(const std::string& newname_) +{ + fDist->RenameParameter(fSigmaName, newname_); + fSigmaName = newname_; +} \ No newline at end of file diff --git a/src/systematic/GaussianConvolution.h b/src/systematic/GaussianConvolution.h new file mode 100644 index 00000000..858ef96a --- /dev/null +++ b/src/systematic/GaussianConvolution.h @@ -0,0 +1,37 @@ +/*******************************************************/ +/* A specialisation of the Convolution class, for the */ +/* specific case of a Gaussian kernel. */ +/*******************************************************/ + +#ifndef __OXSX_GAUSSIANCONVOLUTION__ +#define __OXSX_GAUSSIANCONVOLUTION__ +#include + + +class GaussianConvolution : public Convolution +{ +public: + // Overwrite the allowed constructor + GaussianConvolution(const std::string &name); + // Alternative contructor, if you want a custom CDF cutoff: + // the number of sigmas away from the mean beyond which the CDF is set to 0 or 1 as appropriate. 5 by default. + GaussianConvolution(const std::string &name, double cutoff_); + // Getters/Setters for underlying Gaussian distribution's width + // (FitComponent interface will still work, this is just a quality-of-life feature) + double GetSigma() const; + void SetSigma(double sigma_); + // Get the Gaussian kernel's sigma parameter's name + std::string GetSigmaName() const; + // Rename the Gaussian kernel's sigma parameter's name + // (default: stddevs_0) + void RenameSigma(const std::string& newname_); + +private: + std::string fSigmaName; // tracks name of Gaussian's sigma parameter: used internally so user doesn't need to remember! + + void ConstructSubmatrix(std::vector &column_indices, std::vector &row_indices, + std::vector &vals) const; +}; + + +#endif \ No newline at end of file diff --git a/src/systematic/GaussianSqrtConvolution.cpp b/src/systematic/GaussianSqrtConvolution.cpp new file mode 100644 index 00000000..115a13a9 --- /dev/null +++ b/src/systematic/GaussianSqrtConvolution.cpp @@ -0,0 +1,168 @@ +#include +#include +#include +#include +#include +#include + +GaussianSqrtConvolution::GaussianSqrtConvolution(const std::string &name_) : +Convolution(name_), fSigmaName("grad") +{ + // Set the kernel to be a Gaussian, with square-root scaling + // Gaussian* gauss_tmp = new Gaussian(1, ""); + Gaussian gauss_tmp(1, ""); + + VaryingCDF* smearer = new VaryingCDF(""); + // VaryingCDF smearer(""); + smearer->SetKernel(&gauss_tmp); + + SquareRootScale* sqrt_scaler = new SquareRootScale("sqrt_scale"); + smearer->SetDependance("stddevs_0", sqrt_scaler); + + SetConditionalPDF(smearer); +} + +GaussianSqrtConvolution::GaussianSqrtConvolution(const std::string &name_, double cutoff_) : +Convolution(name_), fSigmaName("grad") +{ + // Set the kernel to be a Gaussian. + Gaussian gauss_tmp(1, ""); + gauss_tmp.SetCdfCutOff(cutoff_); // set a custom CDF cutoff + + VaryingCDF* smearer = new VaryingCDF(""); + // VaryingCDF smearer(""); + smearer->SetKernel(&gauss_tmp); + + SquareRootScale* fSqrtScaler = new SquareRootScale("sqrt_scale"); + smearer->SetDependance("stddevs_0", fSqrtScaler); + + SetConditionalPDF(smearer); +} + +void GaussianSqrtConvolution::ConstructSubmatrix(std::vector &column_indices, std::vector &row_indices, + std::vector &vals) const +{ + /* + * Construct the sub-matrix associated with the convolution in the transformed axes. + * Returns this information, by reference, in the form of column & row indices for non-zero entries of this submatrix, + * and their associated values. This form is what Armadillo prefers for sparse matrices! + */ + // variables storing the axes bin information + std::vector binCentres(fSubMapAxes.GetNDimensions()); + std::vector lowEdges(fSubMapAxes.GetNDimensions()); + std::vector highEdges(fSubMapAxes.GetNDimensions()); + // Pre-allocate memory for the sub-matrix data + // likely to need only a fraction of this memory! + column_indices.reserve(fSubMapAxes.GetNBins() * fSubMapAxes.GetNBins()); + row_indices.reserve(fSubMapAxes.GetNBins() * fSubMapAxes.GetNBins()); + vals.reserve(fSubMapAxes.GetNBins() * fSubMapAxes.GetNBins()); + + // Some tricks we are pulling here to make things go even faster: + // - In most use cases, the binning along the smearing axis is ~equal, + // so the amount of smearing will become translation-invariant! + // Cache calculations of the smearing Integral() (which is expensive) + // in a map object based on the integration endpoints (the output bin edges) + // After the first row, we expect ~all integrals to have already been calculated! + // - For typical smearing kernels, the contribution necessarily goes monotonically down + // for bins further away from the original bin. In other words, we expect this submatrix + // to be quite diagonal-heavy, with no random spikes far away. + // Also, beyond some distance we expect the contribution to go to zero. + // Instead of blithely scanning as usual over the full 2D grid of (origBin, destBin), + // we now start ~along the diagonal and work outwards - first forwards, then backwards. + // Whenever a zero is reached, we immediately exit the loop! + std::map, double> integral_cache; + // Loop over all entries of the sub-matrix to determine their values + for (long long unsigned int origBin = 0; origBin < fSubMapAxes.GetNBins(); origBin++) + { + // get the centre of the bin. Need to offset by this for a convolution + fSubMapAxes.GetBinCentres(origBin, binCentres); + + // loop over the bins it can be smeared into, going forwards + // from the bin after the diagonal first + for (long long unsigned int destBin = origBin+1; destBin < fSubMapAxes.GetNBins(); destBin++) + { + fSubMapAxes.GetBinLowEdges(destBin, lowEdges); + fSubMapAxes.GetBinHighEdges(destBin, highEdges); + + // Calculate destination bin edges along smearing axis, relative + // to origBin's centre + const double xlo = lowEdges.at(0) - binCentres.at(0); + const double xhi = lowEdges.at(0) - binCentres.at(0); + const std::pair edges {xlo, xhi}; + // Has an integral already been calculated for this pair of (relative) bin edges? + const auto it = std::find_if(integral_cache.begin(), integral_cache.end(), [edges](const std::pair, double>& p){ return p.first == edges; }); + double integral = 0.; + if (it == integral_cache.end()) + { + // Nope, need to calculate the integral (and add it to the cache) + integral = fDist->Integral(lowEdges, highEdges, binCentres); + integral_cache[edges] = integral; + } else + { + // Yep, use that result! + integral = it->second; + } + + // Only bother adding to matrix if non-zero! + // No point looking at any further destination bins in this loop + // as we know they'll be zero too + if (integral == 0.) + { + break; + } + column_indices.push_back(origBin); + row_indices.push_back(destBin); + vals.push_back(integral); + } + // Now do the same inner loop, but going backwards, starting from + // the diagonal - this ensures we don't accidentally try negative bins initially + for (long long unsigned int destBin = origBin; destBin < fSubMapAxes.GetNBins(); destBin--) + { + fSubMapAxes.GetBinLowEdges(destBin, lowEdges); + fSubMapAxes.GetBinHighEdges(destBin, highEdges); + + const double xlo = lowEdges.at(0) - binCentres.at(0); + const double xhi = lowEdges.at(0) - binCentres.at(0); + const std::pair edges {xlo, xhi}; + const auto it = std::find_if(integral_cache.begin(), integral_cache.end(), [edges](const std::pair, double>& p){ return p.first == edges; }); + double integral = 0.; + if (it == integral_cache.end()) + { + integral = fDist->Integral(lowEdges, highEdges, binCentres); + integral_cache[edges] = integral; + } else + { + integral = it->second; + } + + if (integral == 0.) + { + break; + } + column_indices.push_back(origBin); + row_indices.push_back(destBin); + vals.push_back(integral); + } + } +} + +double GaussianSqrtConvolution::GetSigma() const +{ + return fDist->GetParameter(fSigmaName); +} + +void GaussianSqrtConvolution::SetSigma(double sigma_) +{ + fDist->SetParameter(fSigmaName, sigma_); +} + +std::string GaussianSqrtConvolution::GetSigmaName() const +{ + return fSigmaName; +} + +void GaussianSqrtConvolution::RenameSigma(const std::string& newname_) +{ + fDist->RenameParameter(fSigmaName, newname_); + fSigmaName = newname_; +} \ No newline at end of file diff --git a/src/systematic/GaussianSqrtConvolution.h b/src/systematic/GaussianSqrtConvolution.h new file mode 100644 index 00000000..5d8efd6c --- /dev/null +++ b/src/systematic/GaussianSqrtConvolution.h @@ -0,0 +1,40 @@ +/****************************************************************/ +/* A specialisation of the Convolution class, for the */ +/* specific case of a Gaussian kernel with a width that */ +/* varies proportional to the square root of the smearing axis. */ +/* Useful in particular for smearing along in energy, say. */ +/****************************************************************/ + +#ifndef __OXSX_GAUSSIANSQRTCONVOLUTION__ +#define __OXSX_GAUSSIANSQRTCONVOLUTION__ +#include + + +class GaussianSqrtConvolution : public Convolution +{ +public: + // Overwrite the allowed constructor + GaussianSqrtConvolution(const std::string &name); + // Alternative contructor, if you want a custom CDF cutoff: + // the number of sigmas away from the mean beyond which the CDF is set to 0 or 1 as appropriate. 5 by default. + GaussianSqrtConvolution(const std::string &name, double cutoff_); + // Getters/Setters for underlying SquareRootScale's width parameter - + // this is the proportionality constant which multiplies the sqrt scaling + // (FitComponent interface will still work, this is just a quality-of-life feature) + double GetSigma() const; + void SetSigma(double sigma_); + // Get the kernel's width parameter's name + std::string GetSigmaName() const; + // Rename the kernel's width parameter's name + // (default: "grad") + void RenameSigma(const std::string& newname_); + +private: + std::string fSigmaName; // tracks name of kernel's width parameter: used internally so user doesn't need to remember! + + void ConstructSubmatrix(std::vector &column_indices, std::vector &row_indices, + std::vector &vals) const; +}; + + +#endif \ No newline at end of file diff --git a/src/systematic/Scale.cpp b/src/systematic/Scale.cpp index a83eb187..d4760395 100644 --- a/src/systematic/Scale.cpp +++ b/src/systematic/Scale.cpp @@ -28,6 +28,7 @@ void Scale::Construct() const BinAxis &scaleAxis = fAxes.GetAxis(fDistObs.GetIndex(scaleAxisName)); // Each (pre-scaled) bin gets a mapping from (post-scaled) bins to non-zero values std::vector> scale_vals(scaleAxis.GetNBins(), std::map{}); + size_t nvals = 0; for (size_t bin = 0; bin < scaleAxis.GetNBins(); bin++) { // First - work out edges of scaled bin interval along axis @@ -43,22 +44,42 @@ void Scale::Construct() if (contribution > 0) { scale_vals[bin][bin_test] = contribution; + nvals++; } } } - // Finally, construct the full response matrix by setting the diagonal + // Finally, construct the full response matrix by setting the non-zero // elements to their appropriate value, using the bin ID mapping to help. fResponse.SetZeros(); + // Set up variables first; pre-allocate memory for the large vectors associated with the full matrix's info + const size_t N = fAxes.GetNBins(); + const size_t n_sub = scaleAxis.GetNBins(); + const size_t n_blocks = N / n_sub; + const size_t size_block = nvals; + std::vector column_indices_bl; + std::vector row_indices_bl; + std::vector vals_bl; + column_indices_bl.reserve(size_block * n_blocks); + row_indices_bl.reserve(size_block * n_blocks); + vals_bl.reserve(size_block * n_blocks); + // Loop over all origin bins for (size_t obs_bin_id = 0; obs_bin_id < fAxes.GetNBins(); obs_bin_id++) { + // ...and also loop over all possible non-zero destination bins const size_t scaleAxisBin = fDistTransBinMapping.at(obs_bin_id); for (const auto &postScaleBinPair : scale_vals.at(scaleAxisBin)) { + // Put info about contribution for this bin into relevant vectors const size_t scaled_bin_id = fMappingDistAndTrans.GetComponent(obs_bin_id, postScaleBinPair.first); - fResponse.SetComponent(scaled_bin_id, obs_bin_id, postScaleBinPair.second); + row_indices_bl.push_back(scaled_bin_id); + column_indices_bl.push_back(obs_bin_id); + vals_bl.push_back(postScaleBinPair.second); } } + // Set elements of full response matrix all in one go using filled vectors + fResponse = SparseMatrix(N, N); + fResponse.SetComponents(row_indices_bl, column_indices_bl, vals_bl); } void Scale::SetScaleFactor(double scaleFactor_) diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt index 22cd1e81..bd7608ec 100644 --- a/test/unit/CMakeLists.txt +++ b/test/unit/CMakeLists.txt @@ -24,6 +24,8 @@ set(tests EventShiftTest.cpp EventSystematicManagerTest.cpp FitParameterTest.cpp + GaussianConvolutionTest.cpp + GaussianSqrtConvolutionTest.cpp GaussianFitterTest.cpp GaussianTest.cpp HistogramGetMultiDSlice.cpp diff --git a/test/unit/GaussianConvolutionTest.cpp b/test/unit/GaussianConvolutionTest.cpp new file mode 100644 index 00000000..3d8c38c7 --- /dev/null +++ b/test/unit/GaussianConvolutionTest.cpp @@ -0,0 +1,272 @@ +#include +#include +#include +#include +#include + +TEST_CASE("Simple GaussianConvolution systematic on 1d PDF") +{ + // First - build base 1D BinnedED object for applying systematic to + AxisCollection axes; + axes.AddAxis(BinAxis("axis0", 0, 4, 4)); + std::vector observables; + observables.push_back("obs0"); + + BinnedED pdf1("pdf1", axes); + pdf1.SetBinContent(1, 10); + pdf1.SetObservables(observables); + + // Now build the GaussianConvolution systematic + double sigma = 0.5; + + GaussianConvolution conv("smear_sys"); + conv.RenameSigma("sigma"); + conv.SetSigma(sigma); + conv.SetAxes(axes); + conv.SetTransformationObs(observables); + conv.SetDistributionObs(observables); + + SECTION("Check GetSigmaName() method") + { + REQUIRE(conv.GetSigmaName() == "sigma"); + } + + SECTION("Check FitComponent interface for GaussianConvolution, 1D") + { + REQUIRE(conv.GetParameterCount() == 1); // one param: Gaussian's sigma + REQUIRE(conv.GetParameter("sigma") == sigma); + REQUIRE(conv.GetName() == "smear_sys"); + + sigma = 0.4; + conv.SetParameter("sigma", sigma); + REQUIRE(conv.GetParameter("sigma") == sigma); + } + + // Add systematic & BinnedED objects to a systematic manager; test smearing + SystematicManager man; + + SECTION("Test 1D GaussianConvolution smearing") + { + man.Add(&conv); + man.AddDist(pdf1, ""); + man.Construct(); + + std::vector pdfs = {pdf1}; + std::vector OrignalPdfs(pdfs); + + man.DistortEDs(OrignalPdfs, pdfs); + + const double mu = 1.5; // centre of bin which has data in it + + std::vector modifiedObs = pdfs.at(0).GetBinContents(); + std::vector correctVals = { + 10. * (gsl_cdf_gaussian_P(1 - mu, sigma) - gsl_cdf_gaussian_P(0 - mu, sigma)), + 10. * (gsl_cdf_gaussian_P(2 - mu, sigma) - gsl_cdf_gaussian_P(1 - mu, sigma)), + 10. * (gsl_cdf_gaussian_P(3 - mu, sigma) - gsl_cdf_gaussian_P(2 - mu, sigma)), + 10. * (gsl_cdf_gaussian_P(4 - mu, sigma) - gsl_cdf_gaussian_P(3 - mu, sigma)), + }; + + REQUIRE(modifiedObs == correctVals); + } +} + +TEST_CASE("Simple GaussianConvolution systematic on 1d PDF, small CDF cutoff") +{ + // First - build base 1D BinnedED object for applying systematic to + AxisCollection axes; + axes.AddAxis(BinAxis("axis0", 0, 4, 4)); + std::vector observables; + observables.push_back("obs0"); + + BinnedED pdf1("pdf1", axes); + pdf1.SetBinContent(1, 10); + pdf1.SetObservables(observables); + + // Now build the GaussianConvolution systematic + double sigma = 0.5; + const double cdfcutoff = 1.; // comedically small cutoff, just to test + GaussianConvolution conv("smear_sys", cdfcutoff); + conv.RenameSigma("sigma"); + conv.SetSigma(sigma); + conv.SetAxes(axes); + conv.SetTransformationObs(observables); + conv.SetDistributionObs(observables); + + SECTION("Check GetSigmaName() method, with custom cutoff in place") + { + REQUIRE(conv.GetSigmaName() == "sigma"); + } + + SECTION("Check FitComponent interface for GaussianConvolution, 1D, with custom cutoff in place") + { + REQUIRE(conv.GetParameterCount() == 1); // one param: Gaussian's sigma + REQUIRE(conv.GetParameter("sigma") == sigma); + REQUIRE(conv.GetName() == "smear_sys"); + + sigma = 0.4; + conv.SetParameter("sigma", sigma); + REQUIRE(conv.GetParameter("sigma") == sigma); + } + + // Add systematic & BinnedED objects to a systematic manager; test smearing + SystematicManager man; + + SECTION("Test 1D GaussianConvolution smearing, with custom cutoff in place") + { + man.Add(&conv); + man.AddDist(pdf1, ""); + man.Construct(); + + std::vector pdfs = {pdf1}; + std::vector OrignalPdfs(pdfs); + + man.DistortEDs(OrignalPdfs, pdfs); + + const double mu = 1.5; // centre of bin which has data in it + + std::vector modifiedObs = pdfs.at(0).GetBinContents(); + // values different now because of the CDF + std::vector correctVals = { + 10. * (gsl_cdf_gaussian_P(1 - mu, sigma) - 0.), + 10. * (gsl_cdf_gaussian_P(2 - mu, sigma) - gsl_cdf_gaussian_P(1 - mu, sigma)), + 10. * (1. - gsl_cdf_gaussian_P(2 - mu, sigma)), + 10. * (0. - 0.), + }; + + REQUIRE(modifiedObs == correctVals); + } +} + +TEST_CASE("Simple GaussianConvolution systematic on 2d PDF") +{ + // First - build base 1D BinnedED object for applying systematic to + AxisCollection axes; + axes.AddAxis(BinAxis("axis0", 0, 4, 4)); + axes.AddAxis(BinAxis("axis1", 0, 2, 2)); + const std::vector observables = {"obs0", "obs1"}; + const std::vector obs_smear = {"obs0"}; + + BinnedED pdf1("pdf1", axes); + pdf1.SetBinContent(axes.FlattenIndices({1, 0}), 10); + pdf1.SetBinContent(axes.FlattenIndices({2, 1}), 20); + pdf1.SetObservables(observables); + + // Now build the Convolution systematic and related objects + double sigma = 0.5; + + GaussianConvolution conv("smear_sys"); + conv.RenameSigma("sigma"); + conv.SetSigma(sigma); + conv.SetAxes(axes); + conv.SetTransformationObs(obs_smear); + conv.SetDistributionObs(observables); + + SECTION("Check FitComponent interface for GaussianConvolution, 2D") + { + REQUIRE(conv.GetParameterCount() == 1); // one param: Gaussian's sigma + REQUIRE(conv.GetParameter("sigma") == sigma); + REQUIRE(conv.GetName() == "smear_sys"); + + sigma = 0.4; + conv.SetParameter("sigma", sigma); + REQUIRE(conv.GetParameter("sigma") == sigma); + } + + // Add systematic & BinnedED objects to a systematic manager; test smearing + SystematicManager man; + + SECTION("Test 2D GaussianConvolution smearing") + { + man.Add(&conv); + man.AddDist(pdf1, ""); + man.Construct(); + + std::vector pdfs = {pdf1}; + std::vector OrignalPdfs(pdfs); + + man.DistortEDs(OrignalPdfs, pdfs); + + const double mu1 = 1.5; // centres of bins which has data in it + const double mu2 = 2.5; // centres of bins which has data in it + + std::vector modifiedObs = pdfs.at(0).GetBinContents(); + std::vector correctVals(pdf1.GetNBins(), 0); + correctVals[axes.FlattenIndices({0, 0})] = 10. * (gsl_cdf_gaussian_P(1 - mu1, sigma) - gsl_cdf_gaussian_P(0 - mu1, sigma)); + correctVals[axes.FlattenIndices({1, 0})] = 10. * (gsl_cdf_gaussian_P(2 - mu1, sigma) - gsl_cdf_gaussian_P(1 - mu1, sigma)); + correctVals[axes.FlattenIndices({2, 0})] = 10. * (gsl_cdf_gaussian_P(3 - mu1, sigma) - gsl_cdf_gaussian_P(2 - mu1, sigma)); + correctVals[axes.FlattenIndices({3, 0})] = 10. * (gsl_cdf_gaussian_P(4 - mu1, sigma) - gsl_cdf_gaussian_P(3 - mu1, sigma)); + correctVals[axes.FlattenIndices({0, 1})] = 20. * (gsl_cdf_gaussian_P(1 - mu2, sigma) - gsl_cdf_gaussian_P(0 - mu2, sigma)); + correctVals[axes.FlattenIndices({1, 1})] = 20. * (gsl_cdf_gaussian_P(2 - mu2, sigma) - gsl_cdf_gaussian_P(1 - mu2, sigma)); + correctVals[axes.FlattenIndices({2, 1})] = 20. * (gsl_cdf_gaussian_P(3 - mu2, sigma) - gsl_cdf_gaussian_P(2 - mu2, sigma)); + correctVals[axes.FlattenIndices({3, 1})] = 20. * (gsl_cdf_gaussian_P(4 - mu2, sigma) - gsl_cdf_gaussian_P(3 - mu2, sigma)); + + REQUIRE(modifiedObs == correctVals); + } +} + +TEST_CASE("Simple GaussianConvolution systematic on 2d PDF, alt axis ordering") +{ + // First - build base 1D BinnedED object for applying systematic to + AxisCollection axes; + axes.AddAxis(BinAxis("axis1", 0, 2, 2)); + axes.AddAxis(BinAxis("axis0", 0, 4, 4)); + const std::vector observables = {"obs1", "obs0"}; + const std::vector obs_smear = {"obs0"}; + + BinnedED pdf1("pdf1", axes); + pdf1.SetBinContent(axes.FlattenIndices({0, 1}), 10); + pdf1.SetBinContent(axes.FlattenIndices({1, 2}), 20); + pdf1.SetObservables(observables); + + // Now build the Convolution systematic and related objects + double sigma = 0.5; + + GaussianConvolution conv("smear_sys"); + conv.RenameSigma("sigma"); + conv.SetSigma(sigma); + conv.SetAxes(axes); + conv.SetTransformationObs(obs_smear); + conv.SetDistributionObs(observables); + + SECTION("Check FitComponent interface for GaussianConvolution, 2D") + { + REQUIRE(conv.GetParameterCount() == 1); // one params: Gaussian's sigma + REQUIRE(conv.GetParameter("sigma") == sigma); + REQUIRE(conv.GetName() == "smear_sys"); + + sigma = 0.4; + conv.SetParameter("sigma", sigma); + REQUIRE(conv.GetParameter("sigma") == sigma); + } + + // Add systematic & BinnedED objects to a systematic manager; test smearing + SystematicManager man; + + SECTION("Test 2D GaussianConvolution smearing") + { + man.Add(&conv); + man.AddDist(pdf1, ""); + man.Construct(); + + std::vector pdfs = {pdf1}; + std::vector OrignalPdfs(pdfs); + + man.DistortEDs(OrignalPdfs, pdfs); + + const double mu1 = 1.5; // centres of bins which has data in it + const double mu2 = 2.5; // centres of bins which has data in it + + std::vector modifiedObs = pdfs.at(0).GetBinContents(); + std::vector correctVals(pdf1.GetNBins(), 0); + correctVals[axes.FlattenIndices({0, 0})] = 10. * (gsl_cdf_gaussian_P(1 - mu1, sigma) - gsl_cdf_gaussian_P(0 - mu1, sigma)); + correctVals[axes.FlattenIndices({0, 1})] = 10. * (gsl_cdf_gaussian_P(2 - mu1, sigma) - gsl_cdf_gaussian_P(1 - mu1, sigma)); + correctVals[axes.FlattenIndices({0, 2})] = 10. * (gsl_cdf_gaussian_P(3 - mu1, sigma) - gsl_cdf_gaussian_P(2 - mu1, sigma)); + correctVals[axes.FlattenIndices({0, 3})] = 10. * (gsl_cdf_gaussian_P(4 - mu1, sigma) - gsl_cdf_gaussian_P(3 - mu1, sigma)); + correctVals[axes.FlattenIndices({1, 0})] = 20. * (gsl_cdf_gaussian_P(1 - mu2, sigma) - gsl_cdf_gaussian_P(0 - mu2, sigma)); + correctVals[axes.FlattenIndices({1, 1})] = 20. * (gsl_cdf_gaussian_P(2 - mu2, sigma) - gsl_cdf_gaussian_P(1 - mu2, sigma)); + correctVals[axes.FlattenIndices({1, 2})] = 20. * (gsl_cdf_gaussian_P(3 - mu2, sigma) - gsl_cdf_gaussian_P(2 - mu2, sigma)); + correctVals[axes.FlattenIndices({1, 3})] = 20. * (gsl_cdf_gaussian_P(4 - mu2, sigma) - gsl_cdf_gaussian_P(3 - mu2, sigma)); + + REQUIRE(modifiedObs == correctVals); + } +} \ No newline at end of file diff --git a/test/unit/GaussianSqrtConvolutionTest.cpp b/test/unit/GaussianSqrtConvolutionTest.cpp new file mode 100644 index 00000000..87110f42 --- /dev/null +++ b/test/unit/GaussianSqrtConvolutionTest.cpp @@ -0,0 +1,274 @@ +#include +#include +#include +#include +#include + +TEST_CASE("Simple GaussianSqrtConvolution systematic on 1d PDF") +{ + // First - build base 1D BinnedED object for applying systematic to + AxisCollection axes; + axes.AddAxis(BinAxis("axis0", 0, 4, 4)); + std::vector observables; + observables.push_back("obs0"); + + BinnedED pdf1("pdf1", axes); + pdf1.SetBinContent(1, 10); + pdf1.SetObservables(observables); + + // Now build the GaussianSqrtConvolution systematic + double sigma = 0.5; + + GaussianSqrtConvolution conv("smear_sys"); + conv.RenameSigma("sigma"); + conv.SetSigma(sigma); + conv.SetAxes(axes); + conv.SetTransformationObs(observables); + conv.SetDistributionObs(observables); + + SECTION("Check GetSigmaName() method") + { + REQUIRE(conv.GetSigmaName() == "sigma"); + } + + SECTION("Check FitComponent interface for GaussianSqrtConvolution, 1D") + { + REQUIRE(conv.GetParameterCount() == 1); // one param: Gaussian's sigma + REQUIRE(conv.GetParameter("sigma") == sigma); + REQUIRE(conv.GetName() == "smear_sys"); + + sigma = 0.4; + conv.SetParameter("sigma", sigma); + REQUIRE(conv.GetParameter("sigma") == sigma); + } + + // Add systematic & BinnedED objects to a systematic manager; test smearing + SystematicManager man; + + SECTION("Test 1D GaussianSqrtConvolution smearing") + { + man.Add(&conv); + man.AddDist(pdf1, ""); + man.Construct(); + + std::vector pdfs = {pdf1}; + std::vector OrignalPdfs(pdfs); + + man.DistortEDs(OrignalPdfs, pdfs); + + const double mu = 1.5; // centre of bin which has data in it + const double s = sigma*std::sqrt(mu); + std::vector modifiedObs = pdfs.at(0).GetBinContents(); + std::vector correctVals = { + 10. * (gsl_cdf_gaussian_P(1 - mu, s) - gsl_cdf_gaussian_P(0 - mu, s)), + 10. * (gsl_cdf_gaussian_P(2 - mu, s) - gsl_cdf_gaussian_P(1 - mu, s)), + 10. * (gsl_cdf_gaussian_P(3 - mu, s) - gsl_cdf_gaussian_P(2 - mu, s)), + 10. * (gsl_cdf_gaussian_P(4 - mu, s) - gsl_cdf_gaussian_P(3 - mu, s)), + }; + + REQUIRE(modifiedObs == correctVals); + } +} + +TEST_CASE("Simple GaussianSqrtConvolution systematic on 1d PDF, small CDF cutoff") +{ + // First - build base 1D BinnedED object for applying systematic to + AxisCollection axes; + axes.AddAxis(BinAxis("axis0", 0, 4, 4)); + std::vector observables; + observables.push_back("obs0"); + + BinnedED pdf1("pdf1", axes); + pdf1.SetBinContent(1, 10); + pdf1.SetObservables(observables); + + // Now build the GaussianSqrtConvolution systematic + double sigma = 0.5; + const double cdfcutoff = 1.; // comedically small cutoff, just to test + GaussianSqrtConvolution conv("smear_sys", cdfcutoff); + conv.RenameSigma("sigma"); + conv.SetSigma(sigma); + conv.SetAxes(axes); + conv.SetTransformationObs(observables); + conv.SetDistributionObs(observables); + + SECTION("Check GetSigmaName() method, with custom cutoff in place") + { + REQUIRE(conv.GetSigmaName() == "sigma"); + } + + SECTION("Check FitComponent interface for GaussianSqrtConvolution, 1D, with custom cutoff in place") + { + REQUIRE(conv.GetParameterCount() == 1); // one param: Gaussian's sigma + REQUIRE(conv.GetParameter("sigma") == sigma); + REQUIRE(conv.GetName() == "smear_sys"); + + sigma = 0.4; + conv.SetParameter("sigma", sigma); + REQUIRE(conv.GetParameter("sigma") == sigma); + } + + // Add systematic & BinnedED objects to a systematic manager; test smearing + SystematicManager man; + + SECTION("Test 1D GaussianSqrtConvolution smearing, with custom cutoff in place") + { + man.Add(&conv); + man.AddDist(pdf1, ""); + man.Construct(); + + std::vector pdfs = {pdf1}; + std::vector OrignalPdfs(pdfs); + + man.DistortEDs(OrignalPdfs, pdfs); + + const double mu = 1.5; // centre of bin which has data in it + const double s = sigma*std::sqrt(mu); + std::vector modifiedObs = pdfs.at(0).GetBinContents(); + // values different now because of the CDF + std::vector correctVals = { + 10. * (gsl_cdf_gaussian_P(1 - mu, s) - 0.), + 10. * (gsl_cdf_gaussian_P(2 - mu, s) - gsl_cdf_gaussian_P(1 - mu, s)), + 10. * (1. - gsl_cdf_gaussian_P(2 - mu, s)), + 10. * (0. - 0.), + }; + + REQUIRE(modifiedObs == correctVals); + } +} + +TEST_CASE("Simple GaussianSqrtConvolution systematic on 2d PDF") +{ + // First - build base 1D BinnedED object for applying systematic to + AxisCollection axes; + axes.AddAxis(BinAxis("axis0", 0, 4, 4)); + axes.AddAxis(BinAxis("axis1", 0, 2, 2)); + const std::vector observables = {"obs0", "obs1"}; + const std::vector obs_smear = {"obs0"}; + + BinnedED pdf1("pdf1", axes); + pdf1.SetBinContent(axes.FlattenIndices({1, 0}), 10); + pdf1.SetBinContent(axes.FlattenIndices({2, 1}), 20); + pdf1.SetObservables(observables); + + // Now build the GaussianSqrtConvolution systematic and related objects + double sigma = 0.5; + + GaussianSqrtConvolution conv("smear_sys"); + conv.RenameSigma("sigma"); + conv.SetSigma(sigma); + conv.SetAxes(axes); + conv.SetTransformationObs(obs_smear); + conv.SetDistributionObs(observables); + + SECTION("Check FitComponent interface for GaussianSqrtConvolution, 2D") + { + REQUIRE(conv.GetParameterCount() == 1); // one param: Gaussian's sigma + REQUIRE(conv.GetParameter("sigma") == sigma); + REQUIRE(conv.GetName() == "smear_sys"); + + sigma = 0.4; + conv.SetParameter("sigma", sigma); + REQUIRE(conv.GetParameter("sigma") == sigma); + } + + // Add systematic & BinnedED objects to a systematic manager; test smearing + SystematicManager man; + + SECTION("Test 2D GaussianSqrtConvolution smearing") + { + man.Add(&conv); + man.AddDist(pdf1, ""); + man.Construct(); + + std::vector pdfs = {pdf1}; + std::vector OrignalPdfs(pdfs); + + man.DistortEDs(OrignalPdfs, pdfs); + + const double mu1 = 1.5; // centres of bins which has data in it + const double mu2 = 2.5; // centres of bins which has data in it + const double s1 = sigma*std::sqrt(mu1); + const double s2 = sigma*std::sqrt(mu2); + std::vector modifiedObs = pdfs.at(0).GetBinContents(); + std::vector correctVals(pdf1.GetNBins(), 0); + correctVals[axes.FlattenIndices({0, 0})] = 10. * (gsl_cdf_gaussian_P(1 - mu1, s1) - gsl_cdf_gaussian_P(0 - mu1, s1)); + correctVals[axes.FlattenIndices({1, 0})] = 10. * (gsl_cdf_gaussian_P(2 - mu1, s1) - gsl_cdf_gaussian_P(1 - mu1, s1)); + correctVals[axes.FlattenIndices({2, 0})] = 10. * (gsl_cdf_gaussian_P(3 - mu1, s1) - gsl_cdf_gaussian_P(2 - mu1, s1)); + correctVals[axes.FlattenIndices({3, 0})] = 10. * (gsl_cdf_gaussian_P(4 - mu1, s1) - gsl_cdf_gaussian_P(3 - mu1, s1)); + correctVals[axes.FlattenIndices({0, 1})] = 20. * (gsl_cdf_gaussian_P(1 - mu2, s2) - gsl_cdf_gaussian_P(0 - mu2, s2)); + correctVals[axes.FlattenIndices({1, 1})] = 20. * (gsl_cdf_gaussian_P(2 - mu2, s2) - gsl_cdf_gaussian_P(1 - mu2, s2)); + correctVals[axes.FlattenIndices({2, 1})] = 20. * (gsl_cdf_gaussian_P(3 - mu2, s2) - gsl_cdf_gaussian_P(2 - mu2, s2)); + correctVals[axes.FlattenIndices({3, 1})] = 20. * (gsl_cdf_gaussian_P(4 - mu2, s2) - gsl_cdf_gaussian_P(3 - mu2, s2)); + + REQUIRE(modifiedObs == correctVals); + } +} + +TEST_CASE("Simple GaussianSqrtConvolution systematic on 2d PDF, alt axis ordering") +{ + // First - build base 1D BinnedED object for applying systematic to + AxisCollection axes; + axes.AddAxis(BinAxis("axis1", 0, 2, 2)); + axes.AddAxis(BinAxis("axis0", 0, 4, 4)); + const std::vector observables = {"obs1", "obs0"}; + const std::vector obs_smear = {"obs0"}; + + BinnedED pdf1("pdf1", axes); + pdf1.SetBinContent(axes.FlattenIndices({0, 1}), 10); + pdf1.SetBinContent(axes.FlattenIndices({1, 2}), 20); + pdf1.SetObservables(observables); + + // Now build the GaussianSqrtConvolution systematic and related objects + double sigma = 0.5; + + GaussianSqrtConvolution conv("smear_sys"); + conv.RenameSigma("sigma"); + conv.SetSigma(sigma); + conv.SetAxes(axes); + conv.SetTransformationObs(obs_smear); + conv.SetDistributionObs(observables); + + SECTION("Check FitComponent interface for GaussianSqrtConvolution, 2D") + { + REQUIRE(conv.GetParameterCount() == 1); // one params: Gaussian's sigma + REQUIRE(conv.GetParameter("sigma") == sigma); + REQUIRE(conv.GetName() == "smear_sys"); + + sigma = 0.4; + conv.SetParameter("sigma", sigma); + REQUIRE(conv.GetParameter("sigma") == sigma); + } + + // Add systematic & BinnedED objects to a systematic manager; test smearing + SystematicManager man; + + SECTION("Test 2D GaussianSqrtConvolution smearing") + { + man.Add(&conv); + man.AddDist(pdf1, ""); + man.Construct(); + + std::vector pdfs = {pdf1}; + std::vector OrignalPdfs(pdfs); + + man.DistortEDs(OrignalPdfs, pdfs); + + const double mu1 = 1.5; // centres of bins which has data in it + const double mu2 = 2.5; // centres of bins which has data in it + const double s1 = sigma*std::sqrt(mu1); + const double s2 = sigma*std::sqrt(mu2); + std::vector modifiedObs = pdfs.at(0).GetBinContents(); + std::vector correctVals(pdf1.GetNBins(), 0); + correctVals[axes.FlattenIndices({0, 0})] = 10. * (gsl_cdf_gaussian_P(1 - mu1, s1) - gsl_cdf_gaussian_P(0 - mu1, s1)); + correctVals[axes.FlattenIndices({0, 1})] = 10. * (gsl_cdf_gaussian_P(2 - mu1, s1) - gsl_cdf_gaussian_P(1 - mu1, s1)); + correctVals[axes.FlattenIndices({0, 2})] = 10. * (gsl_cdf_gaussian_P(3 - mu1, s1) - gsl_cdf_gaussian_P(2 - mu1, s1)); + correctVals[axes.FlattenIndices({0, 3})] = 10. * (gsl_cdf_gaussian_P(4 - mu1, s1) - gsl_cdf_gaussian_P(3 - mu1, s1)); + correctVals[axes.FlattenIndices({1, 0})] = 20. * (gsl_cdf_gaussian_P(1 - mu2, s2) - gsl_cdf_gaussian_P(0 - mu2, s2)); + correctVals[axes.FlattenIndices({1, 1})] = 20. * (gsl_cdf_gaussian_P(2 - mu2, s2) - gsl_cdf_gaussian_P(1 - mu2, s2)); + correctVals[axes.FlattenIndices({1, 2})] = 20. * (gsl_cdf_gaussian_P(3 - mu2, s2) - gsl_cdf_gaussian_P(2 - mu2, s2)); + correctVals[axes.FlattenIndices({1, 3})] = 20. * (gsl_cdf_gaussian_P(4 - mu2, s2) - gsl_cdf_gaussian_P(3 - mu2, s2)); + + REQUIRE(modifiedObs == correctVals); + } +} \ No newline at end of file