From f13877664ac2ee358e449669427f510398e4ebfe Mon Sep 17 00:00:00 2001
From: Andrey Popov <rainfinder@gmail.com>
Date: Thu, 23 Jun 2016 14:14:55 +0200
Subject: [PATCH] Compute systematic variation in plugin LeptonSFWeight

All individual uncertainties are added up in quadrature.
---
 include/mensura/extensions/LeptonSFWeight.hpp |  3 ++-
 modules/extensions/src/LeptonSFWeight.cpp     | 18 ++++++++++++++----
 2 files changed, 16 insertions(+), 5 deletions(-)

diff --git a/include/mensura/extensions/LeptonSFWeight.hpp b/include/mensura/extensions/LeptonSFWeight.hpp
index e885c92..dd972ae 100644
--- a/include/mensura/extensions/LeptonSFWeight.hpp
+++ b/include/mensura/extensions/LeptonSFWeight.hpp
@@ -29,7 +29,8 @@ class LeptonReader;
  * "absEta", "absEtaSC". It no parameter string is found, transverse momentum and (signed)
  * pseudorapidity of the lepton are taken parameters of the histogram.
  * 
- * Currently systematic uncertainty is not evaluated.
+ * Errors of provided scale factors are summed up in quadrature, producing a single systematic
+ * variation.
  */
 class LeptonSFWeight: public EventWeightPlugin
 {
diff --git a/modules/extensions/src/LeptonSFWeight.cpp b/modules/extensions/src/LeptonSFWeight.cpp
index 081629e..eeb4b9d 100644
--- a/modules/extensions/src/LeptonSFWeight.cpp
+++ b/modules/extensions/src/LeptonSFWeight.cpp
@@ -56,8 +56,8 @@ LeptonSFWeight::LeptonSFWeight(Lepton::Flavour targetFlavour_, std::string const
 {
     LoadScaleFactors(srcFileName, histogramNames);
     
-    // The plugin will calculate one weight per event
-    EventWeightPlugin::weights.push_back(0.);
+    // The plugin will calculate one weight per event with one systematic variation
+    EventWeightPlugin::weights = {0., 0., 0.};
 }
 
 
@@ -168,8 +168,11 @@ bool LeptonSFWeight::ProcessEvent()
     auto const &leptons = leptonPlugin->GetLeptons();
     
     
-    // Loop over the leptons and calculate the total scale factor
+    // Loop over the leptons and calculate the total scale factor and its relative uncertainty.
+    //Relative uncertainty of a product of multiple scale factors for, potentially, multiple
+    //leptons is simply a sum of all individual relative uncertainties in quadrature.
     double scaleFactor = 1.;
+    double relUnc2 = 0.;
     
     for (auto const &lepton: leptons)
     {
@@ -177,11 +180,14 @@ bool LeptonSFWeight::ProcessEvent()
         if (lepton.GetFlavour() != targetFlavour)
             continue;
         
+        
         // Loop over components of the scale factor
         for (auto const &sfObject: sfComponents)
         {
             int const bin = sfObject.hist->FindFixBin(sfObject.x(lepton), sfObject.y(lepton));
-            scaleFactor *= sfObject.hist->GetBinContent(bin);
+            double const sf = sfObject.hist->GetBinContent(bin);
+            scaleFactor *= sf;
+            relUnc2 += std::pow(sfObject.hist->GetBinError(bin) / sf, 2);
         }
     }
     
@@ -189,6 +195,10 @@ bool LeptonSFWeight::ProcessEvent()
     // Update the nominal weight
     weights.at(0) = scaleFactor;
     
+    double const unc = scaleFactor * std::sqrt(relUnc2);
+    weights.at(1) = scaleFactor + unc;
+    weights.at(2) = scaleFactor - unc;
+    
     
     return true;
 }