From 70d52090f8fdcc0fd9f55d356b163c136d8b2cbf Mon Sep 17 00:00:00 2001
From: Michael Holzbock <michael.holzbock@cern.ch>
Date: Wed, 15 May 2024 12:37:29 +0200
Subject: [PATCH] Jet pile-up tagging R22 "Phase-2" pre-recommendations

Jet pile-up tagging R22 "Phase-2" pre-recommendations
---
 .../python/JetAnalysisConfig.py               | 45 ++++++++++++++-----
 .../JetJvtEfficiency/FJvtSelectionTool.h      |  4 +-
 .../JetJvtEfficiency/JvtEfficiencyToolBase.h  |  2 +-
 .../JetJvtEfficiency/JvtSelectionToolBase.h   |  2 +-
 .../Root/FJvtSelectionTool.cxx                |  5 ++-
 .../Root/JvtEfficiencyToolBase.cxx            |  3 +-
 .../Root/NNJvtSelectionTool.cxx               |  2 +-
 7 files changed, 45 insertions(+), 18 deletions(-)

diff --git a/PhysicsAnalysis/Algorithms/JetAnalysisAlgorithms/python/JetAnalysisConfig.py b/PhysicsAnalysis/Algorithms/JetAnalysisAlgorithms/python/JetAnalysisConfig.py
index f9f906d8f10e..018e0d251f11 100644
--- a/PhysicsAnalysis/Algorithms/JetAnalysisAlgorithms/python/JetAnalysisConfig.py
+++ b/PhysicsAnalysis/Algorithms/JetAnalysisAlgorithms/python/JetAnalysisConfig.py
@@ -127,6 +127,10 @@ class SmallRJetAnalysisConfig (ConfigBlock) :
             info="whether to run JVT selection. The default is True.")
         self.addOption ('runFJvtSelection', False, type=bool,
             info="whether to run forward JVT selection. The default is False.")
