From 96d5983081af31aa2eff85629afa32e4e6ead3f3 Mon Sep 17 00:00:00 2001
From: Jonathan Burr <jon.burr@cern.ch>
Date: Fri, 28 Aug 2020 21:34:20 +0000
Subject: [PATCH] Apply urgent and simple fixes from MR review

Downgrade INFO message in initialization, fix incorrect operator in operator-= definition, clarify operator[] comments, replace std::pow calls
---
 .../TrigEFMissingET/CMakeLists.txt            |   7 +-
 .../ApproximateTrackToLayerTool.h             |  35 ++
 .../TrigEFMissingET/ExtendTrackToLayerTool.h  |  43 ++
 .../TrigEFMissingET/IExtendTrackToLayerTool.h |  50 ++
 .../TrigEFMissingET/PUClassification.h        |  30 +
 .../TrigEFMissingET/PUSplitGrid.h             |  43 ++
 .../TrigEFMissingET/PeriodicGridBase.h        | 194 +++++++
 .../TrigEFMissingET/PufitGrid.h               | 364 +++++-------
 .../TrigEFMissingET/PufitMultiGrid.h          | 351 ++++++++++++
 .../TrigEFMissingET/PufitMultiGrid.icc        | 313 ++++++++++
 .../TrigEFMissingET/PufitUtils.h              | 436 ++++++++------
 .../TrigEFMissingET/PufitUtils.icc            |  74 +++
 .../python/PUClassification.py                |  17 +
 .../src/ApproximateTrackToLayerTool.cxx       |  27 +
 .../TrigEFMissingET/src/CVFAlg.cxx            | 173 ++++++
 .../TrigEFMissingET/src/CVFAlg.h              |  71 +++
 .../TrigEFMissingET/src/CVFPrepAlg.cxx        |  71 +++
 .../TrigEFMissingET/src/CVFPrepAlg.h          |  62 ++
 .../TrigEFMissingET/src/CellFex.cxx           |  10 +
 .../src/ExtendTrackToLayerTool.cxx            |  66 +++
 .../TrigEFMissingET/src/MHTFex.cxx            |   5 +
 .../TrigEFMissingET/src/MHTPufitFex.cxx       | 204 +++++++
 .../TrigEFMissingET/src/MHTPufitFex.h         | 128 +++++
 .../TrigEFMissingET/src/PFOPrepAlg.cxx        | 137 +++++
 .../TrigEFMissingET/src/PFOPrepAlg.h          | 103 ++++
 .../TrigEFMissingET/src/PFSumFex.cxx          |  10 +
 .../TrigEFMissingET/src/PUSplitGrid.cxx       |  28 +
 .../TrigEFMissingET/src/PUSplitPufitFex.cxx   | 192 +++++++
 .../TrigEFMissingET/src/PUSplitPufitFex.h     | 114 ++++
 .../TrigEFMissingET/src/PeriodicGridBase.cxx  | 157 +++++
 .../TrigEFMissingET/src/PufitGrid.cxx         | 534 +++++++-----------
 .../TrigEFMissingET/src/PufitUtils.cxx        | 483 ++++++++--------
 .../TrigEFMissingET/src/TCFex.cxx             |   5 +
 .../TrigEFMissingET/src/TCPufitFex.cxx        | 147 ++---
 .../TrigEFMissingET/src/TrkMHTFex.cxx         |  23 +-
 .../components/TrigEFMissingET_entries.cxx    |  52 +-
 .../share/ref_RDOtoRDOTrig_mt1_build.ref      |  12 +
 .../share/ref_data_v1Dev_build.ref            |  12 +
 .../TrigEDMConfig/python/TriggerEDMRun3.py    |  12 +
 .../python/HLTMenuConfig/MET/AlgConfigs.py    | 263 +++++++--
 .../python/HLTMenuConfig/MET/ConfigHelpers.py |   8 +-
 .../HLTMenuConfig/MET/METRecoSequences.py     |  96 +++-
 .../python/HLTMenuConfig/Menu/LS2_v1.py       |  39 +-
 .../HLTMenuConfig/Menu/SignatureDicts.py      |  84 +--
 44 files changed, 4109 insertions(+), 1176 deletions(-)
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/ApproximateTrackToLayerTool.h
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/ExtendTrackToLayerTool.h
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/IExtendTrackToLayerTool.h
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PUClassification.h
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PUSplitGrid.h
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PeriodicGridBase.h
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitMultiGrid.h
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitMultiGrid.icc
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitUtils.icc
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/python/PUClassification.py
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/src/ApproximateTrackToLayerTool.cxx
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/src/CVFAlg.cxx
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/src/CVFAlg.h
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/src/CVFPrepAlg.cxx
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/src/CVFPrepAlg.h
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/src/ExtendTrackToLayerTool.cxx
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/src/MHTPufitFex.cxx
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/src/MHTPufitFex.h
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/src/PFOPrepAlg.cxx
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/src/PFOPrepAlg.h
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/src/PUSplitGrid.cxx
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/src/PUSplitPufitFex.cxx
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/src/PUSplitPufitFex.h
 create mode 100644 Trigger/TrigAlgorithms/TrigEFMissingET/src/PeriodicGridBase.cxx

diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/CMakeLists.txt b/Trigger/TrigAlgorithms/TrigEFMissingET/CMakeLists.txt
index 58d595186b02..cbd06b4aad0e 100644
--- a/Trigger/TrigAlgorithms/TrigEFMissingET/CMakeLists.txt
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/CMakeLists.txt
@@ -14,7 +14,12 @@ find_package( Eigen )
 atlas_add_component( TrigEFMissingET
    src/*.cxx src/components/*.cxx
    INCLUDE_DIRS ${EIGEN_INCLUDE_DIRS} ${FASTJET_INCLUDE_DIRS} ${FASTJETCONTRIB_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS} ${TDAQ-COMMON_INCLUDE_DIRS}
-   LINK_LIBRARIES ${EIGEN_LIBRARIES} ${FASTJET_LIBRARIES} ${FASTJETCONTRIB_LIBRARIES} ${ROOT_LIBRARIES} ${TDAQ-COMMON_LIBRARIES} AthContainers AthLinks AthenaBaseComps AthenaMonitoringKernelLib CaloConditions CaloDetDescrLib CaloEvent CaloGeoHelpers CaloIdentifier CaloInterfaceLib CxxUtils FourMomUtils GaudiKernel IRegionSelector Identifier InDetTrackSelectionToolLib JetEDM JetEvent LArCablingLib LArIdentifier LArRecConditions LArRecEvent PathResolver StoreGateLib TrigInterfacesLib TrigMissingEtEvent TrigParticle TrigSteeringEvent TrigT1Interfaces TrigT2CaloCommonLib TrigTimeAlgsLib xAODCaloEvent xAODCore xAODEventInfo xAODJet xAODMuon xAODPFlow xAODTracking )
+   LINK_LIBRARIES ${EIGEN_LIBRARIES} ${FASTJET_LIBRARIES} ${FASTJETCONTRIB_LIBRARIES} ${ROOT_LIBRARIES} ${TDAQ-COMMON_LIBRARIES}
+   AsgTools AthContainers AthLinks AthenaBaseComps AthenaMonitoringKernelLib CaloConditions CaloDetDescrLib CaloEvent
+   CaloGeoHelpers CaloIdentifier CaloInterfaceLib CxxUtils FourMomUtils GaudiKernel IRegionSelector Identifier InDetTrackSelectionToolLib
+   JetEDM JetEvent LArCablingLib LArIdentifier LArRecConditions LArRecEvent PathResolver RecoToolInterfaces StoreGateLib TrigInterfacesLib TrigMissingEtEvent
+   TrigParticle TrigSteeringEvent TrigT1Interfaces TrigT2CaloCommonLib TrigTimeAlgsLib TrkCaloExtension xAODCaloEvent xAODCore xAODEventInfo xAODJet xAODMuon
+   xAODPFlow xAODTracking xAODEventShape TrackVertexAssociationToolLib InDetTrackSelectionToolLib)
 
 # Install files from the package:
 atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} --extend-extensions=ATL900,ATL901 )
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/ApproximateTrackToLayerTool.h b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/ApproximateTrackToLayerTool.h
new file mode 100644
index 000000000000..69e00c336fb1
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/ApproximateTrackToLayerTool.h
@@ -0,0 +1,35 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef TRIGEFMISSINGET_APPROXIMATETRACKTOLAYERTOOL_H
+#define TRIGEFMISSINGET_APPROXIMATETRACKTOLAYERTOOL_H
+
+#include "TrigEFMissingET/IExtendTrackToLayerTool.h"
+#include "AsgTools/AsgTool.h"
+
+class ApproximateTrackToLayerTool : public virtual IExtendTrackToLayerTool, public asg::AsgTool
+{
+public:
+  ASG_TOOL_CLASS(ApproximateTrackToLayerTool, IExtendTrackToLayerTool)
+
+  ApproximateTrackToLayerTool(const std::string &name);
+
+  virtual TrackExtension extendTrack(
+      const EventContext &,
+      const xAOD::TrackParticle &track) const override;
+
+private:
+  // Extrapolation approximation taken from ATL-COM-PHYS-2016-430
+  Gaudi::Property<double> m_trackExtrapolationQuartic{
+      this, "TrackExtrapolationQuarticTerm", 14.6027571,
+      "The quartic term in the track extrapolation"};
+  Gaudi::Property<double> m_trackExtrapolationQuadratic{
+      this, "TrackExtrapolationQuadraticTerm", -44.7818374,
+      "The quadratic term in the track extrapolation"};
+  Gaudi::Property<double> m_trackExtrapolationLinear{
+      this, "TrackExtrapolationQuarticTerm", 540.656643,
+      "The linear term in the track extrapolation"};
+}; //> end class ApproximateTrackToLayerTool
+
+#endif //> !TRIGEFMISSINGET_APPROXIMATETRACKTOLAYERTOOL_H
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/ExtendTrackToLayerTool.h b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/ExtendTrackToLayerTool.h
new file mode 100644
index 000000000000..263f21c2ac96
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/ExtendTrackToLayerTool.h
@@ -0,0 +1,43 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef TRIGEFMISSINGET_EXTENDTRACKTOLAYERTOOL_H
+#define TRIGEFMISSINGET_EXTENDTRACKTOLAYERTOOL_H
+
+#include "TrigEFMissingET/IExtendTrackToLayerTool.h"
+#include "AsgTools/AsgTool.h"
+#include "AsgTools/ToolHandle.h"
+#include "RecoToolInterfaces/IParticleCaloExtensionTool.h"
+#include "CaloGeoHelpers/CaloSampling.h"
+
+/**
+ * @brief Extrapolate tracks using the calo extension tool
+ * @author Jon Burr
+ * 
+ * This is the slowest approach, but probably the most correct...
+ */
+class ExtendTrackToLayerTool : public virtual IExtendTrackToLayerTool, public asg::AsgTool
+{
+public:
+  ASG_TOOL_CLASS(ExtendTrackToLayerTool, IExtendTrackToLayerTool);
+
+  ExtendTrackToLayerTool(const std::string &name);
+
+  virtual StatusCode initialize() override;
+
+  virtual TrackExtension extendTrack(
+    const EventContext& ctx, 
+    const xAOD::TrackParticle &track) const override;
+private:
+  ToolHandle<Trk::IParticleCaloExtensionTool> m_extensionTool{
+      this, "ExtensionTool", "", "The track extension tool"};
+  Gaudi::Property<std::vector<std::string>> m_caloLayerNames{
+      this, "CaloLayerNames", {"EMB2", "EME2"}, "The (ordered) names of layers to which to extend."};
+
+  // Internals
+  /// The calo layer enums
+  std::vector<CaloSampling::CaloSample> m_caloLayers;
+}; //> end class ExtendTrackToLayerTool
+
+#endif //> !TRIGEFMISSINGET_EXTENDTRACKTOLAYERTOOL_H
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/IExtendTrackToLayerTool.h b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/IExtendTrackToLayerTool.h
new file mode 100644
index 000000000000..34fe0e982ee7
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/IExtendTrackToLayerTool.h
@@ -0,0 +1,50 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef TRIGEFMISSINGET_IEXTENDTRACKTOLAYERTOOL_H
+#define TRIGEFMISSINGET_IEXTENDTRACKTOLAYERTOOL_H
+
+#include "AsgTools/IAsgTool.h"
+#include "xAODTracking/TrackParticleContainer.h"
+#include <limits>
+#include <cmath>
+
+class IExtendTrackToLayerTool : public virtual asg::IAsgTool
+{
+  ASG_TOOL_INTERFACE(IExtendTrackToLayerTool)
+public:
+  /// Virtual destructor
+  virtual ~IExtendTrackToLayerTool() {}
+
+  /**
+   * @brief Helper struct to hold track extrapolation information
+   *
+   * Note that the extrapolation could decide that the given track never
+   * reaches the requested layer, in which case phi and eta will be -inf
+   */
+  struct TrackExtension
+  {
+    /// The extrapolated eta
+    double eta{-std::numeric_limits<double>::infinity()};
+    /// The extrapolated phi
+    double phi{-std::numeric_limits<double>::infinity()};
+    /// Whether or not the track crossed the layer
+    inline bool isValid() const
+    {
+      return std::isfinite(phi) && std::isfinite(eta);
+    }
+  };
+
+  /**
+   * @brief Extend the track to a given layer
+   * @param ctx The event context to use
+   * @param track The track to extend
+   * @return A struct containing the new eta/phi coordinates for the track
+   */
+  virtual TrackExtension extendTrack(
+    const EventContext& ctx,
+     const xAOD::TrackParticle &track) const = 0;
+}; //> end class IExtendTrackToLayerTool
+
+#endif //> !TRIGEFMISSINGET_IEXTENDTRACKTOLAYERTOOL_H
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PUClassification.h b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PUClassification.h
new file mode 100644
index 000000000000..460e8f6140dd
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PUClassification.h
@@ -0,0 +1,30 @@
+/*
+ * Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+#ifndef TRIGEFMISSINGET_PUCLASSIFICATION_H
+#define TRIGEFMISSINGET_PUCLASSIFICATION_H
+
+#include <type_traits>
+
+namespace HLT
+{
+  namespace MET
+  {
+    namespace PUClassification
+    {
+      enum PUClassification
+      {
+        NeutralForward = 1 << 0,
+        ChargedPU = 1 << 1,
+        ChargedHS = 1 << 2,
+      };
+      constexpr PUClassification NFcPU = static_cast<PUClassification>(NeutralForward | ChargedPU);
+      constexpr PUClassification NFcHS = static_cast<PUClassification>(NeutralForward | ChargedHS);
+      constexpr PUClassification Charged = static_cast<PUClassification>(ChargedPU | ChargedHS);
+      constexpr PUClassification All = static_cast<PUClassification>(NeutralForward | ChargedPU | ChargedHS);
+    } // namespace PUClassification
+  }   // namespace MET
+} // namespace HLT
+
+#endif //> !TRIGEFMISSINGET_PUCLASSIFICATION_H
\ No newline at end of file
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PUSplitGrid.h b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PUSplitGrid.h
new file mode 100644
index 000000000000..6462fd1f29ee
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PUSplitGrid.h
@@ -0,0 +1,43 @@
+/*
+ * Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+#ifndef TRIGEFMISSINGET_PUSPLITGRID_H
+#define TRIGEFMISSINGET_PUSPLITGRID_H
+
+#include "TrigEFMissingET/PufitMultiGrid.h"
+
+namespace HLT
+{
+  namespace MET
+  {
+    class PUSplitGrid : public PufitMultiGrid<3>
+    {
+    public:
+      /**
+       * @brief Create a new tower grid
+       * @param maxEta The maximum eta range for the grid
+       * @param nEtaTowers The number of eta towers
+       * @param nPhiTowers The number of phi towers
+       * @param displaceEta Whether to displace eta
+       * @param displacePhi Whether to displace phi
+       */
+      PUSplitGrid(
+          double maxEta,
+          std::size_t nEtaTowers,
+          std::size_t nPhiTowers,
+          bool displaceEta = false,
+          bool displacePhi = false);
+
+      /// Construct a grid from the provided parameters
+      PUSplitGrid(const GridParameters &parameters);
+
+      /// Copy constructor
+      PUSplitGrid(const PUSplitGrid &other);
+    }; //> end class PUSplitGrid
+
+    using PUSplitGridSet = PufitMultiGridSet<PUSplitGrid>;
+  } // namespace MET
+} // namespace HLT
+
+#endif //> !TRIGEFMISSINGET_PUSPLITGRID_H
\ No newline at end of file
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PeriodicGridBase.h b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PeriodicGridBase.h
new file mode 100644
index 000000000000..eb207439cbf5
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PeriodicGridBase.h
@@ -0,0 +1,194 @@
+/*
+ * Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+#ifndef TRIGEFMISSINGET_PERIODICGRIDBASE_H
+#define TRIGEFMISSINGET_PERIODICGRIDBASE_H
+
+#include <cstddef>
+#include <utility>
+
+/**
+ * @file TrigEFMissingET/PeriodicGridBase.h
+ * @brief Provide a base class for the grids used in some pufit algorithms
+ * @author Jon bur (jon.burr@cern.ch)
+ */
+
+namespace HLT
+{
+  namespace MET
+  {
+    /// Enum to describe the positioning of the grid
+    enum GridDisplacement
+    {
+      /// The grid is not shifted
+      NoDisplacement = 0,
+      /// The grid is shifted by half a tower width in eta
+      EtaDisplaced = 1,
+      /// The grid is shifted by half a tower width in phi
+      PhiDisplaced = 2,
+      /// The grid is shifted by half a tower width in both eta and phi
+      EtaPhiDisplaced = 3
+    }; //> end enum GridDisplacement
+
+    /**
+   * @brief Parameters describing a grid
+   * 
+   * How the grid divides up the detector is described mainly by the extent in
+   * eta (phi is by definition 2*pi) and the number of divisions in eta and phi.
+   * However, the grids must also support being displaced slightly in eta and
+   * phi, so these members are also required.
+   */
+    struct GridParameters
+    {
+      /// The maximum |eta| value
+      double maxEta{0.};
+      /// The number of divisions along the eta axis
+      std::size_t nEtaTowers{0};
+      /// The number of divisions along the phi axis
+      std::size_t nPhiTowers{0};
+      /// Whether the grid is displaced in eta
+      bool displaceEta{false};
+      /// Whether the grid is displaced in phi
+      bool displacePhi{false};
+      /// Check inequality with other parameter sets
+      bool operator!=(const GridParameters &other) const;
+      bool operator==(const GridParameters &other) const;
+      /// The @see GridDisplacement enum
+      GridDisplacement displacement() const;
+    }; //> end struct GridParameters
+
+    /**
+   * @brief Base class for grids used in some of the pufit algorithms
+   * 
+   * Many versions of the pufit technique require breaking the calorimeter into
+   * a grid in order to measure the variance of energy deposits across the
+   * detector and to identify areas which have a significantly high energy.
+   * 
+   * These grids must also support being displaced slightly in eta and phi to
+   * ensure that high energy deposits are not artificially divided between two
+   * towers.
+   * 
+   * The grids are also periodic in phi - i.e. the tower located at phi is the
+   * same one located at phi + 2pi.
+   */
+    class PeriodicGridBase
+    {
+    public:
+      /**
+     * @brief Base class for towers belonging to the grids
+     */
+      class Tower
+      {
+      public:
+        Tower(std::size_t index);
+        /// The grid which owns this tower
+        virtual const PeriodicGridBase *grid() const = 0;
+
+        /// The central phi coordinate of this tower
+        double towerPhi() const;
+
+        /// The central eta coordinate of this tower
+        double towerEta() const;
+
+        /// The global index of this tower
+        std::size_t index() const;
+
+        /// The eta index of this tower
+        std::size_t etaIndex() const;
+
+        /// The phi index of this tower
+        std::size_t phiIndex() const;
+
+        /// The eta/phi indices of the tower together
+        std::pair<std::size_t, std::size_t> etaPhiIndex() const;
+
+      private:
+        std::size_t m_index;
+      }; //> end class Tower
+
+      /// Construct the grid from its parameters
+      PeriodicGridBase(const GridParameters &parameters);
+
+      /// Construct the grid from its parameters provided directly
+      PeriodicGridBase(
+          double maxEta,
+          std::size_t nEtaTowers,
+          std::size_t nPhiTowers,
+          bool displaceEta = false,
+          bool displacePhi = false);
+
+      /**
+     * @brief Get the index for the given eta, phi values.
+     * @param eta The eta value
+     * @param phi The phi value
+     * @param[out] outOfRange Set to true if outside of the eta range
+     *
+     * If it's out of range the returned index will be nTowers
+     */
+      std::size_t getIndex(
+          double eta,
+          double phi,
+          bool &outOfRange) const;
+
+      /**
+     * @brief Get the eta index for the given value
+     * @param eta The eta value
+     * @param[out] outOfRange Set to true if outside of the eta range
+     *
+     * If it's out of range the returned index will be nEta
+     */
+      std::size_t getEtaIndex(double eta, bool &outOfRange) const;
+
+      /// Get the phi index for the given value
+      std::size_t getPhiIndex(double phi) const;
+
+      /// Convert eta and phi to a global index
+      std::size_t globalIndex(std::size_t iEta, std::size_t iPhi) const;
+
+      /// Convert a global index to an eta/phi index pair
+      std::pair<std::size_t, std::size_t> etaPhiIndex(std::size_t index) const;
+
+      /// Central eta coordinate of the given eta index
+      double centralEta(std::size_t iEta) const;
+      /// Central phi coordinate of the given phi
+      double centralPhi(std::size_t iPhi) const;
+
+      /// The grid parameters
+      const GridParameters &parameters() const;
+
+      /// The maximum eta range for the grid
+      double maxEta() const;
+
+      /// The number of eta bins
+      std::size_t nEtaTowers() const;
+
+      /// The number of phi bins
+      std::size_t nPhiTowers() const;
+
+      /// The number of bins
+      std::size_t nTowers() const;
+
+      /// Whether or not this is displaced in eta
+      bool displaceEta() const;
+
+      /// Whether or not this is displaced in phi
+      bool displacePhi() const;
+
+      /// The grid displacement
+      GridDisplacement displacement() const;
+
+      /// The bin width in eta
+      double etaWidth() const;
+
+      /// The bin width in phi
+      double phiWidth() const;
+
+    private:
+      /// The grid's parameters
+      const GridParameters m_params;
+    }; //> end class PeriodicGridBase
+  }    // namespace MET
+} // namespace HLT
+
+#endif //> !PUFITUTILS_PERIODICGRIDBASE_H
\ No newline at end of file
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitGrid.h b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitGrid.h
index ddb0906a91c1..07f500c6bbb7 100644
--- a/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitGrid.h
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitGrid.h
@@ -16,13 +16,17 @@
 #include <vector>
 #include <utility>
 #include <array>
+#include "TrigEFMissingET/PeriodicGridBase.h"
 #include "TrigEFMissingET/SignedKinematics.h"
 
