Implement matching to generator-level jets

include/PECFwk/PECReader/PECJetMETReader.hpp

@@ -3,6 +3,7 @@
 #include <PECFwk/core/JetMETReader.hpp>
 
 #include <PECFwk/core/LeptonReader.hpp>
+#include <PECFwk/core/GenJetMETReader.hpp>
 #include <PECFwk/PECReader/Candidate.hpp>
 #include <PECFwk/PECReader/Jet.hpp>
 
@@ -20,8 +21,11 @@ class PECInputData;
  * cleaned against tight leptons produced by a LeptonReader with name "Leptons". This behaviour
  * can be configured with method ConfigureLeptonCleaning.
  * 
- * Currently the plugin is not able to reapply JEC, does not perform matching to GEN-level jets,
- * and neglects any systematic variations.
+ * If the name of a plugin that reads generator-level jets is provided with the help of the method
+ * SetGenJetReader, angular matching to them is performed. The maximal allowed angular distance for
+ * matching is set to half of the radius parameter of reconstructed jets.
+ * 
+ * Currently the plugin is not able to reapply JEC and neglects any systematic variations.
  */
 class PECJetMETReader: public JetMETReader
 {
@@ -83,7 +87,10 @@ class PECJetMETReader: public JetMETReader
      */
     virtual double GetJetRadius() const override;
 
-    /// Specifiy desired selection on jets
+    /// Specifies name of the plugin that provides generator-level jets
+    void SetGenJetReader(std::string const name = "GenJetMET");
+
+    /// Specifies desired selection on jets
     void SetSelection(double minPt, double maxAbsEta);
 
 private:
@@ -146,4 +153,14 @@ class PECJetMETReader: public JetMETReader
      * Exploited in jet cleaning.
      */
     double leptonDR2;
+
+    /**
+     * \brief Name of the plugin that produces generator-level jets
+     * 
+     * The name can be empty. In this case no matching to generator-level jets is performed.
+     */
+    std::string genJetPluginName;
+
+    /// Non-owning pointer to a plugin that produces generator-level jets
+    GenJetMETReader const *genJetPlugin;
 };

modules/PECReader/src/PECJetMETReader.cpp

@@ -1,6 +1,7 @@
 #include <PECFwk/PECReader/PECJetMETReader.hpp>
 
 #include <PECFwk/core/Processor.hpp>
+#include <PECFwk/core/PhysicsObjects.hpp>
 #include <PECFwk/core/ROOTLock.hpp>
 #include <PECFwk/PECReader/PECInputData.hpp>
 
@@ -17,7 +18,8 @@ PECJetMETReader::PECJetMETReader(std::string name /*= "JetMET"*/):
     treeName("pecJetMET/JetMET"),
     bfJetPointer(&bfJets), bfMETPointer(&bfMETs),
     minPt(0.), maxAbsEta(std::numeric_limits<double>::infinity()),
-    leptonPluginName("Leptons"), leptonPlugin(nullptr), leptonDR2(0.3 * 0.3)
+    leptonPluginName("Leptons"), leptonPlugin(nullptr), leptonDR2(0.3 * 0.3),
+    genJetPluginName(""), genJetPlugin(nullptr)
 {}
 
 
@@ -29,11 +31,12 @@ PECJetMETReader::PECJetMETReader(PECJetMETReader const &src) noexcept:
     bfJetPointer(&bfJets), bfMETPointer(&bfMETs),
     minPt(src.minPt), maxAbsEta(src.maxAbsEta),
     leptonPluginName(src.leptonPluginName), leptonPlugin(src.leptonPlugin),
-    leptonDR2(src.leptonDR2)
+    leptonDR2(src.leptonDR2),
+    genJetPluginName(src.genJetPluginName), genJetPlugin(src.genJetPlugin)
 {}
 
 
-PECJetMETReader::~PECJetMETReader()
+PECJetMETReader::~PECJetMETReader() noexcept
 {}
 
 
@@ -51,10 +54,13 @@ void PECJetMETReader::BeginRun(Dataset const &)
     inputDataPlugin = dynamic_cast<PECInputData const *>(GetDependencyPlugin(inputDataPluginName));
 
 
