Allow application of T1 MET corrections directly to raw MET

It has been validated with data files that this commit results in the
same values for MET as the subtraction/correction scheme for majority of
events. For few events where the results are different, the difference
is assumed to be caused by the non-trivial jet selection applied by the
PAT MET tool.

include/mensura/extensions/JetMETUpdate.hpp

@@ -26,6 +26,9 @@ class PileUpReader;
  * i.e. with outdated corrections. If stochastic JER corrections have been used in computation of
  * MET, they can only be undone on the average due to their random nature.
  * 
+ * It is also possible to start computation from the raw MET rather than default one. In this case
+ * T1 corrections can be applied directly.
+ * 
  * The source MET can have any additive corrections applied. They will be preserved by this plugin.
  */
 class JetMETUpdate: public JetMETReader
@@ -73,20 +76,27 @@ class JetMETUpdate: public JetMETReader
      * 
      * Arguments are names of instances of JetCorrectorService that evaluate full and L1-only
      * corrections to be applied to compute T1 MET corrections. In addition, corrections used in
-     * the original MET need to be provided to undo outdated T1 corrections. The "full" corrections
+     * the original MET can be provided to undo outdated T1 corrections. The "full" corrections
      * to be applied are not necessarily the same as given to method SetJetCorrection.
      * 
      * For any correction level, empty strings can be given as names of the two corresponding
      * correctors. In this case the corresponding corrections are not applied. This is useful, for
      * instance, when L1-only corrections have not changed, and thus there is no need to evaluate
      * them explicitly since their contributions to MET would cancel out.
+     * 
+     * When computation starts from raw MET (as requested with UseRawMET), it is technically
+     * possible to "undo" MET corrections corresponding to provided original corrections, but this
+     * does not make sense.
      */
     void SetJetCorrectionForMET(std::string const &fullNew, std::string const &l1New,
       std::string const &fullOrig, std::string const &l1Orig);
 
     /// Specifies desired selection on jets
     void SetSelection(double minPt, double maxAbsEta);
 
+    /// Instructs to use raw MET for computation instead of the default one
+    void UseRawMET(bool set = true);
+
 private:
     /**
      * \brief Reads jets and MET from the source JetMETReader and recorrects them
@@ -141,6 +151,9 @@ class JetMETUpdate: public JetMETReader
     /// Minimal transverse momentum for jets to be considered in T1 MET corrections
     double minPtForT1;
 
+    /// Specifies whether computation should start from raw or default MET
+    bool useRawMET;
+
     /// Type of requested systematical variation
     JetCorrectorService::SystType systType;
 

modules/extensions/src/JetMETUpdate.cpp

@@ -17,7 +17,8 @@ JetMETUpdate::JetMETUpdate(std::string const name /*= "JetMET"*/):
     jetCorrForJets(nullptr),
     jetCorrForMETFull(nullptr), jetCorrForMETL1(nullptr),
     jetCorrForMETOrigFull(nullptr), jetCorrForMETOrigL1(nullptr),
-    minPt(0.), maxAbsEta(std::numeric_limits<double>::infinity()), minPtForT1(15.)
+    minPt(0.), maxAbsEta(std::numeric_limits<double>::infinity()), minPtForT1(15.),
+    useRawMET(false)
 {}
 
 
@@ -98,17 +99,6 @@ void JetMETUpdate::SetJetCorrection(std::string const &jetCorrServiceName)
 void JetMETUpdate::SetJetCorrectionForMET(std::string const &fullNew, std::string const &l1New,
   std::string const &fullOrig, std::string const &l1Orig)
 {
-    // Make sure that for each level the two correctors are either both given or both missing
-    if ((fullNew == "") != (fullOrig == "") or (l1New == "") != (l1Orig == ""))
-    {
-        std::ostringstream message;
-        message << "JetMETUpdate[\"" << GetName() << "\"]::SetJetCorrectionForMET: " <<
-          "For each correction level the two correctors must be either provided both or " <<
-          "missing both.";
-        throw std::logic_error(message.str());
-    }
-
-
     // If new and original correctors are the same, drop them since their effect would cancel out
     if (fullNew != fullOrig)
     {
@@ -135,6 +125,12 @@ void JetMETUpdate::SetSelection(double minPt_, double maxAbsEta_)
 }
 
 
+void JetMETUpdate::UseRawMET(bool set /*= true*/)
+{
+    useRawMET = set;
+}
+
+
 bool JetMETUpdate::ProcessEvent()
 {
     jets.clear();
@@ -164,17 +160,24 @@ bool JetMETUpdate::ProcessEvent()
         // Compute the shift in MET due to T1 corrections
         if (jet.Pt() > minPtForT1)
         {
-            // The shift due to different L1 corrections
-            if (jetCorrForMETL1 and jetCorrForMETOrigL1)
-                metShift += srcJet.RawP4() *
-                  (jetCorrForMETL1->Eval(srcJet, rho, systType, systDirection) -
-                  jetCorrForMETOrigL1->Eval(srcJet, rho, systType, systDirection));
+            // Undo applied T1 corrections
+            if (jetCorrForMETOrigL1)
+                metShift -=
+                  srcJet.RawP4() * jetCorrForMETOrigL1->Eval(srcJet, rho, systType, systDirection);
+
+            if (jetCorrForMETOrigFull)
+                metShift +=
+                  srcJet.RawP4() * jetCorrForMETOrigFull->Eval(srcJet, rho, systType, systDirection);
+
+
+            // Apply new T1 corrections
+            if (jetCorrForMETL1)
+                metShift +=
+                  srcJet.RawP4() * jetCorrForMETL1->Eval(srcJet, rho, systType, systDirection);
 
-            // The shift due to different full corrections
-            if (jetCorrForMETFull and jetCorrForMETOrigFull)
-                metShift -= srcJet.RawP4() *
-                  (jetCorrForMETFull->Eval(srcJet, rho, systType, systDirection) -
-                  jetCorrForMETOrigFull->Eval(srcJet, rho, systType, systDirection));
+            if (jetCorrForMETFull)
+                metShift -=
+                  srcJet.RawP4() * jetCorrForMETFull->Eval(srcJet, rho, systType, systDirection);
         }
 
 
@@ -189,7 +192,9 @@ bool JetMETUpdate::ProcessEvent()
 
 
     // Update MET
-    TLorentzVector updatedMET(jetmetPlugin->GetMET().P4() + metShift);
+    TLorentzVector const &startingMET = (useRawMET) ?
+      jetmetPlugin->GetRawMET().P4() : jetmetPlugin->GetMET().P4();
+    TLorentzVector updatedMET(startingMET + metShift);
     met.SetPtEtaPhiM(updatedMET.Pt(), 0., updatedMET.Phi(), 0.);