+        self.addOption ('jvtWP', "FixedEffPt", type=str,
+            info="which Jvt WP to apply. The default is FixedEffPt.")
+        self.addOption ('fJvtWP', "Loose", type=str,
+            info="which fJvt WP to apply. The default is Loose.")
         self.addOption ('runJvtEfficiency', True, type=bool,
             info="whether to calculate the JVT efficiency. The default is True.")
         self.addOption ('runFJvtEfficiency', False, type=bool,
@@ -266,6 +270,14 @@ class SmallRJetAnalysisConfig (ConfigBlock) :
             raise ValueError(
                 "Unsupported input type '{0}' for R=0.4 jets!".format(self.jetInput) )
 
+        if self.jvtWP not in ["FixedEffPt"]:
+            raise ValueError(
+                "Unsupported NNJvt WP '{0}'".format(self.jvtWP) )
+
+        if self.fJvtWP not in ["Loose", "Tight", "Tighter"]:
+            raise ValueError(
+                "Unsupported fJvt WP '{0}'".format(self.fJvtWP) )
+
         if not config.isPhyslite() or self.recalibratePhyslite:
             # Prepare the jet calibration algorithm
             alg = config.createAlgorithm( 'CP::JetCalibrationAlg', 'JetCalibrationAlg'+postfix )
@@ -343,11 +355,12 @@ class SmallRJetAnalysisConfig (ConfigBlock) :
         # Set up the jet efficiency scale factor calculation algorithm
         # Change the truthJetCollection property to AntiKt4TruthWZJets if preferred
         if self.runJvtSelection :
+            assert self.jetInput=="EMPFlow", "NNJvt WPs and SFs only valid for PFlow jets"
             alg = config.createAlgorithm('CP::AsgSelectionAlg', f'JvtSelectionAlg{postfix}')
             config.addPrivateTool('selectionTool', 'CP::NNJvtSelectionTool')
             alg.selectionTool.JetContainer = config.readName(self.containerName)
-            alg.selectionTool.WorkingPoint = "FixedEffPt"
-            alg.selectionTool.MaxPtForJvt = 60e3 if self.jetInput == "EMPFlow" else 120e3
+            alg.selectionTool.WorkingPoint = self.jvtWP
+            alg.selectionTool.MaxPtForJvt = 60e3
             alg.selectionDecoration = "jvt_selection,as_char"
             alg.particles = config.readName(self.containerName)
 
@@ -355,8 +368,12 @@ class SmallRJetAnalysisConfig (ConfigBlock) :
                 alg = config.createAlgorithm( 'CP::JvtEfficiencyAlg', 'JvtEfficiencyAlg'+postfix )
                 config.addPrivateTool( 'efficiencyTool', 'CP::NNJvtEfficiencyTool' )
                 alg.efficiencyTool.JetContainer = config.readName(self.containerName)
-                alg.efficiencyTool.MaxPtForJvt = 60e3 if self.jetInput == "EMPFlow" else 120e3
-                alg.efficiencyTool.WorkingPoint = 'FixedEffPt'
+                alg.efficiencyTool.MaxPtForJvt = 60e3
+                alg.efficiencyTool.WorkingPoint = self.jvtWP
+                if config.geometry() is LHCPeriod.Run2:
+                    alg.efficiencyTool.SFFile = "JetJvtEfficiency/May2024/NNJvtSFFile_Run2_EMPFlow.root"
+                else:
+                    alg.efficiencyTool.SFFile = "JetJvtEfficiency/May2024/NNJvtSFFile_Run3_EMPFlow.root"
                 alg.selection = 'jvt_selection,as_char'
                 alg.scaleFactorDecoration = 'jvt_effSF_%SYS%'
                 alg.outOfValidity = 2
@@ -368,20 +385,23 @@ class SmallRJetAnalysisConfig (ConfigBlock) :
             config.addSelection (self.containerName, 'baselineJvt', 'jvt_selection,as_char', preselection=False)
 
         if self.runFJvtSelection :
+            assert self.jetInput=="EMPFlow", "fJvt WPs and SFs only valid for PFlow jets"
             alg = config.createAlgorithm('CP::AsgSelectionAlg', f'FJvtSelectionAlg{postfix}')
             config.addPrivateTool('selectionTool', 'CP::FJvtSelectionTool')
             alg.selectionTool.JetContainer = config.readName(self.containerName)
-            alg.selectionTool.WorkingPoint = "Loose"
+            alg.selectionTool.WorkingPoint = self.fJvtWP
             alg.selectionDecoration = "fjvt_selection,as_char"
             alg.particles = config.readName(self.containerName)
-            alg = config.createAlgorithm( 'CP::JvtEfficiencyAlg', 'ForwardJvtEfficiencyAlg' )
-            config.addSelection (self.containerName, 'baselineFJvt', 'fjvt_selection,as_char', preselection=False)
 
             if self.runFJvtEfficiency and config.dataType() is not DataType.Data:
                 alg = config.createAlgorithm( 'CP::JvtEfficiencyAlg', 'FJvtEfficiencyAlg'+postfix )
                 config.addPrivateTool( 'efficiencyTool', 'CP::FJvtEfficiencyTool' )
                 alg.efficiencyTool.JetContainer = config.readName(self.containerName)
-                alg.efficiencyTool.WorkingPoint = 'Loose'
+                alg.efficiencyTool.WorkingPoint = self.fJvtWP
+                if config.geometry() is LHCPeriod.Run2:
+                    alg.efficiencyTool.SFFile = "JetJvtEfficiency/May2024/fJvtSFFile_Run2_EMPFlow.root"
+                else:
+                    alg.efficiencyTool.SFFile = "JetJvtEfficiency/May2024/fvtSFFile_Run3_EMPFlow.root"
                 alg.selection = 'fjvt_selection,as_char'
                 alg.scaleFactorDecoration = 'fjvt_effSF_%SYS%'
                 alg.outOfValidity = 2
@@ -389,7 +409,8 @@ class SmallRJetAnalysisConfig (ConfigBlock) :
                 alg.skipBadEfficiency = False
                 alg.jets = config.readName (self.containerName)
                 alg.preselection = config.getPreselection (self.containerName, '')
-                config.addOutputVar (self.containerName, alg.scaleFactorDecoration, 'jvtEfficiency')
+                config.addOutputVar (self.containerName, alg.scaleFactorDecoration, 'fjvtEfficiency')
+            config.addSelection (self.containerName, 'baselineFJvt', 'fjvt_selection,as_char', preselection=False)
 
 
 
@@ -669,6 +690,7 @@ def makeSmallRJetAnalysisConfig( seq, containerName, jetCollection,
                                  jetInput, postfix = None,
                                  runJvtUpdate = None, runNNJvtUpdate = None, runFJvtUpdate = None,
                                  runJvtSelection = None, runFJvtSelection = None,
+                                 jvtWP = None, fJvtWP = None,
                                  runJvtEfficiency = None, runFJvtEfficiency = None,
                                  systematicsModelJES = None, systematicsModelJER = None):
     """Add algorithms for the R=0.4 jets.
@@ -683,6 +705,8 @@ def makeSmallRJetAnalysisConfig( seq, containerName, jetCollection,
         runFJvtUpdate -- Determines whether or not to update forward JVT on the jets
         runJvtSelection -- Determines whether or not to run JVT selection on the jets
         runFJvtSelection -- Determines whether or not to run forward JVT selection on the jets
+        jvtWP -- Defines the NNJvt WP to apply on the jets
+        fJvtWP -- Defines the fJvt WP to apply on the jets
         runJvtEfficiency -- Determines whether or not to calculate the JVT efficiency
         runFJvtEfficiency -- Determines whether or not to calculate the forward JVT efficiency
         systematicsModelJES -- Which NP systematicsModelJES scheme should be used (All, Global, Category, Scenario)
@@ -693,7 +717,6 @@ def makeSmallRJetAnalysisConfig( seq, containerName, jetCollection,
         raise ValueError(
             "Unsupported input type '{0}' for R=0.4 jets!".format(jetInput) )
 
-
     config = SmallRJetAnalysisConfig (containerName, jetCollection, jetInput)
     config.setOptionValue ('postfix', postfix)
     config.setOptionValue ('runJvtUpdate', runJvtUpdate)
@@ -701,6 +724,8 @@ def makeSmallRJetAnalysisConfig( seq, containerName, jetCollection,
     config.setOptionValue ('runFJvtUpdate', runFJvtUpdate)
     config.setOptionValue ('runJvtSelection', runJvtSelection)
     config.setOptionValue ('runFJvtSelection', runFJvtSelection)
+    config.setOptionValue ('jvtWP', jvtWP)
+    config.setOptionValue ('fJvtWP', fJvtWP)
     config.setOptionValue ('runJvtEfficiency', runJvtEfficiency)
     config.setOptionValue ('runFJvtEfficiency', runFJvtEfficiency)
     config.setOptionValue ('systematicsModelJES', systematicsModelJES)
diff --git a/Reconstruction/Jet/JetJvtEfficiency/JetJvtEfficiency/FJvtSelectionTool.h b/Reconstruction/Jet/JetJvtEfficiency/JetJvtEfficiency/FJvtSelectionTool.h
index d780eec7463e..d20037e21f1f 100644
--- a/Reconstruction/Jet/JetJvtEfficiency/JetJvtEfficiency/FJvtSelectionTool.h
+++ b/Reconstruction/Jet/JetJvtEfficiency/JetJvtEfficiency/FJvtSelectionTool.h
@@ -27,11 +27,11 @@ namespace CP {
                 "The working point to use. Set to 'Custom' to manually set the values"};
         Gaudi::Property<float> m_jvtCut{this, "JvtCut", 999, "The JVT selection to make"};
         SG::ReadDecorHandleKey<xAOD::JetContainer> m_jvtMoment{
-                this, "JvtMomentName", "FJvt", "The name of the Jvt moment to use"};
+                this, "JvtMomentName", "DFCommonJets_fJvt", "The name of the Jvt moment to use"};
         SG::ReadDecorHandleKey<xAOD::JetContainer> m_timingMoment{
                 this, "TimingMomentName", "Timing", "The name of the timing moment to use"};
         Gaudi::Property<float> m_timingCut{
-                this, "TimingCut", 10, "Only accept jets with time less than this"};
+                this, "TimingCut", -1, "Only accept jets with time less than this; negative values deactivate timing requirement"};
 
         virtual bool select(const xAOD::IParticle *jet) const override;
 
diff --git a/Reconstruction/Jet/JetJvtEfficiency/JetJvtEfficiency/JvtEfficiencyToolBase.h b/Reconstruction/Jet/JetJvtEfficiency/JetJvtEfficiency/JvtEfficiencyToolBase.h
index 390f89edce6f..96c3583f6daa 100644
--- a/Reconstruction/Jet/JetJvtEfficiency/JetJvtEfficiency/JvtEfficiencyToolBase.h
+++ b/Reconstruction/Jet/JetJvtEfficiency/JetJvtEfficiency/JvtEfficiencyToolBase.h
@@ -51,7 +51,7 @@ namespace CP {
                 this, "DummySFError", 0.1, "The amount by which to vary the dummy SF"};
         // NB: Use a string not a read handle key as this is not written with a write handle key
         Gaudi::Property<std::string> m_jetEtaName{
-                this, "JetEtaName", "DetectorEta", "The name of the jet eta to use."};
+                this, "JetEtaName", "eta", "The name of the jet eta to use."};
         std::unique_ptr<TH2> m_jvtHist;
         std::unique_ptr<TH2> m_effHist;
         bool m_useDummySFs{false};
diff --git a/Reconstruction/Jet/JetJvtEfficiency/JetJvtEfficiency/JvtSelectionToolBase.h b/Reconstruction/Jet/JetJvtEfficiency/JetJvtEfficiency/JvtSelectionToolBase.h
index 9807a7c8c02a..d21b89490f5b 100644
--- a/Reconstruction/Jet/JetJvtEfficiency/JetJvtEfficiency/JvtSelectionToolBase.h
+++ b/Reconstruction/Jet/JetJvtEfficiency/JetJvtEfficiency/JvtSelectionToolBase.h
@@ -34,7 +34,7 @@ namespace CP {
                 this, "MaxEtaForJvt", 2.5, "Accept all jets with |eta| above this"};
         // NB: Use a string not a read handle key as this is not written with a write handle key
         Gaudi::Property<std::string> m_jetEtaName{
-                this, "JetEtaName", "DetectorEta", "The name of the jet eta to use."};
+                this, "JetEtaName", "eta", "The name of the jet eta to use."};
 
         // The template AcceptInfo object
         asg::AcceptInfo m_info;
diff --git a/Reconstruction/Jet/JetJvtEfficiency/Root/FJvtSelectionTool.cxx b/Reconstruction/Jet/JetJvtEfficiency/Root/FJvtSelectionTool.cxx
index 7596d08c6540..c1f8738edb87 100644
--- a/Reconstruction/Jet/JetJvtEfficiency/Root/FJvtSelectionTool.cxx
+++ b/Reconstruction/Jet/JetJvtEfficiency/Root/FJvtSelectionTool.cxx
@@ -2,7 +2,7 @@
 #include "AsgDataHandles/ReadDecorHandle.h"
 
 namespace {
-    const static std::map<std::string, float> workingPoints{{"Loose", 0.5}, {"Tight", 0.4}};
+    const static std::map<std::string, float> workingPoints{{"Loose", 0.5}, {"Tight", 0.4}, {"Tighter", 0.2}};
 }
 
 namespace CP {
@@ -43,7 +43,8 @@ namespace CP {
     }
 
     bool FJvtSelectionTool::select(const xAOD::IParticle *jet) const {
-        return m_jvtAcc(*jet) <= m_jvtCut && m_timingAcc(*jet) <= m_timingCut;
+        // select jet if it passes fJvt requirement and timing cut (if configured)
+        return m_jvtAcc(*jet) <= m_jvtCut && ( m_timingCut > 0 ? std::abs( m_timingAcc(*jet) ) <= m_timingCut : true );
     }
 
 } // namespace CP
diff --git a/Reconstruction/Jet/JetJvtEfficiency/Root/JvtEfficiencyToolBase.cxx b/Reconstruction/Jet/JetJvtEfficiency/Root/JvtEfficiencyToolBase.cxx
index 375138699958..f729359399cd 100644
--- a/Reconstruction/Jet/JetJvtEfficiency/Root/JvtEfficiencyToolBase.cxx
+++ b/Reconstruction/Jet/JetJvtEfficiency/Root/JvtEfficiencyToolBase.cxx
@@ -11,7 +11,7 @@
 namespace {
     bool getBin(const TAxis &axis, float value, int &bin) {
         bin = axis.FindBin(value);
-        return (bin != 0 && bin != axis.GetNbins());
+        return (bin != 0 && bin != axis.GetNbins()+1);
     }
     bool getBinContentAndError(const TH2 &h, float x, float y, float &content, float &error) {
         int xBin{}, yBin{};
@@ -85,6 +85,7 @@ namespace CP {
     StatusCode JvtEfficiencyToolBase::initHists(const std::string &file, const std::string &wp) {
         if (file.empty()) {
             m_useDummySFs = true;
+            ATH_MSG_INFO("No SF file provided, running with dummy SFs of 1 +/- " << m_dummySFError);
             return StatusCode::SUCCESS;
         }
         std::string resolved = PathResolverFindCalibFile(file);
diff --git a/Reconstruction/Jet/JetJvtEfficiency/Root/NNJvtSelectionTool.cxx b/Reconstruction/Jet/JetJvtEfficiency/Root/NNJvtSelectionTool.cxx
index 067b2e9883de..67afe1cde1b7 100644
--- a/Reconstruction/Jet/JetJvtEfficiency/Root/NNJvtSelectionTool.cxx
+++ b/Reconstruction/Jet/JetJvtEfficiency/Root/NNJvtSelectionTool.cxx
@@ -55,6 +55,6 @@ namespace CP {
     }
 
     bool NNJvtSelectionTool::select(const xAOD::IParticle *jet) const {
-        return m_jvtAcc(*jet) > m_cutMap(*jet);
+        return m_jvtAcc(*jet) > m_cutMap(jet->pt(), m_etaAcc(*jet));
     }
 } // namespace CP
-- 
GitLab