-    // Save pointer to plugin that produces leptons
+    // Save pointers to plugins that produce leptons and generator-level jets
     if (leptonPluginName != "")
         leptonPlugin = dynamic_cast<LeptonReader const *>(GetDependencyPlugin(leptonPluginName));
 
+    if (genJetPluginName != "")
+        genJetPlugin = dynamic_cast<GenJetMETReader const *>(GetDependencyPlugin(genJetPluginName));
+
 
     // Set up the tree. Branches with properties that are not currently not used, are disabled
     inputDataPlugin->LoadTree(treeName);
@@ -80,6 +86,12 @@ double PECJetMETReader::GetJetRadius() const
 }
 
 
+void PECJetMETReader::SetGenJetReader(std::string const name /*= "GenJetMET"*/)
+{
+    genJetPluginName = name;
+}
+
+
 void PECJetMETReader::SetSelection(double minPt_, double maxAbsEta_)
 {
     minPt = minPt_;
@@ -152,6 +164,31 @@ bool PECJetMETReader::ProcessEvent()
         jet.SetParentID(j.Flavour());
 
 
+        // Perform matching to generator-level jets if the corresponding reader is available. Choose
+        //the closest jet but require that the angular separation is not larger than half of the
+        //radius parameter of the reconstructed jets
+        if (genJetPlugin)
+        {
+            double minDR2 = std::pow(GetJetRadius() / 2., 2);
+            GenJet const *matchedGenJet = nullptr;
+
+            for (auto const &genJet: genJetPlugin->GetJets())
+            {
+                double const dR2 = std::pow(p4.Eta() - genJet.Eta(), 2) +
+                 std::pow(TVector2::Phi_mpi_pi(p4.Phi() - genJet.Phi()), 2);
+                //^ Do not use TLorentzVector::DeltaR to avoid calculating sqrt
+
+                if (dR2 < minDR2)
+                {
+                    matchedGenJet = &genJet;
+                    minDR2 = dR2;
+                }
+            }
+
+            jet.SetMatchedGenJet(matchedGenJet);
+        }
+
+
         jets.push_back(jet);
     }
 

test/minimal.cpp

@@ -10,6 +10,7 @@
 #include <PECFwk/extensions/MetFilter.hpp>
 
 #include <PECFwk/PECReader/PECGeneratorReader.hpp>
+#include <PECFwk/PECReader/PECGenJetMETReader.hpp>
 #include <PECFwk/PECReader/PECGenParticleReader.hpp>
 #include <PECFwk/PECReader/PECInputData.hpp>
 #include <PECFwk/PECReader/PECJetMETReader.hpp>
@@ -69,9 +70,11 @@ int main()
     processor.RegisterPlugin(BuildPECTriggerFilter(false, triggerRanges));
     processor.RegisterPlugin(new PECLeptonReader);
     processor.RegisterPlugin(new LeptonFilter("LeptonFilter", Lepton::Flavour::Muon, 22., 2.1));
+    processor.RegisterPlugin(new PECGenJetMETReader);
 
     PECJetMETReader *jetReader = new PECJetMETReader;
     jetReader->SetSelection(30., 2.4);
+    jetReader->SetGenJetReader();
     processor.RegisterPlugin(jetReader);
 
     JetFilter *jetFilter = new JetFilter(0., bTagger);
@@ -139,8 +142,20 @@ int main()
         cout << "\nAnalysis jets:\n";
 
         for (auto const &j: jetReader->GetJets())
+        {
             cout << " pt: " << j.Pt() << ", eta: " << j.Eta() << ", b-tag: " <<
               j.BTag(BTagger::Algorithm::CSV) << ", flavour: " << j.GetParentID() << '\n';
+
+            cout << "  pt of matched GEN jet: ";
+            GenJet const *genJet = j.MatchedGenJet();
+
+            if (genJet)
+                cout << genJet->Pt();
+            else
+                cout << "(none)";
+
+            cout << '\n';
+        }
 
         cout << "\nMET: " << jetReader->GetMET().Pt() << '\n';