-namespace HLT { namespace MET {
-  // Forward declare
-  class METComponent;
+namespace HLT
+{
+  namespace MET
+  {
+    // Forward declare
+    class METComponent;
 
-  /**
+    /**
    * @brief Bins energy deposits into a grid.
    *
    * The grid can be displaced in eta and/or phi. In this case the grid will
@@ -30,7 +34,8 @@ namespace HLT { namespace MET {
    * the grid will be treated as periodic in eta, even though the
    * calorimeter isn't really.
    */
-  class PufitGrid {
+    class PufitGrid : public PeriodicGridBase
+    {
     public:
       /**
        * @brief Describes a single element of the grid.
@@ -39,133 +44,103 @@ namespace HLT { namespace MET {
        * which result from the sums of objects going into the towers and the
        * tower coordinates themselves.
        */
-      class Tower {
-        public:
-          /**
+      class Tower : public PeriodicGridBase::Tower
+      {
+      public:
+        /**
            * @brief Create a tower with its parent
            * @param parent The parent grid of this tower
            * @param index The index of this tower in its parent's grid
            */
-          Tower(const PufitGrid* parent, std::size_t index);
+        Tower(const PufitGrid *parent, std::size_t index);
 
-          Tower(const Tower&) = default;
+        Tower(const Tower &) = default;
 
-          /**
+        /**
            * @brief Copy assignment operator
            *
            * This will not copy the other tower's parent. Each tower's parent is
            * fixed from construction! This only copies the energy and masking
            * information over.
            */
-          Tower& operator=(const Tower& other);
-          
-          /// The x-component of this tower's energy
-          double ex() const;
-
-          /// The y-component of this tower's energy
-          double ey() const;
-
-          /// The z-component of this tower's energy
-          double ez() const;
-
-          /// The total sumEt in this tower
-          double sumEt() const;
-
-          /// The total sumE in this tower
-          double sumE() const;
-
-          /// This tower's kinematic phi
-          double phi() const;
-
-          /// This tower's kinematic eta
-          double eta() const;
-
-          /// Whether or not this tower was masked
-          bool masked() const;
-
-          /// Set the mask on this tower
-          void mask(bool value=true);
-
-          /// The central phi of this tower
-          double towerPhi() const;
-
-          /// The central eta of this tower
-          double towerEta() const;
-
-          /// The global index of this tower
-          std::size_t index() const;
-
-          /// The parent grid of this tower
-          const PufitGrid* grid() const;
-
-          /// The eta index of this tower
-          std::size_t etaIndex() const;
-
-          /// The phi index of this tower
-          std::size_t phiIndex() const;
-
-          /// The eta-phi indices of this tower
-          std::pair<std::size_t, std::size_t> etaPhiIndex() const;
-
-          /// Build a kinematics object from this tower
-          SignedKinematics kinematics() const;
-          /// Conversion operator
-          operator SignedKinematics() const;
-
-          /**
-           * @brief Add a signed object to this tower.
-           *
-           * No check is done to make sure that this is the right bin.
-           * The *momentum* of the kinematics will be used. If you do not want
-           * this (i.e. you want mass to be ignored) you should construct the
-           * kinematics without mass.
-           */
-          Tower& operator+=(const SignedKinematics& kin);
-          /**
-           * @brief Remove the energy of a signed object from this tower.
-           *
-           * No check is done to make sure that this is the right bin.
-           * The *momentum* of the kinematics will be used. If you do not want
-           * this (i.e. you want mass to be ignored) you should construct the
-           * kinematics without mass.
-           */
-          Tower& operator-=(const SignedKinematics& kin);
-
-          /**
-           * @brief Add another tower's energies into this one.
-           */
-          Tower& operator+=(const Tower& other);
-          /**
-           * @brief Subtract another tower's energies from this one.
-           */
-          Tower& operator-=(const Tower& other);
-        private:
-          /// The parent grid
-          const PufitGrid* const m_parent;
-
-          /// The index in the parent
-          const std::size_t m_index;
-
-          SignedKinematics m_kin;
-          /// The summed et
-          double m_sumEt{0.};
-          /// The summed energy
-          double m_sumE{0.};
-          /// The mask value
-          bool m_mask{false};
+        Tower &operator=(const Tower &other);
+
+        /// The x-component of this tower's energy
+        double ex() const;
+
+        /// The y-component of this tower's energy
+        double ey() const;
+
+        /// The z-component of this tower's energy
+        double ez() const;
+
+        /// The total sumEt in this tower
+        double sumEt() const;
+
+        /// The total sumE in this tower
+        double sumE() const;
+
+        /// This tower's kinematic phi
+        double phi() const;
+
+        /// This tower's kinematic eta
+        double eta() const;
+
+        /// Whether or not this tower was masked
+        bool masked() const;
+
+        /// Set the mask on this tower
+        void mask(bool value = true);
+
+        /// The parent grid of this tower
+        virtual const PufitGrid *grid() const override;
+
+        /// Build a kinematics object from this tower
+        SignedKinematics kinematics() const;
+        /// Conversion operator
+        operator SignedKinematics() const;
+
+        /**
+         * @brief Add a signed object to this tower.
+         *
+         * No check is done to make sure that this is the right bin.
+         * The *momentum* of the kinematics will be used. If you do not want
+         * this (i.e. you want mass to be ignored) you should construct the
+         * kinematics without mass.
+         */
+        Tower &operator+=(const SignedKinematics &kin);
+        /**
+         * @brief Remove the energy of a signed object from this tower.
+         *
+         * No check is done to make sure that this is the right bin.
+         * The *momentum* of the kinematics will be used. If you do not want
+         * this (i.e. you want mass to be ignored) you should construct the
+         * kinematics without mass.
+         */
+        Tower &operator-=(const SignedKinematics &kin);
+
+        /**
+         * @brief Add another tower's energies into this one.
+         */
+        Tower &operator+=(const Tower &other);
+        /**
+         * @brief Subtract another tower's energies from this one.
+         */
+        Tower &operator-=(const Tower &other);
+
+      private:
+        /// The parent grid
+        const PufitGrid *const m_parent;
+
+        SignedKinematics m_kin;
+        /// The summed et
+        double m_sumEt{0.};
+        /// The summed energy
+        double m_sumE{0.};
+        /// The mask value
+        bool m_mask{false};
       }; //> end class Tower
 
-      /// The grid parameters
-      struct GridParameters {
-        double maxEta;
-        std::size_t nEtaTowers;
-        std::size_t nPhiTowers;
-        bool displaceEta{false};
-        bool displacePhi{false};
-        bool operator!=(const GridParameters& other) const;
-        bool operator==(const GridParameters& other) const;
-      };
-
       /**
        * @brief Create a new tower grid
        * @param maxEta The maximum eta range for the grid
@@ -182,10 +157,10 @@ namespace HLT { namespace MET {
           bool displacePhi = false);
 
       /// Construct a grid from the provided parameters
-      PufitGrid(const GridParameters& parameters);
+      PufitGrid(const GridParameters &parameters);
 
       /// Copy constructor
-      PufitGrid(const PufitGrid& other);
+      PufitGrid(const PufitGrid &other);
 
       /**
        * @brief Assignment operator
@@ -195,7 +170,7 @@ namespace HLT { namespace MET {
        * Take the tower energies/masking from the other grid.
        * This is only allowed between two grids with matching parameters.
        */
-      PufitGrid& operator=(const PufitGrid& other);
+      PufitGrid &operator=(const PufitGrid &other);
 
       /// Reset the internal storage
       void reset();
@@ -209,7 +184,7 @@ namespace HLT { namespace MET {
        * this (i.e. you want mass to be ignored) you should construct the
        * kinematics without mass.
        */
-      PufitGrid& operator+=(const SignedKinematics& kin);
+      PufitGrid &operator+=(const SignedKinematics &kin);
       /**
        * @brief Remove the energy of a signed object from this grid.
        *
@@ -219,43 +194,18 @@ namespace HLT { namespace MET {
        * this (i.e. you want mass to be ignored) you should construct the
        * kinematics without mass.
        */
-      PufitGrid& operator-=(const SignedKinematics& kin);
-
-      /**
-       * @brief Get the tower index for the given eta, phi values.
-       * @param eta The eta value
-       * @param phi The phi value
-       * @param[out] outOfRange Set to true if outside of the eta range
-       *
-       * If it's out of range the returned index will be nTowers
-       */
-      std::size_t getIndex(
-          double eta,
-          double phi,
-          bool& outOfRange) const;
+      PufitGrid &operator-=(const SignedKinematics &kin);
 
-      /**
-       * @brief Get the eta tower index for the given value
-       * @param eta The eta value
-       * @param[out] outOfRange Set to true if outside of the eta range
-       *
-       * If it's out of range the returned index will be nEta
-       */
-      std::size_t getEtaIndex(double eta, bool& outOfRange) const;
+      /// Access stored value by eta/phi index (access is bounds-checked)
+      Tower &operator[](const std::pair<std::size_t, std::size_t> &indices);
+      /// Access stored value by eta/phi index (access is bounds checked)
+      const Tower &operator[](
+          const std::pair<std::size_t, std::size_t> &indices) const;
 
-      /// Get the phi tower index for the given value
-      std::size_t getPhiIndex(double phi) const;
-
-      /// Access stored value by eta/phi index
-      Tower& operator[](const std::pair<std::size_t, std::size_t>& indices);
-      /// Access stored value by eta/phi index
-      const Tower& operator[](
-          const std::pair<std::size_t, std::size_t>& indices) const;
-
-      /// Access stored value by global index number
-      Tower& operator[](std::size_t index);
-      /// Access stored value by global index number
-      const Tower& operator[](std::size_t index) const;
+      /// Access stored value by global index number (access is bounds checked)
+      Tower &operator[](std::size_t index);
+      /// Access stored value by global index number (access is bounds checked)
+      const Tower &operator[](std::size_t index) const;
 
       /// Access by iterator
       std::vector<Tower>::iterator begin();
@@ -266,86 +216,44 @@ namespace HLT { namespace MET {
       /// Iterator end point
       std::vector<Tower>::const_iterator end() const;
 
-      /// Convert eta and phi index to global index
-      std::size_t globalIndex(std::size_t iEta, std::size_t iPhi) const;
-      /// Convert global index to eta and phi indices
-      std::pair<std::size_t, std::size_t> etaPhiIndex(std::size_t index) const;
-
-      /// Central eta coordinate of the given eta index
-      double centralEta(std::size_t iEta) const;
-      /// Central phi coordinate of the given phi tower
-      double centralPhi(std::size_t iPhi) const;
-
-      /// The grid parameters
-      const GridParameters& parameters() const;
-
-      /// The maximum eta range for the grid
-      double maxEta() const;
-
-      /// The number of eta towers
-      std::size_t nEtaTowers() const;
-
-      /// The number of phi towers
-      std::size_t nPhiTowers() const;
-
-      /// The number of towers
-      std::size_t nTowers() const;
-
-      /// Whether or not this is displaced in eta
-      bool displaceEta() const;
-
-      /// Whether or not this is displaced in phi
-      bool displacePhi() const;
-
-      /// The bin width in eta
-      double etaWidth() const;
-
-      /// The bin width in phi
-      double phiWidth() const;
-
       /// Helper enum to describe how to sum over towers
-      enum class SumStrategy {
-        All, //> Sum all towers, regardless of masks
-        Masked, //> Only sum masked towers
+      enum class SumStrategy
+      {
+        All,     //> Sum all towers, regardless of masks
+        Masked,  //> Only sum masked towers
         Unmasked //> Only sum unmasked towers
       };
-      METComponent sum(SumStrategy strategy=SumStrategy::All) const;
+      METComponent sum(SumStrategy strategy = SumStrategy::All) const;
 
       // All of these fail if the grid parameters do not match!
       /// Add a whole grid into this.
-      PufitGrid& operator+=(const PufitGrid& other);
+      PufitGrid &operator+=(const PufitGrid &other);
       /// Subtract a whole grid from this.
-      PufitGrid& operator-=(const PufitGrid& other);
+      PufitGrid &operator-=(const PufitGrid &other);
+
     private:
-      const GridParameters m_params;
-      double m_etaWidth;
-      double m_phiWidth;
       std::vector<Tower> m_towers;
-  }; //> end class PufitGrid
-
-  /// Helper struct to contain a full set of grids
-  struct PufitGridSet {
-    enum GridDisplacement {
-      NoDisplacement = 0,
-      EtaDisplaced = 1,
-      PhiDisplaced = 2,
-      EtaPhiDisplaced = 3
+    }; //> end class PufitGrid
+
+    /// Helper struct to contain a full set of grids
+    struct PufitGridSet
+    {
+      PufitGridSet(double maxEta, std::size_t nEta, std::size_t nPhi);
+      std::array<PufitGrid, 4> grids;
+      /// Add kinematics
+      PufitGridSet &operator+=(const SignedKinematics &kin);
+      /// Subtract kinematics
+      PufitGridSet &operator-=(const SignedKinematics &kin);
+      /// Select a grid
+      PufitGrid &operator[](GridDisplacement displacement);
+      const PufitGrid &operator[](GridDisplacement displacement) const;
     };
-    PufitGridSet(double maxEta, std::size_t nEta, std::size_t nPhi);
-    std::array<PufitGrid, 4> grids;
-    /// Add kinematics
-    PufitGridSet& operator+=(const SignedKinematics& kin);
-    /// Subtract kinematics
-    PufitGridSet& operator-=(const SignedKinematics& kin);
-    /// Select a grid
-    PufitGrid& operator[](GridDisplacement displacement);
-    const PufitGrid& operator[](GridDisplacement displacement) const;
-  };
-
-  /// Elementwise sum
-  PufitGrid operator+(const PufitGrid& lhs, const PufitGrid& rhs);
-  /// Elementwise subtraction
-  PufitGrid operator-(const PufitGrid& lhs, const PufitGrid& rhs);
-} } //> end namespace HLT::MET
+
+    /// Elementwise sum
+    PufitGrid operator+(const PufitGrid &lhs, const PufitGrid &rhs);
+    /// Elementwise subtraction
+    PufitGrid operator-(const PufitGrid &lhs, const PufitGrid &rhs);
+  } // namespace MET
+} // namespace HLT
 
 #endif //> !TRIGEFMISSINGET_PUFITGRID_H
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitMultiGrid.h b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitMultiGrid.h
new file mode 100644
index 000000000000..e3579c36d907
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitMultiGrid.h
@@ -0,0 +1,351 @@
+/*
+ * Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+#ifndef TRIGEFMISSINGET_PUFITMULTIGRID_H
+#define TRIGEFMISSINGET_PUFITMULTIGRID_H
+
+#include "TrigEFMissingET/PeriodicGridBase.h"
+#include "TrigEFMissingET/PufitGrid.h"
+#include <limits>
+#include <type_traits>
+
+namespace HLT
+{
+  namespace MET
+  {
+    /// Compile time check if a number is a power of 2
+    constexpr bool isPow2(std::size_t i)
+    {
+      // An int is a power of 2 if it only contains one 'on' bit
+      // We check for this by shifting an int until its first bit is 1. Then it
+      // is a power of two if that is the only bit set.
+      // We also need a check if this is 0 otherwise the recursion will never
+      // end
+      if (i == 0)
+        // 0 is *not* a power of 2
+        return false;
+      else if (i == 1)
+        // 1 is a power of 2
+        return true;
+      else if (i & 1)
+        // If we're not 1, but the '1' bit is set then this can't be a power of 2
+        return false;
+      else
+        // shifting to the right esse
+        return isPow2(i >> 1);
+    }
+
+    /**
+     * @brief Compile time calculation of the log base 2 of an integer
+     * 
+     * Returns the floor of the log, so intLog2(9) = 3. Log base 2 of 0 is
+     * undefined and here returns the maximum value for std::size_t.
+     */
+    constexpr uint16_t intLog2(std::size_t i, uint16_t tmp = 0)
+    {
+      if (i == 0)
+        return UINT16_MAX;
+      else if (i == 1)
+        return tmp;
+      else
+        return intLog2(i >> 1, tmp + 1);
+    }
+
+    /**
+     * @brief Multiple grids combined into one
+     * 
+     * This class allows representing N different categories in a single grid.
+     */
+    template <std::size_t N>
+    class PufitMultiGrid : public PeriodicGridBase
+    {
+      static_assert(N > 0, "N must be greater than 0");
+      // The C++ standard only guarantees that std::size_t has a bit width of
+      // at least 16
+      static_assert(N < 17, "N must be no greater than 16");
+
+    public:
+      /// Maximum value representable by N bits
+      constexpr static std::size_t All = (1 << N) - 1;
+      /// The number of separate categories in the grid
+      constexpr static std::size_t NCategories = N;
+
+      /**
+       * @brief Tower in the multi-grid
+       * 
+       * Each tower is built from N 'normal' PufitGrid::Towers. Kinematic
+       * quantities can be read from the tower either from any of the
+       * 'subtowers' or from a sum over them, passing in a std::size_t as a bitmask.
+       */
+      class Tower : public PeriodicGridBase::Tower
+      {
+        friend class PufitMulitGrid;
+
+      public:
+        /**
+         * @brief Create a tower with its parent grid
+         * @param parent The parent grid of this tower
+         * @param index The index of this tower in its parent's grid
+         */
+        Tower(PufitMultiGrid *parent, std::size_t index);
+
+        /**
+         * @brief Copy assignment operator
+         *
+         * This will not copy the other tower's parent. Each tower's parent is
+         * fixed from construction! This only copies the energy and masking
+         * information over.
+         */
+        Tower &operator=(const Tower &other);
+
+        /// The x-component of this tower's energy
+        double ex(std::size_t type = All) const;
+
+        /// The y-component of this tower's energy
+        double ey(std::size_t type = All) const;
+
+        /// The z-component of this tower's energy
+        double ez(std::size_t type = All) const;
+
+        /// The total sumEt in this tower
+        double sumEt(std::size_t type = All) const;
+
+        /// The total sumE in this tower
+        double sumE(std::size_t type = All) const;
+
+        /// This tower's kinematic phi
+        double phi(std::size_t type = All) const;
+
+        /// This tower's kinematic eta
+        double eta(std::size_t type = All) const;
+
+        /// Whether or not this tower was masked
+        bool masked() const;
+
+        /// Set the mask on this tower
+        void mask(bool value = true);
+
+        /// The parent grid of this tower
+        const PufitMultiGrid *grid() const override;
+
+        /// Build a kinematics object from this tower
+        SignedKinematics kinematics(std::size_t type = All) const;
+
+        /// Conversion operator, returns the value for the sum of all the
+        /// grids
+        operator SignedKinematics() const;
+
+        /**
+         * @brief Add another tower's energies into this one.
+         */
+        Tower &operator+=(const Tower &other);
+        /**
+         * @brief Subtract another tower's energies from this one.
+         */
+        Tower &operator-=(const Tower &other);
+
+        /**
+         * @brief Get one of the underlying towers
+         *
+         * This is only defined where I is a power of 2
+         */
+        template <std::size_t I,
+                  typename = typename std::enable_if<isPow2(I)>::type,
+                  typename = typename std::enable_if<I <= All>::type>
+        PufitGrid::Tower &get()
+        {
+          return subTower(intLog2(I));
+        }
+
+        /**
+         * @brief Get one of the underlying towers
+         *
+         * This is only defined where I is a power of 2
+         */
+        template <std::size_t I,
+                  typename = typename std::enable_if<isPow2(I)>::type,
+                  typename = typename std::enable_if<I <= All>::type>
+        const PufitGrid::Tower &get() const
+        {
+          return subTower(intLog2(I));
+        }
+
+      private:
+        /// The parent grid
+        PufitMultiGrid *const m_parent;
+
+        /// The mask value
+        bool m_mask{false};
+
+        /// Get a subtower by index
+        PufitGrid::Tower &subTower(std::size_t ii)
+        {
+          return m_parent->m_grids[ii][index()];
+        }
+        /// Get a subtower by index
+        const PufitGrid::Tower &subTower(std::size_t ii) const
+        {
+          return m_parent->m_grids[ii][index()];
+        }
+
+        /// Apply a function to all sub towers
+        void applyToAll(PufitGrid::Tower &(PufitGrid::Tower::*f)(const PufitGrid::Tower &), const PufitMultiGrid::Tower &other)
+        {
+          for (std::size_t ii = 0; ii < N; ++ii)
+            (subTower(ii).*f)(other.subTower(ii));
+        }
+
+        /// Sum over the results of all sub towers whose indices match the
+        /// 'type' mask
+        template <typename T>
+        typename std::decay<T>::type sumOver(int type, T (PufitGrid::Tower::*f)() const) const
+        {
+          typename std::decay<T>::type val{};
+          for (std::size_t ii = 0; ii < N; ++ii)
+            if (1 << ii & type)
+              val += (subTower(ii).*f)();
+          return val;
+        }
+
+      }; //> end class Tower
+
+      /**
+       * @brief Create a new tower grid
+       * @param maxEta The maximum eta range for the grid
+       * @param nEtaTowers The number of eta towers
+       * @param nPhiTowers The number of phi towers
+       * @param displaceEta Whether to displace eta
+       * @param displacePhi Whether to displace phi
+       */
+      PufitMultiGrid(
+          double maxEta,
+          std::size_t nEtaTowers,
+          std::size_t nPhiTowers,
+          bool displaceEta = false,
+          bool displacePhi = false);
+
+      /// Construct a grid from the provided parameters
+      PufitMultiGrid(const GridParameters &parameters);
+
+      /// Copy constructor
+      PufitMultiGrid(const PufitMultiGrid &other);
+
+      /**
+       * @brief Assignment operator
+       * @param other The grid whose towers to take
+       *
+       * Take the tower energies/masking from the other grid.
+       */
+      PufitMultiGrid &operator=(const PufitMultiGrid &other);
+
+      /// Reset the internal storage
+      void reset();
+
+      /// Access stored value by eta/phi index (access is bounds checked)
+      Tower &operator[](const std::pair<std::size_t, std::size_t> &indices);
+      /// Access stored value by eta/phi index (access is bounds checked)
+      const Tower &operator[](const std::pair<std::size_t, std::size_t> &indices) const;
+
+      /// Access stored value by global index number (access is bounds checked)
+      Tower &operator[](std::size_t index);
+      /// Access stored value by global index number (access is bounds checked)
+      const Tower &operator[](std::size_t index) const;
+
+      /// Access by iterator
+      typename std::vector<Tower>::iterator begin();
+      /// Access by iterator
+      typename std::vector<Tower>::const_iterator begin() const;
+      /// Iterator end point
+      typename std::vector<Tower>::iterator end();
+      /// Iterator end point
+      typename std::vector<Tower>::const_iterator end() const;
+
+      // All of these fail if the grid parameters do not match!
+      /// Add a whole grid into this.
+      PufitMultiGrid &operator+=(const PufitMultiGrid &other);
+      /// Subtract a whole grid from this.
+      PufitMultiGrid &operator-=(const PufitMultiGrid &other);
+
+      /**
+       * @brief Get one of the underlying grids
+       *
+       * This is only defined where I is a power of 2
+       */
+      template <std::size_t I,
+                typename = typename std::enable_if<isPow2(I)>::type,
+                typename = typename std::enable_if<I <= All>::type>
+      PufitGrid &get()
+      {
+        return m_grids[intLog2(I)];
+      }
+
+      /**
+       * @brief Get one of the underlying grids
+       *
+       * This is only defined where I is a power of 2
+       */
+      template <std::size_t I,
+                typename = typename std::enable_if<isPow2(I)>::type,
+                typename = typename std::enable_if<I <= All>::type>
+      const PufitGrid &get() const
+      {
+        return m_grids[intLog2(I)];
+      }
+
+      /**
+       * @brief Get one of the underlying grids
+       *
+       * This version is allowed to return a grid where idx is not a power of
+       * 2, in which case the grid will be constructed on the fly from the sum
+       * of its bitmask
+       */
+      PufitGrid get(std::size_t type) const;
+
+    private:
+      std::array<PufitGrid, N> m_grids;
+      std::vector<Tower> m_towers;
+    }; //> end class PufitMultiGrid
+
+    /// Helper struct to forward the SignedKinematics operators nicely
+    template <typename Grid>
+    struct PufitMultiGridSet
+    {
+      PufitMultiGridSet(double maxEta, std::size_t nEta, std::size_t nPhi);
+      std::array<Grid, 4> grids;
+
+      /// Select a grid
+      Grid &operator[](GridDisplacement displacement)
+      {
+        return grids[displacement];
+      }
+
+      /// Select a grid (const)
+      const Grid &operator[](GridDisplacement displacement) const
+      {
+        return grids[displacement];
+      }
+
+      /// Helper struct to forward the SignedKinematics operators nicely
+      template <std::size_t I>
+      struct Element
+      {
+        Element(PufitMultiGridSet &parent);
+        PufitMultiGridSet &parent;
+        Element &operator+=(const SignedKinematics &kin);
+        Element &operator-=(const SignedKinematics &kin);
+      }; //> end struct Element<I>
+
+      template <std::size_t I,
+                typename = typename std::enable_if<isPow2(I)>::type,
+                typename = typename std::enable_if<I <= Grid::All>::type>
+      Element<I> get();
+
+      PufitGridSet get(std::size_t type) const;
+    }; //> end struct PufitMultiGridSet<Grid>
+  }    // namespace MET
+} // namespace HLT
+
+#include "TrigEFMissingET/PufitMultiGrid.icc"
+
+#endif //> !TRIGEFMISSINGET_PUFITMULTIGRID_H
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitMultiGrid.icc b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitMultiGrid.icc
new file mode 100644
index 000000000000..21762f703891
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitMultiGrid.icc
@@ -0,0 +1,313 @@
+/*
+ * Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+#ifndef TRIGEFMISSINGET_PUFITMULTIGRID_ICC
+#define TRIGEFMISSINGET_PUFITMULTIGRID_ICC
+
+namespace HLT
+{
+  namespace MET
+  {
+    namespace detail
+    {
+      /// Initialise an array (thank you stack overflow https://stackoverflow.com/a/41259045/8434292)
+      template <typename T, std::size_t... Is>
+      constexpr std::array<T, sizeof...(Is)>
+      make_array(const T &value, std::index_sequence<Is...>)
+      {
+        return {{(static_cast<void>(Is), value)...}};
+      }
+    } // namespace detail
+
+    template <std::size_t N>
+    PufitMultiGrid<N>::Tower::Tower(PufitMultiGrid *parent, std::size_t index)
+        : PeriodicGridBase::Tower(index),
+          m_parent(parent)
+    {
+    }
+
+    template <std::size_t N>
+    typename PufitMultiGrid<N>::Tower &PufitMultiGrid<N>::Tower::operator=(const Tower &other)
+    {
+      applyToAll(&PufitGrid::Tower::operator=, other);
+      return *this;
+    }
+
+    template <std::size_t N>
+    double PufitMultiGrid<N>::Tower::ex(std::size_t type) const
+    {
+      return sumOver(type, &PufitGrid::Tower::ex);
+    }
+
+    template <std::size_t N>
+    double PufitMultiGrid<N>::Tower::ey(std::size_t type) const
+    {
+      return sumOver(type, &PufitGrid::Tower::ey);
+    }
+
+    template <std::size_t N>
+    double PufitMultiGrid<N>::Tower::ez(std::size_t type) const
+    {
+      return sumOver(type, &PufitGrid::Tower::ez);
+    }
+
+    template <std::size_t N>
+    double PufitMultiGrid<N>::Tower::sumEt(std::size_t type) const
+    {
+      return sumOver(type, &PufitGrid::Tower::sumEt);
+    }
+
+    template <std::size_t N>
+    double PufitMultiGrid<N>::Tower::sumE(std::size_t type) const
+    {
+      return sumOver(type, &PufitGrid::Tower::sumE);
+    }
+
+    template <std::size_t N>
+    double PufitMultiGrid<N>::Tower::phi(std::size_t type) const
+    {
+      return kinematics(type).phi();
+    }
+
+    template <std::size_t N>
+    double PufitMultiGrid<N>::Tower::eta(std::size_t type) const
+    {
+      return kinematics(type).eta();
+    }
+
+    template <std::size_t N>
+    bool PufitMultiGrid<N>::Tower::masked() const { return m_mask; }
+
+    template <std::size_t N>
+    void PufitMultiGrid<N>::Tower::mask(bool value)
+    {
+      m_mask = value;
+      for (std::size_t ii = 0; ii < N; ++ii)
+        subTower(ii).mask(value);
+    }
+
+    template <std::size_t N>
+    const PufitMultiGrid<N> *PufitMultiGrid<N>::Tower::grid() const { return m_parent; }
+
+    template <std::size_t N>
+    SignedKinematics PufitMultiGrid<N>::Tower::kinematics(std::size_t type) const
+    {
+      return sumOver(type, &PufitGrid::Tower::kinematics);
+    }
+
+    template <std::size_t N>
+    PufitMultiGrid<N>::Tower::operator SignedKinematics() const { return kinematics(); }
+
+    template <std::size_t N>
+    typename PufitMultiGrid<N>::Tower &PufitMultiGrid<N>::Tower::operator+=(const Tower &other)
+    {
+      applyToAll(&PufitGrid::Tower::operator+=, other);
+      return *this;
+    }
+
+    template <std::size_t N>
+    typename PufitMultiGrid<N>::Tower &PufitMultiGrid<N>::Tower::operator-=(const Tower &other)
+    {
+      applyToAll(&PufitGrid::Tower::operator-=, other);
+      return *this;
+    }
+
+    template <std::size_t N>
+    PufitMultiGrid<N>::PufitMultiGrid(
+        double maxEta,
+        std::size_t nEtaTowers,
+        std::size_t nPhiTowers,
+        bool displaceEta,
+        bool displacePhi)
+        : PufitMultiGrid(
+              GridParameters{maxEta, nEtaTowers, nPhiTowers, displaceEta, displacePhi})
+    {
+    }
+
+    template <std::size_t N>
+    PufitMultiGrid<N>::PufitMultiGrid(const GridParameters &parameters)
+        : PeriodicGridBase(parameters),
+          m_grids(detail::make_array(PufitGrid(parameters), std::make_index_sequence<N>()))
+    {
+      m_towers.reserve(nTowers());
+      for (std::size_t index = 0; index < nTowers(); ++index)
+        m_towers.emplace_back(this, index);
+    }
+
+    template <std::size_t N>
+    PufitMultiGrid<N>::PufitMultiGrid(const PufitMultiGrid &other)
+        : PufitMultiGrid(other.parameters())
+    {
+      *this = other;
+    }
+
+    template <std::size_t N>
+    PufitMultiGrid<N> &PufitMultiGrid<N>::operator=(const PufitMultiGrid &other)
+    {
+      if (parameters() != other.parameters())
+        throw std::invalid_argument("Grid parameters do not match");
+      std::copy(other.begin(), other.end(), m_towers.begin());
+      std::copy(other.m_grids.begin(), other.m_grids.end(), m_grids.begin());
+      return *this;
+    }
+
+    template <std::size_t N>
+    void PufitMultiGrid<N>::reset()
+    {
+      std::fill(begin(), end(), Tower(nullptr, -1));
+      for (PufitGrid &grid : m_grids)
+        grid.reset();
+    }
+    template <std::size_t N>
+    typename PufitMultiGrid<N>::Tower &PufitMultiGrid<N>::operator[](
+        const std::pair<std::size_t, std::size_t> &indices)
+    {
+      return operator[](globalIndex(indices.first, indices.second));
+    }
+    template <std::size_t N>
+    const typename PufitMultiGrid<N>::Tower &PufitMultiGrid<N>::operator[](
+        const std::pair<std::size_t, std::size_t> &indices) const
+    {
+      return operator[](globalIndex(indices.first, indices.second));
+    }
+
+    template <std::size_t N>
+    typename PufitMultiGrid<N>::Tower &PufitMultiGrid<N>::operator[](std::size_t index)
+    {
+      return m_towers.at(index);
+    }
+    template <std::size_t N>
+    const typename PufitMultiGrid<N>::Tower &PufitMultiGrid<N>::operator[](std::size_t index) const
+    {
+      return m_towers.at(index);
+    }
+
+    template <std::size_t N>
+    typename std::vector<typename PufitMultiGrid<N>::Tower>::iterator PufitMultiGrid<N>::begin()
+    {
+      return m_towers.begin();
+    }
+    template <std::size_t N>
+    typename std::vector<typename PufitMultiGrid<N>::Tower>::const_iterator PufitMultiGrid<N>::begin() const
+    {
+      return m_towers.begin();
+    }
+    template <std::size_t N>
+    typename std::vector<typename PufitMultiGrid<N>::Tower>::iterator PufitMultiGrid<N>::end()
+    {
+      return m_towers.end();
+    }
+    template <std::size_t N>
+    typename std::vector<typename PufitMultiGrid<N>::Tower>::const_iterator PufitMultiGrid<N>::end() const
+    {
+      return m_towers.end();
+    }
+
+    template <std::size_t N>
+    PufitGrid PufitMultiGrid<N>::get(std::size_t type) const
+    {
+      if (isPow2(type))
+        return m_grids[intLog2(type)];
+      PufitGrid val(parameters());
+      for (std::size_t ii = 0; ii < N; ++ii)
+        if (1 << ii & type)
+          val += m_grids[ii];
+      return val;
+    }
+
+    template <typename Grid>
+    PufitMultiGridSet<Grid>::PufitMultiGridSet(
+        double maxEta, std::size_t nEta, std::size_t nPhi)
+        : grids({Grid(maxEta, nEta, nPhi, false, false),
+                 Grid(maxEta, nEta, nPhi, true, false),
+                 Grid(maxEta, nEta, nPhi, false, true),
+                 Grid(maxEta, nEta, nPhi, true, true)})
+    {
+    }
+
+    template <typename Grid>
+    template <std::size_t I>
+    PufitMultiGridSet<Grid>::Element<I>::Element(PufitMultiGridSet &parent)
+        : parent(parent) {}
+
+    template <typename Grid>
+    template <std::size_t I>
+    PufitMultiGridSet<Grid>::Element<I> &PufitMultiGridSet<Grid>::Element<I>::operator+=(
+        const SignedKinematics &kin)
+    {
+      for (Grid &grid : parent.grids)
+        grid.template get<I>() += kin;
+      return *this;
+    }
+
+    template <typename Grid>
+    template <std::size_t I>
+    PufitMultiGridSet<Grid>::Element<I> &PufitMultiGridSet<Grid>::Element<I>::operator-=(
+        const SignedKinematics &kin)
+    {
+      for (Grid &grid : parent.grids)
+        grid.template get<I>() -= kin;
+      return *this;
+    }
+
+    template <typename Grid>
+    template <std::size_t I, typename, typename>
+    PufitMultiGridSet<Grid>::Element<I> PufitMultiGridSet<Grid>::get()
+    {
+      return Element<I>(*this);
+    }
+
+    template <typename Grid>
+    PufitGridSet PufitMultiGridSet<Grid>::get(std::size_t type) const
+    {
+      PufitGridSet gridSet(grids[0].maxEta, grids[0].nEta, grids[0].nPhi);
+      for (std::size_t ii = 0; ii < 4; ++ii)
+        gridSet.grids[ii] = grids[ii].get(type);
+      return gridSet;
+    }
+
+    template <std::size_t N>
+    PufitMultiGrid<N> &PufitMultiGrid<N>::operator+=(const PufitMultiGrid &other)
+    {
+      if (parameters() != other.parameters())
+        throw std::invalid_argument("Grid parameters do not match");
+      auto itr = m_grids.begin();
+      auto otherItr = other.m_grids.begin();
+      for (; itr != end(); ++itr, ++otherItr)
+        *itr += *otherItr;
+      return *this;
+    }
+
+    template <std::size_t N>
+    PufitMultiGrid<N> &PufitMultiGrid<N>::operator-=(const PufitMultiGrid &other)
+    {
+      if (parameters() != other.parameters())
+        throw std::invalid_argument("Grid parameters do not match");
+      auto itr = m_grids.begin();
+      auto otherItr = other.m_grids.begin();
+      for (; itr != end(); ++itr, ++otherItr)
+        *itr -= *otherItr;
+      return *this;
+    }
+
+    template <std::size_t N>
+    PufitMultiGrid<N> operator+(const PufitMultiGrid<N> &lhs, const PufitMultiGrid<N> &rhs)
+    {
+      PufitMultiGrid<N> ret(lhs);
+      ret += rhs;
+      return ret;
+    }
+
+    template <std::size_t N>
+    PufitMultiGrid<N> operator-(const PufitMultiGrid<N> &lhs, const PufitMultiGrid<N> &rhs)
+    {
+      PufitMultiGrid<N> ret(lhs);
+      ret -= rhs;
+      return ret;
+    }
+
+  } // namespace MET
+} // namespace HLT
+
+#endif //> !TRIGEFMISSINGET_PUFITMULTIGRID_ICC
\ No newline at end of file
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitUtils.h b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitUtils.h
index e7d3c3488995..53e7698c83b8 100644
--- a/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitUtils.h
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitUtils.h
@@ -15,195 +15,255 @@
 #include <Eigen/Dense>
 #include "TrigEFMissingET/SignedKinematics.h"
 #include "TrigEFMissingET/PufitGrid.h"
+#include "TrigEFMissingET/PufitMultiGrid.h"
 #include <vector>
 
-namespace HLT { namespace MET { namespace PufitUtils {
-  /// Helper struct to hold the sum over pileup objects and its covariance
-  struct CovarianceSum {
-    /// Default constructor - zero initialize everything.
-    CovarianceSum();
-    /// Construct with an existing sum and matrix
-    CovarianceSum(const Eigen::Vector2d& sum, const Eigen::Matrix2d& covariance);
-    /**
-     * @brief Add a new contribution to the sum
-     * @brief kin The kinematics of the contribution
-     * @brief sigma The resolution of the contribution
-     */
-    CovarianceSum& add(const SignedKinematics& kin, double sigma);
-    /// The sum
-    Eigen::Vector2d sum;
-    /// The covariance matrix
-    Eigen::Matrix2d covariance;
-  };
-
-  /**
-   * @brief Calculate the trimmed mean and variance for a grid of towers
-   * @param grid The grid to calculate
-   * @param trimFraction The fraction of towers to use in the calculation.
-   * @return A pair containing the mean in 'first' and variance in 'second'
-   *
-   * Any masks applied to the grid will be ignored.
-   */
-  std::pair<double, double> trimmedMeanAndVariance(
-      const PufitGrid& grid,
-      double trimFraction);
-
-  /**
-   * @brief Calculate the mean and variance of unmasked towers
-   * @param grid The grid to calculate
-   * @return A pair containing the mean in 'first' and variance in 'second'
-   */
-  std::pair<double, double> unmaskedMeanAndVariance(const PufitGrid& grid);
-
-  /// Select the grid with the highest masked sumEt
-  PufitGridSet::GridDisplacement selectGrid(const PufitGridSet& grids);
-
-  /**
-   * @brief Perform the pile-up fit
-   * @param pileupSum 2-vector sum over the pileup contributions
-   * @param pileupCovariance Covariance matrix of the pileup sum
-   * @param towerExpectations Expected values for each tower based on the pileup
-   * energy density
-   * @param towerVariances Expected values for the variances of each tower based
-   * on variances of background towers
-   * @param correctionDirections The directions of the objects to correct
-   * @param The relative importance of the constraints.
-   * @return A vector containing the magnitudes of the corrections
-   *
-   * This variant provides individual estimates of the expected energy in each
-   * tower and the corresponding variances - these are for cases where the areas
-   * of each object are different.
-   */
-  Eigen::VectorXd pufit(
-      const Eigen::Vector2d& pileupSum,
-      const Eigen::Matrix2d& pileupCovariance,
-      const Eigen::VectorXd& towerExpectations,
-      const Eigen::VectorXd& towerVariances,
-      const Eigen::VectorXd& correctionDirections,
-      double constraintImportance = 1);
-
-  /**
-   * @brief Perform the pile-up fit
-   * @param pileupSum 2-vector sum over the pileup contributions
-   * @param pileupCovariance Covariance matrix of the pileup sum
-   * @param towerExpectations Expected values for each tower based on the pileup
-   * energy density
-   * @param towerVariances Expected values for the variances of each tower based
-   * on variances of background towers
-   * @param cosSin The cosines and sines of the objects to correct
-   * @param The relative importance of the constraints.
-   * @return A vector containing the magnitudes of the corrections
-   *
-   * This variant directly provides the cos-sin matrix for the correction,
-   * potentially allowing many cos/sin calls to be avoided.
-   * This variant provides individual estimates of the expected energy in each
-   * tower and the corresponding variances - these are for cases where the areas
-   * of each object are different.
-   */
-  Eigen::VectorXd pufit(
-      const Eigen::Vector2d& pileupSum,
-      const Eigen::Matrix2d& pileupCovariance,
-      const Eigen::VectorXd& towerExpectations,
-      const Eigen::VectorXd& towerVariances,
-      const Eigen::Matrix<double, 2, Eigen::Dynamic>& cosSin,
-      double constraintImportance = 1);
-
-  /**
-   * @brief Perform the pile-up fit
-   * @param pileupSum 2-vector sum over the pileup contributions
-   * @param pileupCovariance Covariance matrix of the pileup sum
-   * @param towerExpectations Expected values for each tower based on the pileup
-   * energy density
-   * @param towerVariances Expected values for the variances of each tower based
-   * on variances of background towers
-   * @param toCorrect Kinematics of the objects to correct
-   * @param The relative importance of the constraints.
-   * @return The corrections to apply
-   *
-   * Note that this version still returns the corrections that should be
-   * applied, not the corrected objects. The corrections are returned using the
-   * negatives of the value calculated by the fit - i.e. you should add the
-   * returned objects and let them act as negative energy objects (unless the
-   * sign of the correction was negative in which case it will act as a positive
-   * energy object).
-   * This variant provides individual estimates of the expected energy in each
-   * tower and the corresponding variances - these are for cases where the areas
-   * of each object are different.
-   */
-  std::vector<SignedKinematics> pufit(
-      const Eigen::Vector2d& pileupSum,
-      const Eigen::Matrix2d& pileupCovariance,
-      const std::vector<double>& towerExpectations,
-      const std::vector<double>& towerVariances,
-      const std::vector<SignedKinematics>& toCorrect,
-      double constraintImportance = 1);
-
-  /**
-   * @brief Perform the pile-up fit
-   * @param pileupSum 2-vector sum over the pileup contributions
-   * @param pileupCovariance Covariance matrix of the pileup sum
-   * @param towerMean Estimate of the mean pileup tower transverse energy
-   * @param towerVariance Estimate of the variance of pileup tower transverse
-   * energies
-   * @param correctionDirections The directions of the objects to correct
-   * @param The relative importance of the constraints.
-   * @return A vector containing the magnitudes of the corrections
-   */
-  Eigen::VectorXd pufit(
-      const Eigen::Vector2d& pileupSum,
-      const Eigen::Matrix2d& pileupCovariance,
-      double towerMean,
-      double towerVariance,
-      const Eigen::VectorXd& correctionDirections,
-      double constraintImportance = 1);
-
-  /**
-   * @brief Perform the pile-up fit
-   * @param pileupSum 2-vector sum over the pileup contributions
-   * @param pileupCovariance Covariance matrix of the pileup sum
-   * @param towerMean Estimate of the mean pileup tower transverse energy
-   * @param towerVariance Estimate of the variance of pileup tower transverse
-   * energies
-   * @param cosSin The cosines and sines of the objects to correct
-   * @param The relative importance of the constraints.
-   * @return A vector containing the magnitudes of the corrections
-   *
-   * This variant directly provides the cos-sin matrix for the correction,
-   * potentially allowing many cos/sin calls to be avoided.
-   */
-  Eigen::VectorXd pufit(
-      const Eigen::Vector2d& pileupSum,
-      const Eigen::Matrix2d& pileupCovariance,
-      double towerMean,
-      double towerVariance,
-      const Eigen::Matrix<double, 2, Eigen::Dynamic>& cosSin,
-      double constraintImportance = 1);
-
-  /**
-   * @brief Perform the pile-up fit
-   * @param pileupSum 2-vector sum over the pileup contributions
-   * @param pileupCovariance Covariance matrix of the pileup sum
-   * @param towerMean Estimate of the mean pileup tower transverse energy
-   * @param towerVariance Estimate of the variance of pileup tower transverse
-   * energies
-   * @param toCorrect Kinematics of the objects to correct
-   * @param The relative importance of the constraints.
-   * @return The corrections to apply
-   *
-   * Note that this version still returns the corrections that should be
-   * applied, not the corrected objects. The corrections are returned using the
-   * negatives of the value calculated by the fit - i.e. you should add the
-   * returned objects and let them act as negative energy objects (unless the
-   * sign of the correction was negative in which case it will act as a positive
-   * energy object).
-   */
-  std::vector<SignedKinematics> pufit(
-      const Eigen::Vector2d& pileupSum,
-      const Eigen::Matrix2d& pileupCovariance,
-      double towerMean,
-      double towerVariance,
-      const std::vector<SignedKinematics>& toCorrect,
-      double constraintImportance = 1);
-} } } //> end namespace HLT::MET::PufitUtils
+namespace HLT
+{
+    namespace MET
+    {
+        namespace PufitUtils
+        {
+            /// Helper struct to hold the sum over pileup objects and its covariance
+            struct CovarianceSum
+            {
+                /// Default constructor - zero initialize everything.
+                CovarianceSum();
+                /// Construct with an existing sum and matrix
+                CovarianceSum(const Eigen::Vector2d &sum, const Eigen::Matrix2d &covariance);
+                /**
+                 * @brief Add a new contribution to the sum
+                 * @brief kin The kinematics of the contribution
+                 * @brief sigma The resolution of the contribution
+                 */
+                CovarianceSum &add(const SignedKinematics &kin, double sigma);
+                /// The sum
+                Eigen::Vector2d sum;
+                /// The covariance matrix
+                Eigen::Matrix2d covariance;
+            };
+
+            /**
+             * @brief Calculate the trimmed mean and variance for a vector of tower sumEts
+             * @param sorted The sumEts for all towers, sorted in ascending order
+             * @param trimFraction The fraction of towers to use in the calculation
+             * @param[out] mean The calculated mean
+             * @param[out] variance The calculated variance
+             */
+            void trimmedMeanAndVariance(
+                const std::vector<double> &sorted,
+                double trimFraction,
+                double &mean,
+                double &variance);
+
+            /**
+             * @brief Calculate the trimmed mean and variance for a vector of tower sumEts
+             * @param grid The grid to calculate
+             * @param trimFraction The fraction of towers to use in the calculation
+             * @param[out] mean The calculated mean
+             * @param[out] variance The calculated variance
+             */
+            void trimmedMeanAndVariance(
+                const PufitGrid &grid,
+                double trimFraction,
+                double &mean,
+                double &variance);
+
+            /**
+             * @brief Calculate the trimmed mean and variance for a vector of tower sumEts
+             * @param grid The grid to calculate
+             * @param type The int mask of subgrid types
+             * @param trimFraction The fraction of twoers to use in the calculation
+             * @param[out] mean The calculated mean
+             * @param[out] variance The calculated variance
+             */
+            template <std::size_t N>
+            void trimmedMeanAndVariance(
+                const PufitMultiGrid<N> &grid,
+                std::size_t type,
+                double trimFraction,
+                double &mean,
+                double &variance);
+
+            /**
+             * @brief Calculate the mean and variance of unmasked towers
+             * @param grid The grid to calculate
+             * @param[out] mean The calculated mean
+             * @param[out] variance The calculated variance
+             */
+            void unmaskedMeanAndVariance(
+                const PufitGrid &grid, double &mean, double &variance);
+
+            /**
+             * @brief Calculate the mean and variance of unmasked towers
+             * @param grid The grid to calculate
+             * @param type The int mask of subgrid types
+             * @param[out] mean The calculated mean
+             * @param[out] variance The calculated variance
+             */
+            template <std::size_t N>
+            void unmaskedMeanAndVariance(
+                const PufitMultiGrid<N> &grid, int type, double &mean, double &variance);
+
+            /// Select the grid with the highest masked sumEt
+            GridDisplacement selectGrid(const PufitGridSet &grids);
+
+            /// Select the grid with the highest masked sumEt
+            template <typename Grid>
+            GridDisplacement selectGrid(
+                const PufitMultiGridSet<Grid> &grids,
+                std::size_t type);
+
+            /**
+             * @brief Perform the pile-up fit
+             * @param pileupSum 2-vector sum over the pileup contributions
+             * @param pileupCovariance Covariance matrix of the pileup sum
+             * @param towerExpectations Expected values for each tower based on the pileup
+             * energy density
+             * @param towerVariances Expected values for the variances of each tower based
+             * on variances of background towers
+             * @param correctionDirections The directions of the objects to correct
+             * @param The relative importance of the constraints.
+             * @return A vector containing the magnitudes of the corrections
+             *
+             * This variant provides individual estimates of the expected energy in each
+             * tower and the corresponding variances - these are for cases where the areas
+             * of each object are different.
+             */
+            Eigen::VectorXd pufit(
+                const Eigen::Vector2d &pileupSum,
+                const Eigen::Matrix2d &pileupCovariance,
+                const Eigen::VectorXd &towerExpectations,
+                const Eigen::VectorXd &towerVariances,
+                const Eigen::VectorXd &correctionDirections,
+                double constraintImportance = 1);
+
+            /**
+             * @brief Perform the pile-up fit
+             * @param pileupSum 2-vector sum over the pileup contributions
+             * @param pileupCovariance Covariance matrix of the pileup sum
+             * @param towerExpectations Expected values for each tower based on the pileup
+             * energy density
+             * @param towerVariances Expected values for the variances of each tower based
+             * on variances of background towers
+             * @param cosSin The cosines and sines of the objects to correct
+             * @param The relative importance of the constraints.
+             * @return A vector containing the magnitudes of the corrections
+             *
+             * This variant directly provides the cos-sin matrix for the correction,
+             * potentially allowing many cos/sin calls to be avoided.
+             * This variant provides individual estimates of the expected energy in each
+             * tower and the corresponding variances - these are for cases where the areas
+             * of each object are different.
+             */
+            Eigen::VectorXd pufit(
+                const Eigen::Vector2d &pileupSum,
+                const Eigen::Matrix2d &pileupCovariance,
+                const Eigen::VectorXd &towerExpectations,
+                const Eigen::VectorXd &towerVariances,
+                const Eigen::Matrix<double, 2, Eigen::Dynamic> &cosSin,
+                double constraintImportance = 1);
+
+            /**
+             * @brief Perform the pile-up fit
+             * @param pileupSum 2-vector sum over the pileup contributions
+             * @param pileupCovariance Covariance matrix of the pileup sum
+             * @param towerExpectations Expected values for each tower based on the pileup
+             * energy density
+             * @param towerVariances Expected values for the variances of each tower based
+             * on variances of background towers
+             * @param toCorrect Kinematics of the objects to correct
+             * @param The relative importance of the constraints.
+             * @return The corrections to apply
+             *
+             * Note that this version still returns the corrections that should be
+             * applied, not the corrected objects. The corrections are returned using the
+             * negatives of the value calculated by the fit - i.e. you should add the
+             * returned objects and let them act as negative energy objects (unless the
+             * sign of the correction was negative in which case it will act as a positive
+             * energy object).
+             * This variant provides individual estimates of the expected energy in each
+             * tower and the corresponding variances - these are for cases where the areas
+             * of each object are different.
+             */
+            std::vector<SignedKinematics> pufit(
+                const Eigen::Vector2d &pileupSum,
+                const Eigen::Matrix2d &pileupCovariance,
+                const std::vector<double> &towerExpectations,
+                const std::vector<double> &towerVariances,
+                const std::vector<SignedKinematics> &toCorrect,
+                double constraintImportance = 1);
+
+            /**
+             * @brief Perform the pile-up fit
+             * @param pileupSum 2-vector sum over the pileup contributions
+             * @param pileupCovariance Covariance matrix of the pileup sum
+             * @param towerMean Estimate of the mean pileup tower transverse energy
+             * @param towerVariance Estimate of the variance of pileup tower transverse
+             * energies
+             * @param correctionDirections The directions of the objects to correct
+             * @param The relative importance of the constraints.
+             * @return A vector containing the magnitudes of the corrections
+             */
+            Eigen::VectorXd pufit(
+                const Eigen::Vector2d &pileupSum,
+                const Eigen::Matrix2d &pileupCovariance,
+                double towerMean,
+                double towerVariance,
+                const Eigen::VectorXd &correctionDirections,
+                double constraintImportance = 1);
+
+            /**
+             * @brief Perform the pile-up fit
+             * @param pileupSum 2-vector sum over the pileup contributions
+             * @param pileupCovariance Covariance matrix of the pileup sum
+             * @param towerMean Estimate of the mean pileup tower transverse energy
+             * @param towerVariance Estimate of the variance of pileup tower transverse
+             * energies
+             * @param cosSin The cosines and sines of the objects to correct
+             * @param The relative importance of the constraints.
+             * @return A vector containing the magnitudes of the corrections
+             *
+             * This variant directly provides the cos-sin matrix for the correction,
+             * potentially allowing many cos/sin calls to be avoided.
+             */
+            Eigen::VectorXd pufit(
+                const Eigen::Vector2d &pileupSum,
+                const Eigen::Matrix2d &pileupCovariance,
+                double towerMean,
+                double towerVariance,
+                const Eigen::Matrix<double, 2, Eigen::Dynamic> &cosSin,
+                double constraintImportance = 1);
+
+            /**
+             * @brief Perform the pile-up fit
+             * @param pileupSum 2-vector sum over the pileup contributions
+             * @param pileupCovariance Covariance matrix of the pileup sum
+             * @param towerMean Estimate of the mean pileup tower transverse energy
+             * @param towerVariance Estimate of the variance of pileup tower transverse
+             * energies
+             * @param toCorrect Kinematics of the objects to correct
+             * @param The relative importance of the constraints.
+             * @return The corrections to apply
+             *
+             * Note that this version still returns the corrections that should be
+             * applied, not the corrected objects. The corrections are returned using the
+             * negatives of the value calculated by the fit - i.e. you should add the
+             * returned objects and let them act as negative energy objects (unless the
+             * sign of the correction was negative in which case it will act as a positive
+             * energy object).
+             */
+            std::vector<SignedKinematics> pufit(
+                const Eigen::Vector2d &pileupSum,
+                const Eigen::Matrix2d &pileupCovariance,
+                double towerMean,
+                double towerVariance,
+                const std::vector<SignedKinematics> &toCorrect,
+                double constraintImportance = 1);
+        } // namespace PufitUtils
+    }     // namespace MET
+} // namespace HLT
+
+#include "TrigEFMissingET/PufitUtils.icc"
 
 #endif //> !TRIGEFMISSINGET_PUFITUTILS_H
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitUtils.icc b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitUtils.icc
new file mode 100644
index 000000000000..37abecb9b638
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/TrigEFMissingET/PufitUtils.icc
@@ -0,0 +1,74 @@
+/**
+   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+namespace HLT
+{
+  namespace MET
+  {
+    namespace PufitUtils
+    {
+      template <std::size_t N>
+      void trimmedMeanAndVariance(
+          const PufitMultiGrid<N> &grid,
+          std::size_t type,
+          double trimFraction,
+          double &mean,
+          double &variance)
+      {
+        // Construct a list of sorted tower energies
+        std::vector<double> sorted;
+        sorted.reserve(grid.nTowers());
+        for (const typename PufitMultiGrid<N>::Tower &tower : grid)
+        {
+          double sumEt = tower.sumEt(type);
+          sorted.insert(
+              std::lower_bound(sorted.begin(), sorted.end(), sumEt),
+              sumEt);
+        }
+        trimmedMeanAndVariance(sorted, trimFraction, mean, variance);
+      }
+
+      template <std::size_t N>
+      void unmaskedMeanAndVariance(
+          const PufitMultiGrid<N> &grid, std::size_t type, double &mean, double &variance)
+      {
+        double sum = 0;
+        double squaredSum = 0;
+        std::size_t n = 0;
+        for (const typename PufitMultiGrid<N>::Tower &tower : grid)
+        {
+          if (tower.masked())
+            continue;
+          ++n;
+          double sumEt = tower.sumEt(type);
+          sum += sumEt;
+          squaredSum += sumEt * sumEt;
+        }
+        mean = sum / n;
+        variance = squaredSum / n - mean * mean;
+      }
+
+      template <typename Grid>
+      GridDisplacement selectGrid(
+          const PufitMultiGridSet<Grid> &grids, std::size_t type)
+      {
+        GridDisplacement maximum = NoDisplacement;
+        double maxSum = 0;
+        for (std::size_t d = 0; d < 4; ++d)
+        {
+          double sum = 0;
+          for (const typename Grid::Tower &tower : grids[GridDisplacement(d)])
+            if (tower.masked())
+              sum += tower.sumEt(type);
+          if (sum > maxSum)
+          {
+            maximum = GridDisplacement(d);
+            maxSum = sum;
+          }
+        }
+        return maximum;
+      }
+    } // namespace PufitUtils
+  }   // namespace MET
+} // namespace HLT
\ No newline at end of file
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/python/PUClassification.py b/Trigger/TrigAlgorithms/TrigEFMissingET/python/PUClassification.py
new file mode 100644
index 000000000000..ec4c3e8d5f9d
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/python/PUClassification.py
@@ -0,0 +1,17 @@
+#
+#  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+#
+
+""" Python version of the PUClassification enum """
+
+# Neutral/forward components without tracks
+NeutralForward = 1 << 0
+# Tracks not associated to the PV
+ChargedPU = 1 << 1
+# Tracks associated to the PV
+ChargedHS = 1 << 2
+
+NFcPU = NeutralForward | ChargedPU
+NFcHS = NeutralForward | ChargedHS
+Charged = ChargedHS | ChargedPU
+All = NeutralForward | ChargedPU | ChargedHS
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/ApproximateTrackToLayerTool.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/ApproximateTrackToLayerTool.cxx
new file mode 100644
index 000000000000..874daff9fa61
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/ApproximateTrackToLayerTool.cxx
@@ -0,0 +1,27 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "TrigEFMissingET/ApproximateTrackToLayerTool.h"
+
+ApproximateTrackToLayerTool::ApproximateTrackToLayerTool(const std::string &name) : asg::AsgTool(name)
+{
+}
+
+IExtendTrackToLayerTool::TrackExtension ApproximateTrackToLayerTool::extendTrack(
+    const EventContext &,
+    const xAOD::TrackParticle &track) const
+{
+  // The approximation doesn't affect track etas at all
+  double eta = track.eta();
+  double sinDPhi = (m_trackExtrapolationLinear.value() +
+                    m_trackExtrapolationQuadratic.value() * eta*eta +
+                    m_trackExtrapolationQuartic.value() * eta*eta*eta*eta) /
+                   (track.pt() * track.charge());
+
+  if (std::abs(sinDPhi) > 1.0)
+    // The track was extrapolated outside of acceptance
+    return TrackExtension{};
+  else
+    return TrackExtension{track.eta(), track.phi() - std::asin(sinDPhi)};
+}
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/CVFAlg.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/CVFAlg.cxx
new file mode 100644
index 000000000000..b4de435dc06c
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/CVFAlg.cxx
@@ -0,0 +1,173 @@
+/*
+   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+#include "CVFAlg.h"
+#include "StoreGate/DecorKeyHelpers.h"
+#include "StoreGate/ReadHandle.h"
+#include "StoreGate/WriteDecorHandle.h"
+#include "AthContainers/ConstDataVector.h"
+
+// NB: We should try and get something like the faster binned phi approach used
+// in the PFlow construction. Of course, we might *really* want to just find a
+// way of sharing that information? (for the track <-> cluster association)
+namespace
+{
+  bool fastDR(
+      double eta,
+      double phi,
+      const IExtendTrackToLayerTool::TrackExtension &track,
+      double dr,
+      double dr2)
+  {
+    double dEta = std::abs(eta - track.eta);
+    if (dEta > dr)
+      return false;
+    double dPhi = std::abs(phi - track.phi);
+    if (dPhi > dr)
+      return false;
+    return (dEta * dEta + dPhi * dPhi) < dr2;
+  }
+} // namespace
+
+namespace HLT
+{
+  namespace MET
+  {
+    CVFAlg::CVFAlg(const std::string &name, ISvcLocator *pSvcLocator)
+        : AthReentrantAlgorithm(name, pSvcLocator)
+    {
+    }
+
+    StatusCode CVFAlg::initialize()
+    {
+      ATH_MSG_DEBUG("Initializing " << name());
+      CHECK(m_extensionTool.retrieve());
+      if (!m_trackSelTool.empty())
+      {
+        CHECK(m_trackSelTool.retrieve());
+        m_useTrackSelTool = true;
+      }
+      CHECK(m_tvaTool.retrieve());
+      CHECK(m_inputClusterKey.initialize());
+      CHECK(m_inputTrackKey.initialize());
+      CHECK(m_inputVertexKey.initialize());
+      if (m_outputCVFKey.key().find(".") == std::string::npos)
+        m_outputCVFKey = m_inputClusterKey.key() + "." + m_outputCVFKey.key();
+      else if (SG::contKeyFromKey(m_outputCVFKey.key()) != m_inputClusterKey.key())
+      {
+        ATH_MSG_ERROR("CVF key does not match input cluster key!");
+        return StatusCode::FAILURE;
+      }
+      CHECK(m_outputCVFKey.initialize());
+      return StatusCode::SUCCESS;
+    }
+
+    StatusCode CVFAlg::execute(const EventContext &ctx) const
+    {
+      using TrackExtension = IExtendTrackToLayerTool::TrackExtension;
+
+      auto clusters = SG::makeHandle(m_inputClusterKey, ctx);
+      if (!clusters.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_inputClusterKey);
+        return StatusCode::FAILURE;
+      }
+      auto tracks = SG::makeHandle(m_inputTrackKey, ctx);
+      if (!tracks.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_inputTrackKey);
+        return StatusCode::FAILURE;
+      }
+      auto vertexHandle = SG::makeHandle(m_inputVertexKey, ctx);
+      if (!vertexHandle.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_inputVertexKey);
+        return StatusCode::FAILURE;
+      }
+      auto decCVF = SG::makeHandle<float>(m_outputCVFKey, ctx);
+
+      // Find the primary vertex
+      const xAOD::Vertex *priVtx = nullptr;
+      std::vector<const xAOD::Vertex *> vertices;
+      vertices.reserve(vertexHandle->size());
+      for (const xAOD::Vertex *vtx : *vertexHandle)
+      {
+        if (!priVtx && vtx->vertexType() == xAOD::VxType::PriVtx)
+          priVtx = vtx;
+        vertices.push_back(vtx);
+      }
+      if (!priVtx)
+      {
+        // If we fail to find a primary vertex then we give every cluster a
+        // CVF value of -1
+        for (const xAOD::CaloCluster *iclus : *clusters)
+          decCVF(*iclus) = -1;
+        return StatusCode::SUCCESS;
+      }
+
+      const xAOD::TrackParticleContainer *filteredTracks = nullptr;
+      ConstDataVector<xAOD::TrackParticleContainer> viewTracks(SG::VIEW_ELEMENTS);
+      if (m_useTrackSelTool)
+      {
+        for (const xAOD::TrackParticle *track : *tracks)
+          if (m_trackSelTool->accept(*track))
+            viewTracks.push_back(track);
+        filteredTracks = viewTracks.asDataVector();
+      }
+      else
+        filteredTracks = tracks.ptr();
+
+      // To speed things up later, check the maximum positive and negative range of
+      // the track etas. This means we can avoid too many loops and dR comparisons
+      // later on.
+      double minTrkEta = std::numeric_limits<double>::infinity();
+      double maxTrkEta = -std::numeric_limits<double>::infinity();
+      ATH_MSG_DEBUG("Building extrapolation");
+      std::vector<std::tuple<bool, const xAOD::TrackParticle *, TrackExtension>> extensions;
+      extensions.reserve(filteredTracks->size());
+      for (const xAOD::TrackParticle *track : *filteredTracks)
+      {
+        TrackExtension ext = m_extensionTool->extendTrack(ctx, *track);
+        if (!ext.isValid())
+          continue;
+        minTrkEta = std::min(minTrkEta, ext.eta);
+        maxTrkEta = std::max(maxTrkEta, ext.eta);
+        bool isPV;
+        if (m_useCompatible)
+          isPV = m_tvaTool->isCompatible(*track, *priVtx);
+        else
+          isPV = m_tvaTool->getUniqueMatchVertex(*track, vertices) == priVtx;
+        extensions.push_back(std::make_tuple(isPV, track, ext));
+      }
+
+      double DR2 = m_clusterMatchDR.value() * m_clusterMatchDR.value();
+      for (const xAOD::CaloCluster *iclus : *clusters)
+      {
+        double eta = iclus->eta();
+        if (iclus->e() ||
+            eta - m_clusterMatchDR > maxTrkEta ||
+            eta + m_clusterMatchDR < minTrkEta)
+        {
+          decCVF(*iclus) = -1;
+          continue;
+        }
+
+        double phi = iclus->phi();
+        double totalSum = 0;
+        double PVSum = 0;
+        for (const std::tuple<bool, const xAOD::TrackParticle *, TrackExtension> &ext : extensions)
+        {
+          // match the cluster to the track
+          if (!fastDR(eta, phi, std::get<2>(ext), m_clusterMatchDR, DR2))
+            continue;
+          totalSum += std::get<1>(ext)->pt();
+          if (std::get<0>(ext))
+            PVSum += std::get<1>(ext)->pt();
+        }
+        decCVF(*iclus) = totalSum == 0 ? -1 : PVSum / totalSum;
+      }
+      return StatusCode::SUCCESS;
+    }
+  } // namespace MET
+} // namespace HLT
\ No newline at end of file
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/CVFAlg.h b/Trigger/TrigAlgorithms/TrigEFMissingET/src/CVFAlg.h
new file mode 100644
index 000000000000..b9d03ba922f5
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/CVFAlg.h
@@ -0,0 +1,71 @@
+/*
+   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+/******************************************************************************
+ * @package Trigger/TrigAlgorithms/TrigEFMissingET
+ * @class   CVFAlg
+ * @brief   Algorithm class to calculate CVF
+ * @author  Jon Burr
+ *****************************************************************************/
+
+#ifndef TRIGEFMISSINGET_CVFALG_H
+#define TRIGEFMISSINGET_CVFALG_H
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "AsgTools/ToolHandle.h"
+#include "TrigEFMissingET/IExtendTrackToLayerTool.h"
+#include "InDetTrackSelectionTool/IInDetTrackSelectionTool.h"
+#include "TrackVertexAssociationTool/ITrackVertexAssociationTool.h"
+#include "StoreGate/ReadHandleKey.h"
+#include "StoreGate/WriteDecorHandleKey.h"
+#include "xAODCaloEvent/CaloClusterContainer.h"
+#include "xAODTracking/TrackParticleContainer.h"
+#include "xAODTracking/VertexContainer.h"
+#include "Gaudi/Property.h"
+
+namespace HLT
+{
+  namespace MET
+  {
+    class CVFAlg : public ::AthReentrantAlgorithm
+    {
+    public:
+      /// Constructor
+      CVFAlg(const std::string &name, ISvcLocator *pSvcLocator);
+
+      /// Initialise the algorithm
+      virtual StatusCode initialize() override;
+
+      /// Run the algorithm
+      virtual StatusCode execute(const EventContext &ctx) const override;
+
+    private:
+      ToolHandle<IExtendTrackToLayerTool> m_extensionTool{
+          this, "ExtensionTool", "", "The extension tool"};
+      ToolHandle<InDet::IInDetTrackSelectionTool> m_trackSelTool{
+          this, "TrackSelectionTool", "", "An optional track selection tool to filter tracks"};
+      ToolHandle<CP::ITrackVertexAssociationTool> m_tvaTool{
+          this, "TVATool", "", "Track -> vertex associationTool"};
+      SG::ReadHandleKey<xAOD::CaloClusterContainer> m_inputClusterKey{
+          this, "InputClusterKey", "", "Input cluster container"};
+      SG::ReadHandleKey<xAOD::TrackParticleContainer> m_inputTrackKey{
+          this, "InputTrackKey", "", "Input track container"};
+      SG::ReadHandleKey<xAOD::VertexContainer> m_inputVertexKey{
+          this, "InputVertexKey", "", "Input vertex container"};
+      SG::WriteDecorHandleKey<xAOD::CaloClusterContainer> m_outputCVFKey{
+          this, "OutputCVFKey", "", "Output CVF name"};
+      Gaudi::Property<double> m_clusterMatchDR{
+          this, "ClusterMatchDR", 0.1, "Max DR to match clusters to extrapolated tracks"};
+      Gaudi::Property<bool> m_useCompatible{
+          this, "UseCompatible", true,
+          "Use the track -> vertex compatibility over unique matching."};
+
+      // Internal
+      /// Is a track selection tool being used?
+      bool m_useTrackSelTool{false};
+    }; //> end class CVFAlg
+  }    // namespace MET
+} // namespace HLT
+
+#endif //> !TRIGEFMISSINGET_CVFALG_H
\ No newline at end of file
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/CVFPrepAlg.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/CVFPrepAlg.cxx
new file mode 100644
index 000000000000..d102e43783a3
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/CVFPrepAlg.cxx
@@ -0,0 +1,71 @@
+/*
+   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+#include "CVFPrepAlg.h"
+#include "StoreGate/DecorKeyHelpers.h"
+#include "StoreGate/ReadHandle.h"
+#include "StoreGate/ReadDecorHandle.h"
+#include "StoreGate/WriteDecorHandle.h"
+#include "TrigEFMissingET/PUClassification.h"
+
+namespace HLT
+{
+  namespace MET
+  {
+    CVFPrepAlg::CVFPrepAlg(const std::string &name, ISvcLocator *pSvcLocator)
+        : AthReentrantAlgorithm(name, pSvcLocator)
+    {
+    }
+
+    StatusCode CVFPrepAlg::initialize()
+    {
+      CHECK(m_inputClusterKey.initialize());
+
+      if (m_inputCVFKey.key().find(".") == std::string::npos)
+        m_inputCVFKey = m_inputClusterKey.key() + "." + m_inputCVFKey.key();
+      else if (SG::contKeyFromKey(m_inputCVFKey.key()) != m_inputClusterKey.key())
+      {
+        ATH_MSG_ERROR("CVF key does not match input cluster container key!");
+        return StatusCode::FAILURE;
+      }
+      CHECK(m_inputCVFKey.initialize());
+
+      if (m_outputCategoryKey.key().find(".") == std::string::npos)
+        m_outputCategoryKey = m_inputClusterKey.key() + "." + m_outputCategoryKey.key();
+      else if (SG::contKeyFromKey(m_outputCategoryKey.key()) != m_inputClusterKey.key())
+      {
+        ATH_MSG_ERROR("Output category key does not match input cluster container key!");
+        return StatusCode::FAILURE;
+      }
+      CHECK(m_outputCategoryKey.initialize());
+
+      return StatusCode::SUCCESS;
+    }
+
+    StatusCode CVFPrepAlg::execute(const EventContext &ctx) const
+    {
+      auto clusters = SG::makeHandle(m_inputClusterKey, ctx);
+      if (!clusters.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_inputClusterKey);
+        return StatusCode::FAILURE;
+      }
+      auto accCVF = SG::makeHandle<float>(m_inputCVFKey, ctx);
+      auto decCategory = SG::makeHandle<int>(m_outputCategoryKey, ctx);
+
+      for (const xAOD::CaloCluster *iclus : *clusters)
+      {
+        float CVF = accCVF(*iclus);
+        if (CVF < 0)
+          decCategory(*iclus) = PUClassification::NeutralForward;
+        else if (CVF < m_cvfThreshold)
+          decCategory(*iclus) = PUClassification::ChargedPU;
+        else
+          decCategory(*iclus) = PUClassification::ChargedHS;
+      }
+      return StatusCode::SUCCESS;
+    }
+
+  } // namespace MET
+} // namespace HLT
\ No newline at end of file
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/CVFPrepAlg.h b/Trigger/TrigAlgorithms/TrigEFMissingET/src/CVFPrepAlg.h
new file mode 100644
index 000000000000..1ed8a7332529
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/CVFPrepAlg.h
@@ -0,0 +1,62 @@
+/*
+   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+/******************************************************************************
+ * @package Trigger/TrigAlgorithms/TrigEFMissingET
+ * @class   CVFPrepAlg
+ *
+ * @brief Algorithm to prepare clusters for the PUSplit pufit alg
+ * @author Jon Burr
+ *****************************************************************************/
+
+#ifndef TRIGEFMISSINGET_CVFPREPALG_H
+#define TRIGEFMISSINGET_CVFPREPALG_H
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "StoreGate/ReadHandleKey.h"
+#include "StoreGate/ReadDecorHandleKey.h"
+#include "StoreGate/WriteDecorHandleKey.h"
+#include "xAODCaloEvent/CaloClusterContainer.h"
+#include "Gaudi/Property.h"
+
+namespace HLT
+{
+  namespace MET
+  {
+    /**************************************************************************
+     * @class PFOPrepAlg
+     * 
+     * Class to prepare the clusters for the PUSplit pufit alg.
+     * 
+     * This just uses the CVF value per cluster to add the necessary
+     * classification decoration
+     **************************************************************************/
+    class CVFPrepAlg : public AthReentrantAlgorithm
+    {
+    public:
+      /// Constructor
+      CVFPrepAlg(const std::string &name, ISvcLocator *pSvcLocator);
+
+      virtual StatusCode initialize() override;
+
+      virtual StatusCode execute(const EventContext &ctx) const override;
+
+    private:
+      /// The input container
+      SG::ReadHandleKey<xAOD::CaloClusterContainer> m_inputClusterKey{
+          this, "InputClusterKey", "", "Input cluster container"};
+      /// The input CVF key
+      SG::ReadDecorHandleKey<xAOD::CaloClusterContainer> m_inputCVFKey{
+          this, "InputCVFKey", "", "Input CVF key"};
+      /// The output classifcation
+      SG::WriteDecorHandleKey<xAOD::CaloClusterContainer> m_outputCategoryKey{
+          this, "OutputCategoryKey", "", "Output category name"};
+      /// The CVF threshold
+      Gaudi::Property<float> m_cvfThreshold{
+          this, "CVFThreshold", 0.5, "The CVF threshold to apply"};
+    }; //> end class CVFPrepAlg
+  }    // namespace MET
+} // namespace HLT
+
+#endif //> !TRIGEFMISSINGET_CVFPREPALG_H
\ No newline at end of file
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/CellFex.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/CellFex.cxx
index 5d5295d690c2..f0881be1ad7b 100644
--- a/Trigger/TrigAlgorithms/TrigEFMissingET/src/CellFex.cxx
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/CellFex.cxx
@@ -58,8 +58,18 @@ namespace HLT { namespace MET {
   {
     // Retrieve the inputs
     auto cells = SG::makeHandle(m_cellsKey, context);
+      if (!cells.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_cellsKey);
+        return StatusCode::FAILURE;
+      }
     // NB - there's no makeHandle overload for ReadCondHandle
     SG::ReadCondHandle noiseCDO(m_noiseCDOKey, context);
+      if (!noiseCDO.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_noiseCDOKey);
+        return StatusCode::FAILURE;
+      }
 
     // Prepare the individual components
     std::array<METComponent, N_SAMPLINGS> sums;
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/ExtendTrackToLayerTool.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/ExtendTrackToLayerTool.cxx
new file mode 100644
index 000000000000..e1d2d4af07f6
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/ExtendTrackToLayerTool.cxx
@@ -0,0 +1,66 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "TrigEFMissingET/ExtendTrackToLayerTool.h"
+#include "TrkCaloExtension/CaloExtensionHelpers.h"
+
+ExtendTrackToLayerTool::ExtendTrackToLayerTool(const std::string &name) : asg::AsgTool(name)
+{
+}
+
+StatusCode ExtendTrackToLayerTool::initialize()
+{
+  ATH_MSG_DEBUG("Initialize " << name());
+  CHECK(m_extensionTool.retrieve());
+
+  if (m_caloLayerNames.value().size() == 0)
+  {
+    ATH_MSG_ERROR("No layers requested, tool would not do anything!");
+    return StatusCode::FAILURE;
+  }
+
+  std::map<std::string, CaloSampling::CaloSample> samplings;
+  for (unsigned int ii = 0; ii < CaloSampling::getNumberOfSamplings(); ++ii)
+    samplings[CaloSampling::getSamplingName(ii)] = static_cast<CaloSampling::CaloSample>(ii);
+  m_caloLayers.reserve(m_caloLayerNames.value().size());
+  for (const std::string &name : m_caloLayerNames.value())
+  {
+    auto itr = samplings.find(name);
+    if (itr == samplings.end())
+    {
+      ATH_MSG_ERROR("Unknown sampling '" << name << "' requested!");
+      return StatusCode::FAILURE;
+    }
+  }
+  return StatusCode::SUCCESS;
+}
+
+IExtendTrackToLayerTool::TrackExtension ExtendTrackToLayerTool::extendTrack(
+    const EventContext &ctx,
+    const xAOD::TrackParticle &track) const
+{
+  std::unique_ptr<Trk::CaloExtension> extension = m_extensionTool->caloExtension(
+      ctx, track);
+  if (!extension)
+    // This means that track wasn't extrapolated into the calorimeter
+    return TrackExtension{};
+  // Get the eta-phi entry positions of the intersections of this track with the
+  // calorimeter
+  CaloExtensionHelpers::EtaPhiHashLookupVector result;
+  CaloExtensionHelpers::entryEtaPhiHashLookupVector(*extension, result);
+  // Go through the layers we requested in order
+  bool isCrossed;
+  double eta;
+  double phi;
+  for (CaloSampling::CaloSample layer : m_caloLayers)
+  {
+    std::tie(isCrossed, eta, phi) = result[layer];
+    // Return the first that is crossed
+    if (isCrossed)
+      return TrackExtension{eta, phi};
+  }
+  // If we got here, it crossed no layer that we're interested in, return an
+  // invalid extension
+  return TrackExtension{};
+}
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/MHTFex.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/MHTFex.cxx
index a46bc079d6ff..2465aeb14212 100644
--- a/Trigger/TrigAlgorithms/TrigEFMissingET/src/MHTFex.cxx
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/MHTFex.cxx
@@ -32,6 +32,11 @@ namespace HLT { namespace MET {
   {
     // Retrieve the inputs
     auto jets = SG::makeHandle(m_jetKey, context);
+    if (!jets.isValid())
+    {
+      ATH_MSG_ERROR("Failed to retrieve " << m_jetKey);
+      return StatusCode::FAILURE;
+    }
 
     // Prepare the output values
     std::array<METComponent, 4> mhtSums;
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/MHTPufitFex.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/MHTPufitFex.cxx
new file mode 100644
index 000000000000..c16c94f88f07
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/MHTPufitFex.cxx
@@ -0,0 +1,204 @@
+/*
+   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+/******************************************************************************
+ * @package Trigger/TrigAlgorithms/TrigEFMissingET
+ * @file MHTPufitFex.cxx
+ *
+ * Implementation of mhtpufit fex class
+ * @author Jon Burr
+ *****************************************************************************/
+
+#include "MHTPufitFex.h"
+#include "StoreGate/ReadHandle.h"
+#include "StoreGate/ReadDecorHandle.h"
+#include "StoreGate/DecorKeyHelpers.h"
+#include "TrigEFMissingET/METComponent.h"
+#include "TrigEFMissingET/PufitUtils.h"
+#include "TrigEFMissingET/PufitGrid.h"
+#include "xAODBase/IParticleHelpers.h"
+
+namespace
+{
+  const static SG::AuxElement::ConstAccessor<float> accArea("ActiveArea4vec_pt");
+  const static SG::AuxElement::ConstAccessor<float> accDetectorEta("DetectorEta");
+} // namespace
+
+namespace HLT
+{
+  namespace MET
+  {
+    double MHTPufitFex::getSigma(const SignedKinematics &kin) const
+    {
+      return m_caloNoise * m_caloNoise + kin.absPt() * m_caloStoch * m_caloStoch;
+    }
+
+    MHTPufitFex::MHTPufitFex(const std::string &name, ISvcLocator *pSvcLocator)
+        : FexBase(name, pSvcLocator)
+    {
+    }
+
+    StatusCode MHTPufitFex::initialize()
+    {
+      CHECK(m_inputJetsKey.initialize());
+      if (m_inputJvtKey.key().find(".") == std::string::npos)
+        m_inputJvtKey = m_inputJetsKey.key() + "." + m_inputJvtKey.key();
+      else if (SG::contKeyFromKey(m_inputJvtKey.key()) != m_inputJvtKey.key())
+      {
+        ATH_MSG_ERROR("Input JVT key does not match the input jet key!");
+        return StatusCode::FAILURE;
+      }
+      CHECK(m_inputJvtKey.initialize());
+
+      CHECK(m_inputKey.initialize());
+
+      CHECK(m_rhoKey.initialize(m_jetCalibIncludesAreaSub));
+
+      return initializeBase({"UncorrJets"});
+    }
+
+    StatusCode MHTPufitFex::fillMET(
+        xAOD::TrigMissingET &met,
+        const EventContext &context,
+        MonGroupBuilder &) const
+    {
+      using namespace PufitUtils;
+      // Retrieve the input
+      auto inputJets = SG::makeHandle(m_inputJetsKey, context);
+      if (!inputJets.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_inputJetsKey);
+        return StatusCode::FAILURE;
+      }
+      auto jvtAcc = SG::makeHandle<float>(m_inputJvtKey, context);
+      auto inputs = SG::makeHandle(m_inputKey, context);
+      if (!inputs.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_inputKey);
+        return StatusCode::FAILURE;
+      }
+
+      double rho = 0;
+      if (m_jetCalibIncludesAreaSub)
+      {
+        auto rhoHandle = SG::makeHandle(m_rhoKey, context);
+        if (!rhoHandle->getDensity(xAOD::EventShape::Density, rho))
+        {
+          ATH_MSG_ERROR("EventShape density not available!");
+          return StatusCode::FAILURE;
+        }
+      }
+
+      // Create the grid
+      PufitGrid grid(m_maxEta, m_nEtaBins, m_nPhiBins);
+      CovarianceSum pileupSum;
+      METComponent sum;
+      METComponent uncorrSum;
+
+      // Keep track of which clusters/PFOs we have decided are included in a HS jet
+      std::set<const xAOD::IParticle *> hardScatterInputs;
+      std::vector<SignedKinematics> hardScatterJets;
+      std::vector<float> jetAreas;
+      for (const xAOD::Jet *ijet : *inputJets)
+      {
+        float eta = m_useDetectorEta ? accDetectorEta(*ijet) : ijet->eta();
+        if (std::abs(eta) > m_centralEta)
+        {
+          if (ijet->pt() < m_forwardPt)
+            continue;
+        }
+        else
+        {
+          if (ijet->pt() < m_minPt)
+            continue;
+          else if (ijet->pt() < m_maxPt && jvtAcc(*ijet) < m_jvtCut)
+            continue;
+        }
+        // If we're here then the jet has passed the selection
+        hardScatterJets.push_back(*ijet);
+        jetAreas.push_back(accArea(*ijet));
+        for (const auto &link : ijet->constituentLinks())
+        {
+          if (link.isValid())
+          {
+            const xAOD::IParticle *constituent = *link;
+            if (const xAOD::IParticle *original = xAOD::getOriginalObject(*constituent))
+              hardScatterInputs.insert(original);
+            else
+              hardScatterInputs.insert(constituent);
+          }
+          else
+            ATH_MSG_WARNING("Invalid constituent link!");
+        }
+      } //> end loop over jets
+
+      if (hardScatterJets.size() > 0)
+      {
+        for (const xAOD::IParticle *ipart : *inputs)
+        {
+          const xAOD::IParticle *original = xAOD::getOriginalObject(*ipart);
+          if (!original)
+            original = ipart;
+          bool isHS = hardScatterInputs.count(original);
+          SignedKinematics kin = SignedKinematics::fromEnergyEtaPhi(
+              ipart->e(), ipart->eta(), ipart->phi());
+          if (!isHS)
+            pileupSum.add(kin, getSigma(kin));
+
+          // Figure out which tower this particle belongs in
+          bool outOfRange = false;
+          std::size_t idx = grid.getIndex(ipart->eta(), ipart->phi(), outOfRange);
+          if (!outOfRange)
+          {
+            PufitGrid::Tower &tower = grid[idx];
+            tower += kin;
+            if (isHS)
+              // mask a tower that has *any* clusters from HS jets in it
+              tower.mask();
+          }
+        }
+
+        // Calculate the expected PU contributions to each jet
+        double towerMean = 0;
+        double towerVariance = 0;
+        unmaskedMeanAndVariance(grid, towerMean, towerVariance);
+        float towerArea = grid.etaWidth() * grid.phiWidth();
+        std::vector<double> means;
+        std::vector<double> variances;
+        means.reserve(jetAreas.size());
+        variances.reserve(jetAreas.size());
+        for (float area : jetAreas)
+        {
+          means.push_back(towerMean * area / towerArea);
+          variances.push_back(towerVariance * area / towerArea);
+        }
+
+        std::vector<SignedKinematics> corrections = pufit(
+            pileupSum.sum,
+            pileupSum.covariance,
+            means,
+            variances,
+            hardScatterJets,
+            m_constraintImportance);
+
+        for (std::size_t ii = 0; ii < hardScatterJets.size(); ++ii)
+        {
+          SignedKinematics kin = hardScatterJets.at(ii);
+          uncorrSum += kin;
+          if (m_jetCalibIncludesAreaSub)
+            // Add back the energy that was subtracted by the rho sub part of the
+            // calibration otherwise pufit will overcount this effect
+            kin += SignedKinematics::fromEnergyEtaPhi(
+                jetAreas.at(ii) * rho, kin.eta(), kin.phi());
+          sum += kin;
+          sum += corrections.at(ii);
+        }
+      }
+
+      uncorrSum.fillMETComponent(0, met);
+      sum.fillMET(met);
+
+      return StatusCode::SUCCESS;
+    }
+  } // namespace MET
+} // namespace HLT
\ No newline at end of file
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/MHTPufitFex.h b/Trigger/TrigAlgorithms/TrigEFMissingET/src/MHTPufitFex.h
new file mode 100644
index 000000000000..ef918c7d88e0
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/MHTPufitFex.h
@@ -0,0 +1,128 @@
+/*
+   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+/******************************************************************************
+ * @package Trigger/TrigAlgorithms/TrigEFMissingET
+ * @class   MHTPufitFex
+ *
+ * @brief Fex class for the mhtpufit algorithm, using the pufit procedure to
+ * correct the MET produced by a sum over jets
+ * @author Jon Burr
+ *****************************************************************************/
+
+#ifndef TRIGEFMISSINGET_MHTPUFITFEX_H
+#define TRIGEFMISSINGET_MHTPUFITFEX_H
+
+#include "FexBase.h"
+#include "TrigEFMissingET/SignedKinematics.h"
+#include "xAODJet/JetContainer.h"
+#include "xAODBase/IParticleContainer.h"
+#include "xAODEventShape/EventShape.h"
+#include "StoreGate/ReadDecorHandleKey.h"
+#include "Gaudi/Property.h"
+#include "GaudiKernel/SystemOfUnits.h"
+
+namespace HLT
+{
+  namespace MET
+  {
+    /**************************************************************************
+     * @class MHTPufitFex
+     * 
+     * Class to calculate MET from a sum over jets using a correction derived
+     * with the pufit technique
+     *************************************************************************/
+    class MHTPufitFex : public FexBase
+    {
+    public:
+      /// Constructor
+      MHTPufitFex(const std::string &name, ISvcLocator *pSvcLocator);
+
+      /// Initialize the fex
+      virtual StatusCode initialize() override;
+
+    private:
+      /************************************************************************
+       * Properties
+       ***********************************************************************/
+      /// Input objects
+      SG::ReadHandleKey<xAOD::JetContainer> m_inputJetsKey{
+          this, "InputJetsName", "", "The input jet container"};
+      /// The input JVT decoration
+      SG::ReadDecorHandleKey<xAOD::JetContainer> m_inputJvtKey{
+          this, "InputJvtName", "Jvt", "The input JVT name"};
+      /// The input clusters or PFOs
+      SG::ReadHandleKey<xAOD::IParticleContainer> m_inputKey{
+          this, "InputName", "", "The input clusters or PFOs"};
+      Gaudi::Property<bool> m_jetCalibIncludesAreaSub{
+          this, "JetCalibIncludesAreaSub", true,
+          "Whether the calibration applied to the jets includes area subtraction"};
+      SG::ReadHandleKey<xAOD::EventShape> m_rhoKey{
+          this, "JetEventShapeName", "",
+          "The name of the event shape container for the area correction"};
+      /// The sigma threshold
+      Gaudi::Property<float> m_nSigma{
+          this, "NSigma", 5, "Set the threshold at mean + NSigma*variance"};
+      /// The eta range of the grid
+      Gaudi::Property<float> m_maxEta{
+          this, "MaxEta", 5, "The maximum eta range"};
+      /// The number of bins in eta
+      Gaudi::Property<std::size_t> m_nEtaBins{
+          this, "NEtaBins", 14, "The number of eta bins"};
+      Gaudi::Property<std::size_t> m_nPhiBins{
+          this, "NPhiBins", 8, "The number of phi bins"};
+      /// The trimming fraction
+      Gaudi::Property<float> m_trimFraction{
+          this, "TrimFraction", 0.9,
+          "The fraction of bins to use when calculating the mean and variance"};
+      /// The coefficient of the noise term in the calo resolution estimate
+      Gaudi::Property<float> m_caloNoise{
+          this, "CaloNoise", 50,
+          "The coefficient of the noise term in the calorimeter resolution estimate [MeV]"};
+      /// The coefficient of the stochastic term in the calo resolution estimate
+      Gaudi::Property<float> m_caloStoch{
+          this, "CaloStochastic", 15.81,
+          "The coefficient of the stochastic term in the calorimeter resolution estimate [MeV^1/2]"};
+      /// The relative constraint importance
+      Gaudi::Property<float> m_constraintImportance{
+          this, "ConstraintImportance", 1,
+          "The relative importance of the two constraints to the fit"};
+      // Properties relating to jet selection
+      Gaudi::Property<float> m_jvtCut{
+          this, "JvtCut", 0.59, "The JVT selection in the central region"};
+      Gaudi::Property<float> m_minPt{
+          this, "MinPt", 20*Gaudi::Units::GeV, "The minimum pT (in the central region)"};
+      Gaudi::Property<float> m_maxPt{
+          this, "MaxPt", 120*Gaudi::Units::GeV,
+          "The maximum pT (in the central region), above which the JVT selection is not applied"};
+      Gaudi::Property<float> m_forwardPt{
+          this, "ForwardPt", 30*Gaudi::Units::GeV, "The minimum pt in the forward region"};
+      Gaudi::Property<float> m_centralEta{
+          this, "CentralEta", 2.4, "The boundary between the central and border regions"};
+      Gaudi::Property<bool> m_useDetectorEta{
+          this, "UseDetectorEta", true, "Whether to use the 'DectectorEta' value to select central/forward jets"};
+      /************************************************************************
+       * Internal functions
+       ***********************************************************************/
+      /**
+       * @brief Calculate and fill the output MET value
+       * @param met The object to fill
+       * @param context The event context
+       * @param monitors[out] Any extra monitors to fill
+       */
+      virtual StatusCode fillMET(
+          xAOD::TrigMissingET &met,
+          const EventContext &context,
+          MonGroupBuilder &monitors) const override;
+
+      /**
+       * @brief Calculate the estimate on the variance of a tower
+       * @param kin The kinematics of the tower
+       */
+      double getSigma(const SignedKinematics &kin) const;
+    }; //> end class MHTPufitFex
+  }    // namespace MET
+} // namespace HLT
+
+#endif //> !TRIGEFMISSINGET_MHTPUFITFEX_H
\ No newline at end of file
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/PFOPrepAlg.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PFOPrepAlg.cxx
new file mode 100644
index 000000000000..b70b6c13e700
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PFOPrepAlg.cxx
@@ -0,0 +1,137 @@
+/*
+   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+#include "PFOPrepAlg.h"
+#include "StoreGate/ReadHandle.h"
+#include "StoreGate/WriteHandle.h"
+#include "StoreGate/DecorKeyHelpers.h"
+#include "TrigEFMissingET/PUClassification.h"
+#include "AthContainers/ConstDataVector.h"
+
+/// Anonymous namespace
+namespace
+{
+  const static SG::AuxElement::ConstAccessor<char> accPVMatched("matchedToPV");
+}
+
+namespace HLT
+{
+  namespace MET
+  {
+    PFOPrepAlg::PFOPrepAlg(const std::string &name, ISvcLocator *pSvcLocator)
+        : AthReentrantAlgorithm(name, pSvcLocator)
+    {
+    }
+
+    StatusCode PFOPrepAlg::initialize()
+    {
+      m_manualTVA = !m_tvaTool.empty();
+      CHECK(m_inputNeutralKey.initialize());
+      CHECK(m_inputChargedKey.initialize());
+      CHECK(m_outputKey.initialize());
+      if (m_outputCategoryKey.key().find(".") == std::string::npos)
+        m_outputCategoryKey = m_outputKey.key() + "." + m_outputCategoryKey.key();
+      else if (SG::contKeyFromKey(m_outputCategoryKey.key()) != m_outputKey.key())
+      {
+        ATH_MSG_ERROR("Category key does not match output container key!");
+        return StatusCode::FAILURE;
+      }
+      CHECK(m_outputCategoryKey.initialize());
+      m_decCategory.emplace(SG::decorKeyFromKey(m_outputCategoryKey.key()));
+
+      if (m_manualTVA)
+      {
+        CHECK(m_tvaTool.retrieve());
+        if (!m_trackSelTool.empty())
+          CHECK(m_trackSelTool.retrieve());
+      }
+      CHECK(m_inputVertexKey.initialize(m_manualTVA));
+      return StatusCode::SUCCESS;
+    }
+
+    StatusCode PFOPrepAlg::execute(const EventContext &ctx) const
+    {
+      auto charged = SG::makeHandle(m_inputChargedKey, ctx);
+      if (!charged.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_inputChargedKey);
+        return StatusCode::FAILURE;
+      }
+      auto neutral = SG::makeHandle(m_inputNeutralKey, ctx);
+      if (!neutral.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_inputNeutralKey);
+        return StatusCode::FAILURE;
+      }
+      const xAOD::Vertex *priVtx = nullptr;
+      std::vector<const xAOD::Vertex *> vertices;
+      if (m_manualTVA)
+      {
+        auto vertexHandle = SG::makeHandle(m_inputVertexKey, ctx);
+        if (!vertexHandle.isValid())
+        {
+          ATH_MSG_ERROR("Failed to retrieve " << m_inputVertexKey);
+          return StatusCode::FAILURE;
+        }
+        // Find the primary vertex and build up the vector for the TVA tool
+        vertices.reserve(vertexHandle->size());
+        for (const xAOD::Vertex *ivtx : *vertexHandle)
+        {
+          if (!priVtx && ivtx->vertexType() == xAOD::VxType::PriVtx)
+            priVtx = ivtx;
+          vertices.push_back(ivtx);
+        }
+      }
+      auto outputHandle = SG::makeHandle(m_outputKey, ctx);
+      auto &decCategory = *m_decCategory;
+
+      auto combined = std::make_unique<ConstDataVector<xAOD::PFOContainer>>(SG::VIEW_ELEMENTS);
+      combined->reserve(charged->size() + neutral->size());
+
+      for (const xAOD::PFO *ipfo : *charged)
+      {
+        combined->push_back(ipfo);
+        if (m_manualTVA)
+        {
+          if (!priVtx)
+            // If we could not identify a primary vertex, classify as neutral
+            // Note that this means in order to distinguish cPFOs from nPFOs the 'isCharged'
+            // function should still be used
+            decCategory(*ipfo) = PUClassification::NeutralForward;
+          else
+          {
+            const xAOD::TrackParticle *itrk = ipfo->track(0);
+            if (itrk && (m_trackSelTool.empty() || m_trackSelTool->accept(*itrk)))
+            {
+              bool isPVMatched;
+              if (m_useCompatible)
+                isPVMatched = m_tvaTool->isCompatible(*itrk, *priVtx);
+              else
+                isPVMatched = m_tvaTool->getUniqueMatchVertex(*itrk, vertices) == priVtx;
+              if (isPVMatched)
+                decCategory(*ipfo) = PUClassification::ChargedHS;
+              else
+                decCategory(*ipfo) = PUClassification::ChargedPU;
+            }
+          }
+        }
+        else
+        {
+          if (accPVMatched(*ipfo))
+            decCategory(*ipfo) = PUClassification::ChargedHS;
+          else
+            decCategory(*ipfo) = PUClassification::ChargedPU;
+        }
+      }
+      for (const xAOD::PFO *ipfo : *neutral)
+      {
+        combined->push_back(ipfo);
+        decCategory(*ipfo) = PUClassification::NeutralForward;
+      }
+
+      CHECK(outputHandle.put(ctx, std::move(combined)) != nullptr);
+      return StatusCode::SUCCESS;
+    }
+  } // namespace MET
+} // namespace HLT
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/PFOPrepAlg.h b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PFOPrepAlg.h
new file mode 100644
index 000000000000..7093a71e13ba
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PFOPrepAlg.h
@@ -0,0 +1,103 @@
+/*
+   Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+/******************************************************************************
+ * @package Trigger/TrigAlgorithms/TrigEFMissingET
+ * @class   PFOPrepAlg
+ *
+ * @brief Algorithm to prepare PFOs for the PUSplit pufit alg
+ * @author Jon Burr
+ *****************************************************************************/
+
+#ifndef TRIGEFMISSINGET_PFOPREPALG_H
+#define TRIGEFMISSINGET_PFOPREPALG_H 1
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "StoreGate/ReadHandleKey.h"
+#include "StoreGate/WriteDecorHandleKey.h"
+#include "StoreGate/WriteHandleKey.h"
+#include "xAODPFlow/PFOContainer.h"
+#include "xAODTracking/VertexContainer.h"
+#include "AsgTools/ToolHandle.h"
+#include "TrackVertexAssociationTool/ITrackVertexAssociationTool.h"
+#include "InDetTrackSelectionTool/IInDetTrackSelectionTool.h"
+#include <utility>
+
+namespace HLT
+{
+  namespace MET
+  {
+    /**************************************************************************
+     * @class PFOPrepAlg
+     * 
+     * Class to prepare the PFOs for the PUSplit pufit alg.
+     * 
+     * The PUSplit alg expects a single input container, which has a decoration
+     * classifying each object according to the @see PUClassification enum. This
+     * algorithm reads in the input charged and neutral PFO containers and
+     * combines them.
+     * 
+     * The classification either uses the 'matchedToPV' decoration from the cPFO
+     * creation, or a TVA tool, (optional) track selection tool and vertex
+     * container to determine the classification.
+     **************************************************************************/
+    class PFOPrepAlg : public AthReentrantAlgorithm
+    {
+      // Alias optional to make it clear how we're using it
+      template <typename T>
+      using deferred_t = std::optional<T>;
+
+    public:
+      /// Constructor
+      PFOPrepAlg(const std::string &name, ISvcLocator *pSvcLocator);
+
+      virtual StatusCode initialize() override;
+      /// Run the algorithm
+      virtual StatusCode execute(const EventContext &context) const override;
+
+    private:
+      /// The input neutral container
+      SG::ReadHandleKey<xAOD::PFOContainer> m_inputNeutralKey{
+          this, "InputNeutralKey", "", "Input neutral container"};
+      /// The input charged container
+      SG::ReadHandleKey<xAOD::PFOContainer> m_inputChargedKey{
+          this, "InputChargedKey", "", "Input charged container"};
+      /// The input vertex container
+      SG::ReadHandleKey<xAOD::VertexContainer> m_inputVertexKey{
+          this, "InputVertexKey", "",
+          "Input vertex container (only needed if TVA tool is supplied."};
+      /// The output combined container
+      SG::WriteHandleKey<xAOD::PFOContainer> m_outputKey{
+          this, "OutputKey", "", "Output combined container"};
+      /// The output classification
+      SG::WriteDecorHandleKey<xAOD::PFOContainer> m_outputCategoryKey{
+          this, "OutputCategoryKey", "", "Output category name"};
+      /// Track -> vertex association
+      ToolHandle<CP::ITrackVertexAssociationTool> m_tvaTool{
+          this, "TVATool", "",
+          "Track -> vertex association tool for the cPFOS. "
+          "If not provided the matchedToPV decoration is used."};
+      /// Track selection tool
+      ToolHandle<InDet::IInDetTrackSelectionTool> m_trackSelTool{
+          this, "TrackSelTool", "",
+          "The track selection tool (optional). Only used if a TVA tool is provided"};
+      /// Choose between the unique (many-to-one) track->vertex association and
+      /// just checking compatibility with the PV (many-to-many).
+      /// Studies with offline objects show a fairly clear preference for the
+      /// 'compatible' approach so that is the default
+      Gaudi::Property<bool> m_useCompatible{
+          this, "UseCompatible", true,
+          "Use the track -> vertex compatibility over unique matching. "
+          "This only takes effect if a TVA tool is supplied"};
+
+      /// Whether to perform track -> vertex matching manually.
+      bool m_manualTVA{false};
+
+      // Internal
+      deferred_t<SG::AuxElement::Decorator<int>> m_decCategory;
+    };
+  } // namespace MET
+} // namespace HLT
+
+#endif //> !TRIGEFMISSINGET_PFOPREPALG
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/PFSumFex.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PFSumFex.cxx
index d884b86edb6a..7003f196d20e 100644
--- a/Trigger/TrigAlgorithms/TrigEFMissingET/src/PFSumFex.cxx
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PFSumFex.cxx
@@ -33,7 +33,17 @@ namespace HLT { namespace MET {
       MonGroupBuilder&) const
   {
     auto charged = SG::makeHandle(m_chargedPFOKey, context);
+      if (!charged.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_chargedPFOKey);
+        return StatusCode::FAILURE;
+      }
     auto neutral = SG::makeHandle(m_neutralPFOKey, context);
+      if (!neutral.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_neutralPFOKey);
+        return StatusCode::FAILURE;
+      }
 
     std::array<METComponent, 3> pfoSums;
     for (const xAOD::PFO* ipfo : *neutral)
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/PUSplitGrid.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PUSplitGrid.cxx
new file mode 100644
index 000000000000..6a7b945df987
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PUSplitGrid.cxx
@@ -0,0 +1,28 @@
+/*
+ * Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+#include "TrigEFMissingET/PUSplitGrid.h"
+
+namespace HLT
+{
+  namespace MET
+  {
+    PUSplitGrid::PUSplitGrid(
+        double maxEta,
+        std::size_t nEtaTowers,
+        std::size_t nPhiTowers,
+        bool displaceEta,
+        bool displacePhi) : PUSplitGrid(GridParameters{maxEta, nEtaTowers, nPhiTowers, displaceEta, displacePhi})
+    {
+    }
+
+    PUSplitGrid::PUSplitGrid(const GridParameters &parameters) : PufitMultiGrid<3>(parameters)
+    {
+    }
+
+    PUSplitGrid::PUSplitGrid(const PUSplitGrid &other) : PufitMultiGrid<3>(other)
+    {
+    }
+  } // namespace MET
+} // namespace HLT
\ No newline at end of file
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/PUSplitPufitFex.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PUSplitPufitFex.cxx
new file mode 100644
index 000000000000..8d9887d03674
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PUSplitPufitFex.cxx
@@ -0,0 +1,192 @@
+/*
+ * Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+/******************************************************************************
+ * @package Trigger/TrigAlgorithms/TrigEFMissingET
+ * @file PUSplitPufitFex.cxx
+ *
+ * Implementation of pileup-split pufit fex class
+ * @author Jon Burr
+ *****************************************************************************/
+
+#include "PUSplitPufitFex.h"
+#include "TrigEFMissingET/METComponent.h"
+#include "TrigEFMissingET/PUSplitGrid.h"
+#include "TrigEFMissingET/PufitUtils.h"
+#include "StoreGate/DecorKeyHelpers.h"
+#include "StoreGate/ReadDecorHandle.h"
+
+namespace HLT
+{
+  namespace MET
+  {
+    double PUSplitPufitFex::getSigma(const SignedKinematics &kin) const
+    {
+      return m_caloNoise * m_caloNoise + kin.absPt() * m_caloStoch * m_caloStoch;
+    }
+
+    PUSplitPufitFex::PUSplitPufitFex(const std::string &name, ISvcLocator *pSvcLocator)
+        : FexBase(name, pSvcLocator)
+    {
+    }
+
+    StatusCode PUSplitPufitFex::initialize()
+    {
+      CHECK(m_inputKey.initialize());
+
+      // Update the decor keys if necessary
+      if (m_inputCategoryKey.key().find(".") == std::string::npos)
+        m_inputCategoryKey = m_inputKey.key() + "." + m_inputCategoryKey.key();
+      else if (SG::contKeyFromKey(m_inputCategoryKey.key()) != m_inputKey.key())
+      {
+        ATH_MSG_ERROR("Input category key does not match the input key!");
+        return StatusCode::FAILURE;
+      }
+      CHECK(m_inputCategoryKey.initialize());
+      return initializeBase({"NeutralForward", "ChargedHS", "ChargedPU", "UncorrSelNF"});
+    }
+
+    StatusCode PUSplitPufitFex::fillMET(
+        xAOD::TrigMissingET &met,
+        const EventContext &context,
+        MonGroupBuilder &) const
+    {
+      // Retrieve the input
+      auto input = SG::makeHandle(m_inputKey, context);
+      if (!input.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_inputKey);
+        return StatusCode::FAILURE;
+      }
+      auto categoryAcc = SG::makeHandle<int>(m_inputCategoryKey, context);
+
+      // Create the gridset
+      PUSplitGridSet gridset(m_maxEta, m_nEtaBins, m_nPhiBins);
+
+      // Fill the grids from the input objects
+      for (const xAOD::IParticle *ipart : *input)
+      {
+        PUClassification::PUClassification type = PUClassification::PUClassification(categoryAcc(*ipart));
+        SignedKinematics kin = SignedKinematics::fromEnergyEtaPhi(
+            ipart->e(), ipart->eta(), ipart->phi());
+        switch (type)
+        {
+        case PUClassification::NeutralForward:
+          gridset.get<PUClassification::NeutralForward>() += kin;
+          break;
+        case PUClassification::ChargedHS:
+          gridset.get<PUClassification::ChargedHS>() += kin;
+          break;
+        case PUClassification::ChargedPU:
+          gridset.get<PUClassification::ChargedPU>() += kin;
+          break;
+        default:
+          ATH_MSG_ERROR("Invalid PU category: " << type);
+          return StatusCode::FAILURE;
+        }
+      }
+
+      // Calculate the mean and variance from the 'nominal' grid
+      double mean;
+      double variance;
+      PufitUtils::trimmedMeanAndVariance(
+          gridset[NoDisplacement], m_neutralThresholdMode, m_trimFraction, mean, variance);
+      // Calculate the threshold
+      double threshold = mean + m_nSigma * sqrt(variance);
+
+      // Mask towers above the threshold
+      std::size_t count = 0;
+      for (PUSplitGrid &grid : gridset.grids)
+        for (PUSplitGrid::Tower &tower : grid)
+          if (tower.sumEt(m_neutralThresholdMode) > threshold)
+          {
+            tower.mask(true);
+            ++count;
+          }
+
+      // Select the right grid
+      GridDisplacement displacement = PufitUtils::selectGrid(
+          gridset, m_neutralThresholdMode);
+      const PUSplitGrid &grid = gridset[displacement];
+
+      // Build the background sum, the masked tower kinematics (used for their
+      // directions), the list of the expected pileup values in each masked tower
+      // and the variances on each of those expected values
+      PufitUtils::CovarianceSum pileupSum;
+      std::vector<SignedKinematics> masked;
+      masked.reserve(count);
+      std::vector<double> means;
+      means.reserve(count);
+      std::vector<double> variances(count, variance);
+
+      for (const PUSplitGrid::Tower &tower : grid)
+      {
+        const SignedKinematics &cPUKin = tower.get<PUClassification::ChargedPU>();
+        if (!tower.masked())
+        {
+          // If the tower is masked, then add its neutral component to the background sum
+          const SignedKinematics &kin = tower.get<PUClassification::NeutralForward>();
+          pileupSum.add(kin, getSigma(kin));
+        }
+        else
+        {
+          masked.push_back(tower.kinematics(PUClassification::NFcHS));
+          // Add the expected pileup contribution. Optionally, exclude the cPU
+          // component from this as this is already accounted
+          if (m_subtractCPUFromMean)
+            means.push_back(mean - cPUKin.pt());
+          else
+            means.push_back(mean);
+        }
+        // cPU elements always count as background
+        // (There's a possible TODO here - maybe this should only be done when subtractCPUFromMean is true?)
+        pileupSum.add(cPUKin, getSigma(cPUKin));
+      }
+
+      if (msgLvl(MSG::VERBOSE))
+      {
+        ATH_MSG_VERBOSE("Pileup sum = " << pileupSum.sum);
+        ATH_MSG_VERBOSE("Pileup covariance = " << pileupSum.covariance);
+        ATH_MSG_VERBOSE("Mean background energy = " << mean);
+        ATH_MSG_VERBOSE("Background energy variance = " << variance);
+        ATH_MSG_VERBOSE("Constraint importance = " << m_constraintImportance);
+        ATH_MSG_VERBOSE("Masked tower directions (sin, cos): ");
+        for (const SignedKinematics &kin : masked)
+          ATH_MSG_VERBOSE("  (" << kin.sinPhi() << ", " << kin.cosPhi() << ")");
+      }
+      ATH_MSG_DEBUG("Calculate corrections");
+
+      std::vector<SignedKinematics> corrections = PufitUtils::pufit(
+          pileupSum.sum,
+          pileupSum.covariance,
+          means,
+          variances,
+          masked,
+          m_constraintImportance);
+
+      // Fill components
+      grid.get<PUClassification::NeutralForward>().sum(PufitGrid::SumStrategy::All).fillMETComponent(0, met);
+      grid.get<PUClassification::ChargedHS>().sum(PufitGrid::SumStrategy::All).fillMETComponent(1, met);
+      grid.get<PUClassification::ChargedPU>().sum(PufitGrid::SumStrategy::All).fillMETComponent(2, met);
+      grid.get<PUClassification::NeutralForward>().sum(PufitGrid::SumStrategy::Masked).fillMETComponent(3, met);
+
+      // Now build the final sum
+      METComponent sum;
+      // Add the NF components of all masked towers
+      sum += grid.get<PUClassification::NeutralForward>().sum(PufitGrid::SumStrategy::Masked);
+      // Add the cHS components of all towers
+      sum += grid.get<PUClassification::ChargedHS>().sum(PufitGrid::SumStrategy::All);
+      ATH_MSG_DEBUG("Before corrections: " << sum);
+      // Add the corrections
+      for (const SignedKinematics &kin : corrections)
+      {
+        ATH_MSG_DEBUG("Correction: (px, py, pz, et) = (" << kin.px() << ", " << kin.py() << ", " << kin.pz() << ", " << kin.et() << ")");
+        sum += kin;
+      }
+      ATH_MSG_DEBUG("Final MET: " << sum);
+      sum.fillMET(met);
+
+      return StatusCode::SUCCESS;
+    }
+  } // namespace MET
+} // namespace HLT
\ No newline at end of file
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/PUSplitPufitFex.h b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PUSplitPufitFex.h
new file mode 100644
index 000000000000..081e4dc10d19
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PUSplitPufitFex.h
@@ -0,0 +1,114 @@
+/*
+ * Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+/******************************************************************************
+ * @package Trigger/TrigAlgorithms/TrigEFMissingET
+ * @class   PUSplitPufitFex
+ *
+ * @brief Fex class for the pufit algorithm, where inputs are split into 
+ * pile-up categories
+ * @author Jon Burr
+ *****************************************************************************/
+
+#ifndef TRIGEFMISSINGET_PUSPLITPUFITFEX_H
+#define TRIGEFMISSINGET_PUSPLITPUFITFEX_H
+
+#include "FexBase.h"
+#include "TrigEFMissingET/SignedKinematics.h"
+#include "xAODBase/IParticleContainer.h"
+#include "StoreGate/ReadDecorHandleKey.h"
+#include "TrigEFMissingET/PUClassification.h"
+#include "Gaudi/Property.h"
+
+namespace HLT
+{
+  namespace MET
+  {
+    /**************************************************************************
+     * @class PUSplitPufitFex
+     * 
+     * Class to calculate pufit MET using pu-split inputs
+     * 
+     * Uses a generalised pufit technique. The inputs are classified using the
+     * @see PUClassification enum and the algorithm is able to treat them
+     * differently
+     *************************************************************************/
+    class PUSplitPufitFex : public FexBase
+    {
+    public:
+      /// Constructor
+      PUSplitPufitFex(const std::string &name, ISvcLocator *pSvcLocator);
+
+      /// Initialize the fex
+      virtual StatusCode initialize() override;
+
+    private:
+      /************************************************************************
+       * Properties
+       ***********************************************************************/
+      /// Input objects
+      SG::ReadHandleKey<xAOD::IParticleContainer> m_inputKey{
+          this, "InputName", "", "The input particle collection"};
+      SG::ReadDecorHandleKey<xAOD::IParticleContainer> m_inputCategoryKey{
+          this, "InputCategoryName", "PUCategory",
+          "The name of PU classification auxdata"};
+      /// The sigma threshold
+      Gaudi::Property<float> m_nSigma{
+          this, "NSigma", 5, "Set the threshold at mean + NSigma*variance"};
+      /// The eta range of the grid
+      Gaudi::Property<float> m_maxEta{
+          this, "MaxEta", 5, "The maximum eta range"};
+      /// The number of bins in eta
+      Gaudi::Property<std::size_t> m_nEtaBins{
+          this, "NEtaBins", 14, "The number of eta bins"};
+      Gaudi::Property<std::size_t> m_nPhiBins{
+          this, "NPhiBins", 8, "The number of phi bins"};
+      /// The trimming fraction
+      Gaudi::Property<float> m_trimFraction{
+          this, "TrimFraction", 0.9,
+          "The fraction of bins to use when calculating the mean and variance"};
+      /// The coefficient of the noise term in the calo resolution estimate
+      Gaudi::Property<float> m_caloNoise{
+          this, "CaloNoise", 50,
+          "The coefficient of the noise term in the calorimeter resolution estimate [MeV]"};
+      /// The coefficient of the stochastic term in the calo resolution estimate
+      Gaudi::Property<float> m_caloStoch{
+          this, "CaloStochastic", 15.81,
+          "The coefficient of the stochastic term in the calorimeter resolution estimate [MeV^1/2]"};
+      /// The relative constraint importance
+      Gaudi::Property<float> m_constraintImportance{
+          this, "ConstraintImportance", 1,
+          "The relative importance of the two constraints to the fit"};
+      /// The neutral threshold mode
+      Gaudi::Property<std::size_t> m_neutralThresholdMode{
+          this, "NeutralThresholdMode", PUClassification::All,
+          "Which towers to use to calculate mean/variance and the masking threshold"};
+      /// Whether to remove the cPU component from the tower expectations
+      Gaudi::Property<bool> m_subtractCPUFromMean{
+          this, "SubtractCPUFromMean", false,
+          "Whether to remove the cPU component from the expected pileup value"};
+      /************************************************************************
+       * Internal functions
+       ***********************************************************************/
+      /**
+       * @brief Calculate and fill the output MET value
+       * @param met The object to fill
+       * @param context The event context
+       * @param monitors[out] Any extra monitors to fill
+       */
+      virtual StatusCode fillMET(
+          xAOD::TrigMissingET &met,
+          const EventContext &context,
+          MonGroupBuilder &monitors) const override;
+
+      /**
+       * @brief Calculate the estimate on the variance of a tower
+       * @param kin The kinematics of the tower
+       */
+      double getSigma(const SignedKinematics &kin) const;
+    }; //> end class PUSplitPufixFex
+  }    // namespace MET
+} // namespace HLT
+
+#endif //> !TRIGEFMISSINGET_PUSPLITPUFITFEX_H
\ No newline at end of file
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/PeriodicGridBase.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PeriodicGridBase.cxx
new file mode 100644
index 000000000000..c8c06dd2c8dd
--- /dev/null
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PeriodicGridBase.cxx
@@ -0,0 +1,157 @@
+/*
+ * Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+#include "TrigEFMissingET/PeriodicGridBase.h"
+#include <cmath>
+
+namespace
+{
+  constexpr double pi = 3.141592653589793238462643383279502884;
+}
+
+namespace HLT
+{
+  namespace MET
+  {
+    bool GridParameters::operator!=(const GridParameters &other) const
+    {
+      return maxEta != other.maxEta ||
+             nEtaTowers != other.nEtaTowers ||
+             nPhiTowers != other.nPhiTowers ||
+             displaceEta != other.displaceEta ||
+             displacePhi != other.displacePhi;
+    }
+
+    bool GridParameters::operator==(const GridParameters &other) const
+    {
+      return !(*this != other);
+    }
+
+    GridDisplacement GridParameters::displacement() const
+    {
+      return static_cast<GridDisplacement>(
+          (displaceEta ? EtaDisplaced : NoDisplacement) |
+          (displacePhi ? PhiDisplaced : NoDisplacement));
+    }
+
+    PeriodicGridBase::Tower::Tower(std::size_t index) : m_index(index) {}
+
+    double PeriodicGridBase::Tower::towerPhi() const
+    {
+      return grid()->centralPhi(m_index);
+    }
+
+    double PeriodicGridBase::Tower::towerEta() const
+    {
+      return grid()->centralEta(m_index);
+    }
+
+    std::size_t PeriodicGridBase::Tower::index() const
+    {
+      return m_index;
+    }
+
+    std::size_t PeriodicGridBase::Tower::etaIndex() const
+    {
+      return etaPhiIndex().first;
+    }
+
+    std::size_t PeriodicGridBase::Tower::phiIndex() const
+    {
+      return etaPhiIndex().second;
+    }
+
+    std::pair<std::size_t, std::size_t> PeriodicGridBase::Tower::etaPhiIndex() const
+    {
+      return grid()->etaPhiIndex(m_index);
+    }
+
+    PeriodicGridBase::PeriodicGridBase(const GridParameters &parameters)
+        : m_params(parameters)
+    {
+    }
+
+    PeriodicGridBase::PeriodicGridBase(
+        double maxEta,
+        std::size_t nEtaTowers,
+        std::size_t nPhiTowers,
+        bool displaceEta,
+        bool displacePhi)
+        : PeriodicGridBase(GridParameters{maxEta, nEtaTowers, nPhiTowers, displaceEta, displacePhi})
+    {
+    }
+
+    std::size_t PeriodicGridBase::getIndex(double eta, double phi, bool &outOfRange) const
+    {
+      std::size_t etaIndex = getEtaIndex(eta, outOfRange);
+      if (outOfRange)
+        return nTowers();
+      return globalIndex(etaIndex, getPhiIndex(phi));
+    }
+
+    std::size_t PeriodicGridBase::getEtaIndex(double eta, bool &outOfRange) const
+    {
+      if (std::abs(eta) >= maxEta())
+      {
+        outOfRange = true;
+        return nEtaTowers();
+      }
+      outOfRange = false;
+      if (displaceEta())
+      {
+        // Apply the displacement by adding the displacement to the input
+        // coordinate. This shifts the grid in the negative direction
+        eta += etaWidth() / 2;
+        // If necessary apply eta periodicity here
+        if (eta >= maxEta())
+          eta -= 2 * maxEta();
+      }
+      return (eta + maxEta()) / etaWidth();
+    }
+
+    std::size_t PeriodicGridBase::getPhiIndex(double phi) const
+    {
+      if (displacePhi())
+        // Apply the displacement by adding the displacement to the input
+        // coordinate. This shifts the grid in the negative direction
+        phi += phiWidth() / 2;
+      // Apply periodicity
+      phi = std::fmod(phi, 2 * pi);
+      if (phi < 0)
+        phi += 2 * pi;
+      return phi / phiWidth();
+    }
+
+    std::size_t PeriodicGridBase::globalIndex(std::size_t iEta, std::size_t iPhi) const
+    {
+      return iEta * nPhiTowers() + iPhi;
+    }
+    std::pair<std::size_t, std::size_t> PeriodicGridBase::etaPhiIndex(std::size_t index) const
+    {
+      return std::make_pair(index / nPhiTowers(), index % nPhiTowers());
+    }
+
+    double PeriodicGridBase::centralEta(std::size_t iEta) const
+    {
+      return -maxEta() + etaWidth() * (iEta + 0.5) + (displaceEta() ? etaWidth() / 2 : 0);
+    }
+    double PeriodicGridBase::centralPhi(std::size_t iPhi) const
+    {
+      return phiWidth() * (iPhi + 0.5) + (displacePhi() ? phiWidth() / 2 : 0);
+    }
+
+    const GridParameters &PeriodicGridBase::parameters() const { return m_params; }
+    double PeriodicGridBase::maxEta() const { return parameters().maxEta; }
+    std::size_t PeriodicGridBase::nEtaTowers() const { return parameters().nEtaTowers; }
+    std::size_t PeriodicGridBase::nPhiTowers() const { return parameters().nPhiTowers; }
+    std::size_t PeriodicGridBase::nTowers() const { return nEtaTowers() * nPhiTowers(); }
+    bool PeriodicGridBase::displaceEta() const { return parameters().displaceEta; }
+    bool PeriodicGridBase::displacePhi() const { return parameters().displacePhi; }
+    GridDisplacement PeriodicGridBase::displacement() const
+    {
+      return parameters().displacement();
+    }
+    double PeriodicGridBase::etaWidth() const { return 2 * maxEta() / nEtaTowers(); }
+    double PeriodicGridBase::phiWidth() const { return 2 * pi / nPhiTowers(); }
+  } // namespace MET
+} // namespace HLT
\ No newline at end of file
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/PufitGrid.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PufitGrid.cxx
index db70958f5458..367b9f04da29 100644
--- a/Trigger/TrigAlgorithms/TrigEFMissingET/src/PufitGrid.cxx
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PufitGrid.cxx
@@ -7,350 +7,252 @@
 #include <algorithm>
 #include <cmath>
 
-
-namespace HLT { namespace MET {
-
-  PufitGrid::Tower::Tower(const PufitGrid* parent, std::size_t index) :
-    m_parent(parent),
-    m_index(index)
+namespace HLT
+{
+  namespace MET
   {
-  }
-
-  PufitGrid::Tower& PufitGrid::Tower::operator=(const Tower& other)
-  {
-    m_kin = other.kinematics();
-    m_sumEt = other.sumEt();
-    m_sumE = other.sumE();
-    m_mask = other.masked();
-    return *this;
-  }
-
-  double PufitGrid::Tower::ex() const { return m_kin.px(); }
-  double PufitGrid::Tower::ey() const { return m_kin.py(); }
-  double PufitGrid::Tower::ez() const { return m_kin.pz(); }
-  double PufitGrid::Tower::sumEt() const { return m_sumEt; }
-  double PufitGrid::Tower::sumE() const { return m_sumE; }
-  double PufitGrid::Tower::phi() const { return kinematics().phi(); }
-  double PufitGrid::Tower::eta() const { return kinematics().eta(); }
-  bool PufitGrid::Tower::masked() const { return m_mask; }
-  void PufitGrid::Tower::mask(bool value) { m_mask = value; }
 
-  double PufitGrid::Tower::towerPhi() const {
-    return grid()->centralPhi(phiIndex() );
-  }
-  double PufitGrid::Tower::towerEta() const {
-    return grid()->centralEta(etaIndex() );
-  }
-  std::size_t PufitGrid::Tower::index() const { return m_index; }
-  std::size_t PufitGrid::Tower::etaIndex() const { 
-    return etaPhiIndex().first;
-  }
-  std::size_t PufitGrid::Tower::phiIndex() const {
-    return etaPhiIndex().second;
-  }
-  std::pair<std::size_t, std::size_t> PufitGrid::Tower::etaPhiIndex() const {
-    return grid()->etaPhiIndex(index() );
-  }
-  const PufitGrid* PufitGrid::Tower::grid() const { return m_parent; }
-
-  SignedKinematics PufitGrid::Tower::kinematics() const {
-    return m_kin;
-  }
-  PufitGrid::Tower::operator SignedKinematics() const { return kinematics(); }
-
-  PufitGrid::Tower& PufitGrid::Tower::operator+=(const SignedKinematics& kin)
-  {
-    m_kin += kin;
-    m_sumEt += kin.pt();
-    m_sumE += kin.energy();
-    return *this;
-  }
-  PufitGrid::Tower& PufitGrid::Tower::operator-=(const SignedKinematics& kin)
-  {
-    m_kin -= kin;
-    m_sumEt -= kin.pt();
-    m_sumE -= kin.energy();
-    return *this;
-  }
-  PufitGrid::Tower& PufitGrid::Tower::operator+=(const Tower& other)
-  {
-    m_kin += other.kinematics();
-    m_sumEt += other.sumEt();
-    m_sumE += other.sumE();
-    return *this;
-  }
-  PufitGrid::Tower& PufitGrid::Tower::operator-=(const Tower& other)
-  {
-    m_kin -= other.kinematics();
-    m_sumEt -= other.sumEt();
-    m_sumE -= other.sumE();
-    return *this;
-  }
-
-  bool PufitGrid::GridParameters::operator!=(const GridParameters& other) const
-  {
-    return maxEta != other.maxEta ||
-      nEtaTowers != other.nEtaTowers ||
-      nPhiTowers != other.nPhiTowers ||
-      displaceEta != other.displaceEta ||
-      displacePhi != other.displacePhi;
-  }
-  bool PufitGrid::GridParameters::operator==(const GridParameters& other) const
-  {
-    return !(*this != other);
-  }
-
-  PufitGrid::PufitGrid(
-      double maxEta,
-      std::size_t nEtaTowers,
-      std::size_t nPhiTowers,
-      bool displaceEta,
-      bool displacePhi) :
-    PufitGrid(
-        GridParameters{maxEta, nEtaTowers, nPhiTowers, displaceEta, displacePhi})
-  {}
+    PufitGrid::Tower::Tower(const PufitGrid *parent, std::size_t index) : PeriodicGridBase::Tower(index),
+                                                                          m_parent(parent)
+    {
+    }
 
-  PufitGrid::PufitGrid(const GridParameters& parameters) :
-    m_params(parameters),
-    m_etaWidth(2*maxEta()/nEtaTowers() ),
-    m_phiWidth(2*maxEta()/nPhiTowers() )
-  {
-    m_towers.reserve(nTowers());
-    for (std::size_t index = 0; index < nTowers(); ++index)
-      m_towers.emplace_back(this, index);
-  }
+    PufitGrid::Tower &PufitGrid::Tower::operator=(const Tower &other)
+    {
+      m_kin = other.kinematics();
+      m_sumEt = other.sumEt();
+      m_sumE = other.sumE();
+      m_mask = other.masked();
+      return *this;
+    }
 
-  PufitGrid::PufitGrid(const PufitGrid& other) :
-    PufitGrid(other.parameters() )
-  {
-    *this = other;
-  }
+    double PufitGrid::Tower::ex() const { return m_kin.px(); }
+    double PufitGrid::Tower::ey() const { return m_kin.py(); }
+    double PufitGrid::Tower::ez() const { return m_kin.pz(); }
+    double PufitGrid::Tower::sumEt() const { return m_sumEt; }
+    double PufitGrid::Tower::sumE() const { return m_sumE; }
+    double PufitGrid::Tower::phi() const { return kinematics().phi(); }
+    double PufitGrid::Tower::eta() const { return kinematics().eta(); }
+    bool PufitGrid::Tower::masked() const { return m_mask; }
+    void PufitGrid::Tower::mask(bool value) { m_mask = value; }
 
-  PufitGrid& PufitGrid::operator=(const PufitGrid& other)
-  {
-    if (parameters() != other.parameters() )
-      throw std::invalid_argument("Grid parameters do not match!");
-    std::copy(other.begin(), other.end(), m_towers.begin() );
-    return *this;
-  }
+    const PufitGrid *PufitGrid::Tower::grid() const { return m_parent; }
 
-  void PufitGrid::reset()
-  {
-    // Here we construct a standalone tower to copy in 0 values to the tower
-    // kinematics. This kind of object is safe here but in general it's not a
-    // good idea to have one.
-    std::fill(begin(), end(), Tower(nullptr, -1) );
-  }
+    SignedKinematics PufitGrid::Tower::kinematics() const
+    {
+      return m_kin;
+    }
+    PufitGrid::Tower::operator SignedKinematics() const { return kinematics(); }
 
-  PufitGrid& PufitGrid::operator+=(const SignedKinematics& kin)
-  {
-    // Find the right tower to add to
-    bool outOfRange = false;
-    std::size_t index = getIndex(kin.eta(), kin.phi(), outOfRange);
-    if (!outOfRange)
-      operator[](index) += kin;
-    return *this;
-  }
-  PufitGrid& PufitGrid::operator-=(const SignedKinematics& kin)
-  {
-    // Find the right tower to subtract from
-    bool outOfRange = false;
-    std::size_t index = getIndex(kin.eta(), kin.phi(), outOfRange);
-    if (!outOfRange)
-      operator[](index) -= kin;
-    return *this;
-  }
+    PufitGrid::Tower &PufitGrid::Tower::operator+=(const SignedKinematics &kin)
+    {
+      m_kin += kin;
+      m_sumEt += kin.pt();
+      m_sumE += kin.energy();
+      return *this;
+    }
+    PufitGrid::Tower &PufitGrid::Tower::operator-=(const SignedKinematics &kin)
+    {
+      m_kin -= kin;
+      m_sumEt -= kin.pt();
+      m_sumE -= kin.energy();
+      return *this;
+    }
+    PufitGrid::Tower &PufitGrid::Tower::operator+=(const Tower &other)
+    {
+      m_kin += other.kinematics();
+      m_sumEt += other.sumEt();
+      m_sumE += other.sumE();
+      return *this;
+    }
+    PufitGrid::Tower &PufitGrid::Tower::operator-=(const Tower &other)
+    {
+      m_kin -= other.kinematics();
+      m_sumEt -= other.sumEt();
+      m_sumE -= other.sumE();
+      return *this;
+    }
 
-  std::size_t PufitGrid::getIndex(
-      double eta,
-      double phi,
-      bool& outOfRange) const
-  {
-    std::size_t etaIndex = getEtaIndex(eta, outOfRange);
-    if (outOfRange)
-      return nTowers();
-    std::size_t phiIndex = getPhiIndex(phi);
-    return globalIndex(etaIndex, phiIndex);
-  }
+    PufitGrid::PufitGrid(
+        double maxEta,
+        std::size_t nEtaTowers,
+        std::size_t nPhiTowers,
+        bool displaceEta,
+        bool displacePhi)
+        : PufitGrid(GridParameters{maxEta, nEtaTowers, nPhiTowers, displaceEta, displacePhi})
+    {
+    }
 
-  std::size_t PufitGrid::getEtaIndex(double eta, bool& outOfRange) const
-  {
-    if (fabs(eta) >= maxEta() ) {
-      outOfRange = true;
-      return nEtaTowers();
-    }
-    outOfRange = false;
-    if (displaceEta() ) {
-      // Apply the displacement by adding the displacement to the input
-      // coordinate. This shifts the grid in the negative direction
-      eta += etaWidth()/2;
-      // If necessary apply eta periodicity here
-      if (eta > + maxEta() )
-        eta -= 2 * maxEta();
-    }
-    return (eta + maxEta() ) / etaWidth();
-  }
+    PufitGrid::PufitGrid(const GridParameters &parameters)
+        : PeriodicGridBase(parameters)
+    {
+      m_towers.reserve(nTowers());
+      for (std::size_t index = 0; index < nTowers(); ++index)
+        m_towers.emplace_back(this, index);
+    }
 
-  std::size_t PufitGrid::getPhiIndex(double phi) const
-  {
-    if (displacePhi() )
-      // Apply the displacement by adding the displacement to the input
-      // coordinate. This shifts the grid in the negative direction
-      phi += phiWidth()/2;
-    // Apply periodicity
-    phi = fmod(phi, 2*TMath::Pi() );
-    if (phi < 0)
-      phi += 2*TMath::Pi();
-    return phi / phiWidth();
-  }
+    PufitGrid::PufitGrid(const PufitGrid &other) : PufitGrid(other.parameters())
+    {
+      *this = other;
+    }
 
-  PufitGrid::Tower& PufitGrid::operator[](
-      const std::pair<std::size_t, std::size_t>& indices)
-  {
-    return operator[](globalIndex(indices.first, indices.second) );
-  }
-  const PufitGrid::Tower& PufitGrid::operator[](
-      const std::pair<std::size_t, std::size_t>& indices) const
-  {
-    return operator[](globalIndex(indices.first, indices.second) );
-  }
+    PufitGrid &PufitGrid::operator=(const PufitGrid &other)
+    {
+      if (parameters() != other.parameters())
+        throw std::invalid_argument("Grid parameters do not match!");
+      std::copy(other.begin(), other.end(), m_towers.begin());
+      return *this;
+    }
 
-  PufitGrid::Tower& PufitGrid::operator[](std::size_t index)
-  {
-    return m_towers.at(index);
-  }
-  const PufitGrid::Tower& PufitGrid::operator[](std::size_t index) const
-  {
-    return m_towers.at(index);
-  }
+    void PufitGrid::reset()
+    {
+      // Here we construct a standalone tower to copy in 0 values to the tower
+      // kinematics. This kind of object is safe here but in general it's not a
+      // good idea to have one.
+      std::fill(begin(), end(), Tower(nullptr, -1));
+    }
 
-  std::vector<PufitGrid::Tower>::iterator PufitGrid::begin()
-  {
-    return m_towers.begin();
-  }
-  std::vector<PufitGrid::Tower>::const_iterator PufitGrid::begin() const
-  {
-    return m_towers.begin();
-  }
-  std::vector<PufitGrid::Tower>::iterator PufitGrid::end()
-  {
-    return m_towers.end();
-  }
-  std::vector<PufitGrid::Tower>::const_iterator PufitGrid::end() const
-  {
-    return m_towers.end();
-  }
+    PufitGrid &PufitGrid::operator+=(const SignedKinematics &kin)
+    {
+      // Find the right tower to add to
+      bool outOfRange = false;
+      std::size_t index = getIndex(kin.eta(), kin.phi(), outOfRange);
+      if (!outOfRange)
+        operator[](index) += kin;
+      return *this;
+    }
+    PufitGrid &PufitGrid::operator-=(const SignedKinematics &kin)
+    {
+      // Find the right tower to subtract from
+      bool outOfRange = false;
+      std::size_t index = getIndex(kin.eta(), kin.phi(), outOfRange);
+      if (!outOfRange)
+        operator[](index) -= kin;
+      return *this;
+    }
 
-  std::size_t PufitGrid::globalIndex(std::size_t iEta, std::size_t iPhi) const
-  {
-    return iEta*nPhiTowers() + iPhi;
-  }
-  std::pair<std::size_t, std::size_t> PufitGrid::etaPhiIndex(std::size_t index) const
-  {
-    return std::make_pair(index/nPhiTowers(), index%nPhiTowers() );
-  }
+    PufitGrid::Tower &PufitGrid::operator[](
+        const std::pair<std::size_t, std::size_t> &indices)
+    {
+      return operator[](globalIndex(indices.first, indices.second));
+    }
+    const PufitGrid::Tower &PufitGrid::operator[](
+        const std::pair<std::size_t, std::size_t> &indices) const
+    {
+      return operator[](globalIndex(indices.first, indices.second));
+    }
 
-  double PufitGrid::centralEta(std::size_t iEta) const
-  {
-    return -maxEta() + etaWidth()*(iEta+0.5) + (displaceEta() ? etaWidth()/2 : 0);
-  }
-  double PufitGrid::centralPhi(std::size_t iPhi) const
-  {
-    return phiWidth()*(iPhi+0.5) + (displacePhi() ? phiWidth()/2 : 0);
-  }
+    PufitGrid::Tower &PufitGrid::operator[](std::size_t index)
+    {
+      return m_towers.at(index);
+    }
+    const PufitGrid::Tower &PufitGrid::operator[](std::size_t index) const
+    {
+      return m_towers.at(index);
+    }
 
-  const PufitGrid::GridParameters& PufitGrid::parameters() const { return m_params; }
-  double PufitGrid::maxEta() const { return parameters().maxEta; }
-  std::size_t PufitGrid::nEtaTowers() const { return parameters().nEtaTowers; }
-  std::size_t PufitGrid::nPhiTowers() const { return parameters().nPhiTowers; }
-  std::size_t PufitGrid::nTowers() const { return nEtaTowers()*nPhiTowers(); }
-  bool PufitGrid::displaceEta() const { return parameters().displaceEta; }
-  bool PufitGrid::displacePhi() const { return parameters().displacePhi; }
-  double PufitGrid::etaWidth() const { return 2*maxEta() / nEtaTowers(); }
-  double PufitGrid::phiWidth() const { return 2*TMath::Pi() / nPhiTowers(); }
+    std::vector<PufitGrid::Tower>::iterator PufitGrid::begin()
+    {
+      return m_towers.begin();
+    }
+    std::vector<PufitGrid::Tower>::const_iterator PufitGrid::begin() const
+    {
+      return m_towers.begin();
+    }
+    std::vector<PufitGrid::Tower>::iterator PufitGrid::end()
+    {
+      return m_towers.end();
+    }
+    std::vector<PufitGrid::Tower>::const_iterator PufitGrid::end() const
+    {
+      return m_towers.end();
+    }
 
-  METComponent PufitGrid::sum(SumStrategy strategy) const
-  {
-    METComponent total;
-    switch (strategy) {
+    METComponent PufitGrid::sum(SumStrategy strategy) const
+    {
+      METComponent total;
+      switch (strategy)
+      {
       case SumStrategy::All:
-        for (const Tower& tower : *this) total += tower;
+        for (const Tower &tower : *this)
+          total += tower;
         break;
       case SumStrategy::Masked:
-        for (const Tower& tower : *this)
-          if (tower.masked() )
+        for (const Tower &tower : *this)
+          if (tower.masked())
             total += tower;
         break;
       case SumStrategy::Unmasked:
-        for (const Tower& tower : *this)
-          if (!tower.masked() )
+        for (const Tower &tower : *this)
+          if (!tower.masked())
             total += tower;
         break;
+      }
+      return total;
     }
-    return total;
-  }
 
-  PufitGrid& PufitGrid::operator+=(const PufitGrid& other)
-  {
-    if (parameters() != other.parameters() )
-      throw std::invalid_argument("Grid parameters do not match");
-    auto itr = begin();
-    auto otherItr = other.begin();
-    for (; itr != end(); ++itr, ++otherItr)
-      *itr += *otherItr;
-    return *this;
-  }
-  PufitGrid& PufitGrid::operator-=(const PufitGrid& other)
-  {
-    if (parameters() != other.parameters() )
-      throw std::invalid_argument("Grid parameters do not match");
-    auto itr = begin();
-    auto otherItr = other.begin();
-    for (; itr != end(); ++itr, ++otherItr)
-      *itr -= *otherItr;
-    return *this;
-  }
+    PufitGrid &PufitGrid::operator+=(const PufitGrid &other)
+    {
+      if (parameters() != other.parameters())
+        throw std::invalid_argument("Grid parameters do not match");
+      auto itr = begin();
+      auto otherItr = other.begin();
+      for (; itr != end(); ++itr, ++otherItr)
+        *itr += *otherItr;
+      return *this;
+    }
+    PufitGrid &PufitGrid::operator-=(const PufitGrid &other)
+    {
+      if (parameters() != other.parameters())
+        throw std::invalid_argument("Grid parameters do not match");
+      auto itr = begin();
+      auto otherItr = other.begin();
+      for (; itr != end(); ++itr, ++otherItr)
+        *itr -= *otherItr;
+      return *this;
+    }
 
-  PufitGridSet::PufitGridSet(double maxEta, std::size_t nEta, std::size_t nPhi) :
-    grids({
-      PufitGrid(maxEta, nEta, nPhi, false, false),
-      PufitGrid(maxEta, nEta, nPhi, true,  false),
-      PufitGrid(maxEta, nEta, nPhi, false, true),
-      PufitGrid(maxEta, nEta, nPhi, true,  true)})
-  {}
+    PufitGridSet::PufitGridSet(double maxEta, std::size_t nEta, std::size_t nPhi)
+        : grids({PufitGrid(maxEta, nEta, nPhi, false, false),
+                 PufitGrid(maxEta, nEta, nPhi, true, false),
+                 PufitGrid(maxEta, nEta, nPhi, false, true),
+                 PufitGrid(maxEta, nEta, nPhi, true, true)})
+    {
+    }
 
-  PufitGridSet& PufitGridSet::operator+=(const SignedKinematics& kin)
-  {
-    for (PufitGrid& grid : grids)
-      grid += kin;
-    return *this;
-  }
-  PufitGridSet& PufitGridSet::operator-=(const SignedKinematics& kin)
-  {
-    for (PufitGrid& grid : grids)
-      grid += kin;
-    return *this;
-  }
-  PufitGrid& PufitGridSet::operator[](GridDisplacement displacement)
-  {
-    return grids[displacement];
-  }
-  const PufitGrid& PufitGridSet::operator[](GridDisplacement displacement) const
-  {
-    return grids[displacement];
-  }
+    PufitGridSet &PufitGridSet::operator+=(const SignedKinematics &kin)
+    {
+      for (PufitGrid &grid : grids)
+        grid += kin;
+      return *this;
+    }
+    PufitGridSet &PufitGridSet::operator-=(const SignedKinematics &kin)
+    {
+      for (PufitGrid &grid : grids)
+        grid -= kin;
+      return *this;
+    }
+    PufitGrid &PufitGridSet::operator[](GridDisplacement displacement)
+    {
+      return grids[displacement];
+    }
+    const PufitGrid &PufitGridSet::operator[](GridDisplacement displacement) const
+    {
+      return grids[displacement];
+    }
 
-  PufitGrid operator+(const PufitGrid& lhs, const PufitGrid& rhs)
-  {
-    PufitGrid ret(lhs);
-    ret += rhs;
-    return ret;
-  }
-  PufitGrid operator-(const PufitGrid& lhs, const PufitGrid& rhs)
-  {
-    PufitGrid ret(lhs);
-    ret -= rhs;
-    return ret;
-  }
+    PufitGrid operator+(const PufitGrid &lhs, const PufitGrid &rhs)
+    {
+      PufitGrid ret(lhs);
+      ret += rhs;
+      return ret;
+    }
+    PufitGrid operator-(const PufitGrid &lhs, const PufitGrid &rhs)
+    {
+      PufitGrid ret(lhs);
+      ret -= rhs;
+      return ret;
+    }
 
-} } //> end namespace HLT::MET
+  } // namespace MET
+} // namespace HLT
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/PufitUtils.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PufitUtils.cxx
index eb104b30bd2f..23b054156966 100644
--- a/Trigger/TrigAlgorithms/TrigEFMissingET/src/PufitUtils.cxx
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/PufitUtils.cxx
@@ -6,254 +6,271 @@
 #include <algorithm>
 #include <numeric>
 
+namespace HLT
+{
+  namespace MET
+  {
+    namespace PufitUtils
+    {
 
-namespace HLT { namespace MET { namespace PufitUtils {
-
-  CovarianceSum::CovarianceSum() :
-    sum(Eigen::Vector2d::Zero() ),
-    covariance(Eigen::Matrix2d::Zero() )
-  {}
+      CovarianceSum::CovarianceSum() : sum(Eigen::Vector2d::Zero()),
+                                       covariance(Eigen::Matrix2d::Zero())
+      {
+      }
 
-  CovarianceSum::CovarianceSum(
-      const Eigen::Vector2d& sum,
-      const Eigen::Matrix2d& covariance) :
-    sum(sum), covariance(covariance) {}
+      CovarianceSum::CovarianceSum(
+          const Eigen::Vector2d &sum,
+          const Eigen::Matrix2d &covariance) : sum(sum), covariance(covariance) {}
 
-  CovarianceSum& CovarianceSum::add(
-      const SignedKinematics& kin,
-      double sigma)
-  {
-    Eigen::Vector2d cosSin(kin.cosPhi(), kin.sinPhi() );
-    sum += Eigen::Vector2d(kin.px(), kin.py() );
-    covariance += sigma*cosSin*cosSin.transpose();
-    return *this;
-  }
+      CovarianceSum &CovarianceSum::add(
+          const SignedKinematics &kin,
+          double sigma)
+      {
+        Eigen::Vector2d cosSin(kin.cosPhi(), kin.sinPhi());
+        sum += Eigen::Vector2d(kin.px(), kin.py());
+        covariance += sigma * cosSin * cosSin.transpose();
+        return *this;
+      }
 
-  std::pair<double, double> trimmedMeanAndVariance(
-      const PufitGrid& grid,
-      double trimFraction)
-  {
-    // Construct a list of sorted tower energies
-    std::vector<double> sorted;
-    sorted.reserve(grid.nTowers() );
-    for (const PufitGrid::Tower& tower : grid)
-      sorted.insert(
-          std::lower_bound(sorted.begin(), sorted.end(), tower.sumEt() ),
-          tower.sumEt() );
+      void trimmedMeanAndVariance(
+          const std::vector<double> &sorted,
+          double trimFraction,
+          double &mean,
+          double &variance)
+      {
+        // The number to trim from each side
+        std::size_t nTrim = sorted.size() * (1 - trimFraction) / 2;
 
-    // The number to trim from each side
-    std::size_t nTrim = grid.nTowers() * (1-trimFraction)/2;
-    //
-    double mean = std::accumulate(sorted.begin()+nTrim, sorted.end()-nTrim, 0.);
-    mean /= (grid.nTowers() - 2*nTrim);
+        mean = std::accumulate(sorted.begin() + nTrim, sorted.end() - nTrim, 0.);
+        mean /= (sorted.size() - 2 * nTrim);
 
-    // Now get the variance. The corresponding accumulate call here would be
-    // rather ugly so just write the loop
-    double variance = 0;
-    for (std::size_t ii = 0; ii < grid.nTowers()-nTrim; ++ii) {
-      double val = sorted.at(ii)-mean;
-      variance += val*val * (ii < nTrim ? 2 : 1);
-    }
-    variance /= grid.nTowers();
-    return std::make_pair(mean, variance);
-  }
+        // Now get the variance
+        variance = 0;
+        for (std::size_t ii = 0; ii < sorted.size() - nTrim; ++ii)
+        {
+          double val = sorted.at(ii) - mean;
+          variance += val * val * (ii < nTrim ? 2 : 1);
+        }
+        variance /= sorted.size();
+      }
 
-  std::pair<double, double> unmaskedMeanAndVariance(const PufitGrid& grid)
-  {
-    double sum = 0;
-    double squaredSum = 0;
-    std::size_t n = 0;
-    for (const PufitGrid::Tower& tower : grid) {
-      if (tower.masked() )
-        continue;
-      ++n;
-      sum += tower.sumEt();
-      squaredSum += tower.sumEt()*tower.sumEt();
-    }
-    double mean = sum/n;
-    // Note that this could result in catastrophic cancellation in some cases.
-    // However the amount of precision present in a double is around 10^15 so we
-    // can afford to lose some significant figures
-    double variance = squaredSum/n - mean*mean;
-    return std::make_pair(mean, variance);
-  }
+      void trimmedMeanAndVariance(
+          const PufitGrid &grid,
+          double trimFraction,
+          double &mean,
+          double &variance)
+      {
+        std::vector<double> sorted;
+        sorted.reserve(grid.nTowers());
+        for (const PufitGrid::Tower &tower : grid)
+          sorted.insert(
+              std::lower_bound(sorted.begin(), sorted.end(), tower.sumEt()),
+              tower.sumEt());
+        trimmedMeanAndVariance(sorted, trimFraction, mean, variance);
+      }
 
-  PufitGridSet::GridDisplacement selectGrid(const PufitGridSet& grids)
-  {
-    PufitGridSet::GridDisplacement maximum = PufitGridSet::NoDisplacement;
-    double maxSum = 0;
-    for (std::size_t d = 0; d < 4; ++d) {
-      double sum = 0;
-      for (const PufitGrid::Tower& tower : grids[PufitGridSet::GridDisplacement(d)])
-        if (tower.masked() )
+      void unmaskedMeanAndVariance(
+          const PufitGrid &grid,
+          double &mean,
+          double &variance)
+      {
+        double sum = 0;
+        double squaredSum = 0;
+        std::size_t n = 0;
+        for (const PufitGrid::Tower &tower : grid)
+        {
+          if (tower.masked())
+            continue;
+          ++n;
           sum += tower.sumEt();
-      if (sum > maxSum) {
-        maximum = PufitGridSet::GridDisplacement(d);
-        maxSum = sum;
+          squaredSum += tower.sumEt() * tower.sumEt();
+        }
+        mean = sum / n;
+        // Note that this could result in catastrophic cancellation in some cases.
+        // However the amount of precision present in a double is around 10^15 so we
+        // can afford to lose some significant figures
+        variance = squaredSum / n - mean * mean;
       }
-    }
-    return maximum;
-  }
-
-  Eigen::VectorXd pufit(
-      const Eigen::Vector2d& pileupSum,
-      const Eigen::Matrix2d& pileupCovariance,
-      const Eigen::VectorXd& towerExpectations,
-      const Eigen::VectorXd& towerVariances,
-      const Eigen::VectorXd& correctionDirections,
-      double constraintImportance)
-  {
-    // Keep track of the number of corrections to derive
-    std::size_t nCorr = correctionDirections.size();
-    if (nCorr == 0)
-      // No need to do anything to derive corrections to nothing
-      return {};
-    // Turn the directions into cos/sin pairs
-    Eigen::Matrix<double, 2, Eigen::Dynamic> cosSin(2, nCorr);
-    cosSin.row(0) = correctionDirections.array().cos();
-    cosSin.row(1) = correctionDirections.array().sin();
-    return pufit(
-        pileupSum,
-        pileupCovariance,
-        towerExpectations,
-        towerVariances,
-        cosSin,
-        constraintImportance);
-  }
 
-  Eigen::VectorXd pufit(
-      const Eigen::Vector2d& pileupSum,
-      const Eigen::Matrix2d& pileupCovariance,
-      double towerMean,
-      double towerVariance,
-      const Eigen::VectorXd& correctionDirections,
-      double constraintImportance)
-  {
-    std::size_t nCorr = correctionDirections.size();
-    return pufit(
-        pileupSum,
-        pileupCovariance,
-        Eigen::VectorXd::Constant(nCorr, towerMean),
-        Eigen::VectorXd::Constant(nCorr, towerVariance),
-        correctionDirections,
-        constraintImportance);
-  }
+      GridDisplacement selectGrid(const PufitGridSet &grids)
+      {
+        GridDisplacement maximum = NoDisplacement;
+        double maxSum = 0;
+        for (std::size_t d = 0; d < 4; ++d)
+        {
+          double sum = 0;
+          for (const PufitGrid::Tower &tower : grids[GridDisplacement(d)])
+            if (tower.masked())
+              sum += tower.sumEt();
+          if (sum > maxSum)
+          {
+            maximum = GridDisplacement(d);
+            maxSum = sum;
+          }
+        }
+        return maximum;
+      }
 
-  Eigen::VectorXd pufit(
-      const Eigen::Vector2d& pileupSum,
-      const Eigen::Matrix2d& pileupCovariance,
-      const Eigen::VectorXd& towerExpectations,
-      const Eigen::VectorXd& towerVariances,
-      const Eigen::Matrix<double, 2, Eigen::Dynamic>& cosSin,
-      double constraintImportance)
-  {
-    // Keep track of the number of corrections to derive
-    std::size_t nCorr = cosSin.cols();
-    if (nCorr == 0)
-      // No need to do anything to derive corrections to nothing
-      return {};
-    // Eigen internally uses somewhat complicated 'expression templates' to hold
-    // the mathematics of a matrix calculation without actually doing the loop.
-    // This allows it to avoid doing more loops than it actually needs to. It
-    // will do this by default unless you force it to evaluate an expression by
-    // binding it to a vector matrix type. Therefore we use 'auto' to store the
-    // individual expression pieces, rather than forcing an early evaluation.
-    auto constants =
-      // Constant part derived from the uniformity constraint
-      towerExpectations.cwiseQuotient(towerVariances) - 
-      // Constant part derived from the pileup balance constraint
-      constraintImportance * cosSin.transpose()*pileupCovariance.inverse()*pileupSum;
-    Eigen::MatrixXd diagonal = towerVariances.cwiseInverse().asDiagonal();
-    auto coeffMatrix = 
-      // Matrix part derived from the uniformity constraint
-      diagonal + 
-      // Matrix part derived from the pileup balance constraint
-      constraintImportance * (cosSin.transpose()*pileupCovariance.inverse()*cosSin);
-    // Now return the actual corrections
-    return coeffMatrix.inverse() * constants;
+      Eigen::VectorXd pufit(
+          const Eigen::Vector2d &pileupSum,
+          const Eigen::Matrix2d &pileupCovariance,
+          const Eigen::VectorXd &towerExpectations,
+          const Eigen::VectorXd &towerVariances,
+          const Eigen::VectorXd &correctionDirections,
+          double constraintImportance)
+      {
+        // Keep track of the number of corrections to derive
+        std::size_t nCorr = correctionDirections.size();
+        if (nCorr == 0)
+          // No need to do anything to derive corrections to nothing
+          return {};
+        // Turn the directions into cos/sin pairs
+        Eigen::Matrix<double, 2, Eigen::Dynamic> cosSin(2, nCorr);
+        cosSin.row(0) = correctionDirections.array().cos();
+        cosSin.row(1) = correctionDirections.array().sin();
+        return pufit(
+            pileupSum,
+            pileupCovariance,
+            towerExpectations,
+            towerVariances,
+            cosSin,
+            constraintImportance);
+      }
 
-  }
+      Eigen::VectorXd pufit(
+          const Eigen::Vector2d &pileupSum,
+          const Eigen::Matrix2d &pileupCovariance,
+          double towerMean,
+          double towerVariance,
+          const Eigen::VectorXd &correctionDirections,
+          double constraintImportance)
+      {
+        std::size_t nCorr = correctionDirections.size();
+        return pufit(
+            pileupSum,
+            pileupCovariance,
+            Eigen::VectorXd::Constant(nCorr, towerMean),
+            Eigen::VectorXd::Constant(nCorr, towerVariance),
+            correctionDirections,
+            constraintImportance);
+      }
 
-  Eigen::VectorXd pufit(
-      const Eigen::Vector2d& pileupSum,
-      const Eigen::Matrix2d& pileupCovariance,
-      double towerMean,
-      double towerVariance,
-      const Eigen::Matrix<double, 2, Eigen::Dynamic>& cosSin,
-      double constraintImportance)
-  {
-    std::size_t nCorr = cosSin.cols();
-    return pufit(
-        pileupSum,
-        pileupCovariance,
-        Eigen::VectorXd::Constant(nCorr, towerMean),
-        Eigen::VectorXd::Constant(nCorr, towerVariance),
-        cosSin,
-        constraintImportance);
-  }
+      Eigen::VectorXd pufit(
+          const Eigen::Vector2d &pileupSum,
+          const Eigen::Matrix2d &pileupCovariance,
+          const Eigen::VectorXd &towerExpectations,
+          const Eigen::VectorXd &towerVariances,
+          const Eigen::Matrix<double, 2, Eigen::Dynamic> &cosSin,
+          double constraintImportance)
+      {
+        // Keep track of the number of corrections to derive
+        std::size_t nCorr = cosSin.cols();
+        if (nCorr == 0)
+          // No need to do anything to derive corrections to nothing
+          return {};
+        // Eigen internally uses somewhat complicated 'expression templates' to hold
+        // the mathematics of a matrix calculation without actually doing the loop.
+        // This allows it to avoid doing more loops than it actually needs to. It
+        // will do this by default unless you force it to evaluate an expression by
+        // binding it to a vector matrix type. Therefore we use 'auto' to store the
+        // individual expression pieces, rather than forcing an early evaluation.
+        auto constants =
+            // Constant part derived from the uniformity constraint
+            towerExpectations.cwiseQuotient(towerVariances) -
+            // Constant part derived from the pileup balance constraint
+            constraintImportance * cosSin.transpose() * pileupCovariance.inverse() * pileupSum;
+        Eigen::MatrixXd diagonal = towerVariances.cwiseInverse().asDiagonal();
+        auto coeffMatrix =
+            // Matrix part derived from the uniformity constraint
+            diagonal +
+            // Matrix part derived from the pileup balance constraint
+            constraintImportance * (cosSin.transpose() * pileupCovariance.inverse() * cosSin);
+        // Now return the actual corrections
+        return coeffMatrix.inverse() * constants;
+      }
 
-  std::vector<SignedKinematics> pufit(
-      const Eigen::Vector2d& pileupSum,
-      const Eigen::Matrix2d& pileupCovariance,
-      const std::vector<double>& towerExpectations,
-      const std::vector<double>& towerVariances,
-      const std::vector<SignedKinematics>& toCorrect,
-      double constraintImportance)
-  {
-    // Keep track of the number of corrections to derive
-    std::size_t nCorr = toCorrect.size();
-    if (nCorr == 0)
-      // No need to do anything to derive corrections to nothing
-      return {};
+      Eigen::VectorXd pufit(
+          const Eigen::Vector2d &pileupSum,
+          const Eigen::Matrix2d &pileupCovariance,
+          double towerMean,
+          double towerVariance,
+          const Eigen::Matrix<double, 2, Eigen::Dynamic> &cosSin,
+          double constraintImportance)
+      {
+        std::size_t nCorr = cosSin.cols();
+        return pufit(
+            pileupSum,
+            pileupCovariance,
+            Eigen::VectorXd::Constant(nCorr, towerMean),
+            Eigen::VectorXd::Constant(nCorr, towerVariance),
+            cosSin,
+            constraintImportance);
+      }
 
-    // Get the cos-sin matrix
-    Eigen::Matrix<double, 2, Eigen::Dynamic> cosSin(2, nCorr);
-    Eigen::VectorXd vecTowerExpectations(nCorr);
-    Eigen::VectorXd vecTowerVariances(nCorr);
-    for (std::size_t ii = 0; ii < nCorr; ++ii) {
-      cosSin(0, ii) = toCorrect.at(ii).cosPhi();
-      cosSin(1, ii) = toCorrect.at(ii).sinPhi();
-      vecTowerExpectations(ii) = towerExpectations.at(ii);
-      vecTowerVariances(ii) = towerVariances.at(ii);
-    }
-    Eigen::VectorXd corrections = pufit(
-        pileupSum,
-        pileupCovariance,
-        vecTowerExpectations,
-        vecTowerVariances,
-        cosSin,
-        constraintImportance);
-    // Now convert these into kinematics
-    std::vector<SignedKinematics> kinCorrections;
-    kinCorrections.reserve(nCorr);
-    for (std::size_t ii = 0; ii < nCorr; ++ii)
-      kinCorrections.push_back(SignedKinematics::fromEtEtaPhi(
-            -corrections.coeff(ii),
-            toCorrect.at(ii).eta(),
-            toCorrect.at(ii).phi() ) );
-    return kinCorrections;
+      std::vector<SignedKinematics> pufit(
+          const Eigen::Vector2d &pileupSum,
+          const Eigen::Matrix2d &pileupCovariance,
+          const std::vector<double> &towerExpectations,
+          const std::vector<double> &towerVariances,
+          const std::vector<SignedKinematics> &toCorrect,
+          double constraintImportance)
+      {
+        // Keep track of the number of corrections to derive
+        std::size_t nCorr = toCorrect.size();
+        if (nCorr == 0)
+          // No need to do anything to derive corrections to nothing
+          return {};
 
-  }
-  std::vector<SignedKinematics> pufit(
-      const Eigen::Vector2d& pileupSum,
-      const Eigen::Matrix2d& pileupCovariance,
-      double towerMean,
-      double towerVariance,
-      const std::vector<SignedKinematics>& toCorrect,
-      double constraintImportance)
-  {
-    // Keep track of the number of corrections to derive
-    std::size_t nCorr = toCorrect.size();
-    if (nCorr == 0)
-      // No need to do anything to derive corrections to nothing
-      return {};
+        // Get the cos-sin matrix
+        Eigen::Matrix<double, 2, Eigen::Dynamic> cosSin(2, nCorr);
+        Eigen::VectorXd vecTowerExpectations(nCorr);
+        Eigen::VectorXd vecTowerVariances(nCorr);
+        for (std::size_t ii = 0; ii < nCorr; ++ii)
+        {
+          cosSin(0, ii) = toCorrect.at(ii).cosPhi();
+          cosSin(1, ii) = toCorrect.at(ii).sinPhi();
+          vecTowerExpectations(ii) = towerExpectations.at(ii);
+          vecTowerVariances(ii) = towerVariances.at(ii);
+        }
+        Eigen::VectorXd corrections = pufit(
+            pileupSum,
+            pileupCovariance,
+            vecTowerExpectations,
+            vecTowerVariances,
+            cosSin,
+            constraintImportance);
+        // Now convert these into kinematics
+        std::vector<SignedKinematics> kinCorrections;
+        kinCorrections.reserve(nCorr);
+        for (std::size_t ii = 0; ii < nCorr; ++ii)
+          kinCorrections.push_back(SignedKinematics::fromEtEtaPhi(
+              -corrections.coeff(ii),
+              toCorrect.at(ii).eta(),
+              toCorrect.at(ii).phi()));
+        return kinCorrections;
+      }
+      std::vector<SignedKinematics> pufit(
+          const Eigen::Vector2d &pileupSum,
+          const Eigen::Matrix2d &pileupCovariance,
+          double towerMean,
+          double towerVariance,
+          const std::vector<SignedKinematics> &toCorrect,
+          double constraintImportance)
+      {
+        // Keep track of the number of corrections to derive
+        std::size_t nCorr = toCorrect.size();
+        if (nCorr == 0)
+          // No need to do anything to derive corrections to nothing
+          return {};
 
-    return pufit(
-        pileupSum,
-        pileupCovariance,
-        std::vector<double>(nCorr, towerMean),
-        std::vector<double>(nCorr, towerVariance),
-        toCorrect,
-        constraintImportance);
-  }
-} } } //> end namespace HLT::MET::PufitUtils
+        return pufit(
+            pileupSum,
+            pileupCovariance,
+            std::vector<double>(nCorr, towerMean),
+            std::vector<double>(nCorr, towerVariance),
+            toCorrect,
+            constraintImportance);
+      }
+    } // namespace PufitUtils
+  }   // namespace MET
+} // namespace HLT
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/TCFex.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/TCFex.cxx
index 564c141e3179..c5f4ae5e221e 100644
--- a/Trigger/TrigAlgorithms/TrigEFMissingET/src/TCFex.cxx
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/TCFex.cxx
@@ -39,6 +39,11 @@ namespace HLT { namespace MET {
   {
     // Retrieve the inputs
     auto clusters = SG::makeHandle(m_clusterKey, context);
+    if (!clusters.isValid())
+    {
+      ATH_MSG_ERROR("Failed to retrieve " << m_clusterKey);
+      return StatusCode::FAILURE;
+    }
 
     auto state = m_useUncalibrated ? 
       xAOD::CaloCluster::UNCALIBRATED : 
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/TCPufitFex.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/TCPufitFex.cxx
index c6c4b1e0e0b7..13c4596ef003 100644
--- a/Trigger/TrigAlgorithms/TrigEFMissingET/src/TCPufitFex.cxx
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/TCPufitFex.cxx
@@ -15,80 +15,93 @@
 #include "TrigEFMissingET/PufitGrid.h"
 #include "TrigEFMissingET/PufitUtils.h"
 
-namespace HLT { namespace MET {
-  TCPufitFex::TCPufitFex(const std::string& name, ISvcLocator* pSvcLocator) :
-    FexBase(name, pSvcLocator)
-  {}
-
-  StatusCode TCPufitFex::initialize()
+namespace HLT
+{
+  namespace MET
   {
-    CHECK(m_clusterKey.initialize() );
-    return initializeBase({"AllTowers", "UncorrSelTowers"});
-  }
+    TCPufitFex::TCPufitFex(const std::string &name, ISvcLocator *pSvcLocator) : FexBase(name, pSvcLocator)
+    {
+    }
 
-  StatusCode TCPufitFex::fillMET(
-      xAOD::TrigMissingET& met,
-      const EventContext& context,
-      MonGroupBuilder&) const
-  {
-    // Retrieve the inputs
-    auto clusters = SG::makeHandle(m_clusterKey, context);
+    StatusCode TCPufitFex::initialize()
+    {
+      CHECK(m_clusterKey.initialize());
+      return initializeBase({"AllTowers", "UncorrSelTowers"});
+    }
 
-    // Create the pufit grids
-    PufitGridSet gridset(m_maxEta, m_nEtaBins, m_nPhiBins);
+    StatusCode TCPufitFex::fillMET(
+        xAOD::TrigMissingET &met,
+        const EventContext &context,
+        MonGroupBuilder &) const
+    {
+      // Retrieve the inputs
+      auto clusters = SG::makeHandle(m_clusterKey, context);
+      if (!clusters.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_clusterKey);
+        return StatusCode::FAILURE;
+      }
+
+      // Create the pufit grids
+      PufitGridSet gridset(m_maxEta, m_nEtaBins, m_nPhiBins);
 
-    // Fill the grids from the input clusters
-    for (const xAOD::CaloCluster* iclus : *clusters)
-      gridset += SignedKinematics::fromEnergyEtaPhi(
-          iclus->e(), iclus->eta(), iclus->phi() );
+      // Fill the grids from the input clusters
+      for (const xAOD::CaloCluster *iclus : *clusters)
+        gridset += SignedKinematics::fromEnergyEtaPhi(
+            iclus->e(), iclus->eta(), iclus->phi());
 
-    // We calculate the mean and variance from the 'nominal' grid
-    std::pair<double, double> mav = PufitUtils::trimmedMeanAndVariance(
-        gridset[PufitGridSet::NoDisplacement], m_trimFraction);
-    // Calculate the threshold
-    float threshold = mav.first + m_nSigma * sqrt(mav.second);
+      // We calculate the mean and variance from the 'nominal' grid
+      double mean;
+      double variance;
+      PufitUtils::trimmedMeanAndVariance(
+          gridset[NoDisplacement], m_trimFraction, mean, variance);
+      // Calculate the threshold
+      float threshold = mean + m_nSigma * sqrt(variance);
 
-    // Apply the masks
-    for (PufitGrid& grid : gridset.grids)
-      for (PufitGrid::Tower& tower : grid)
-        tower.mask(tower.sumEt() > threshold);
+      // Apply the masks
+      for (PufitGrid &grid : gridset.grids)
+        for (PufitGrid::Tower &tower : grid)
+          tower.mask(tower.sumEt() > threshold);
 
-    // Select the right grid
-    PufitGridSet::GridDisplacement idx = PufitUtils::selectGrid(gridset);
-    const PufitGrid& grid = gridset[idx];
+      // Select the right grid
+      GridDisplacement idx = PufitUtils::selectGrid(gridset);
+      const PufitGrid &grid = gridset[idx];
 
-    // Build the pileup sum and prepare the vector of corrections
-    PufitUtils::CovarianceSum pileupSum;
-    std::vector<SignedKinematics> masked;
-    for (const PufitGrid::Tower& tower : grid) {
-      if (!tower.masked() ) {
-        float sigma =
-          m_caloNoise*m_caloNoise +
-          tower.kinematics().absPt()*m_caloStoch*m_caloStoch;
-        pileupSum.add(tower, sigma);
+      // Build the pileup sum and prepare the vector of corrections
+      PufitUtils::CovarianceSum pileupSum;
+      std::vector<SignedKinematics> masked;
+      for (const PufitGrid::Tower &tower : grid)
+      {
+        if (!tower.masked())
+        {
+          float sigma =
+              m_caloNoise * m_caloNoise +
+              tower.kinematics().absPt() * m_caloStoch * m_caloStoch;
+          pileupSum.add(tower, sigma);
+        }
+        else
+          masked.push_back(tower);
       }
-      else
-        masked.push_back(tower);
-    }
-    // Perform the fit
-    std::vector<SignedKinematics> corrections = PufitUtils::pufit(
-        pileupSum.sum,
-        pileupSum.covariance,
-        mav.first,
-        mav.second,
-        masked,
-        m_constraintImportance);
+      // Perform the fit
+      std::vector<SignedKinematics> corrections = PufitUtils::pufit(
+          pileupSum.sum,
+          pileupSum.covariance,
+          mean,
+          variance,
+          masked,
+          m_constraintImportance);
 
-    // Save the sum over all towers to the corresponding component
-    grid.sum(PufitGrid::SumStrategy::All).fillMETComponent(0, met);
-    // Get the sum over the selected towers
-    METComponent sum = grid.sum(PufitGrid::SumStrategy::Masked);
-    // Save this uncorrected sum to the corresponding component
-    sum.fillMETComponent(1, met);
-    // Apply the corrections
-    for (const SignedKinematics& kin : corrections)
-      sum += kin;
-    sum.fillMET(met);
-    return StatusCode::SUCCESS;
-  }
-} } //> end namespace HLT::MET
+      // Save the sum over all towers to the corresponding component
+      grid.sum(PufitGrid::SumStrategy::All).fillMETComponent(0, met);
+      // Get the sum over the selected towers
+      METComponent sum = grid.sum(PufitGrid::SumStrategy::Masked);
+      // Save this uncorrected sum to the corresponding component
+      sum.fillMETComponent(1, met);
+      // Apply the corrections
+      for (const SignedKinematics &kin : corrections)
+        sum += kin;
+      sum.fillMET(met);
+      return StatusCode::SUCCESS;
+    }
+  } // namespace MET
+} // namespace HLT
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/TrkMHTFex.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/TrkMHTFex.cxx
index c55a4d2a19a2..e19b1d6ea92e 100644
--- a/Trigger/TrigAlgorithms/TrigEFMissingET/src/TrkMHTFex.cxx
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/TrkMHTFex.cxx
@@ -39,14 +39,14 @@ namespace HLT { namespace MET {
     CHECK( m_trackSelTool.retrieve() );
 
     // Update the decor keys if necessary
-    if (m_jvtKey.key().find(".") == std::string::npos) {
+      if (m_jvtKey.key().find(".") == std::string::npos)
       m_jvtKey = m_jetKey.key() + "." + m_jvtKey.key();
-      CHECK( m_jvtKey.initialize() );
-    }
-    else if (SG::contKeyFromKey(m_jvtKey.key() ) != m_jetKey.key() ) {
+      else if (SG::contKeyFromKey(m_jvtKey.key()) != m_jetKey.key())
+      {
       ATH_MSG_ERROR("JVT container key does not match jet key!");
       return StatusCode::FAILURE;
     }
+      CHECK(m_jvtKey.initialize());
     m_trackGA.emplace(m_trackGAName);
 
     return initializeBase(
@@ -60,8 +60,23 @@ namespace HLT { namespace MET {
   {
     // Retrieve the inputs
     auto jets = SG::makeHandle(m_jetKey, context);
+      if (!jets.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_jetKey);
+        return StatusCode::FAILURE;
+      }
     auto tracks = SG::makeHandle(m_trackKey, context);
+      if (!tracks.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_trackKey);
+        return StatusCode::FAILURE;
+      }
     auto vertices = SG::makeHandle(m_vertexKey, context);
+      if (!vertices.isValid())
+      {
+        ATH_MSG_ERROR("Failed to retrieve " << m_vertexKey);
+        return StatusCode::FAILURE;
+      }
     auto tva = SG::makeHandle(m_tvaKey, context);
     auto jvtAcc = SG::makeHandle<float>(m_jvtKey, context);
     auto& trackLinksAcc = *m_trackGA;
diff --git a/Trigger/TrigAlgorithms/TrigEFMissingET/src/components/TrigEFMissingET_entries.cxx b/Trigger/TrigAlgorithms/TrigEFMissingET/src/components/TrigEFMissingET_entries.cxx
index 96a1ea112c8c..54b92397e67f 100644
--- a/Trigger/TrigAlgorithms/TrigEFMissingET/src/components/TrigEFMissingET_entries.cxx
+++ b/Trigger/TrigAlgorithms/TrigEFMissingET/src/components/TrigEFMissingET_entries.cxx
@@ -18,23 +18,37 @@
 #include "../TCFex.h"
 #include "../TCPufitFex.h"
 #include "../PFSumFex.h"
+#include "../PUSplitPufitFex.h"
+#include "../CVFAlg.h"
+#include "../CVFPrepAlg.h"
+#include "../PFOPrepAlg.h"
+#include "TrigEFMissingET/ApproximateTrackToLayerTool.h"
+#include "TrigEFMissingET/ExtendTrackToLayerTool.h"
+#include "../MHTPufitFex.h"
 
-DECLARE_COMPONENT( EFMissingET )
-DECLARE_COMPONENT( EFMissingETBaseTool )
-DECLARE_COMPONENT( EFMissingETFromCells )
-DECLARE_COMPONENT( EFMissingETFromClusters )
-DECLARE_COMPONENT( EFMissingETFromClustersPS )
-DECLARE_COMPONENT( EFMissingETFromClustersPUC )
-DECLARE_COMPONENT( EFMissingETFromFEBHeader )
-DECLARE_COMPONENT( EFMissingETFromJets )
-DECLARE_COMPONENT( EFMissingETFromTrackAndJets )
-DECLARE_COMPONENT( EFMissingETFromClustersTracksPUC )
-DECLARE_COMPONENT( EFMissingETFromTrackAndClusters )
-DECLARE_COMPONENT( EFMissingETFlags )
-DECLARE_COMPONENT( EFMissingETFromHelper )
-DECLARE_COMPONENT( HLT::MET::TrkMHTFex )
-DECLARE_COMPONENT( HLT::MET::CellFex )
-DECLARE_COMPONENT( HLT::MET::MHTFex )
-DECLARE_COMPONENT( HLT::MET::TCFex )
-DECLARE_COMPONENT( HLT::MET::TCPufitFex )
-DECLARE_COMPONENT( HLT::MET::PFSumFex )
+DECLARE_COMPONENT(EFMissingET)
+DECLARE_COMPONENT(EFMissingETBaseTool)
+DECLARE_COMPONENT(EFMissingETFromCells)
+DECLARE_COMPONENT(EFMissingETFromClusters)
+DECLARE_COMPONENT(EFMissingETFromClustersPS)
+DECLARE_COMPONENT(EFMissingETFromClustersPUC)
+DECLARE_COMPONENT(EFMissingETFromFEBHeader)
+DECLARE_COMPONENT(EFMissingETFromJets)
+DECLARE_COMPONENT(EFMissingETFromTrackAndJets)
+DECLARE_COMPONENT(EFMissingETFromClustersTracksPUC)
+DECLARE_COMPONENT(EFMissingETFromTrackAndClusters)
+DECLARE_COMPONENT(EFMissingETFlags)
+DECLARE_COMPONENT(EFMissingETFromHelper)
+DECLARE_COMPONENT(HLT::MET::TrkMHTFex)
+DECLARE_COMPONENT(HLT::MET::CellFex)
+DECLARE_COMPONENT(HLT::MET::MHTFex)
+DECLARE_COMPONENT(HLT::MET::TCFex)
+DECLARE_COMPONENT(HLT::MET::TCPufitFex)
+DECLARE_COMPONENT(HLT::MET::PFSumFex)
+DECLARE_COMPONENT(HLT::MET::PUSplitPufitFex)
+DECLARE_COMPONENT(HLT::MET::CVFAlg)
+DECLARE_COMPONENT(HLT::MET::CVFPrepAlg)
+DECLARE_COMPONENT(HLT::MET::PFOPrepAlg)
+DECLARE_COMPONENT(ApproximateTrackToLayerTool)
+DECLARE_COMPONENT(ExtendTrackToLayerTool)
+DECLARE_COMPONENT(HLT::MET::MHTPufitFex)
\ No newline at end of file
diff --git a/Trigger/TrigValidation/TrigAnalysisTest/share/ref_RDOtoRDOTrig_mt1_build.ref b/Trigger/TrigValidation/TrigAnalysisTest/share/ref_RDOtoRDOTrig_mt1_build.ref
index d357d2034bba..18b41e35da0d 100644
--- a/Trigger/TrigValidation/TrigAnalysisTest/share/ref_RDOtoRDOTrig_mt1_build.ref
+++ b/Trigger/TrigValidation/TrigAnalysisTest/share/ref_RDOtoRDOTrig_mt1_build.ref
@@ -772,9 +772,21 @@ TrigSignatureMoniMT                                 INFO -- #1649696554 Features
 TrigSignatureMoniMT                                 INFO HLT_xe30_cell_xe30_tcpufit_L1XE10 #3768353779
 TrigSignatureMoniMT                                 INFO -- #3768353779 Events         19         19         0          0          0          0          0          0          0          0          0          0          14         -          14
 TrigSignatureMoniMT                                 INFO -- #3768353779 Features                             0          0          0          0          0          0          0          0          0          0          14         -
+TrigSignatureMoniMT                                 INFO HLT_xe30_cvfpufit_L1XE10 #3860749499
+TrigSignatureMoniMT                                 INFO -- #3860749499 Events         19         19         0          0          0          0          0          0          0          0          0          0          14         -          14
+TrigSignatureMoniMT                                 INFO -- #3860749499 Features                             0          0          0          0          0          0          0          0          0          0          14         -
 TrigSignatureMoniMT                                 INFO HLT_xe30_mht_L1XE10 #3626903018
 TrigSignatureMoniMT                                 INFO -- #3626903018 Events         19         19         0          0          0          0          0          0          0          0          0          0          19         -          19
 TrigSignatureMoniMT                                 INFO -- #3626903018 Features                             0          0          0          0          0          0          0          0          0          0          19         -
+TrigSignatureMoniMT                                 INFO HLT_xe30_mhtpufit_em_subjesgscIS_L1XE10 #689201557
+TrigSignatureMoniMT                                 INFO -- #689201557 Events          19         19         0          0          0          0          0          0          0          0          0          0          15         -          15
+TrigSignatureMoniMT                                 INFO -- #689201557 Features                              0          0          0          0          0          0          0          0          0          0          15         -
+TrigSignatureMoniMT                                 INFO HLT_xe30_mhtpufit_pf_subjesgscIS_L1XE10 #1886909707
+TrigSignatureMoniMT                                 INFO -- #1886909707 Events         19         19         0          0          0          0          0          0          0          0          0          0          14         -          14
+TrigSignatureMoniMT                                 INFO -- #1886909707 Features                             0          0          0          0          0          0          0          0          0          0          14         -
+TrigSignatureMoniMT                                 INFO HLT_xe30_pfopufit_L1XE10 #2252641537
+TrigSignatureMoniMT                                 INFO -- #2252641537 Events         19         19         0          0          0          0          0          0          0          0          0          0          14         -          14
+TrigSignatureMoniMT                                 INFO -- #2252641537 Features                             0          0          0          0          0          0          0          0          0          0          14         -
 TrigSignatureMoniMT                                 INFO HLT_xe30_pfsum_L1XE10 #998713382
 TrigSignatureMoniMT                                 INFO -- #998713382 Events          19         19         0          0          0          0          0          0          0          0          0          0          16         -          16
 TrigSignatureMoniMT                                 INFO -- #998713382 Features                              0          0          0          0          0          0          0          0          0          0          16         -
diff --git a/Trigger/TrigValidation/TriggerTest/share/ref_data_v1Dev_build.ref b/Trigger/TrigValidation/TriggerTest/share/ref_data_v1Dev_build.ref
index 9629f6ef14ad..5f3fabea7793 100644
--- a/Trigger/TrigValidation/TriggerTest/share/ref_data_v1Dev_build.ref
+++ b/Trigger/TrigValidation/TriggerTest/share/ref_data_v1Dev_build.ref
@@ -772,9 +772,21 @@ TrigSignatureMoniMT                                 INFO -- #1649696554 Features
 TrigSignatureMoniMT                                 INFO HLT_xe30_cell_xe30_tcpufit_L1XE10 #3768353779
 TrigSignatureMoniMT                                 INFO -- #3768353779 Events         20         20         0          0          0          0          0          0          0          0          0          0          3          -          3          
 TrigSignatureMoniMT                                 INFO -- #3768353779 Features                             0          0          0          0          0          0          0          0          0          0          3          -          
+TrigSignatureMoniMT                                 INFO HLT_xe30_cvfpufit_L1XE10 #3860749499
+TrigSignatureMoniMT                                 INFO -- #3860749499 Events         20         20         0          0          0          0          0          0          0          0          0          0          6          -          6          
+TrigSignatureMoniMT                                 INFO -- #3860749499 Features                             0          0          0          0          0          0          0          0          0          0          6          -          
 TrigSignatureMoniMT                                 INFO HLT_xe30_mht_L1XE10 #3626903018
 TrigSignatureMoniMT                                 INFO -- #3626903018 Events         20         20         0          0          0          0          0          0          0          0          0          0          14         -          14         
 TrigSignatureMoniMT                                 INFO -- #3626903018 Features                             0          0          0          0          0          0          0          0          0          0          14         -          
+TrigSignatureMoniMT                                 INFO HLT_xe30_mhtpufit_em_subjesgscIS_L1XE10 #689201557
+TrigSignatureMoniMT                                 INFO -- #689201557 Events          20         20         0          0          0          0          0          0          0          0          0          0          4          -          4          
+TrigSignatureMoniMT                                 INFO -- #689201557 Features                              0          0          0          0          0          0          0          0          0          0          4          -          
+TrigSignatureMoniMT                                 INFO HLT_xe30_mhtpufit_pf_subjesgscIS_L1XE10 #1886909707
+TrigSignatureMoniMT                                 INFO -- #1886909707 Events         20         20         0          0          0          0          0          0          0          0          0          0          3          -          3          
+TrigSignatureMoniMT                                 INFO -- #1886909707 Features                             0          0          0          0          0          0          0          0          0          0          3          -          
+TrigSignatureMoniMT                                 INFO HLT_xe30_pfopufit_L1XE10 #2252641537
+TrigSignatureMoniMT                                 INFO -- #2252641537 Events         20         20         0          0          0          0          0          0          0          0          0          0          4          -          4          
+TrigSignatureMoniMT                                 INFO -- #2252641537 Features                             0          0          0          0          0          0          0          0          0          0          4          -          
 TrigSignatureMoniMT                                 INFO HLT_xe30_pfsum_L1XE10 #998713382
 TrigSignatureMoniMT                                 INFO -- #998713382 Events          20         20         0          0          0          0          0          0          0          0          0          0          5          -          5          
 TrigSignatureMoniMT                                 INFO -- #998713382 Features                              0          0          0          0          0          0          0          0          0          0          5          -          
diff --git a/Trigger/TriggerCommon/TrigEDMConfig/python/TriggerEDMRun3.py b/Trigger/TriggerCommon/TrigEDMConfig/python/TriggerEDMRun3.py
index 6e9164536f59..f38de5d4a1f1 100644
--- a/Trigger/TriggerCommon/TrigEDMConfig/python/TriggerEDMRun3.py
+++ b/Trigger/TriggerCommon/TrigEDMConfig/python/TriggerEDMRun3.py
@@ -351,6 +351,18 @@ TriggerHLTListRun3 = [
     ('xAOD::TrigMissingETContainer#HLT_MET_pfsum',                         'BS ESD AODFULL AODSLIM AODVERYSLIM', 'MET'),
     ('xAOD::TrigMissingETAuxContainer#HLT_MET_pfsumAux.',                  'BS ESD AODFULL AODSLIM AODVERYSLIM', 'MET'),
 
+    ('xAOD::TrigMissingETContainer#HLT_MET_pfopufit',                      'BS ESD AODFULL AODSLIM AODVERYSLIM', 'MET'),
+    ('xAOD::TrigMissingETAuxContainer#HLT_MET_pfopufitAux.',               'BS ESD AODFULL AODSLIM AODVERYSLIM', 'MET'),
+
+    ('xAOD::TrigMissingETContainer#HLT_MET_cvfpufit',                      'BS ESD AODFULL AODSLIM AODVERYSLIM', 'MET'),
+    ('xAOD::TrigMissingETAuxContainer#HLT_MET_cvfpufitAux.',               'BS ESD AODFULL AODSLIM AODVERYSLIM', 'MET'),
+
+    ('xAOD::TrigMissingETContainer#HLT_MET_mhtpufit_pf_subjesgscIS',       'BS ESD AODFULL AODSLIM AODVERYSLIM', 'MET'),
+    ('xAOD::TrigMissingETAuxContainer#HLT_MET_mhtpufit_pf_subjesgscISAux.', 'BS ESD AODFULL AODSLIM AODVERYSLIM', 'MET'),
+
+    ('xAOD::TrigMissingETContainer#HLT_MET_mhtpufit_em_subjesgscIS',       'BS ESD AODFULL AODSLIM AODVERYSLIM', 'MET'),
+    ('xAOD::TrigMissingETAuxContainer#HLT_MET_mhtpufit_em_subjesgscISAux.', 'BS ESD AODFULL AODSLIM AODVERYSLIM', 'MET'),
+
     ('xAOD::CaloClusterContainer#HLT_TopoCaloClustersFS',                  'BS ESD AODFULL AODSLIM AODVERYSLIM', 'MET'),
     ('xAOD::CaloClusterTrigAuxContainer#HLT_TopoCaloClustersFSAux.nCells', 'BS ESD AODFULL AODSLIM AODVERYSLIM', 'MET'),
 
diff --git a/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/MET/AlgConfigs.py b/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/MET/AlgConfigs.py
index e5653ddbab2b..0240e8e15b2b 100644
--- a/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/MET/AlgConfigs.py
+++ b/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/MET/AlgConfigs.py
@@ -6,10 +6,13 @@ from .ConfigHelpers import AlgConfig, jetRecoDictForMET
 from ..Menu.MenuComponents import RecoFragmentsPool
 from ..Menu.SignatureDicts import METChainParts
 import GaudiKernel.SystemOfUnits as Units
+import TrigEFMissingET.PUClassification as PUClassification
 
 from AthenaCommon.Logging import logging
+
 log = logging.getLogger(__name__)
 
+
 def test_configs():
     """ Make sure that all algorithms defined in the METChainParts have
     configurations
@@ -25,26 +28,30 @@ def test_configs():
         else:
             unknown_algs.append(alg)
     assert len(unknown_algs) == 0, (
-             "The following EFrecoAlgs do not have AlgConfig classes: "
-             "{}".format(unknown_algs) )
-    
+        "The following EFrecoAlgs do not have AlgConfig classes: "
+        "{}".format(unknown_algs)
+    )
+
 
 class CellConfig(AlgConfig):
     @classmethod
     def algType(cls):
         return "cell"
-    
+
     def __init__(self, **recoDict):
         super(CellConfig, self).__init__(**recoDict)
-        from TriggerMenuMT.HLTMenuConfig.CommonSequences.CaloSequenceSetup import cellRecoSequence
+        from TriggerMenuMT.HLTMenuConfig.CommonSequences.CaloSequenceSetup import (
+            cellRecoSequence,
+        )
         from TrigEFMissingET.TrigEFMissingETConf import HLT__MET__CellFex
+
         cellMakerSeq, cellName = RecoFragmentsPool.retrieve(
-                cellRecoSequence, flags=None, RoIs=self.inputMaker.RoIs)
+            cellRecoSequence, flags=None, RoIs=self.inputMaker.RoIs
+        )
 
         self.inputs = [cellMakerSeq]
-        self.fexAlg = self._make_fex_alg(
-                HLT__MET__CellFex,
-                CellName = cellName)
+        self.fexAlg = self._make_fex_alg(HLT__MET__CellFex, CellName=cellName)
+
 
 class TCConfig(AlgConfig):
     @classmethod
@@ -54,20 +61,24 @@ class TCConfig(AlgConfig):
     def __init__(self, calib, **recoDict):
         super(TCConfig, self).__init__(calib=calib, **recoDict)
         from TriggerMenuMT.HLTMenuConfig.CommonSequences.CaloSequenceSetup import (
-                caloClusterRecoSequence, LCCaloClusterRecoSequence)
+            caloClusterRecoSequence,
+            LCCaloClusterRecoSequence,
+        )
         from TrigEFMissingET.TrigEFMissingETConf import HLT__MET__TCFex
+
         RoIs = self.inputMaker.RoIs
-        if calib == "em": 
+        if calib == "em":
             tcSeq, clusterName = RecoFragmentsPool.retrieve(
-                    caloClusterRecoSequence, flags = None, RoIs=RoIs)
+                caloClusterRecoSequence, flags=None, RoIs=RoIs
+            )
         elif calib == "lcw":
             tcSeq, clusterName = RecoFragmentsPool.retrieve(
-                    LCCaloClusterRecoSequence, flag = None, RoIs=RoIs)
+                LCCaloClusterRecoSequence, flag=None, RoIs=RoIs
+            )
 
         self.inputs = [tcSeq]
-        self.fexAlg = self._make_fex_alg(
-                HLT__MET__TCFex,
-                ClusterName = clusterName)
+        self.fexAlg = self._make_fex_alg(HLT__MET__TCFex, ClusterName=clusterName)
+
 
 class TCPufitConfig(AlgConfig):
     @classmethod
@@ -77,20 +88,24 @@ class TCPufitConfig(AlgConfig):
     def __init__(self, calib, **recoDict):
         super(TCPufitConfig, self).__init__(calib=calib, **recoDict)
         from TriggerMenuMT.HLTMenuConfig.CommonSequences.CaloSequenceSetup import (
-                caloClusterRecoSequence, LCCaloClusterRecoSequence)
+            caloClusterRecoSequence,
+            LCCaloClusterRecoSequence,
+        )
         from TrigEFMissingET.TrigEFMissingETConf import HLT__MET__TCPufitFex
+
         RoIs = self.inputMaker.RoIs
-        if calib == "em": 
+        if calib == "em":
             tcSeq, clusterName = RecoFragmentsPool.retrieve(
-                    caloClusterRecoSequence, flags=None, RoIs=RoIs)
+                caloClusterRecoSequence, flags=None, RoIs=RoIs
+            )
         elif calib == "lcw":
             tcSeq, clusterName = RecoFragmentsPool.retrieve(
-                    LCCaloClusterRecoSequence, flags=None, RoIs=RoIs)
+                LCCaloClusterRecoSequence, flags=None, RoIs=RoIs
+            )
 
         self.inputs = [tcSeq]
-        self.fexAlg = self._make_fex_alg(
-                HLT__MET__TCPufitFex,
-                ClusterName = clusterName)
+        self.fexAlg = self._make_fex_alg(HLT__MET__TCPufitFex, ClusterName=clusterName)
+
 
 class MHTConfig(AlgConfig):
     @classmethod
@@ -101,17 +116,18 @@ class MHTConfig(AlgConfig):
         super(MHTConfig, self).__init__(**recoDict)
         from ..Jet.JetRecoSequences import jetRecoSequence
         from TrigEFMissingET.TrigEFMissingETConf import HLT__MET__MHTFex
+
         jetRecoDict = jetRecoDictForMET(**recoDict)
         # TODO - right now jet calibration is hardcoded to EM
         jetRecoDict["calib"] = "em"
         jetRecoDict["jetCalib"] = "subjes"
         jetSeq, jetName = RecoFragmentsPool.retrieve(
-                jetRecoSequence, None, dataSource="data", **jetRecoDict)
+            jetRecoSequence, None, dataSource="data", **jetRecoDict
+        )
 
         self.inputs = [jetSeq]
-        self.fexAlg = self._make_fex_alg(
-                HLT__MET__MHTFex,
-                JetName = jetName)
+        self.fexAlg = self._make_fex_alg(HLT__MET__MHTFex, JetName=jetName)
+
 
 class TrkMHTConfig(AlgConfig):
     @classmethod
@@ -122,11 +138,13 @@ class TrkMHTConfig(AlgConfig):
         super(TrkMHTConfig, self).__init__(**recoDict)
         from ..Jet.JetRecoSequences import jetRecoSequence
         from TrigEFMissingET.TrigEFMissingETConf import HLT__MET__TrkMHTFex
+
         jetRecoDict = jetRecoDictForMET(trkopt="ftf", **recoDict)
         # TODO - right now jet calibration is hardcoded to EM
         jetRecoDict["calib"] = "em"
         jetSeq, jetName = RecoFragmentsPool.retrieve(
-                jetRecoSequence, None, dataSource="data", **jetRecoDict)
+            jetRecoSequence, None, dataSource="data", **jetRecoDict
+        )
 
         # These are the names set by the upstream algorithms. Unfortunately
         # these aren't passed to us - we just have to 'know' them
@@ -137,16 +155,18 @@ class TrkMHTConfig(AlgConfig):
 
         self.inputs = [jetSeq]
         self.fexAlg = self._make_fex_alg(
-                HLT__MET__TrkMHTFex,
-                JetName       = jetName,
-                TrackName     = tracks,
-                VertexName    = vertices,
-                TVAName       = tva,
-                TrackLinkName = track_links)
-        self.fexAlg.TrackSelTool.CutLevel         = "Loose"
-        self.fexAlg.TrackSelTool.maxZ0SinTheta    = 1.5
+            HLT__MET__TrkMHTFex,
+            JetName=jetName,
+            TrackName=tracks,
+            VertexName=vertices,
+            TVAName=tva,
+            TrackLinkName=track_links,
+        )
+        self.fexAlg.TrackSelTool.CutLevel = "Loose"
+        self.fexAlg.TrackSelTool.maxZ0SinTheta = 1.5
         self.fexAlg.TrackSelTool.maxD0overSigmaD0 = 3
-        self.fexAlg.TrackSelTool.minPt            = 1 * Units.GeV
+        self.fexAlg.TrackSelTool.minPt = 1 * Units.GeV
+
 
 class PFSumConfig(AlgConfig):
     @classmethod
@@ -156,28 +176,155 @@ class PFSumConfig(AlgConfig):
     def __init__(self, **recoDict):
         super(PFSumConfig, self).__init__(**recoDict)
 
-        from TriggerMenuMT.HLTMenuConfig.CommonSequences.CaloSequenceSetup import (
-                caloClusterRecoSequence)
-        from eflowRec.PFHLTSequence import PFHLTSequence
-        from ..Jet.JetRecoConfiguration import defineJetConstit
         from TrigEFMissingET.TrigEFMissingETConf import HLT__MET__PFSumFex
-        jetRecoDict = jetRecoDictForMET(trkopt="ftf", **recoDict)
-        jetRecoDict["calib"] = "em"
-        jetRecoDict["dataType"] = "pf"
+        from .METRecoSequences import pfoRecoSequence
+
+        self.inputs, pfoPrefix = RecoFragmentsPool.retrieve(
+            pfoRecoSequence, flags=None, RoIs=self.inputMaker.RoIs
+        )
+        self.fexAlg = self._make_fex_alg(
+            HLT__MET__PFSumFex,
+            NeutralPFOName=pfoPrefix + "CHSNeutralParticleFlowObjects",
+            ChargedPFOName=pfoPrefix + "CHSChargedParticleFlowObjects",
+        )
+
+
+class PFOPufitConfig(AlgConfig):
+    @classmethod
+    def algType(cls):
+        return "pfopufit"
+
+    def __init__(self, **recoDict):
+        super(PFOPufitConfig, self).__init__(**recoDict)
+
+        from TrigEFMissingET.TrigEFMissingETConf import (
+            HLT__MET__PUSplitPufitFex,
+            HLT__MET__PFOPrepAlg,
+        )
+        from .METRecoSequences import pfoRecoSequence
+
+        pfoInputs, pfoPrefix = RecoFragmentsPool.retrieve(
+            pfoRecoSequence, flags=None, RoIs=self.inputMaker.RoIs
+        )
+        # NB for this, we might be slightly misusing the 'flags' parameter in
+        # the reco fragments pool. Here, we let it just pass the name parameter
+        # through to the underlying alg config class parameter
+        prepAlg = RecoFragmentsPool.retrieve(
+            HLT__MET__PFOPrepAlg,
+            f"{pfoPrefix}PFOPufitPrepAlg",
+            InputNeutralKey=f"{pfoPrefix}CHSNeutralParticleFlowObjects",
+            InputChargedKey=f"{pfoPrefix}CHSChargedParticleFlowObjects",
+            OutputKey=f"{pfoPrefix}METTrigCombinedParticleFlowObjects",
+            OutputCategoryKey="PUClassification",
+        )
+        self.inputs = pfoInputs + [prepAlg]
+        # TODO - make the neutral threshold mode settable in the chain name?
+        self.fexAlg = self._make_fex_alg(
+            HLT__MET__PUSplitPufitFex,
+            InputName=prepAlg.OutputKey,
+            InputCategoryName=prepAlg.OutputCategoryKey,
+            NeutralThresholdMode=PUClassification.NeutralForward,
+        )
+
+
+class CVFPufitConfig(AlgConfig):
+    @classmethod
+    def algType(cls):
+        return "cvfpufit"
+
+    def __init__(self, **recoDict):
+        super().__init__(**recoDict)
+        from .METRecoSequences import cvfClusterSequence
+        from TrigEFMissingET.TrigEFMissingETConf import (
+            HLT__MET__CVFPrepAlg,
+            HLT__MET__PUSplitPufitFex,
+        )
 
         RoIs = self.inputMaker.RoIs
-        tcSeq, clusterName = RecoFragmentsPool.retrieve(
-                caloClusterRecoSequence, flags=None, RoIs=RoIs)
-        pfseq, pfoPrefix = RecoFragmentsPool.retrieve(
-                PFHLTSequence, None, clustersin = clusterName, tracktype="ftf")
-        constit = defineJetConstit(jetRecoDict, pfoPrefix=pfoPrefix)
-        from JetRecConfig.ConstModHelpers import getConstitModAlg
-        constit_mod_seq = getConstitModAlg(
-                constit, "HLT",
-                tvaKey="JetTrackVtxAssoc_{trkopt}".format(**jetRecoDict),
-                vtxKey="HLT_IDVertex_FS")
-        self.inputs = [tcSeq, pfseq, constit_mod_seq]
+        calib = recoDict["calib"]
+        inputs, clusterName, cvfName = RecoFragmentsPool.retrieve(
+            cvfClusterSequence, flags=None, RoIs=RoIs, **recoDict
+        )
+
+        prepAlg = RecoFragmentsPool.retrieve(
+            HLT__MET__CVFPrepAlg,
+            f"{calib}ClusterCVFPrepAlg",
+            InputClusterKey=clusterName,
+            InputCVFKey=cvfName,
+            OutputCategoryKey="PUClassification",
+        )
+
+        self.inputs = inputs + [prepAlg]
+        self.fexAlg = self._make_fex_alg(
+            HLT__MET__PUSplitPufitFex,
+            InputName=clusterName,
+            InputCategoryName=prepAlg.OutputCategoryKey,
+            NeutralThresholdMode=PUClassification.All,
+        )
+
+
+class MHTPufitConfig(AlgConfig):
+    @classmethod
+    def algType(cls):
+        return "mhtpufit"
+
+    def __init__(self, **recoDict):
+        super().__init__(**recoDict)
+        from ..Jet.JetRecoSequences import jetRecoSequence
+        from ..Jet.JetRecoConfiguration import defineJets
+        from TriggerMenuMT.HLTMenuConfig.CommonSequences.CaloSequenceSetup import (
+            caloClusterRecoSequence,
+        )
+        from TrigEFMissingET.TrigEFMissingETConf import HLT__MET__MHTPufitFex
+
+        jetRecoDict = jetRecoDictForMET(trkopt="ftf", **recoDict)
+        # If this is PFlow then set the calib type to "em"
+        if recoDict["jetDataType"] == "pf":
+            jetRecoDict["calib"] = "em"
+        jetSeq, jetName = RecoFragmentsPool.retrieve(
+            jetRecoSequence, flags=None, dataSource="data", **jetRecoDict
+        )
+
+        # We need to get the input name that the jet sequence used
+        _, clusterName = RecoFragmentsPool.retrieve(
+            caloClusterRecoSequence, flags=None, RoIs=self.inputMaker.RoIs
+        )
+        if jetRecoDict["dataType"] == "pf":
+            from eflowRec.PFHLTSequence import PFHLTSequence
+
+            _, pfoPrefix = RecoFragmentsPool.retrieve(
+                PFHLTSequence,
+                flags=None,
+                clustersin=clusterName,
+                tracktype=jetRecoDict["trkopt"],
+            )
+            jetDef = defineJets(jetRecoDict, pfoPrefix=pfoPrefix)
+        elif jetRecoDict["dataType"] == "tc":
+            jetDef = defineJets(jetRecoDict, clustersKey=clusterName)
+        else:
+            raise ValueError(
+                "Unexpected jetDataType {}".format(jetRecoDict["dataType"])
+            )
+        inputName = jetDef.inputdef.inputname
+        calibHasAreaSub = "sub" in jetRecoDict["jetCalib"]
+        if calibHasAreaSub:
+            from JetRecConfig.JetRecConfig import getEventShapeAlg, getConstitPJGAlg
+
+            evtShapeAlg = getEventShapeAlg(
+                jetDef.inputdef,
+                getConstitPJGAlg(jetDef.inputdef).OutputContainer,
+                "HLT_",
+            )
+            rhoKey = evtShapeAlg.EventDensityTool.OutputContainer
+        else:
+            rhoKey = ""
+
+        self.inputs = [jetSeq]
         self.fexAlg = self._make_fex_alg(
-                HLT__MET__PFSumFex,
-                NeutralPFOName = pfoPrefix+"CHSNeutralParticleFlowObjects",
-                ChargedPFOName = pfoPrefix+"CHSChargedParticleFlowObjects")
+            HLT__MET__MHTPufitFex,
+            InputJetsName=jetName,
+            InputName=inputName,
+            JetCalibIncludesAreaSub=calibHasAreaSub,
+            JetEventShapeName=rhoKey,
+        )
+
diff --git a/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/MET/ConfigHelpers.py b/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/MET/ConfigHelpers.py
index 64dea08ca274..34aeda220639 100644
--- a/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/MET/ConfigHelpers.py
+++ b/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/MET/ConfigHelpers.py
@@ -11,7 +11,7 @@ import six
 
 # The keys from the MET chain dict that directly affect reconstruction
 # The order here is important as it also controls the dict -> string conversion
-recoKeys = ["EFrecoAlg", "calib", "addInfo"]
+recoKeys = ["EFrecoAlg", "calib", "jetDataType", "jetCalib", "addInfo"]
 
 def extractMETRecoDict(chainDict, fillDefaults=True):
     """ Extract the keys relevant to reconstruction from a provided dictionary
@@ -40,7 +40,11 @@ def metRecoDictToString(recoDict, skipDefaults=True):
 def jetRecoDictForMET(**recoDict):
     """ Get a jet reco dict that's usable for the MET slice """
     keys = ["recoAlg", "dataType", "calib", "jetCalib", "trkopt", "cleaning"]
-    return {k: recoDict.get(k, JetChainParts_Default[k]) for k in keys}
+    jrd = {k: recoDict.get(k, JetChainParts_Default[k]) for k in keys}
+    if "jetDataType" in recoDict:
+        # Allow for the renaming dataType -> jetDataType
+        jrd["dataType"] = recoDict["jetDataType"]
+    return jrd
 
 
 
diff --git a/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/MET/METRecoSequences.py b/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/MET/METRecoSequences.py
index f5614c708710..4bfaa2be9f78 100644
--- a/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/MET/METRecoSequences.py
+++ b/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/MET/METRecoSequences.py
@@ -3,16 +3,108 @@
 #
 
 
-from .ConfigHelpers import AlgConfig
+from .ConfigHelpers import AlgConfig, jetRecoDictForMET
 from AthenaCommon.Logging import logging
+from TriggerMenuMT.HLTMenuConfig.CommonSequences.CaloSequenceSetup import (
+    caloClusterRecoSequence,
+    LCCaloClusterRecoSequence,
+)
+from eflowRec.PFHLTSequence import PFHLTSequence
+from ..Jet.JetRecoConfiguration import defineJetConstit
+from ..Jet.JetTrackingConfig import JetTrackingSequence
+from JetRecConfig.ConstModHelpers import getConstitModAlg
+from ..Menu.MenuComponents import RecoFragmentsPool
+from TrigEFMissingET.TrigEFMissingETConf import (
+    HLT__MET__CVFAlg,
+    ApproximateTrackToLayerTool,
+)
+from InDetTrackSelectionTool.InDetTrackSelectionToolConf import (
+    InDet__InDetTrackSelectionTool,
+)
+from TrackVertexAssociationTool.TrackVertexAssociationToolConf import (
+    CP__TrackVertexAssociationTool,
+)
+
+
 log = logging.getLogger(__name__)
 
-                    
 
 def metAthSequence(dummyFlags, **recoDict):
     conf = AlgConfig.fromRecoDict(**recoDict)
     return conf.athSequence, conf.inputMaker, conf.outputKey
 
+
 def metRecoSequence(dummyFlags, **recoDict):
     conf = AlgConfig.fromRecoDict(**recoDict)
     return conf.recoSequence, conf.outputKey
+
+
+def pfoRecoSequence(dummyFlags, RoIs, **recoDict):
+    """ Get PFOs as inputs for the MET algorithms
+    
+    Returns the list of input sequences and the pfo prefix
+    """
+
+    tcSeq, clusterName = RecoFragmentsPool.retrieve(
+        caloClusterRecoSequence, flags=None, RoIs=RoIs
+    )
+    pfSeq, pfoPrefix = RecoFragmentsPool.retrieve(
+        PFHLTSequence, flags=None, clustersin=clusterName, tracktype="ftf"
+    )
+    # The jet constituent modifier sequence here is to apply the correct weights
+    # and decorate the PV matching decoration
+    jetRecoDict = jetRecoDictForMET(trkopt="ftf", **recoDict)
+    jetRecoDict["calib"] = "em"
+    jetRecoDict["dataType"] = "pf"
+    constit = defineJetConstit(jetRecoDict, pfoPrefix=pfoPrefix)
+    constit_mod_seq = getConstitModAlg(
+        constit,
+        "HLT",
+        tvaKey="JetTrackVtxAssoc_{trkopt}".format(**jetRecoDict),
+        vtxKey="HLT_IDVertex_FS",
+    )
+    return [tcSeq, pfSeq, constit_mod_seq], pfoPrefix
+
+
+def cvfClusterSequence(dummyFlags, RoIs, **recoDict):
+    """ Get the clusters with the CVF calculated
+
+    Returns the list of input sequences, the cluster container name and the CVF
+    decoration name
+    """
+    calib = recoDict["calib"]
+    if calib == "em":
+        tcSeq, clusterName = RecoFragmentsPool.retrieve(
+            caloClusterRecoSequence, flags=None, RoIs=RoIs
+        )
+    elif calib == "lcw":
+        tcSeq, clusterName = RecoFragmentsPool.retrieve(
+            LCCaloClusterRecoSequence, flags=None, RoIs=RoIs
+        )
+    else:
+        raise ValueError(f"Unsupported calib state {calib} requested!")
+
+    trkopt = "ftf"
+    trackSeq, trackColls = RecoFragmentsPool.retrieve(
+        JetTrackingSequence, flags=None, trkopt=trkopt, RoIs=RoIs
+    )
+
+    # These are the names set by the upstream algorithms. Unfortunately
+    # these aren't passed to us - we just have to 'know' them
+    tracks = "HLT_IDTrack_FS_FTF"
+    vertices = "HLT_IDVertex_FS"
+
+    cvfAlg = RecoFragmentsPool.retrieve(
+        HLT__MET__CVFAlg,
+        f"{calib}{trkopt}ClusterCVFAlg",
+        InputClusterKey=clusterName,
+        InputTrackKey=tracks,
+        InputVertexKey=vertices,
+        OutputCVFKey="CVF",
+        TrackSelectionTool=InDet__InDetTrackSelectionTool(CutLevel="TightPrimary"),
+        TVATool=CP__TrackVertexAssociationTool(
+            WorkingPoint="Loose", VertexContainer=vertices
+        ),
+        ExtensionTool=ApproximateTrackToLayerTool(),
+    )
+    return [tcSeq, trackSeq, cvfAlg], clusterName, "CVF"
diff --git a/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/Menu/LS2_v1.py b/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/Menu/LS2_v1.py
index 15829ab323f7..ab3082aae6f1 100644
--- a/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/Menu/LS2_v1.py
+++ b/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/Menu/LS2_v1.py
@@ -33,13 +33,13 @@ def setupMenu():
     TriggerFlags.MuonSlice.signatures = TriggerFlags.MuonSlice.signatures() + [
         #ATR-19985
         ChainProp(name='HLT_mu6_idperf_L1MU6', groups=SingleMuonGroup),
-        ChainProp(name='HLT_mu24_idperf_L1MU20', groups=SingleMuonGroup),        
-        ChainProp(name='HLT_mu6_mu6noL1_L1MU6', l1SeedThresholds=['MU6','FSNOSEED'], mergingStrategy='serial', groups=MultiMuonGroup),        
-
+        ChainProp(name='HLT_mu24_idperf_L1MU20', groups=SingleMuonGroup),
+        ChainProp(name='HLT_mu6_mu6noL1_L1MU6', l1SeedThresholds=['MU6','FSNOSEED'], mergingStrategy='serial', groups=MultiMuonGroup),
+    
         #test chains
         ChainProp(name='HLT_mu6fast_L1MU6', groups=SingleMuonGroup),
         ChainProp(name='HLT_mu6Comb_L1MU6', groups=SingleMuonGroup),
-        ChainProp(name='HLT_mu6_L1MU6',     groups=SingleMuonGroup),  
+        ChainProp(name='HLT_mu6_L1MU6',     groups=SingleMuonGroup),
 
         ChainProp(name='HLT_mu20_ivar_L1MU6',      groups=SingleMuonGroup),
         ChainProp(name='HLT_mu6_ivarmedium_L1MU6', groups=SingleMuonGroup),
@@ -54,11 +54,11 @@ def setupMenu():
 
         # this is for test only
         ChainProp(name='HLT_2mu6_Dr_L12MU4',  groups=MultiMuonGroup),
-         #  ChainProp(name='HLT_mu6_Dr_mu4_Dr_L12MU4', l1SeedThresholds=['MU4']*2, groups=MultiMuonGroup),
+        #  ChainProp(name='HLT_mu6_Dr_mu4_Dr_L12MU4', l1SeedThresholds=['MU4']*2, groups=MultiMuonGroup),
         # ATR-20049
         ChainProp(name='HLT_mu6_mu4_L12MU4',  l1SeedThresholds=['MU4']*2, groups=MultiMuonGroup),
         ChainProp(name='HLT_2mu6Comb_L12MU6', l1SeedThresholds=['MU6'],   groups=MultiMuonGroup),
-
+        
         # in planned primary as an option
         ChainProp(name='HLT_mu24_ivarmedium_L1MU20', groups=SingleMuonGroup),
 
@@ -71,18 +71,18 @@ def setupMenu():
         ChainProp(name='HLT_mu60_L1MU20', groups=SingleMuonGroup),
         ChainProp(name='HLT_mu24_mu10noL1_L1MU20', l1SeedThresholds=['MU20','FSNOSEED'], mergingStrategy='serial', groups=MultiMuonGroup),
         ChainProp(name='HLT_mu26_mu8noL1_L1MU20', l1SeedThresholds=['MU20','FSNOSEED'], mergingStrategy='serial', groups=MultiMuonGroup),
-        ChainProp(name='HLT_mu26_mu10noL1_L1MU20', l1SeedThresholds=['MU20','FSNOSEED'], mergingStrategy='serial', groups=MultiMuonGroup),
+        ChainProp(name='HLT_mu26_mu10noL1_L1MU20', l1SeedThresholds=['MU20', 'FSNOSEED'], mergingStrategy='serial', groups=MultiMuonGroup),
         ChainProp(name='HLT_mu28_mu8noL1_L1MU20', l1SeedThresholds=['MU20','FSNOSEED'], mergingStrategy='serial', groups=MultiMuonGroup),
         ChainProp(name='HLT_mu22_2mu4noL1_L1MU20', l1SeedThresholds=['MU20','FSNOSEED'], mergingStrategy='serial', groups=MultiMuonGroup),
         ChainProp(name='HLT_mu24_2mu4noL1_L1MU20', l1SeedThresholds=['MU20','FSNOSEED'], mergingStrategy='serial', groups=MultiMuonGroup),
         ChainProp(name="HLT_mu10_L1MU10", groups=SingleMuonGroup),
         ChainProp(name='HLT_2mu4_L12MU4',  groups=MultiMuonGroup),
 
-       # ATR-19452
+        # ATR-19452
         ChainProp(name='HLT_2mu4_muonqual_L12MU4',  groups=MultiMuonGroup),
         ChainProp(name='HLT_2mu6_muonqual_L12MU6',  groups=MultiMuonGroup),
 
-       # ATR-20650
+        # ATR-20650
         ChainProp(name='HLT_mu0_muoncalib_L1MU4_EMPTY', groups=SingleMuonGroup),
         ChainProp(name='HLT_mu0_muoncalib_L1MU20',      groups=SingleMuonGroup),
 
@@ -93,7 +93,7 @@ def setupMenu():
         ChainProp(name="HLT_mu26_L1MU20", groups=SingleMuonGroup),
 
 
-     ]
+    ]
 
     TriggerFlags.EgammaSlice.signatures = TriggerFlags.EgammaSlice.signatures() + [
         # ElectronChains----------
@@ -133,7 +133,7 @@ def setupMenu():
         ChainProp(name='HLT_e20_lhmedium_e15_lhmedium_Zee_L12EM3', groups=MultiElectronGroup),    
         # for moving to PhysicsP1, ATR-21242
         # ChainProp(name='HLT_2e17_etcut_L12EM15VHI', stream=[PhysicsStream], groups=MultiElectronGroup),
-
+   
         # PhotonChains------------
         # these are to debug photon working points should be removed in production
         ChainProp(name='HLT_g5_etcut_L1EM3', groups=SinglePhotonGroup),
@@ -152,7 +152,7 @@ def setupMenu():
 
         # ATR-19360
         ChainProp(name='HLT_g5_etcut_LArPEB_L1EM3',stream=['LArCells'], groups=SinglePhotonGroup),
-
+    
         # for moving to PhysicsP1, ATR-21242
         ChainProp(name='HLT_g140_etcut_L1EM22VHI', groups=SinglePhotonGroup),
     ]
@@ -163,12 +163,15 @@ def setupMenu():
         ChainProp(name='HLT_xe30_tcpufit_L1XE10', groups=SingleMETGroup),
         ChainProp(name='HLT_xe30_trkmht_L1XE10', groups=SingleMETGroup),
         ChainProp(name='HLT_xe30_pfsum_L1XE10', groups=SingleMETGroup),
-
+        ChainProp(name='HLT_xe30_pfopufit_L1XE10', groups=SingleMETGroup),
+        ChainProp(name='HLT_xe30_cvfpufit_L1XE10', groups=SingleMETGroup),
+        ChainProp(name='HLT_xe30_mhtpufit_em_subjesgscIS_L1XE10', groups=SingleMETGroup),
+        ChainProp(name='HLT_xe30_mhtpufit_pf_subjesgscIS_L1XE10', groups=SingleMETGroup),
         ChainProp(name='HLT_xe110_tc_em_L1XE50', groups=SingleMETGroup),
         ChainProp(name='HLT_xe110_mht_L1XE50', groups=SingleMETGroup),
         ChainProp(name='HLT_xe110_tcpufit_L1XE50', groups=SingleMETGroup),
         # MultiMET Chain
-        ChainProp(name='HLT_xe30_cell_xe30_tcpufit_L1XE10',l1SeedThresholds=['XE10']*2, groups=MultiMETGroup), #must be FS seeded
+        ChainProp(name='HLT_xe30_cell_xe30_tcpufit_L1XE10', l1SeedThresholds=['XE10']*2, groups=MultiMETGroup), #must be FS seeded
     ]
 
     TriggerFlags.JetSlice.signatures = TriggerFlags.JetSlice.signatures() + [
@@ -268,15 +271,15 @@ def setupMenu():
         ChainProp(name='HLT_e7_lhmedium_mu24_L1MU20',l1SeedThresholds=['EM3','MU20'],  mergingStrategy='serial', stream=[PhysicsStream], groups=MultiElectronGroup),
         # Test photon-muon chain (isolation is there to have different number of steps)
         ChainProp(name='HLT_g25_medium_mu24_ivarmedium_L1MU20',l1SeedThresholds=['EM15VH','MU20'], mergingStrategy='serial', stream=[PhysicsStream], groups=MultiElectronGroup),
-
+    
         # electron + photon stay in the same step - these need to be parallel merged!
         ChainProp(name='HLT_e3_etcut1step_g5_etcut_L12EM3',l1SeedThresholds=['EM3','EM3'], mergingStrategy='parallel', stream=[PhysicsStream], groups=MultiElectronGroup),
-
+  
         # Test chains for muon + jet/MET merging/aligning
         ChainProp(name='HLT_mu6fast_xe30_mht_L1XE10', l1SeedThresholds=['MU6','XE10'], mergingStrategy='serial', stream=[PhysicsStream], groups=SingleMETGroup),
         ChainProp(name='HLT_mu6fast_j45_nojcalib_L1J20', l1SeedThresholds=['MU6','J20'], mergingStrategy='serial', stream=[PhysicsStream], groups=SingleMETGroup),
-
-
+    
+    
     ]
     TriggerFlags.HeavyIonSlice.signatures  = TriggerFlags.HeavyIonSlice.signatures() + []
     TriggerFlags.BeamspotSlice.signatures  = TriggerFlags.BeamspotSlice.signatures() + [
diff --git a/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/Menu/SignatureDicts.py b/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/Menu/SignatureDicts.py
index eb42124b4c2a..22495dc13a3c 100644
--- a/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/Menu/SignatureDicts.py
+++ b/Trigger/TriggerCommon/TriggerMenuMT/python/HLTMenuConfig/Menu/SignatureDicts.py
@@ -32,13 +32,13 @@ SliceIDDict = {
 
 AllowedSignatures = ["jet", "bjet", "ht",
                      "electron", "photon", "egamma",
-                     "muon", 
+                     "muon",
                      "met",
-                     "tau", 
-                     "minbias", 
-                     "heavyion", 
-                     "cosmic", 
-                     "calibration", "streaming", "monitoring", 'eb']        
+                     "tau",
+                     "minbias",
+                     "heavyion",
+                     "cosmic",
+                     "calibration", "streaming", "monitoring", 'eb']
 
 #==========================================================
 # ---- Generic Template for all chains          ----
@@ -70,17 +70,17 @@ TestChainParts = {
     'trigType'       : ['TestChain'],
     'threshold'      : '',
     'addInfo'        : [''],
-    }
+}
 
 # ---- Test Dictionary of default Values ----
 TestChainParts_Default = {
     'signature'      : ['Test'],
     'L1threshold'    : '',
-    'multiplicity'   : '',    
+    'multiplicity'   : '',
     'trigType'       : '',
     'threshold'      : '',
     'addInfo'        : [],
-    }
+}
 
 #==========================================================
 # Jet
@@ -106,7 +106,7 @@ JetChainParts = {
     'calib'        : ['em', 'lcw'],
     'jetCalib'     : ['jes', 'subjes', 'subjesIS', 'subjesgscIS', 'subresjesgscIS', 'nojcalib'],
     'scan'         : ['FS',],
-    'addInfo'      : ['perf'],    
+    'addInfo'      : ['perf'],
 
     'TLA'          : [],
     'dataScouting' : [],
@@ -114,7 +114,7 @@ JetChainParts = {
     'topo'         : AllowedTopos_jet,
 
     'bTag'         : ['boffperf'  ,
-                      'bmv2c2040' , 'bmv2c2050' , 'bmv2c2060' , 'bmv2c2070' , 'bmv2c2077' , 'bmv2c2085' , 
+                      'bmv2c2040' , 'bmv2c2050' , 'bmv2c2060' , 'bmv2c2070' , 'bmv2c2077' , 'bmv2c2085' ,
                       'bmv2c1040' , 'bmv2c1050' , 'bmv2c1060' , 'bmv2c1070' , 'bmv2c1077' , 'bmv2c1085' ,
                       'bhmv2c1040', 'bhmv2c1050', 'bhmv2c1060', 'bhmv2c1070', 'bhmv2c1077', 'bhmv2c1085'],
     'bTracking'    : [],
@@ -145,9 +145,9 @@ JetChainParts_Default = {
     'calib'        :'em',
     'jetCalib'     :'subjesIS',
     'scan'         :'FS',
-    'addInfo'      : [],    
-    'TLA'          : '',    
-    'topo'         : [],  
+    'addInfo'      : [],
+    'TLA'          : '',
+    'topo'         : [],
     'bTag'         : '',
     'bTracking'    : '',
     'bConfig'      : [],
@@ -156,7 +156,7 @@ JetChainParts_Default = {
     'trkopt'       : 'notrk',
     'hypoScenario' : 'simple',
     'smc'          : 'nosmc',
-    }
+}
 
 # ---- bJet Dictionary of default Values that are different to the ones for normal jet chains ----
 bJetChainParts_Default = {
@@ -179,7 +179,7 @@ HTChainParts_Default['trigType']     = 'ht'
 HTChainParts_Default['extra']     = ''
 
 #==========================================================
-# Muon 
+# Muon
 #==========================================================
 AllowedTopos_mu = []
 
@@ -188,7 +188,7 @@ MuonChainParts = {
     'signature'      : ['Muon'],
     'L1threshold'    : '',
     'chainPartName'  : [],
-    'multiplicity'   : '',    
+    'multiplicity'   : '',
     'trigType'       : ['mu'],
     'etaRange'       : ['0eta2550','0eta105'],
     'threshold'      : '',
@@ -199,23 +199,23 @@ MuonChainParts = {
     'addInfo'        : ['1step','idperf','3layersEC','cosmic',"muonqual"],
     'topo'           : AllowedTopos_mu,
     'flavour'        : [],
-    }
+}
 # ---- MuonDictinary of default Values ----
 MuonChainParts_Default = {
     'signature'      : ['Muon'],
     'L1threshold'    : '',
-    'multiplicity'   : '',    
+    'multiplicity'   : '',
     'trigType'       : ['mu'],
-    'etaRange'       : '0eta250',    
+    'etaRange'       : '0eta250',
     'threshold'      : '',
-    'extra'          : '',    
-    'IDinfo'         : '',    
-    'isoInfo'        : '',    
+    'extra'          : '',
+    'IDinfo'         : '',
+    'isoInfo'        : '',
     'addInfo'        : [],
     'invMassInfo'    : '',
     'topo'           : [],
     'flavour'        : '',
-    }
+}
 
 #==========================================================
 # Bphysics
@@ -246,7 +246,7 @@ TauChainParts = {
     'preselection' : ['track', 'tracktwo', 'tracktwoEF', 'tracktwoMVA', 'tracktwoEFmvaTES', 'ptonly', ],
     'selection'    : ['medium1', 'verylooseRNN', 'looseRNN', 'mediumRNN', 'tightRNN', 'perf', 'idperf'],
     'multiplicity' : '',
-    'trigType'     : ['tau'],   
+    'trigType'     : ['tau'],
     'trkInfo'      : '',
     'extra'        : '',
     'recoAlg'      : '',
@@ -262,7 +262,7 @@ TauChainParts_Default = {
     'preselection' : 'tracktwo',
     'selection'    : 'medium1',
     'multiplicity' : '',
-    'trigType'     : ['tau'],   
+    'trigType'     : ['tau'],
     'trkInfo'      : [],
     'extra'        : '',
     'recoAlg'      : '',
@@ -283,15 +283,17 @@ METChainParts = {
     'threshold'    : '',
     'multiplicity' : '',
     'topo'         : AllowedTopos_xe,
-    'trigType'     : ['xe'],   
+    'trigType'     : ['xe'],
     'extra'        : ['noL1'],
-    'calib'        : ['lcw','em'],    
+    'calib'        : ['lcw','em'],
+    'jetCalib'     : JetChainParts['jetCalib'],
     'L2recoAlg'    : [],
-    'EFrecoAlg'    : ['cell', 'tc', 'tcpufit', 'mht', 'trkmht', 'pfsum'],
+    'EFrecoAlg'    : ['cell', 'tc', 'tcpufit', 'mht', 'trkmht', 'pfsum', 'cvfpufit', 'pfopufit', 'mhtpufit'],
+    'jetDataType'  : JetChainParts['dataType'],
     'L2muonCorr'   : [],
     'EFmuonCorr'   : [],
     'addInfo'      : ['FStracks'],
-    }
+}
 # ---- MetDictinary of default Values ----
 METChainParts_Default = {
     'signature'      : ['MET'],
@@ -300,20 +302,22 @@ METChainParts_Default = {
     'threshold'      : '',
     'extra'          : '',
     'calib'          : 'lcw',
+    'jetCalib'       : JetChainParts_Default['jetCalib'],
     'L2recoAlg'      : '',
     'EFrecoAlg'      : '',
     'L2muonCorr'     : '',
     'EFmuonCorr'     : '',
     'addInfo'        : '',
-    }
+    'jetDataType'    : 'tc',
+}
 
 #==========================================================
 # XS
 #==========================================================
 # ---- xs Dictinary of all allowed Values ----
-XSChainParts = METChainParts 
+XSChainParts = METChainParts
 XSChainParts['signature'] = ['XS']
-XSChainParts['trigType']  = ['xs'],   
+XSChainParts['trigType']  = ['xs']
 
 # ---- xs Dictinary of default Values ----
 XSChainParts_Default = METChainParts_Default
@@ -324,9 +328,9 @@ XSChainParts_Default['trigType']  = ['xs']
 # TE
 #==========================================================
 # ---- te Dictinary of all allowed Values ----
-TEChainParts = METChainParts 
+TEChainParts = METChainParts
 TEChainParts['signature'] = ['TE']
-TEChainParts['trigType']  = ['te'],   
+TEChainParts['trigType']  = ['te']
 
 # ---- te Dictinary of default Values ----
 TEChainParts_Default = METChainParts_Default
@@ -343,7 +347,7 @@ ElectronChainParts = {
     'chainPartName'  : '',
     'L1threshold'    : '',
     'extra'          : '',
-    'multiplicity'   : '',    
+    'multiplicity'   : '',
     'trigType'       : ['e'],
     'threshold'      : '',
     'etaRange'       : [],
@@ -354,11 +358,11 @@ ElectronChainParts = {
     'lhInfo'         : [],
     'L2IDAlg'        : ['noringer'],
     'addInfo'        : [ 'etcut', 'etcut1step',"v2","v3"],
-    }
+}
 # ---- Egamma Dictinary of default Values ----
 ElectronChainParts_Default = {
     'signature'      : ['Electron'],
-    'multiplicity'   : '',    
+    'multiplicity'   : '',
     'L1threshold'         : '',
     'trigType'       : '',
     'threshold'      : '',
@@ -369,14 +373,14 @@ ElectronChainParts_Default = {
     'isoInfo'        : '',
     'reccalibInfo'   : '',
     'trkInfo'        : '',
-    'caloInfo'       : '',   
+    'caloInfo'       : '',
     'lhInfo'         : '',
     'L2IDAlg'        : '',
     'hypoInfo'       : '',
     'recoAlg'        : '',
     'FSinfo'         : '',
     'addInfo'        : [],
-    }
+}
 
 #==========================================================
 # Photon chains
-- 
GitLab