New GaussianConvolution & GaussianSqrtConvolution classes; Update to Scale systematic

.gitignore

@@ -24,4 +24,5 @@ gsl/
 cmake-build-debug/
 cmake-build/
 *.sif
-cmake-build-forprof/
+cmake-build-forprof/
+cmake-build-valgrind/

examples/MCMC.cpp

@@ -24,14 +24,17 @@ in these sorts of circumstances, which is valuable for us to know.
 #include <VaryingCDF.h>
 #include <SquareRootScale.h>
 #include <Convolution.h>
+#include <GaussianConvolution.h>
+#include <GaussianSqrtConvolution.h>
 // global consts
 constexpr double N_SIGNAL = 100.; // Truth signal rate in dataset
 constexpr double N_BACK = 300.; // Truth background rate in dataset
 constexpr double SIGMA_BACK = 50.; // Constraint uncertainty on background normalisation
 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<BinnedED> load_mc(const AxisCollection& ax, const std::vector<std::string>& 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<std::string>({"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<std::string>({"energy", "r3"})));
+    rsmear->SetTransformationObs(ObsSet("r3"));
+    rsmear->Construct();
+
+    lh_function.AddSystematic(rsmear);
+    lh_function.SetConstraint("radial_smear_sigma", 0., SIGMA_RSMEAR);
 }
 
 std::pair<double, double> 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,

src/function/Gaussian.cpp

@@ -36,7 +36,7 @@ Gaussian::Gaussian(const std::vector<double> &mean_,
     Initialise(mean_, stdDev_, name_);
 }
 
-Gaussian::Gaussian(const Gaussian &copy_) : fFitter(this, copy_.GetMeanNames(), copy_.GetStDevNames())
+Gaussian::Gaussian(const Gaussian &copy_) : fFitter(this, copy_.GetMeanNames(), copy_.GetStDevNames(), copy_.GetHideMeanParameters())
 {
     fMeans = copy_.fMeans;
     fStdDevs = copy_.fStdDevs;
@@ -53,7 +53,7 @@ Gaussian::operator=(const Gaussian &copy_)
     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 //
 /////////////////

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;

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<std::string> &meanNames_, const std::vector<std::string> &stdDevNames_)
+GaussianFitter::GaussianFitter(Gaussian *gaus, const std::vector<std::string> &meanNames_, const std::vector<std::string> &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<std::string>::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<std::string>::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<double> 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<std::string> 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<std::string>
@@ -145,7 +151,7 @@ GaussianFitter::GetParameterNames() const
     std::set<std::string> 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;

src/function/GaussianFitter.h

@@ -12,7 +12,7 @@ class GaussianFitter
 {
 public:
     GaussianFitter(Gaussian *gaus, const size_t &nDims);
-    GaussianFitter(Gaussian *gaus, const std::vector<std::string> &, const std::vector<std::string> &);
+    GaussianFitter(Gaussian *gaus, const std::vector<std::string> &, const std::vector<std::string> &, 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<std::string> GetMeanNames() const;
     std::vector<std::string> 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<std::string> fMeansNames;
     std::vector<std::string> fStdDevsNames;
+    bool fHideMeans;
 };
 #endif

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

src/systematic/Convolution.cpp

@@ -40,6 +40,7 @@ void Convolution::ConstructSubmatrix(std::vector<long long unsigned int> &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<long long unsigned int> &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!
      *

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: