diff --git a/Event/xAOD/xAODJet/Root/JetContainerInfo.cxx b/Event/xAOD/xAODJet/Root/JetContainerInfo.cxx
index 2ad0f7af044f569fcd8e11806d05e051959597c8..f226a7727409c6a24c712f4dae065094cf731f02 100644
--- a/Event/xAOD/xAODJet/Root/JetContainerInfo.cxx
+++ b/Event/xAOD/xAODJet/Root/JetContainerInfo.cxx
@@ -110,6 +110,7 @@ namespace xAOD {
       { "LCPFlow",               LCPFlow },
       { "EMPFlow",               EMPFlow },
       { "EMCPFlow",              EMCPFlow },
+      { "EMPFlowFE",             EMPFlowFE },
       { "TrackCaloCluster",      TrackCaloCluster },
       { "EMTopoOriginSK",        EMTopoOriginSK },
       { "EMTopoOriginCS",        EMTopoOriginCS },
@@ -160,6 +161,7 @@ namespace xAOD {
       { LCPFlow,               "LCPFlow" },
       { EMPFlow,               "EMPFlow" },
       { EMCPFlow,              "EMCPFlow" },
+      { EMPFlowFE,             "EMPFlowFE" },
       { TrackCaloCluster,      "TrackCaloCluster" },
       { EMTopoOriginSK,        "EMTopoOriginSK" },
       { EMTopoOriginCS,        "EMTopoOriginCS" },
diff --git a/Event/xAOD/xAODJet/xAODJet/JetContainerInfo.h b/Event/xAOD/xAODJet/xAODJet/JetContainerInfo.h
index 45d48a4fabe9cd7616ecf10a59001d62325ad552..a6aa9c8940c2582fcdec6634c4cb00ff67ddb64f 100644
--- a/Event/xAOD/xAODJet/xAODJet/JetContainerInfo.h
+++ b/Event/xAOD/xAODJet/xAODJet/JetContainerInfo.h
@@ -95,6 +95,7 @@ namespace xAOD {
       HI,
       HIClusters,
       Other = 100,
+      EMPFlowFE = 200, // Temporary, until xAOD::PFO is phased out and replaced with xAOD::FlowElement
       Uncategorized= 1000
     };
 
diff --git a/Event/xAOD/xAODPFlow/Root/FlowElement_v1.cxx b/Event/xAOD/xAODPFlow/Root/FlowElement_v1.cxx
index 2beecd79cb68459ce5f039ee769a974313736f4f..6bf9fe6f72c4a126721881311d23502b25c257f8 100644
--- a/Event/xAOD/xAODPFlow/Root/FlowElement_v1.cxx
+++ b/Event/xAOD/xAODPFlow/Root/FlowElement_v1.cxx
@@ -48,6 +48,9 @@ namespace xAOD {
 
   }
 
+  void FlowElement_v1::setP4(FourMom_t p4){
+    setP4(p4.Pt(), p4.Eta(), p4.Phi(), p4.M());
+  }
 
   AUXSTORE_PRIMITIVE_SETTER_AND_GETTER(FlowElement_v1, float, charge, setCharge) 
   bool  FlowElement_v1::isCharged() const { return !bool( signalType()& Neutral );}
diff --git a/Event/xAOD/xAODPFlow/xAODPFlow/versions/FlowElement_v1.h b/Event/xAOD/xAODPFlow/xAODPFlow/versions/FlowElement_v1.h
index 79b43c8094e91e6c19fc26363c809414eadbd020..cfac70c55344898755d5c487133cacfab3a41614 100644
--- a/Event/xAOD/xAODPFlow/xAODPFlow/versions/FlowElement_v1.h
+++ b/Event/xAOD/xAODPFlow/xAODPFlow/versions/FlowElement_v1.h
@@ -76,6 +76,7 @@ namespace xAOD {
     virtual Type::ObjectType type() const override ;
     
     void setP4(float pt, float eta, float phi, float m) ;
+    void setP4(FourMom_t p4);
     ///@}
 
     // *************************************************
diff --git a/PhysicsAnalysis/AnalysisCommon/ThinningUtils/src/ThinNegativeEnergyNeutralPFOsAlg.cxx b/PhysicsAnalysis/AnalysisCommon/ThinningUtils/src/ThinNegativeEnergyNeutralPFOsAlg.cxx
index 77a29ea721ffb62e1bb4111adbb2b5e0a1462911..42caf7ee5f6dac9db05fc098eb0b329559745283 100644
--- a/PhysicsAnalysis/AnalysisCommon/ThinningUtils/src/ThinNegativeEnergyNeutralPFOsAlg.cxx
+++ b/PhysicsAnalysis/AnalysisCommon/ThinningUtils/src/ThinNegativeEnergyNeutralPFOsAlg.cxx
@@ -65,7 +65,8 @@ StatusCode ThinNegativeEnergyNeutralPFOsAlg::initialize()
       ATH_MSG_ERROR( "StreamName property has not been initialized." );
       return StatusCode::FAILURE;
     }
-    ATH_CHECK( m_neutralPFOsKey.initialize (m_streamName, m_doThinning) );
+    ATH_CHECK( m_neutralPFOsKey.initialize (m_streamName, m_doThinning && !m_neutralPFOsKey.key().empty()) );
+    ATH_CHECK( m_neutralPFOsFEKey.initialize (m_streamName, m_doThinning && !m_neutralPFOsFEKey.key().empty()) );
     
     // Initialize the counters to zero
     m_nEventsProcessed = 0;
@@ -101,25 +102,43 @@ StatusCode ThinNegativeEnergyNeutralPFOsAlg::execute()
 
     const EventContext& ctx = Gaudi::Hive::currentContext();
     
-    // Retrieve the container
-    SG::ThinningHandle<xAOD::PFOContainer> neutralPFOs (m_neutralPFOsKey, ctx);
-
-    // Set up masks
-    std::vector<bool> mask;
-    int nNeutralPFOs = neutralPFOs->size();
-    m_nNeutralPFOsProcessed += nNeutralPFOs;
-    mask.assign(nNeutralPFOs,false);
-    
-    // Loop over NeutralPFOs and update mask
-    for (int i=0; i<nNeutralPFOs; ++i) {
-        const xAOD::PFO* neutralPFO = (*neutralPFOs)[i];
-        // Retain postive energy neutral PFOs
-        if (neutralPFO->ptEM()>0.0) {mask[i] = true;}
-        else {++m_nNeutralPFOsThinned;}
+    if(!m_neutralPFOsFEKey.key().empty()){
+
+        // Retrieve the container
+        SG::ThinningHandle<xAOD::PFOContainer> neutralPFOs (m_neutralPFOsKey, ctx);
+
+        // Set up masks
+        std::vector<bool> mask;
+        int nNeutralPFOs = neutralPFOs->size();
+        m_nNeutralPFOsProcessed += nNeutralPFOs;
+        mask.assign(nNeutralPFOs,false);
+
+        // Loop over NeutralPFOs and update mask
+        for (int i=0; i<nNeutralPFOs; ++i) {
+            const xAOD::PFO* neutralPFO = (*neutralPFOs)[i];
+            // Retain postive energy neutral PFOs
+            if (neutralPFO->ptEM()>0.0) {mask[i] = true;}
+            else {++m_nNeutralPFOsThinned;}
+        }
+
+        // Apply masks to thinning service
+        neutralPFOs.keep (mask);
     }
 
-    // Apply masks to thinning service
-    neutralPFOs.keep (mask);
+    if(!m_neutralPFOsFEKey.key().empty()){
+        SG::ThinningHandle<xAOD::FlowElementContainer> neutralFEs (m_neutralPFOsFEKey, ctx);
+        std::vector<bool> mask;
+        int nNeutralFEs = neutralFEs->size();
+        mask.assign(nNeutralFEs, false);
+
+        for(int i=0; i<nNeutralFEs; i++){
+            const xAOD::FlowElement* neutralFE = (*neutralFEs)[i];
+            // TODO: Is this OK for LC-scale PFOs?
+            //       Can't access EM-scale momentum without some link to EM-scale FlowElements.
+            if(neutralFE->pt() > 0.0) mask[i] = true;
+        }
+        neutralFEs.keep(mask);
+    }
     
     return StatusCode::SUCCESS;
 }
diff --git a/PhysicsAnalysis/AnalysisCommon/ThinningUtils/src/ThinNegativeEnergyNeutralPFOsAlg.h b/PhysicsAnalysis/AnalysisCommon/ThinningUtils/src/ThinNegativeEnergyNeutralPFOsAlg.h
index 59ffbc43e716a01300b1b465cb8f46bb86b29a00..5ab917823e3c14e4704725f2c8de96ada6fa5218 100644
--- a/PhysicsAnalysis/AnalysisCommon/ThinningUtils/src/ThinNegativeEnergyNeutralPFOsAlg.h
+++ b/PhysicsAnalysis/AnalysisCommon/ThinningUtils/src/ThinNegativeEnergyNeutralPFOsAlg.h
@@ -23,6 +23,7 @@
 #include "StoreGate/ThinningHandleKey.h"
 
 #include "xAODPFlow/PFOContainer.h"
+#include "xAODPFlow/FlowElementContainer.h"
 
 class ThinNegativeEnergyNeutralPFOsAlg
 : public ::AthAlgorithm
@@ -54,7 +55,9 @@ private:
    
     /// Names of the containers to thin
     SG::ThinningHandleKey<xAOD::PFOContainer> m_neutralPFOsKey
-    { this, "NeutralPFOsKey", "JetETMissNeutralParticleFlowObjects", "StoreGate key for the PFOContainer to be thinned" };
+    { this, "NeutralPFOsKey", "JetETMissNeutralParticleFlowObjects", "StoreGate key for the PFOContainer to be thinned (if any)" };
+    SG::ThinningHandleKey<xAOD::FlowElementContainer> m_neutralPFOsFEKey
+    { this, "NeutralPFOsFEKey", "", "StoreGate key for the FlowElementContainer to be thinned (if any)" };
  
     /// Counters
     unsigned long m_nEventsProcessed;
diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkCore/python/ContainersOnTheFly.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkCore/python/ContainersOnTheFly.py
index ea1cbc28ff4da4e37117e2e732f3ca23ab16b44d..f6076b4dd75c412debbb183a643ead53d2a7c69b 100644
--- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkCore/python/ContainersOnTheFly.py
+++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkCore/python/ContainersOnTheFly.py
@@ -68,6 +68,8 @@ ContainersOnTheFly = [
                         ["AntiKt4EMTopoCSSKJetsAux","xAOD::JetAuxContainer"],
                         ["AntiKt4EMPFlowCSSKJets","xAOD::JetContainer"],
                         ["AntiKt4EMPFlowCSSKJetsAux","xAOD::JetAuxContainer"],
+                        ["AntiKt4EMPFlowFEJets","xAOD::JetContainer"],
+                        ["AntiKt4EMPFlowFEJetsAux","xAOD::JetContainerAux"],
                         ["AntiKt4EMTopoJets_BTagging201810","xAOD::JetContainer"],
                         ["AntiKt4EMTopoJets_BTagging201810Aux","xAOD::ShallowAuxContainer"],
                         ["AntiKt4EMPFlowJets_BTagging201903","xAOD::JetContainer"],
diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkCore/python/FullListOfSmartContainers.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkCore/python/FullListOfSmartContainers.py
index 9752dbbbccf3df70f258207c7f24b739fe8ac533..c24f368067682201c40becbae8fe4e2cb01f23e6 100644
--- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkCore/python/FullListOfSmartContainers.py
+++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkCore/python/FullListOfSmartContainers.py
@@ -46,6 +46,7 @@ FullListOfSmartContainers = [
    "AntiKt10LCTopoTrimmedPtFrac5SmallR20ExKt2GASubJets",
    "AntiKt10LCTopoTrimmedPtFrac5SmallR20ExKt3GASubJets",
    "AntiKt4EMPFlowJets",
+   "AntiKt4EMPFlowFEJets",
    "AntiKt4EMPFlowJets_BTagging201903",
    "AntiKt4EMPFlowJets_BTagging201810",
    "AntiKt4EMTopoJets_BTagging201810",
diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkCore/python/SlimmingHelper.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkCore/python/SlimmingHelper.py
index 38393d5a747485cba1df5b7a0cca311431dae09c..a0dcba819a6a0b94812342885c931070aecccaec 100644
--- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkCore/python/SlimmingHelper.py
+++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkCore/python/SlimmingHelper.py
@@ -499,8 +499,10 @@ class SlimmingHelper:
                         items.extend(AntiKt10UFOCSSKRecursiveSoftDropBeta100Zcut5NinfJetsCPContent)
                 elif collectionName=="AntiKt4EMPFlowJets":
                         from DerivationFrameworkJetEtMiss.AntiKt4EMPFlowJetsCPContent import AntiKt4EMPFlowJetsCPContent
-                        #from DerivationFrameworkCore.AntiKt4EMPFlowJetsCPContent import AntiKt4EMPFlowJetsCPContent
                         items.extend(AntiKt4EMPFlowJetsCPContent)
+                elif collectionName=="AntiKt4EMPFlowFEJets":
+                        from DerivationFrameworkJetEtMiss.AntiKt4EMPFlowFEJetsCPContent import AntiKt4EMPFlowFEJetsCPContent
+                        items.extend(AntiKt4EMPFlowFEJetsCPContent)
                 elif collectionName=="AntiKt4EMPFlowJets_BTagging201810":
                         if "AntiKt4EMPFlowJets_BTagging201810" not in self.AppendToDictionary:
                                 self.AppendToDictionary["AntiKt4EMPFlowJets_BTagging201810"]='xAOD::JetContainer'
diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkJetEtMiss/python/AntiKt4EMPFlowFEJetsCPContent.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkJetEtMiss/python/AntiKt4EMPFlowFEJetsCPContent.py
new file mode 100644
index 0000000000000000000000000000000000000000..bc27e6867cea4215539e4f4d4fb3329c5a873067
--- /dev/null
+++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkJetEtMiss/python/AntiKt4EMPFlowFEJetsCPContent.py
@@ -0,0 +1,10 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+AntiKt4EMPFlowFEJetsCPContent = [
+"Kt4EMPFlowFEEventShape",
+"Kt4EMPFlowFEEventShapeAux.Density",
+"AntiKt4EMPFlowFEJets",
+"AntiKt4EMPFlowFEJetsAux.pt.eta.phi.m.JetConstitScaleMomentum_pt.JetConstitScaleMomentum_eta.JetConstitScaleMomentum_phi.JetConstitScaleMomentum_m.NumTrkPt500.SumPtTrkPt500.NumChargedPFOPt500.SumPtChargedPFOPt500.EnergyPerSampling.ActiveArea4vec_eta.ActiveArea4vec_m.ActiveArea4vec_phi.ActiveArea4vec_pt.DetectorEta.DetectorY.FracSamplingMax.FracSamplingMaxIndex.GhostTrack.Jvt.JVFCorr.JvtRpt.NumTrkPt1000.NumChargedPFOPt1000.TrackWidthPt1000.ChargedPFOWidthPt1000.GhostMuonSegmentCount.PartonTruthLabelID.ConeTruthLabelID.HadronConeExclExtendedTruthLabelID.HadronConeExclTruthLabelID.TrueFlavor.Timing",
+"PrimaryVertices",
+"PrimaryVerticesAux.vertexType" 
+]
diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkJetEtMiss/python/ExtendedJetCommon.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkJetEtMiss/python/ExtendedJetCommon.py
index 90fd39d0e95472342afa49bcb93903d01c9bfab2..608b87894bb5b1fc2a2d845491679332ecfdc902 100644
--- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkJetEtMiss/python/ExtendedJetCommon.py
+++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkJetEtMiss/python/ExtendedJetCommon.py
@@ -122,6 +122,10 @@ def addAntiKt10TruthWZJets(sequence,outputlist):
     if DerivationFrameworkIsMonteCarlo:
         addStandardJets("AntiKt", 1.0, "TruthWZ", ptmin=40000, mods="truth_ungroomed", algseq=sequence, outputGroup=outputlist)
 
+def addAntiKt4EMPFlowJetsFE(sequence, outputlist):
+    addCHSPFlowObjectsFE()
+    addStandardJets("AntiKt", 0.4, "EMPFlowFE", ptmin=10000, ptminFilter=15000, mods="pflow_ungroomed", algseq=sequence, outputGroup=outputlist)
+
 ##################################################################  
 
 def replaceAODReducedJets(jetlist,sequence,outputlist):
@@ -395,7 +399,7 @@ def addJetPtAssociation(jetalg, truthjetalg, sequence, algname):
     if hasattr(ToolSvc,jetptassociationtoolname):
         jetaugtool.JetPtAssociationTool = getattr(ToolSvc,jetptassociationtoolname)
     else:
-        jetptassociationtool = CfgMgr.JetPtAssociationTool(jetptassociationtoolname, InputContainer=truthjetalg, AssociationName="GhostTruth")
+        jetptassociationtool = CfgMgr.JetPtAssociationTool(jetptassociationtoolname, JetContainer=jetalg+'Jets', MatchingJetContainer=truthjetalg, AssociationName="GhostTruth")
         ToolSvc += jetptassociationtool
         jetaugtool.JetPtAssociationTool = jetptassociationtool
 
@@ -480,7 +484,7 @@ def addQGTaggerTool(jetalg, sequence, algname, truthjetalg=None ):
         if hasattr(ToolSvc,jetptassociationtoolname):
             jetaugtool.JetPtAssociationTool = getattr(ToolSvc,jetptassociationtoolname)
         else:
-            jetptassociationtool = CfgMgr.JetPtAssociationTool(jetptassociationtoolname, InputContainer=truthjetalg, AssociationName="GhostTruth")
+            jetptassociationtool = CfgMgr.JetPtAssociationTool(jetptassociationtoolname, JetContainer=jetalg+'Jets', MatchingJetContainer=truthjetalg, AssociationName="GhostTruth")
             ToolSvc += jetptassociationtool
             jetaugtool.JetPtAssociationTool = jetptassociationtool
 
@@ -693,6 +697,8 @@ def addConstModJets(jetalg,radius,inputtype,constmods,sequence,outputlist,custom
         constit = JetConstit( xAODType.CaloCluster, ["LC","Origin"])
     elif inputtype == "EMPFlow":
         constit = JetConstit( xAODType.ParticleFlow )
+    elif inputtype == "EMPFlowFE":
+        constit = JetConstit( xAODType.FlowElement )
 
     constit.modifiers += constmods
 
@@ -773,6 +779,26 @@ def addCHSPFlowObjects():
             job.jetalg.Tools.append(jtm.jetconstitCHSPFlow)
             extjetlog.info("Added CHS PFlow sequence to \'jetalg\'")
             extjetlog.info(job.jetalg.Tools)
+
+def addCHSPFlowObjectsFE():
+    # Only act if the collection does not already exist
+    from RecExConfig.AutoConfiguration import IsInInputFile
+    if not IsInInputFile("xAOD::FlowElementContainer","CHSFlowElements"):
+        # Check that an alg doing this has not already been inserted
+        from AthenaCommon.AlgSequence import AlgSequence
+        job = AlgSequence()
+        from JetRec.JetRecStandard import jtm
+        if not hasattr(job,"jetalgCHSPFlowFE") and not hasattr(jtm,"jetconstitCHSPFlowFE"):
+            from JetRec.JetRecConf import JetToolRunner
+            jtm += JetToolRunner("jetconstitCHSPFlowFE",
+                                 EventShapeTools=[],
+                                 Tools=[jtm.JetConstitSeq_PFlowCHS_FE])
+            # Add this tool runner to the JetAlgorithm instance "jetalg"
+            # which runs all preparatory tools
+            # This was added by JetCommon
+            job.jetalg.Tools.append(jtm.jetconstitCHSPFlowFE)
+            extjetlog.info("Added CHS PFlow (FlowElement) sequence to \'jetalg\'")
+            extjetlog.info(job.jetalg.Tools)
 ##################################################################
 applyJetCalibration_xAODColl("AntiKt4EMTopo")
 updateJVT_xAODColl("AntiKt4EMTopo")
diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkJetEtMiss/python/JetCommon.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkJetEtMiss/python/JetCommon.py
index 08edf7e56256a71b21cf42d66d6f65b9e6b30597..cca3d4b75189fa213237e34e2e3e6bc14232c76e 100644
--- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkJetEtMiss/python/JetCommon.py
+++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkJetEtMiss/python/JetCommon.py
@@ -143,11 +143,12 @@ def reCreatePseudoJets(jetalg, rsize, inputtype, variableRMassScale=-1.0, variab
             return [jtm.tools[tmpName]]
 
         # no container exist. simply build a new one.
-        if inputtype=="LCTopo" or inputtype=="EMTopo" or inputtype == "EMPFlow" or inputtype == "EMCPFlow":
+        if inputtype=="LCTopo" or inputtype=="EMTopo" or inputtype == "EMPFlow" or inputtype == "EMCPFlow" or inputtype == "EMPFlowFE":
             defaultmods = {"EMTopo":"emtopo_ungroomed",
                            "LCTopo":"lctopo_ungroomed",
                            "EMPFlow":"pflow_ungroomed",
                            "EMCPFlow":"pflow_ungroomed",
+                           "EMPFlowFE":"pflow_ungroomed",
                            "Truth":"truth_ungroomed",
                            "TruthWZ":"truth_ungroomed",
                            "PV0Track":"track_ungroomed"}
@@ -191,6 +192,8 @@ def reCreatePseudoJets(jetalg, rsize, inputtype, variableRMassScale=-1.0, variab
             constit = JetConstit( xAODType.CaloCluster, ["LC","Origin"])
         elif inputtype == "EMPFlow":
             constit = JetConstit( xAODType.ParticleFlow )
+        elif inputtype == "EMPFlowFE":
+            constit = JetConstit( xAODType.FlowElement )
 
         constit.modifiers += constmods
 
@@ -441,6 +444,7 @@ def addStandardJets(jetalg, rsize, inputtype, ptmin=0., ptminFilter=0.,
                        "LCTopo":"lctopo_ungroomed",
                        "EMPFlow":"pflow_ungroomed",
                        "EMCPFlow":"pflow_ungroomed",
+                       "EMPFlowFE":"pflow_ungroomed",
                        "Truth":"truth_ungroomed",
                        "TruthWZ":"truth_ungroomed",
                        "PV0Track":"track_ungroomed",
@@ -459,7 +463,7 @@ def addStandardJets(jetalg, rsize, inputtype, ptmin=0., ptminFilter=0.,
             finderArgs['overwrite']=True
     
         # map the input to the jtm code for PseudoJetGetter
-        getterMap = dict( LCTopo = 'lctopo', EMTopo = 'emtopo', EMPFlow = 'empflow', EMCPFlow = 'emcpflow', 
+        getterMap = dict( LCTopo = 'lctopo', EMTopo = 'emtopo', EMPFlow = 'empflow', EMCPFlow = 'emcpflow', EMPFlowFE = 'empflowfe',
                           Truth = 'truth',  TruthWZ = 'truthwz', TruthDressedWZ = 'truthdressedwz', TruthCharged = 'truthcharged',
                           PV0Track='pv0track')
         # create the finder for the temporary collection.
@@ -483,7 +487,7 @@ def addStandardJets(jetalg, rsize, inputtype, ptmin=0., ptminFilter=0.,
         alg = JetAlgorithm(algname, Tools = pretools+[finderTool])
         dfjetlog.info( "Added "+algname+" to sequence "+algseq.name() )
         algseq += alg
-        DFJetAlgs[algname] = alg;
+        DFJetAlgs[algname] = alg
 
 ################################################################## 
 # Schedule the adding of BCID info
diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloEnergies.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloEnergies.h
index 73a499e5b1517a817d3d1a033d255b8577d92d17..bbab1cc914ad645b78ee813666b5fe89a17053f7 100644
--- a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloEnergies.h
+++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetCaloEnergies.h
@@ -30,6 +30,7 @@ public:
 protected:
   void fillEperSamplingCluster(const xAOD::Jet &jet, std::vector<float> & ePerSampling ) const ;
   void fillEperSamplingPFO(const xAOD::Jet &jet, std::vector<float> & ePerSampling ) const ;
+  void fillEperSamplingFE(const xAOD::Jet &jet, std::vector<float> & ePerSampling ) const ;
 
 private:
   Gaudi::Property<std::string> m_jetContainerName{this, "JetContainer", "", "SG key for the input jet container"};
diff --git a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetPtAssociationTool.h b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetPtAssociationTool.h
index 275d72f46cb640cf54a4836682794d5d97559dc4..b180c7507533083a8b9a2d42c8ffcb48cf88a348 100644
--- a/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetPtAssociationTool.h
+++ b/Reconstruction/Jet/JetMomentTools/JetMomentTools/JetPtAssociationTool.h
@@ -65,7 +65,8 @@ private:  // data
 
   /// Properties.
   Gaudi::Property<std::string> m_aname{this, "AssociationName", "", "Key for association vector"};
-  SG::ReadHandleKey<xAOD::JetContainer> m_jetContainerName{this, "JetContainer", "", "Input jet container"};
+  Gaudi::Property<std::string> m_jetContainerName{this, "JetContainer", "", "Input jet container"};
+  SG::ReadHandleKey<xAOD::JetContainer> m_matchingContainerName{this, "MatchingJetContainer", "", "Jet container to match to"};
 
   SG::WriteDecorHandleKey<xAOD::JetContainer> m_assocFracKey{this, "AssociationFractionName", "AssociationFraction", "SG key for output AssociationFraction decoration"};
   SG::WriteDecorHandleKey<xAOD::JetContainer> m_assocLinkKey{this, "AssociationLinkName", "AssociationLink", "SG key for output AssociationLink decoration"};
diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetCaloEnergies.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetCaloEnergies.cxx
index 4a6ae0627097623de345051add25e23983c7c95e..4393e2949fbc28a9d7c03b23191bf02e7c9465ee 100644
--- a/Reconstruction/Jet/JetMomentTools/Root/JetCaloEnergies.cxx
+++ b/Reconstruction/Jet/JetMomentTools/Root/JetCaloEnergies.cxx
@@ -7,9 +7,11 @@
 #include "xAODJet/JetAccessorMap.h"
 
 #include "JetUtils/JetCaloQualityUtils.h"
+#include "PFlowUtils/FEHelpers.h"
 
 #include "xAODCaloEvent/CaloCluster.h"
 #include "xAODPFlow/PFO.h"
+#include "xAODPFlow/FlowElement.h"
 
 #include <vector>
 
@@ -62,13 +64,17 @@ StatusCode JetCaloEnergies::decorate(const xAOD::JetContainer& jets) const {
     if ( ctype  == xAOD::Type::CaloCluster ) {
       ATH_MSG_VERBOSE("  Constituents are calo clusters.");
       fillEperSamplingCluster(*jet, ePerSampling);
-    
+
     } else if (ctype  == xAOD::Type::ParticleFlow) {
       ATH_MSG_VERBOSE("  Constituents are pflow objects.");
       fillEperSamplingPFO(*jet, ePerSampling);
 
+    } else if (ctype  == xAOD::Type::FlowElement) {
+      ATH_MSG_VERBOSE("  Constituents are FlowElements.");
+      fillEperSamplingFE(*jet, ePerSampling);
+
     }else {
-      ATH_MSG_VERBOSE("Constituents are not calo clusters nor pflow objects.");
+      ATH_MSG_VERBOSE("Constituents are not CaloClusters, PFOs, or FlowElements.");
     }
 
   }
@@ -97,10 +103,10 @@ void JetCaloEnergies::fillEperSamplingCluster(const xAOD::Jet& jet, std::vector<
   psFracHandle(jet) = jet::JetCaloQualityUtils::presamplerFraction( &jet );
 }
 
-#define FillESamplingPFO( LAYERNAME )					\
-  float E_##LAYERNAME = 0.0;						\
+#define FillESamplingPFO( LAYERNAME )                                        \
+  float E_##LAYERNAME = 0.0;                                                \
   if (constit->attribute(xAOD::PFODetails::eflowRec_LAYERENERGY_##LAYERNAME, E_##LAYERNAME)) { \
-    ePerSampling[CaloSampling::LAYERNAME] += E_##LAYERNAME;		\
+    ePerSampling[CaloSampling::LAYERNAME] += E_##LAYERNAME;                \
   }
   
 void JetCaloEnergies::fillEperSamplingPFO(const xAOD::Jet & jet, std::vector<float> & ePerSampling ) const {
@@ -114,50 +120,50 @@ void JetCaloEnergies::fillEperSamplingPFO(const xAOD::Jet & jet, std::vector<flo
     if (jet.rawConstituent(i)->type()==xAOD::Type::ParticleFlow){
       const xAOD::PFO* constit = static_cast<const xAOD::PFO*>(jet.rawConstituent(i));
       if ( fabs(constit->charge())>FLT_MIN ){
-	eTot += constit->track(0)->e();
+        eTot += constit->track(0)->e();
       } else {
-	eTot += constit->e();
-	FillESamplingPFO(PreSamplerB);
-	FillESamplingPFO(EMB1);
-	FillESamplingPFO(EMB2);
-	FillESamplingPFO(EMB3);
-
-	FillESamplingPFO(PreSamplerE);
-	FillESamplingPFO(EME1);
-	FillESamplingPFO(EME2);
-	FillESamplingPFO(EME3);
-
-	FillESamplingPFO(HEC0);
-	FillESamplingPFO(HEC1);
-	FillESamplingPFO(HEC2);
-	FillESamplingPFO(HEC3);
-
-	FillESamplingPFO(TileBar0);
-	FillESamplingPFO(TileBar1);
-	FillESamplingPFO(TileBar2);
-
-	FillESamplingPFO(TileGap1);
-	FillESamplingPFO(TileGap2);
-	FillESamplingPFO(TileGap3);
-
-	FillESamplingPFO(TileExt0);
-	FillESamplingPFO(TileExt1);
-	FillESamplingPFO(TileExt2);
-
-	FillESamplingPFO(FCAL0);
-	FillESamplingPFO(FCAL1);
-	FillESamplingPFO(FCAL2);
-
-	FillESamplingPFO(MINIFCAL0);
-	FillESamplingPFO(MINIFCAL1);
-	FillESamplingPFO(MINIFCAL2);
-	FillESamplingPFO(MINIFCAL3);	
-
-	emTot += ( E_PreSamplerB+E_EMB1+E_EMB2+E_EMB3+
-		  E_PreSamplerE+E_EME1+E_EME2+E_EME3+
-		  E_FCAL0 );
-	hecTot += ( E_HEC0+E_HEC1+E_HEC2+E_HEC3 );
-	
+        eTot += constit->e();
+        FillESamplingPFO(PreSamplerB);
+        FillESamplingPFO(EMB1);
+        FillESamplingPFO(EMB2);
+        FillESamplingPFO(EMB3);
+
+        FillESamplingPFO(PreSamplerE);
+        FillESamplingPFO(EME1);
+        FillESamplingPFO(EME2);
+        FillESamplingPFO(EME3);
+
+        FillESamplingPFO(HEC0);
+        FillESamplingPFO(HEC1);
+        FillESamplingPFO(HEC2);
+        FillESamplingPFO(HEC3);
+
+        FillESamplingPFO(TileBar0);
+        FillESamplingPFO(TileBar1);
+        FillESamplingPFO(TileBar2);
+
+        FillESamplingPFO(TileGap1);
+        FillESamplingPFO(TileGap2);
+        FillESamplingPFO(TileGap3);
+
+        FillESamplingPFO(TileExt0);
+        FillESamplingPFO(TileExt1);
+        FillESamplingPFO(TileExt2);
+
+        FillESamplingPFO(FCAL0);
+        FillESamplingPFO(FCAL1);
+        FillESamplingPFO(FCAL2);
+
+        FillESamplingPFO(MINIFCAL0);
+        FillESamplingPFO(MINIFCAL1);
+        FillESamplingPFO(MINIFCAL2);
+        FillESamplingPFO(MINIFCAL3);        
+
+        emTot += ( E_PreSamplerB+E_EMB1+E_EMB2+E_EMB3+
+                  E_PreSamplerE+E_EME1+E_EME2+E_EME3+
+                  E_FCAL0 );
+        hecTot += ( E_HEC0+E_HEC1+E_HEC2+E_HEC3 );
+        
       }//only consider neutral PFO
     } else {
       ATH_MSG_WARNING("Tried to call fillEperSamplingPFlow with a jet constituent that is not a PFO!");
@@ -186,4 +192,44 @@ void JetCaloEnergies::fillEperSamplingPFO(const xAOD::Jet & jet, std::vector<flo
   
 }
 
-//**********************************************************************
+void JetCaloEnergies::fillEperSamplingFE(const xAOD::Jet& jet, std::vector<float> & ePerSampling ) const {
+  float emTot = 0.;
+  float hecTot = 0.;
+  float psTot = 0.;
+  float eTot = 0.;
+  size_t numConstit = jet.numConstituents();    
+  for ( size_t i=0; i<numConstit; i++ ) {
+    if(jet.rawConstituent(i)->type()!=xAOD::Type::FlowElement) {
+      ATH_MSG_WARNING("Tried to call fillEperSamplingFE with a jet constituent that is not a FlowElement!");
+      continue;
+    }
+    const xAOD::FlowElement* constit = static_cast<const xAOD::FlowElement*>(jet.rawConstituent(i));
+    if(constit->isCharged()){
+      eTot += constit->chargedObject(0)->e();
+    }
+    else{
+      eTot += constit->e();
+      std::vector<float> constitEPerSampling = FEHelpers::getEnergiesPerSampling(*constit);
+      for ( size_t s = CaloSampling::PreSamplerB; s < CaloSampling::Unknown; s++ ) {
+        ePerSampling[s] += constitEPerSampling[s];
+      }
+      emTot += (constitEPerSampling[CaloSampling::PreSamplerB] + constitEPerSampling[CaloSampling::EMB1]
+              + constitEPerSampling[CaloSampling::EMB2]        + constitEPerSampling[CaloSampling::EMB3]
+              + constitEPerSampling[CaloSampling::PreSamplerE] + constitEPerSampling[CaloSampling::EME1]
+              + constitEPerSampling[CaloSampling::EME2]        + constitEPerSampling[CaloSampling::EME3]
+              + constitEPerSampling[CaloSampling::FCAL0]);
+
+      hecTot += (constitEPerSampling[CaloSampling::HEC0] + constitEPerSampling[CaloSampling::HEC1]
+               + constitEPerSampling[CaloSampling::HEC2] + constitEPerSampling[CaloSampling::HEC3]);
+
+      psTot += (constitEPerSampling[CaloSampling::PreSamplerB] + constitEPerSampling[CaloSampling::PreSamplerE]);
+    }
+  }
+  SG::WriteDecorHandle<xAOD::JetContainer, float> emFracHandle(m_emFracKey);
+  SG::WriteDecorHandle<xAOD::JetContainer, float> hecFracHandle(m_hecFracKey);
+  SG::WriteDecorHandle<xAOD::JetContainer, float> psFracHandle(m_psFracKey);
+  
+  emFracHandle(jet)  = eTot != 0. ? emTot/eTot  : 0.;
+  hecFracHandle(jet) = eTot != 0. ? hecTot/eTot : 0.;
+  psFracHandle(jet)  = eTot != 0. ? psTot/eTot  : 0.;
+}
\ No newline at end of file
diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetECPSFractionTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetECPSFractionTool.cxx
index d0658d332f87a2123860d70eef83b1aa6473edd2..97323b4b73b6af7e9b6c16b565a9e49b56090d42 100644
--- a/Reconstruction/Jet/JetMomentTools/Root/JetECPSFractionTool.cxx
+++ b/Reconstruction/Jet/JetMomentTools/Root/JetECPSFractionTool.cxx
@@ -8,11 +8,11 @@
 #include "JetMomentTools/JetECPSFractionTool.h"
 #include "xAODCaloEvent/CaloCluster.h"
 #include "xAODPFlow/PFO.h"
+#include "xAODPFlow/FlowElement.h"
 #include "CaloGeoHelpers/CaloSampling.h"
 
 using xAOD::IParticle;
 using xAOD::CaloCluster;
-using xAOD::PFO;
 using xAOD::JetConstituent;
 using xAOD::JetFourMom_t;
 
@@ -62,11 +62,16 @@ double JetECPSFractionTool::energyFraction(const xAOD::Jet& jet) const {
       continue;
     }
     const IParticle* ppar = pcon->rawConstituent();
+    if(ppar == nullptr){
+      ATH_MSG_WARNING("Constituent has no associated raw constituent.");
+      continue;
+    }
     double conEnergy = pcon->e();
+
     // If constituent is pflow, then look to see if it is cluster-based.
-    const PFO* ppfl = dynamic_cast<const PFO*>(ppar);
     const CaloCluster* pclu = nullptr;
-    if ( ppfl != nullptr ) {
+    if ( ppar->type() == xAOD::Type::ParticleFlow){
+      const xAOD::PFO* ppfl = static_cast<const xAOD::PFO*>(ppar);
       if ( !ppfl->isCharged() ) {
         pclu = ppfl->cluster(0);  // Assume PFO has exactly one cluster.
         if ( pclu != nullptr ) ATH_MSG_VERBOSE("  Constituent is a PFO pointing to a cluster");
@@ -74,8 +79,22 @@ double JetECPSFractionTool::energyFraction(const xAOD::Jet& jet) const {
       } else {
         ATH_MSG_VERBOSE("  Constituent is a PFO pointing to a track");
       }
-    // Otherwise check if constituent is cluster.
-    } else {
+    }
+    else if ( ppar->type() == xAOD::Type::FlowElement){
+      const xAOD::FlowElement* pfe = static_cast<const xAOD::FlowElement*>(ppar);
+      if(!(pfe->signalType() & xAOD::FlowElement::PFlow)){
+        ATH_MSG_WARNING("  Constituent is a FlowElement but not a PFO. This isn't supported; skipping.");
+        continue;
+      }
+      if ( !pfe->isCharged() ) {
+        pclu = dynamic_cast<const CaloCluster*>(pfe->otherObject(0));  // Assume PFO has exactly one cluster.
+        if ( pclu != nullptr ) ATH_MSG_VERBOSE("  Constituent is a PFO pointing to a cluster");
+        else ATH_MSG_WARNING("  Constituent is a non-charge PFO but no cluster is found.");
+      }
+      else ATH_MSG_VERBOSE("  Constituent is a charged PFO");
+    }
+    else {
+      // Not a PFO. Check if constituent is a cluster.
       pclu = dynamic_cast<const CaloCluster*>(ppar);
       if ( pclu != nullptr ) ATH_MSG_VERBOSE("  Constituent is a cluster");
     }
diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetPtAssociationTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetPtAssociationTool.cxx
index a30e24d34bdd5f1219dd0825d9042ac8d36eecce..f7e6cca9ae07b6d63ede782bfba1c37354ad4e2d 100644
--- a/Reconstruction/Jet/JetMomentTools/Root/JetPtAssociationTool.cxx
+++ b/Reconstruction/Jet/JetMomentTools/Root/JetPtAssociationTool.cxx
@@ -31,12 +31,14 @@ StatusCode JetPtAssociationTool::initialize() {
     ATH_MSG_ERROR("JetPtAssociationTool needs to have its input jet container configured!");
     return StatusCode::FAILURE;
   }
-  m_assocFracKey = m_jetContainerName.key() + "." + m_aname + m_assocFracKey.key();
-  m_assocLinkKey = m_jetContainerName.key() + "." + m_aname + m_assocLinkKey.key();
+  m_assocFracKey = m_jetContainerName + "." + m_aname + m_assocFracKey.key();
+  m_assocLinkKey = m_jetContainerName + "." + m_aname + m_assocLinkKey.key();
 
   ATH_CHECK(m_assocFracKey.initialize());
   ATH_CHECK(m_assocLinkKey.initialize());
 
+  ATH_CHECK(m_matchingContainerName.initialize());
+
   return StatusCode::SUCCESS;
 
 }
@@ -59,10 +61,10 @@ StatusCode JetPtAssociationTool::decorate(const xAOD::JetContainer& jets) const
     }
 
     // Retrieve the container of jets to be matched.
-    SG::ReadHandle<xAOD::JetContainer> handle (m_jetContainerName);
+    SG::ReadHandle<xAOD::JetContainer> handle (m_matchingContainerName);
     if (!handle.isValid()){
       ATH_MSG_WARNING("Matching jet container not found: "
-                      << m_jetContainerName); 
+                      << m_matchingContainerName); 
       return StatusCode::FAILURE;
     }
 
diff --git a/Reconstruction/Jet/JetMomentTools/Root/JetTrackMomentsTool.cxx b/Reconstruction/Jet/JetMomentTools/Root/JetTrackMomentsTool.cxx
index 8b953bf3638d1a8206b8fe5aff6ffd089f67b263..c5ac221efcd8cb406456559546c76a8cdcc18375 100644
--- a/Reconstruction/Jet/JetMomentTools/Root/JetTrackMomentsTool.cxx
+++ b/Reconstruction/Jet/JetMomentTools/Root/JetTrackMomentsTool.cxx
@@ -7,6 +7,7 @@
 #include <sstream>
 #include "JetUtils/JetDistances.h"
 #include "xAODPFlow/PFO.h"
+#include "xAODPFlow/FlowElement.h"
 #include "AsgDataHandles/WriteDecorHandle.h"
 
 JetTrackMomentsTool::JetTrackMomentsTool(const std::string& name)
@@ -105,6 +106,20 @@ StatusCode JetTrackMomentsTool::decorate(const xAOD::JetContainer& jets) const {
         }// We have a charged PFO
       }// Loop on jet constituents
     }// This jet is made from xAOD::PFO, so we do calculate the pflow moments
+    else if (ctype  == xAOD::Type::FlowElement) {
+      size_t numConstit = jet->numConstituents();
+      for ( size_t i=0; i<numConstit; i++ ) {
+        const xAOD::FlowElement* constit = dynamic_cast<const xAOD::FlowElement*>(jet->rawConstituent(i));
+        if(constit != nullptr && (constit->signalType() & xAOD::FlowElement::PFlow)){
+          isPFlowJet = true;
+          if (constit->isCharged()){
+            const xAOD::TrackParticle *thisTrack = dynamic_cast<const xAOD::TrackParticle*>(constit->chargedObject(0));//PFO should have only 1 track
+            if(thisTrack != nullptr) pflowTracks.push_back(thisTrack);
+            else ATH_MSG_WARNING("Charged PFO had no associated TrackParticle");
+          }// We have a charged PFO
+        }// The FlowElement is a PFO
+      }// Loop on jet constituents
+    }// This jet is made from xAOD::FlowElement, so we calculate the pflow moments if they're PFOs
 
     // For each track cut, compute and set the associated moments
     for (size_t iCut = 0; iCut < m_minTrackPt.size(); ++iCut) {
diff --git a/Reconstruction/Jet/JetMomentTools/python/JetMomentToolsConfig.py b/Reconstruction/Jet/JetMomentTools/python/JetMomentToolsConfig.py
index 225ff9f519ce6793e44a5d3f63a418c27cb7a1e4..dfab14517e76287729c34ebab17be64b252d27bd 100644
--- a/Reconstruction/Jet/JetMomentTools/python/JetMomentToolsConfig.py
+++ b/Reconstruction/Jet/JetMomentTools/python/JetMomentToolsConfig.py
@@ -34,7 +34,7 @@ def getEMScaleMomTool(jetdef):
     if jetdef.inputdef.basetype==xAODType.CaloCluster:
         builtFromEMClusters = jetdef.inputdef.inputname in ["CaloCalTopoClusters","HLT_CaloTopoClustersFS"] and jetdef.inputdef.modifiers==["EM"]
         useUncalibConstits = not builtFromEMClusters
-    elif jetdef.inputdef.basetype==xAODType.ParticleFlow:
+    elif (jetdef.inputdef.basetype==xAODType.ParticleFlow or jetdef.inputdef.basetype==xAODType.FlowElement):
         useUncalibConstits = True
     else:
         raise ValueError("EM scale momentum not defined for input type {}".format(jetdef.inputdef.basetype))
diff --git a/Reconstruction/Jet/JetRec/JetRec/PseudoJetGetter.h b/Reconstruction/Jet/JetRec/JetRec/PseudoJetGetter.h
index eebbc0592fb27c0bb2f328b90f7ec180b4ef7d86..35a9d447cb14e97abdee68f75b7f33c9ef3ca012 100644
--- a/Reconstruction/Jet/JetRec/JetRec/PseudoJetGetter.h
+++ b/Reconstruction/Jet/JetRec/JetRec/PseudoJetGetter.h
@@ -22,6 +22,7 @@
 #ifndef GENERATIONBASE
 #include "xAODCaloEvent/CaloCluster.h"
 #include "xAODPFlow/PFO.h"
+#include "xAODPFlow/FlowElement.h"
 #endif
 
 namespace PseudoJetGetter {
@@ -116,6 +117,15 @@ namespace PseudoJetGetter {
     }
 
     bool operator()(const xAOD::IParticle* ip){
+
+      if(ip->type() == xAOD::Type::FlowElement){
+        const xAOD::FlowElement* pfo = dynamic_cast<const xAOD::FlowElement*>(ip);
+        if( pfo->isCharged() ){
+          const static SG::AuxElement::ConstAccessor<char> PVMatchedAcc("matchedToPV");
+	        return !PVMatchedAcc(*pfo);
+        }
+        return skipNegativeEnergy && pfo->e()<FLT_MIN;
+      }
     
       const xAOD::PFO* pfo = dynamic_cast<const xAOD::PFO*>(ip);
     
diff --git a/Reconstruction/Jet/JetRec/python/JetAlgorithm.py b/Reconstruction/Jet/JetRec/python/JetAlgorithm.py
index f2764d774b7f3486e920509c5f005e5f70a970b3..e65b1effc910a12f409f68eab31d4eb0848866ff 100644
--- a/Reconstruction/Jet/JetRec/python/JetAlgorithm.py
+++ b/Reconstruction/Jet/JetRec/python/JetAlgorithm.py
@@ -60,6 +60,9 @@ def addJetRecoToAlgSequence(job =None, useTruth =None, eventShapeTools =None,
     "empflow"  : ("EMPFlowEventShape",  jtm.empflowget),
   }
 
+  if jetFlags.usePFlowFE():
+    evsDict["empflow_fe"] = ("EMPFlowFEEventShape",  jtm.empflowget_fe)
+
   if jetFlags.useTracks():
     if jetFlags.useVertices():
       evsDict["emtopo"] = ("EMTopoOriginEventShape",   jtm.emoriginget)
@@ -97,6 +100,7 @@ def addJetRecoToAlgSequence(job =None, useTruth =None, eventShapeTools =None,
 
   ## if jetFlags.useCells():
   ##   ctools += [jtm.missingcells] commented out : incompatible with trigger : ATR-9696
+
   if jetFlags.useTracks:
     ctools += [jtm.tracksel, jtm.trackselloose_trackjets]
     if jetFlags.useVertices:
@@ -144,6 +148,19 @@ def addJetRecoToAlgSequence(job =None, useTruth =None, eventShapeTools =None,
             StreamName = 'StreamAOD'
             )
           postalgs.append(CHSnPFOsThinAlg)
+    if jetFlags.usePFlowFE() and not IsInInputFile("xAOD::FlowElementContainer","CHSFlowElements"):
+      if not hasattr(job,"jetalgCHSPFlowFE"):
+        ctools += [jtm.JetConstitSeq_PFlowCHS_FE]
+        if thinneg:
+          from ThinningUtils.ThinningUtilsConf import ThinNegativeEnergyNeutralPFOsAlg
+          CHSnFEsThinAlg = ThinNegativeEnergyNeutralPFOsAlg(
+            "ThinNegativeEnergyCHSNeutralFEsAlg",
+            NeutralPFOsKey="", # don't do the usual PFOs here
+            NeutralPFOsFEKey = "CHSNeutralFlowElements",
+            ThinNegativeEnergyNeutralPFOs = True,
+            StreamName = 'StreamAOD'
+            )
+          postalgs.append(CHSnFEsThinAlg)
 
   from JetRec.JetRecConf import JetToolRunner
   from JetRec.JetRecConf import JetAlgorithm
diff --git a/Reconstruction/Jet/JetRec/python/JetRecFlags.py b/Reconstruction/Jet/JetRec/python/JetRecFlags.py
index da956e946dc57c959e77487132587af9e81f8dfa..25952538c1391e87b8e187a06f6f72e78a687ade 100644
--- a/Reconstruction/Jet/JetRec/python/JetRecFlags.py
+++ b/Reconstruction/Jet/JetRec/python/JetRecFlags.py
@@ -23,6 +23,7 @@
 #   useVertices - Toggles whether PFlow jet reconstruction makes use of vertex information.
 #   useMuonSegmentss - Muon segemnt association is enabled
 #   usePFlow - PFlow jets and associations are enabled
+#   usePFlowFE - PFlow jets using FlowElements are enabled
 #   useInDetTrackSelection - The inner detector track selection
 #     tool is used. This requires track propagator exist.
 #   jetAODList - The list of jet collections to be written out
@@ -109,6 +110,13 @@ class usePFlow(JobProperty):
   allowedTypes = ['bool']  # type
   StoredValue  = True      # default value
 
+class usePFlowFE(JobProperty):
+  """ If true, pflow objects are present as FlowElements and used in jet reconstruction.
+  """
+  statusOn     = True
+  allowedTypes = ['bool']  # type
+  StoredValue  = False     # default value
+
 class eventShapeTools(JobProperty):
   """ List of event shape tools that should be called to calculate rho.
       Allowed values are "emtopo", "lctopo", "emorig", "lcorig", "empflow", "emcpflow", "lcpflow".
@@ -251,6 +259,7 @@ jobproperties.JetRecFlags.add_JobProperty(useVertices)
 jobproperties.JetRecFlags.add_JobProperty(useInDetTrackSelection)
 jobproperties.JetRecFlags.add_JobProperty(useMuonSegments)
 jobproperties.JetRecFlags.add_JobProperty(usePFlow)
+jobproperties.JetRecFlags.add_JobProperty(usePFlowFE)
 jobproperties.JetRecFlags.add_JobProperty(eventShapeTools)
 jobproperties.JetRecFlags.add_JobProperty(jetAODList)
 jobproperties.JetRecFlags.add_JobProperty(useCells)
diff --git a/Reconstruction/Jet/JetRec/python/JetRecStandard.py b/Reconstruction/Jet/JetRec/python/JetRecStandard.py
index d5c93c1caa6633b3165d96edc23326fddc9ac0ac..2773bc56b02f1af4086e9d1f0512253594b822e1 100644
--- a/Reconstruction/Jet/JetRec/python/JetRecStandard.py
+++ b/Reconstruction/Jet/JetRec/python/JetRecStandard.py
@@ -112,6 +112,8 @@ if jetFlags.eventShapeTools() == None:
     jetFlags.eventShapeTools += ['emtopo', 'lctopo']
   if jetFlags.usePFlow():
     jetFlags.eventShapeTools += ['empflow']
+  if jetFlags.usePFlowFE():
+    jetFlags.eventShapeTools += ['empflow_fe']
 
 # Import the jet tool manager.
 from JetRec.JetRecStandardToolManager import jtm
diff --git a/Reconstruction/Jet/JetRec/python/JetRecStandardToolManager.py b/Reconstruction/Jet/JetRec/python/JetRecStandardToolManager.py
index 5abe0c060d0429f5f1f10fe6725e593beef2b86c..971154b1de1448b329c7b0d889f4682aac71e295 100644
--- a/Reconstruction/Jet/JetRec/python/JetRecStandardToolManager.py
+++ b/Reconstruction/Jet/JetRec/python/JetRecStandardToolManager.py
@@ -99,6 +99,7 @@ def filterout(skiptoolnames, tools):
 
 # Pseudojet getters
 empfgetters =  [jtm.empflowget]
+empfgetters_fe = [jtm.empflowget_fe]
 
 trackgetters = [jtm.trackget]
 # Add track ghosts
@@ -114,11 +115,13 @@ if jetFlags.useTracks():
   emgetters += [jtm.gtrackget]
   lcgetters += [jtm.gtrackget]
   empfgetters += [jtm.gtrackget]
+  empfgetters_fe += [jtm.gtrackget]
 
 if jetFlags.useMuonSegments():
   emgetters += [jtm.gmusegget]
   lcgetters += [jtm.gmusegget]
   empfgetters  += [jtm.gmusegget]
+  empfgetters_fe  += [jtm.gmusegget]
 # Add jet ghosts.
 if 1:
   for gettername in jetFlags.additionalTopoGetters():
@@ -133,6 +136,7 @@ if jetFlags.useTruth():
   emgetters += [jtm.gtruthget]
   lcgetters += [jtm.gtruthget]
   empfgetters += [jtm.gtruthget]
+  empfgetters_fe += [jtm.gtruthget]
   # Add truth cone matching and truth flavor ghosts.
   flavorgetters = []
   for ptype in jetFlags.truthFlavorTags():
@@ -143,6 +147,7 @@ if jetFlags.useTruth():
   truthwzgetters += flavorgetters
   trackgetters   += flavorgetters
   empfgetters    += flavorgetters
+  empfgetters_fe += flavorgetters
 # Add track jet ghosts.
 if jetFlags.useTracks():
   trackjetgetters = []
@@ -152,12 +157,14 @@ if jetFlags.useTracks():
   emgetters += trackjetgetters
   lcgetters += trackjetgetters
   empfgetters += trackjetgetters
+  empfgetters_fe += trackjetgetters
 
 
 # Add getter lists to jtm indexed by input type name.
 jtm.gettersMap["emtopo"]    = list(emgetters)
 jtm.gettersMap["lctopo"]    = list(lcgetters)
 jtm.gettersMap["empflow"]   = list(empfgetters)
+jtm.gettersMap["empflowfe"] = list(empfgetters_fe)
 jtm.gettersMap["track"]     = list(trackgetters)
 jtm.gettersMap["pv0track"]  = list(trackgetters)
 if jetFlags.useTruth():
@@ -167,6 +174,7 @@ if jetFlags.useTruth():
 jtm.gettersMap["emtopo_reduced"]  = filterout(["gakt2trackget","gakt4trackget"],emgetters)
 jtm.gettersMap["lctopo_reduced"]  = filterout(["gakt2trackget","gakt4trackget"],lcgetters)
 jtm.gettersMap["empflow_reduced"] = filterout(["gakt2trackget","gakt4trackget"],empfgetters)
+jtm.gettersMap["empflowfe_reduced"] = filterout(["gakt2trackget","gakt4trackget"],empfgetters_fe)
 
 
 #########################################################
diff --git a/Reconstruction/Jet/JetRec/python/JetRecStandardTools.py b/Reconstruction/Jet/JetRec/python/JetRecStandardTools.py
index 8e764915b017108d9f30fed2a138f5138b387244..954540ed1059119ad26602d40014d243ce6914a5 100644
--- a/Reconstruction/Jet/JetRec/python/JetRecStandardTools.py
+++ b/Reconstruction/Jet/JetRec/python/JetRecStandardTools.py
@@ -332,15 +332,29 @@ ctm.add( CorrectPFOTool("CorrectPFOTool",
                         ),
          alias = 'correctPFO' )
 
+ctm.add( CorrectPFOTool("CorrectPFOTool_FE",
+                        WeightPFOTool = jtm.pflowweighter,
+                        InputIsEM = True,
+                        CalibratePFO = False,
+                        UseChargedWeights = True,
+                        InputType = xAODType.FlowElement
+                        ),
+         alias = 'correctPFO_FE' )
+
 # this removes (weights momenta to 0) charged PFOs from non-hard-scatter vertices
 ctm.add( ChargedHadronSubtractionTool("CHSTool", InputType = xAODType.ParticleFlow),
          alias = 'chsPFO' )
 
+ctm.add( ChargedHadronSubtractionTool("CHSTool_FE", InputType = xAODType.FlowElement),
+         alias = 'chsPFO_FE' )
+
 # Options to disable dependence on primary vertex container
 # for PFO corrections (e.g. when running cosmics)
 if not (jetFlags.useTracks and jetFlags.useVertices):
   ctm.modifiersMap['correctPFO'].CorrectNeutral=False
+  ctm.modifiersMap['correctPFO_FE'].CorrectNeutral=False
   ctm.modifiersMap['chsPFO'].IgnoreVertex=True
+  ctm.modifiersMap['chsPFO_FE'].IgnoreVertex=True
 
 # Run the above tools to modify PFO
 jtm += ctm.buildConstitModifSequence( "JetConstitSeq_PFlowCHS",
@@ -348,6 +362,12 @@ jtm += ctm.buildConstitModifSequence( "JetConstitSeq_PFlowCHS",
                                       OutputContainer = "CHS",  #"ParticleFlowObjects" will be appended later
                                       modList = ['correctPFO', 'chsPFO'] )
 
+jtm += ctm.buildConstitModifSequence( "JetConstitSeq_PFlowCHS_FE",
+                                      InputContainer = "JetETMiss",
+                                      OutputContainer = "CHS",  #"FlowElements" will be appended later
+                                      modList = ['correctPFO_FE', 'chsPFO_FE'],
+                                      InputType = xAODType.FlowElement )
+
 # EM-scale pflow.
 jtm += PseudoJetAlgorithm(
   "empflowget",
@@ -357,6 +377,14 @@ jtm += PseudoJetAlgorithm(
   SkipNegativeEnergy = True,
 )
 
+jtm += PseudoJetAlgorithm(
+  "empflowget_fe",
+  Label = "EMPFlowFE",
+  InputContainer = "CHSFlowElements",
+  OutputContainer = "PseudoJetEMPFlowFE",
+  SkipNegativeEnergy = True,
+)
+
 # AntiKt2 track jets.
 jtm += PseudoJetAlgorithm(
   "gakt2trackget", # give a unique name
diff --git a/Reconstruction/Jet/JetRec/python/JetToolSupport.py b/Reconstruction/Jet/JetRec/python/JetToolSupport.py
index a3f8defa6bdd4bbbb1c73353b66a6929219dc35e..16864e40eedd70ac80e9d87c981e39de31e63482 100644
--- a/Reconstruction/Jet/JetRec/python/JetToolSupport.py
+++ b/Reconstruction/Jet/JetRec/python/JetToolSupport.py
@@ -246,7 +246,7 @@ class JetToolManager:
         tname = mod + "_" + salg + srad
         if not tname in self.tools:
           from JetMomentTools.JetMomentToolsConf import JetPtAssociationTool
-          self += JetPtAssociationTool(tname, InputContainer=cname, AssociationName="GhostTruth")
+          self += JetPtAssociationTool(tname, JetContainer=output, MatchingJetContainer=cname, AssociationName="GhostTruth")
         outmods += [self.tools[tname]]
       # trackassoc - Does track jet association replacing the input name with "Track"
       elif mod == "trackassoc":
@@ -265,7 +265,7 @@ class JetToolManager:
           tname = mod + "_" + salg + srad
           if not tname in self.tools:
             from JetMomentTools.JetMomentToolsConf import JetPtAssociationTool
-            self += JetPtAssociationTool(tname, InputContainer=cname, AssociationName="GhostTrack")
+            self += JetPtAssociationTool(tname, JetContainer=output, MatchingJetContainer=cname, AssociationName="GhostTrack")
           outmods += [self.tools[tname]]
       # jetfilter - Filter to remove jets with pT < self.ptminFilter
       elif mod == "jetfilter":
diff --git a/Reconstruction/Jet/JetRecTools/JetRecTools/ChargedHadronSubtractionTool.h b/Reconstruction/Jet/JetRecTools/JetRecTools/ChargedHadronSubtractionTool.h
index f99efdf2d9e442dc71257b5173148fd2d18e596c..fd335f874706fd8368f3985f642961354437a1ba 100644
--- a/Reconstruction/Jet/JetRecTools/JetRecTools/ChargedHadronSubtractionTool.h
+++ b/Reconstruction/Jet/JetRecTools/JetRecTools/ChargedHadronSubtractionTool.h
@@ -22,6 +22,7 @@
 #include "JetEDM/TrackVertexAssociation.h"
 #include "xAODTracking/VertexContainer.h" 
 #include "xAODPFlow/PFOContainer.h"
+#include "xAODPFlow/FlowElementContainer.h"
 
 class ChargedHadronSubtractionTool : public JetConstituentModifierBase{
   ASG_TOOL_CLASS(ChargedHadronSubtractionTool, IJetConstituentModifier)
@@ -38,6 +39,7 @@ class ChargedHadronSubtractionTool : public JetConstituentModifierBase{
   StatusCode process_impl(xAOD::IParticleContainer* cont) const; 
   // Type-specific operation
   StatusCode matchToPrimaryVertex(xAOD::PFOContainer& cont) const;
+  StatusCode matchToPrimaryVertex(xAOD::FlowElementContainer& cont) const;
 
   const xAOD::Vertex* getPrimaryVertex() const;
   bool m_useTrackToVertexTool;
diff --git a/Reconstruction/Jet/JetRecTools/JetRecTools/CorrectPFOTool.h b/Reconstruction/Jet/JetRecTools/JetRecTools/CorrectPFOTool.h
index a4bd5718dbd8200bb62d71a4394fb99f8069a128..9397e7366f6a9f8ebd3544d53f372f0ab7c74f73 100644
--- a/Reconstruction/Jet/JetRecTools/JetRecTools/CorrectPFOTool.h
+++ b/Reconstruction/Jet/JetRecTools/JetRecTools/CorrectPFOTool.h
@@ -22,6 +22,7 @@
 #include "JetRecTools/JetConstituentModifierBase.h"
 #include "xAODTracking/VertexContainer.h" 
 #include "xAODPFlow/PFOContainer.h"
+#include "xAODPFlow/FlowElementContainer.h"
 #include "AsgTools/ToolHandle.h"
 #include "PFlowUtils/IWeightPFOTool.h"
 
@@ -41,10 +42,14 @@ class CorrectPFOTool : public JetConstituentModifierBase{
   StatusCode process_impl(xAOD::IParticleContainer* cont) const override final;
   // Type-specific operation
   StatusCode correctPFO(xAOD::PFOContainer& cont) const;
+  StatusCode correctPFO(xAOD::FlowElementContainer& cont) const;
+
 
   const xAOD::Vertex* getPrimaryVertex() const;
   StatusCode applyNeutralCorrection(xAOD::PFO& pfo, const xAOD::Vertex& vtx) const;
+  StatusCode applyNeutralCorrection(xAOD::FlowElement& pfo, const xAOD::Vertex& vtx) const;
   StatusCode applyChargedCorrection(xAOD::PFO& pfo) const;
+  StatusCode applyChargedCorrection(xAOD::FlowElement& pfo) const;
 
   bool m_inputIsEM;   /// If true EM clusters are used for neutral PFOs.
   bool m_calibrate;   /// If true, EM PFOs are calibrated to LC.
diff --git a/Reconstruction/Jet/JetRecTools/JetRecTools/JetConstituentModSequence.h b/Reconstruction/Jet/JetRecTools/JetRecTools/JetConstituentModSequence.h
index d9bd8670609ee82fee74c491e2e32bdf36d602b0..e2517af01afdce44f1f637a184edaaa9eff46aa1 100644
--- a/Reconstruction/Jet/JetRecTools/JetRecTools/JetConstituentModSequence.h
+++ b/Reconstruction/Jet/JetRecTools/JetRecTools/JetConstituentModSequence.h
@@ -25,6 +25,7 @@
 #include "xAODPFlow/TrackCaloClusterContainer.h"
 #include "xAODTracking/TrackParticleContainer.h"
 #include "xAODPFlow/PFOContainer.h"
+#include "xAODPFlow/FlowElementContainer.h"
 
 #include "AthenaMonitoringKernel/GenericMonitoringTool.h"
 
@@ -73,7 +74,16 @@ protected:
 
   SG::WriteHandleKey<xAOD::PFOContainer> m_outAllPFOKey{this, "OutAllPFOKey", "", "WriteHandleKey for all modified PFlow Objects"};
 
+  SG::ReadHandleKey<xAOD::FlowElementContainer> m_inChargedFEKey{this, "InChargedFEKey", "", "ReadHandleKey for modified Charged FlowElements"};
+  SG::WriteHandleKey<xAOD::FlowElementContainer> m_outChargedFEKey{this, "OutChargedFEKey", "", "WriteHandleKey for modified Charged FlowElements"};
+
+  SG::ReadHandleKey<xAOD::FlowElementContainer> m_inNeutralFEKey{this, "InNeutralFEKey", "", "ReadHandleKey for modified Neutral FlowElements"};
+  SG::WriteHandleKey<xAOD::FlowElementContainer> m_outNeutralFEKey{this, "OutNeutralFEKey", "", "WriteHandleKey for modified Neutral FlowElements"};
+
+  SG::WriteHandleKey<xAOD::FlowElementContainer> m_outAllFEKey{this, "OutAllFEKey", "", "WriteHandleKey for all modified FlowElements"};
+
   StatusCode copyModRecordPFO() const;
+  StatusCode copyModRecordFE() const;
 
   /// helper function to cast, shallow copy and record a container.
 
diff --git a/Reconstruction/Jet/JetRecTools/JetRecTools/PuppiWeightTool.h b/Reconstruction/Jet/JetRecTools/JetRecTools/PuppiWeightTool.h
index 55a740965bf2e7e541dfcdc084245109147372d7..47b5931afc9510653ac82d05112c61b233a8fdea 100644
--- a/Reconstruction/Jet/JetRecTools/JetRecTools/PuppiWeightTool.h
+++ b/Reconstruction/Jet/JetRecTools/JetRecTools/PuppiWeightTool.h
@@ -11,6 +11,7 @@
 #include "xAODCaloEvent/CaloClusterContainer.h"
 #include "xAODTracking/VertexContainer.h"
 #include "xAODPFlow/PFOContainer.h"
+#include "xAODPFlow/FlowElementContainer.h"
 #include "JetRecTools/Puppi.h"
 #include <string>
 
@@ -25,7 +26,8 @@ class PuppiWeightTool: public JetConstituentModifierBase {
   StatusCode finalize() override final;
 
   StatusCode process_impl(xAOD::IParticleContainer* cont) const override final;
-  StatusCode applyPuppiWeights(xAOD::PFOContainer* cont) const; 
+  StatusCode applyPuppiWeights(xAOD::PFOContainer* cont) const;
+  StatusCode applyPuppiWeights(xAOD::FlowElementContainer* cont) const;
 
  private:
 
diff --git a/Reconstruction/Jet/JetRecTools/Root/ChargedHadronSubtractionTool.cxx b/Reconstruction/Jet/JetRecTools/Root/ChargedHadronSubtractionTool.cxx
index ca226604232013e7f3173cc5bd246580e4395393..22138763d0c29d125e83cb4cb168ba0a93422ac5 100644
--- a/Reconstruction/Jet/JetRecTools/Root/ChargedHadronSubtractionTool.cxx
+++ b/Reconstruction/Jet/JetRecTools/Root/ChargedHadronSubtractionTool.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 #include "JetRecTools/ChargedHadronSubtractionTool.h"
 
@@ -22,9 +22,9 @@ ChargedHadronSubtractionTool::ChargedHadronSubtractionTool(const std::string& na
 }
 
 StatusCode ChargedHadronSubtractionTool::initialize() {
-  if(m_inputType!=xAOD::Type::ParticleFlow) {
+  if(m_inputType!=xAOD::Type::ParticleFlow && m_inputType!=xAOD::Type::FlowElement) {
     ATH_MSG_ERROR("ChargedHadronSubtractionTool requires PFO inputs. It cannot operate on objects of type "
-		  << m_inputType);
+                  << m_inputType);
     return StatusCode::FAILURE;
   }
   // Only initialise the DataHandle we will use, to avoid superfluous dependencies
@@ -38,8 +38,21 @@ StatusCode ChargedHadronSubtractionTool::initialize() {
 }
 
 StatusCode ChargedHadronSubtractionTool::process_impl(xAOD::IParticleContainer* cont) const {
+  
   // Type-checking happens in the JetConstituentModifierBase class
   // so it is safe just to static_cast
+  if(m_inputType==xAOD::Type::FlowElement){
+    xAOD::FlowElementContainer* feCont = static_cast<xAOD::FlowElementContainer*>(cont);
+
+    // Check that these FlowElements are actually PFOs; nothing else is supported
+    if(!feCont->empty() && !(feCont->front()->signalType() & xAOD::FlowElement::PFlow)){
+      ATH_MSG_ERROR("This tool correctly recieved a FlowElement, but it wasn't a PFO!"
+                 << " Signal type was " << feCont->front()->signalType());
+      return StatusCode::FAILURE;
+    }
+
+    return matchToPrimaryVertex(*feCont);
+  }
   xAOD::PFOContainer* pfoCont = static_cast<xAOD::PFOContainer*> (cont);
   return matchToPrimaryVertex(*pfoCont);
 }
@@ -91,17 +104,17 @@ StatusCode ChargedHadronSubtractionTool::matchToPrimaryVertex(xAOD::PFOContainer
     if(m_useTrackToVertexTool) {
       auto handle = SG::makeHandle(m_trkVtxAssoc_key);
       if(!handle.isValid()){
-	ATH_MSG_ERROR("Can't retrieve TrackVertexAssociation : "<< m_trkVtxAssoc_key.key()); 
-	return StatusCode::FAILURE;
+        ATH_MSG_ERROR("Can't retrieve TrackVertexAssociation : "<< m_trkVtxAssoc_key.key()); 
+        return StatusCode::FAILURE;
       }
       trkVtxAssoc = handle.cptr();
     } else {
       vtx = getPrimaryVertex();
       if(vtx==nullptr) {
-	ATH_MSG_ERROR("Primary vertex container was empty or no valid vertex found!");
-	return StatusCode::FAILURE;
+        ATH_MSG_ERROR("Primary vertex container was empty or no valid vertex found!");
+        return StatusCode::FAILURE;
       } else if (vtx->vertexType()==xAOD::VxType::NoVtx) {
-	ATH_MSG_VERBOSE("No genuine primary vertex found. Will consider all PFOs matched.");
+        ATH_MSG_VERBOSE("No genuine primary vertex found. Will consider all PFOs matched.");
       }
     }
   }
@@ -118,22 +131,85 @@ StatusCode ChargedHadronSubtractionTool::matchToPrimaryVertex(xAOD::PFOContainer
     } else {
       const xAOD::TrackParticle* ptrk = ppfo->track(0);
       if(ptrk==nullptr) {
-	ATH_MSG_WARNING("Charged PFO with index " << ppfo->index() << " has no ID track!");
-	continue;
+        ATH_MSG_WARNING("Charged PFO with index " << ppfo->index() << " has no ID track!");
+        continue;
+      }
+      if(trkVtxAssoc) { // Use TrackVertexAssociation
+        const xAOD::Vertex* thisTracksVertex = trkVtxAssoc->associatedVertex(ptrk);
+        matchedToPrimaryVertex = (xAOD::VxType::PriVtx == thisTracksVertex->vertexType());
+      } else { // Use Primary Vertex
+        if(vtx->vertexType()==xAOD::VxType::NoVtx) { // No reconstructed vertices
+          matchedToPrimaryVertex = true; // simply match all cPFOs in this case
+        } else { // Had a good reconstructed vertex.
+          // vtz.z() provides z of that vertex w.r.t the center of the beamspot (z = 0).
+          // Thus we correct the track z0 to be w.r.t z = 0
+          float z0 = ptrk->z0() + ptrk->vz() - vtx->z();
+          float theta = ptrk->theta();
+          matchedToPrimaryVertex = ( fabs(z0*sin(theta)) < 2.0 );
+        }
+      } // TVA vs PV decision
+    }
+    PVMatchedAcc(*ppfo) = matchedToPrimaryVertex;
+  }
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode ChargedHadronSubtractionTool::matchToPrimaryVertex(xAOD::FlowElementContainer& cont) const {
+  const static SG::AuxElement::Accessor<char> PVMatchedAcc("matchedToPV");
+
+  // Use only one of TVA or PV
+  const jet::TrackVertexAssociation* trkVtxAssoc = nullptr;
+  const xAOD::Vertex* vtx = nullptr;
+  if(!m_ignoreVertex) {
+    // In cosmics, there's no PV container so we need to avoid attempting
+    // to retrieve anything related to it
+    if(m_useTrackToVertexTool) {
+      auto handle = SG::makeHandle(m_trkVtxAssoc_key);
+      if(!handle.isValid()){
+        ATH_MSG_ERROR("Can't retrieve TrackVertexAssociation : "<< m_trkVtxAssoc_key.key()); 
+        return StatusCode::FAILURE;
+      }
+      trkVtxAssoc = handle.cptr();
+    } else {
+      vtx = getPrimaryVertex();
+      if(vtx==nullptr) {
+        ATH_MSG_ERROR("Primary vertex container was empty or no valid vertex found!");
+        return StatusCode::FAILURE;
+      } else if (vtx->vertexType()==xAOD::VxType::NoVtx) {
+        ATH_MSG_VERBOSE("No genuine primary vertex found. Will consider all PFOs matched.");
+      }
+    }
+  }
+
+  for ( xAOD::FlowElement* ppfo : cont ) {
+    // Ignore neutral PFOs
+    if(!ppfo->isCharged()) continue;
+
+    bool matchedToPrimaryVertex = false;
+    if(m_ignoreVertex) {
+      // If we don't use vertex information, don't bother computing the decision
+      // Just pass every cPFO -- there shouldn't be many in cosmics!
+      matchedToPrimaryVertex = true;
+    } else {
+      const xAOD::TrackParticle* ptrk = dynamic_cast<const xAOD::TrackParticle*>(ppfo->chargedObject(0));
+      if(ptrk==nullptr) {
+        ATH_MSG_WARNING("Charged PFO with index " << ppfo->index() << " has no ID track!");
+        continue;
       }
       if(trkVtxAssoc) { // Use TrackVertexAssociation
-	const xAOD::Vertex* thisTracksVertex = trkVtxAssoc->associatedVertex(ptrk);
-	matchedToPrimaryVertex = (xAOD::VxType::PriVtx == thisTracksVertex->vertexType());
+        const xAOD::Vertex* thisTracksVertex = trkVtxAssoc->associatedVertex(ptrk);
+        matchedToPrimaryVertex = (xAOD::VxType::PriVtx == thisTracksVertex->vertexType());
       } else { // Use Primary Vertex
-	if(vtx->vertexType()==xAOD::VxType::NoVtx) { // No reconstructed vertices
-	  matchedToPrimaryVertex = true; // simply match all cPFOs in this case
-	} else { // Had a good reconstructed vertex.
-	  // vtz.z() provides z of that vertex w.r.t the center of the beamspot (z = 0).
-	  // Thus we correct the track z0 to be w.r.t z = 0
-	  float z0 = ptrk->z0() + ptrk->vz() - vtx->z();
-	  float theta = ptrk->theta();
-	  matchedToPrimaryVertex = ( fabs(z0*sin(theta)) < 2.0 );
-	}
+        if(vtx->vertexType()==xAOD::VxType::NoVtx) { // No reconstructed vertices
+          matchedToPrimaryVertex = true; // simply match all cPFOs in this case
+        } else { // Had a good reconstructed vertex.
+          // vtz.z() provides z of that vertex w.r.t the center of the beamspot (z = 0).
+          // Thus we correct the track z0 to be w.r.t z = 0
+          float z0 = ptrk->z0() + ptrk->vz() - vtx->z();
+          float theta = ptrk->theta();
+          matchedToPrimaryVertex = ( fabs(z0*sin(theta)) < 2.0 );
+        }
       } // TVA vs PV decision
     }
     PVMatchedAcc(*ppfo) = matchedToPrimaryVertex;
diff --git a/Reconstruction/Jet/JetRecTools/Root/ConstitTimeCutTool.cxx b/Reconstruction/Jet/JetRecTools/Root/ConstitTimeCutTool.cxx
index e3cbee3fda6622add3186c232095bd6bd766e6fd..8aba5151e96e3a4bbcddbf0b5c3308cf42613f6d 100644
--- a/Reconstruction/Jet/JetRecTools/Root/ConstitTimeCutTool.cxx
+++ b/Reconstruction/Jet/JetRecTools/Root/ConstitTimeCutTool.cxx
@@ -6,7 +6,9 @@
 #include "JetRecTools/ConstitTimeCutTool.h"
 #include "xAODCaloEvent/CaloClusterContainer.h"
 #include "xAODPFlow/PFOContainer.h"
+#include "xAODPFlow/FlowElementContainer.h"
 #include "xAODPFlow/PFODefs.h"
+#include "PFlowUtils/FEHelpers.h"
 
 using namespace std;
 
@@ -30,9 +32,9 @@ StatusCode ConstitTimeCutTool::initialize() {
     }
 
   } else {
-    if(m_inputType!=xAOD::Type::CaloCluster) {
+    if(m_inputType!=xAOD::Type::CaloCluster && m_inputType!= xAOD::Type::FlowElement) {
       ATH_MSG_ERROR("Incompatible configuration: ConstitTimeCutTool is not specialised for inputs of type "
-		    << m_inputType);
+                    << m_inputType);
       return StatusCode::FAILURE;
     }
   }
@@ -59,21 +61,41 @@ StatusCode ConstitTimeCutTool::process_impl(xAOD::IParticleContainer* cont) cons
     {
       xAOD::PFOContainer* pfos = static_cast<xAOD::PFOContainer*> (cont);
       for(xAOD::PFO* pfo : *pfos) {
-	if(fabs(pfo->charge())<FLT_MIN || m_applyToChargedPFO) { // only apply to neutrals if m_applyToChargedPFO is false. If m_applyToChargedPFO is true then apply to all POs.
-	  float time(0.);
-	  float quality(0.);
-	  float lambda_center(0.);
-	  // Only apply cut if retrieval succeeded, else warn
-	  if(pfo->attribute(xAOD::PFODetails::eflowRec_TIMING,time) &&
+        if(fabs(pfo->charge())<FLT_MIN || m_applyToChargedPFO) { // only apply to neutrals if m_applyToChargedPFO is false. If m_applyToChargedPFO is true then apply to all POs.
+          float time(0.);
+          float quality(0.);
+          float lambda_center(0.);
+          // Only apply cut if retrieval succeeded, else warn
+          if(pfo->attribute(xAOD::PFODetails::eflowRec_TIMING,time) &&
              pfo->attribute(xAOD::PFODetails::eflowRec_AVG_LAR_Q,quality) &&
              pfo->attribute(xAOD::PFODetails::eflowRec_CENTER_LAMBDA,lambda_center)
             ) {
             //quality is on [0,2^16-1] scale
-	    ATH_CHECK( applyTimingCut(pfo, time, quality/65535, lambda_center) );
-	  } else {
-	    ATH_MSG_WARNING("Failed to retrieve the PFO informations necessary for timing cut at PFO #" << pfo->index());
-	  }
-	}
+            ATH_CHECK( applyTimingCut(pfo, time, quality/65535, lambda_center) );
+          } else {
+            ATH_MSG_WARNING("Failed to retrieve the PFO informations necessary for timing cut at PFO #" << pfo->index());
+          }
+        }
+      }
+    }
+    break;
+  case xAOD::Type::FlowElement:
+    {
+      xAOD::FlowElementContainer* fes = static_cast<xAOD::FlowElementContainer*>(cont);
+      if(!fes->empty() && !(fes->front()->signalType() & xAOD::FlowElement::PFlow)){
+        ATH_MSG_ERROR("ConstitTimeCutTool received FlowElements that aren't PFOs, this isn't supported!");
+        return StatusCode::FAILURE;
+      }
+      for(xAOD::FlowElement* fe : *fes){
+        if(!fe->isCharged() || m_applyToChargedPFO){
+          const static SG::AuxElement::ConstAccessor<float> acc_timing("TIMING");
+          const static SG::AuxElement::ConstAccessor<float> acc_larq("AVG_LAR_Q");
+          const static SG::AuxElement::ConstAccessor<float> acc_clambda("CENTER_LAMBDA");
+          float time = acc_timing(*fe);
+          float quality = acc_larq(*fe);
+          float lambda_center = acc_clambda(*fe);
+          ATH_CHECK( applyTimingCut(fe, time, quality/65535, lambda_center) );
+        }
       }
     }
     break;
@@ -82,11 +104,9 @@ StatusCode ConstitTimeCutTool::process_impl(xAOD::IParticleContainer* cont) cons
     ATH_MSG_ERROR("No specialisation for object type " << m_inputType);
     return StatusCode::FAILURE;
   }
-
   return StatusCode::SUCCESS;
 }
 
-
 StatusCode ConstitTimeCutTool::applyTimingCut(xAOD::IParticle* part, const float& time, const float& quality, const float& lambda_center) const {
   if(abs( part->eta() ) < m_etaMax){
 
@@ -101,4 +121,4 @@ StatusCode ConstitTimeCutTool::applyTimingCut(xAOD::IParticle* part, const float
   }
 
   return StatusCode::SUCCESS;
-}
+}
\ No newline at end of file
diff --git a/Reconstruction/Jet/JetRecTools/Root/ConstituentSubtractorTool.cxx b/Reconstruction/Jet/JetRecTools/Root/ConstituentSubtractorTool.cxx
index 63f69e7f8abb82d3ab58055f37caff8de2c69f53..37f16de33d8304e72aab00782f6d03edf06a0cd3 100644
--- a/Reconstruction/Jet/JetRecTools/Root/ConstituentSubtractorTool.cxx
+++ b/Reconstruction/Jet/JetRecTools/Root/ConstituentSubtractorTool.cxx
@@ -14,6 +14,7 @@
 
 #include "xAODPFlow/PFO.h"
 #include "xAODPFlow/TrackCaloCluster.h"
+#include "xAODPFlow/FlowElement.h"
 
 using namespace fastjet;
 ConstituentSubtractorTool::ConstituentSubtractorTool(const std::string & name): JetConstituentModifierBase(name) {
@@ -73,6 +74,12 @@ StatusCode ConstituentSubtractorTool::process_impl(xAOD::IParticleContainer* con
       xAOD::TrackCaloCluster* tcc = static_cast<xAOD::TrackCaloCluster*>(part);
       accept &= (tcc->taste()!= 0);
     }
+    if(m_inputType==xAOD::Type::FlowElement){
+      xAOD::FlowElement* fe = static_cast<xAOD::FlowElement*>(part);
+      if(((fe->signalType() & xAOD::FlowElement::PFlow) && m_ignoreChargedPFOs) ||
+          (fe->signalType() & xAOD::FlowElement::TCC))
+        accept &= !(fe->isCharged());
+    }
     // Reject object if outside maximum eta range
     accept &= fabs(part->eta()) <= m_maxEta;
 
diff --git a/Reconstruction/Jet/JetRecTools/Root/CorrectPFOTool.cxx b/Reconstruction/Jet/JetRecTools/Root/CorrectPFOTool.cxx
index 5ba17ec0765cd839d2626e576944dbb2047c4d9d..941339c0d10f53f978a00d592542f3cadfd5d34d 100644
--- a/Reconstruction/Jet/JetRecTools/Root/CorrectPFOTool.cxx
+++ b/Reconstruction/Jet/JetRecTools/Root/CorrectPFOTool.cxx
@@ -5,6 +5,7 @@
 // CorrectPFOTool.cxx
 
 #include "JetRecTools/CorrectPFOTool.h"
+#include "PFlowUtils/FEHelpers.h"
 #include <cmath>
 
 CorrectPFOTool::CorrectPFOTool(const std::string &name):
@@ -30,9 +31,9 @@ StatusCode CorrectPFOTool::initialize() {
     ATH_MSG_ERROR("This tool is configured to do nothing!");
     return StatusCode::FAILURE;
   }
-  if(m_inputType!=xAOD::Type::ParticleFlow) {
+  if(m_inputType!=xAOD::Type::ParticleFlow && m_inputType!=xAOD::Type::FlowElement) {
     ATH_MSG_ERROR("CorrectPFOTool requires PFO inputs. It cannot operate on objects of type "
-		  << m_inputType);
+                  << m_inputType);
     return StatusCode::FAILURE;
   }
   if(m_useChargedWeights) {
@@ -45,6 +46,14 @@ StatusCode CorrectPFOTool::initialize() {
 StatusCode CorrectPFOTool::process_impl(xAOD::IParticleContainer* cont) const {
   // Type-checking happens in the JetConstituentModifierBase class
   // so it is safe just to static_cast
+  if(m_inputType == xAOD::Type::FlowElement){
+    xAOD::FlowElementContainer* feCont = static_cast<xAOD::FlowElementContainer*>(cont);
+    if(!feCont->empty() && !(feCont->front()->signalType() & xAOD::FlowElement::PFlow)){
+      ATH_MSG_ERROR("CorrectPFOTool received FlowElements that aren't PFOs");
+      return StatusCode::FAILURE;
+    }
+    return correctPFO(*feCont);
+  }
   xAOD::PFOContainer* pfoCont = static_cast<xAOD::PFOContainer*> (cont);
   return correctPFO(*pfoCont);
 }
@@ -101,11 +110,40 @@ StatusCode CorrectPFOTool::correctPFO(xAOD::PFOContainer& cont) const {
 
     if ( fabs(ppfo->charge())<FLT_MIN) { // Neutral PFOs
       if(m_correctneutral) {
-	ATH_CHECK( applyNeutralCorrection(*ppfo, *vtx) );
+        ATH_CHECK( applyNeutralCorrection(*ppfo, *vtx) );
       }
     } else { // Charged PFOs
       if(m_correctcharged) {
-	ATH_CHECK( applyChargedCorrection(*ppfo) );
+        ATH_CHECK( applyChargedCorrection(*ppfo) );
+      }      
+    }
+  } // PFO loop
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode CorrectPFOTool::correctPFO(xAOD::FlowElementContainer& cont) const { 
+
+  const xAOD::Vertex* vtx = nullptr;
+  if(m_correctneutral) {
+    vtx = getPrimaryVertex();
+    if(vtx==nullptr) {
+      ATH_MSG_ERROR("Primary vertex container was empty or no valid vertex found!");
+      return StatusCode::FAILURE;
+    } else if (vtx->vertexType()==xAOD::VxType::NoVtx) {
+      ATH_MSG_VERBOSE("No genuine primary vertex found. Will not apply origin correction");
+    }
+  }
+
+  for ( xAOD::FlowElement* ppfo : cont ) {
+
+    if ( !ppfo->isCharged()) { // Neutral PFOs
+      if(m_correctneutral) {
+        ATH_CHECK( applyNeutralCorrection(*ppfo, *vtx) );
+      }
+    } else { // Charged PFOs
+      if(m_correctcharged) {
+        ATH_CHECK( applyChargedCorrection(*ppfo) );
       }      
     }
   } // PFO loop
@@ -120,18 +158,29 @@ StatusCode CorrectPFOTool::applyNeutralCorrection(xAOD::PFO& pfo, const xAOD::Ve
     if ( !m_inputIsEM || m_calibrate ) { // Use LC four-vector
       // Only correct if we really reconstructed a vertex
       if(vtx.vertexType()==xAOD::VxType::PriVtx) {
-	pfo.setP4(pfo.GetVertexCorrectedFourVec(vtx));
+        pfo.setP4(pfo.GetVertexCorrectedFourVec(vtx));
       }
     } else { // Use EM four-vector
       // Only apply origin correction if we really reconstructed a vertex
       if(vtx.vertexType()==xAOD::VxType::PriVtx) {
-	pfo.setP4(pfo.GetVertexCorrectedEMFourVec(vtx));
+        pfo.setP4(pfo.GetVertexCorrectedEMFourVec(vtx));
       } else {pfo.setP4(pfo.p4EM());} // Just set EM 4-vec
     }
   }
   return StatusCode::SUCCESS;
 }
 
+StatusCode CorrectPFOTool::applyNeutralCorrection(xAOD::FlowElement& pfo, const xAOD::Vertex& vtx) const {
+  if (pfo.e() < FLT_MIN) { //This is necessary to avoid changing sign of pT for pT<0 PFO
+    pfo.setP4(0,0,0,0);
+  }
+  // Only apply origin correction if we really reconstructed a vertex
+  else if(vtx.vertexType() == xAOD::VxType::PriVtx) {
+    pfo.setP4(FEHelpers::getVertexCorrectedFourVec(pfo, vtx));
+  }
+  return StatusCode::SUCCESS;
+}
+
 StatusCode CorrectPFOTool::applyChargedCorrection(xAOD::PFO& pfo) const {
   if (m_useChargedWeights) {
     float weight = 0.0;
@@ -142,3 +191,12 @@ StatusCode CorrectPFOTool::applyChargedCorrection(xAOD::PFO& pfo) const {
   return StatusCode::SUCCESS;
 }
 
+StatusCode CorrectPFOTool::applyChargedCorrection(xAOD::FlowElement& pfo) const {
+  if (m_useChargedWeights) {
+    float weight = 0.0;
+    ATH_CHECK( m_weightPFOTool->fillWeight( pfo, weight ) );
+    ATH_MSG_VERBOSE("Fill pseudojet for CPFO with weighted pt: " << pfo.pt()*weight);
+    pfo.setP4(pfo.p4()*weight);
+  }//if should use charged PFO weighting scheme
+  return StatusCode::SUCCESS;
+}
\ No newline at end of file
diff --git a/Reconstruction/Jet/JetRecTools/Root/JetConstituentModSequence.cxx b/Reconstruction/Jet/JetRecTools/Root/JetConstituentModSequence.cxx
index 8368124a12e6074f67b09629e109b51b0d9b2f79..29733bfec1130915fb3202d3110f884d0ca8c67b 100644
--- a/Reconstruction/Jet/JetRecTools/Root/JetConstituentModSequence.cxx
+++ b/Reconstruction/Jet/JetRecTools/Root/JetConstituentModSequence.cxx
@@ -20,6 +20,9 @@
 #include "xAODPFlow/TrackCaloCluster.h"
 #include "xAODPFlow/TrackCaloClusterContainer.h"
 #include "xAODPFlow/TrackCaloClusterAuxContainer.h"
+#include "xAODPFlow/FlowElement.h"
+#include "xAODPFlow/FlowElementContainer.h"
+#include "xAODPFlow/FlowElementAuxContainer.h"
 #include "AthenaMonitoringKernel/Monitored.h"
 
 JetConstituentModSequence::JetConstituentModSequence(const std::string &name):
@@ -46,6 +49,7 @@ StatusCode JetConstituentModSequence::initialize() {
   ATH_CHECK( m_modifiers.retrieve() );
 
   ATH_CHECK( m_monTool.retrieve( DisableTool{m_monTool.empty()} ) );
+
   
   // Set and initialise DataHandleKeys only for the correct input type
   // Die if the input type is unsupported
@@ -83,6 +87,24 @@ StatusCode JetConstituentModSequence::initialize() {
     ATH_CHECK(m_inTCCKey.initialize());
     ATH_CHECK(m_outTCCKey.initialize());
     break;
+  case xAOD::Type::FlowElement:
+    {
+      // TODO: This assumes a PFlow-style neutral and charged collection.
+      //       More general FlowElements (e.g. CaloClusters) may necessitate a rework here later.
+      m_inChargedFEKey = m_inputContainer + "ChargedFlowElements";
+      m_inNeutralFEKey = m_inputContainer + "NeutralFlowElements";
+
+      m_outChargedFEKey = m_outputContainer+"ChargedFlowElements";
+      m_outNeutralFEKey = m_outputContainer+"NeutralFlowElements";
+      m_outAllFEKey = m_outputContainer+"FlowElements";
+
+      ATH_CHECK(m_inChargedFEKey.initialize());
+      ATH_CHECK(m_inNeutralFEKey.initialize());
+      ATH_CHECK(m_outChargedFEKey.initialize());
+      ATH_CHECK(m_outNeutralFEKey.initialize());
+      ATH_CHECK(m_outAllFEKey.initialize());
+      break;
+    }
   default:
     ATH_MSG_ERROR(" Unsupported input type "<< m_inputType );
     return StatusCode::FAILURE;
@@ -93,48 +115,57 @@ StatusCode JetConstituentModSequence::initialize() {
   
 int JetConstituentModSequence::execute() const {
 
-  // Define monitored quantities							
-  auto t_exec     = Monitored::Timer<std::chrono::milliseconds>( "TIME_execute"  );	
-  auto t_subtract = Monitored::Timer<std::chrono::milliseconds>( "TIME_subtract" );	
-  // Explicitly start/stop the timer around the subtraction tool calls			
+  // Define monitored quantities
+  auto t_exec     = Monitored::Timer<std::chrono::milliseconds>( "TIME_execute"  );
+  auto t_subtract = Monitored::Timer<std::chrono::milliseconds>( "TIME_subtract" );
+  // Explicitly start/stop the timer around the subtraction tool calls
   t_subtract.start();
 
   // Create the shallow copy according to the input type
   switch(m_inputType){
 
-  case xAOD::Type::CaloCluster: {
-    if (m_trigger){
-      auto clustersin = dynamic_cast<const xAOD::CaloClusterContainer*>(m_trigInputConstits);
-      if(clustersin==nullptr) {
-	ATH_MSG_ERROR("Failed to cast trigInputConstits to CaloCluster");
-	return(3);
+    case xAOD::Type::CaloCluster: {
+      if (m_trigger){
+        auto clustersin = dynamic_cast<const xAOD::CaloClusterContainer*>(m_trigInputConstits);
+        if(clustersin==nullptr) {
+          ATH_MSG_ERROR("Failed to cast trigInputConstits to CaloCluster");
+          return(3);
+        }
+        auto sc  = copyModForTrigger<xAOD::CaloClusterContainer,xAOD::CaloClusterAuxContainer,xAOD::CaloCluster>(*clustersin);
+        if(!sc.isSuccess()) return 1;
+      } else {
+        auto sc  = copyModRecord(m_inClusterKey, 
+                                m_outClusterKey);
+        if(!sc.isSuccess()) return 1;
       }
-      auto sc  = copyModForTrigger<xAOD::CaloClusterContainer,xAOD::CaloClusterAuxContainer,xAOD::CaloCluster>(*clustersin);
-      if(!sc.isSuccess()){return 1;}
-    } else {
-      auto sc  = copyModRecord(m_inClusterKey, 
-			       m_outClusterKey);
-      if(!sc.isSuccess()){return 1;}
+      break;
     }
-    break; }
 
-  case xAOD::Type::ParticleFlow: {
-    if (m_trigger){return 2;}
-    auto sc = copyModRecordPFO();
-    if(!sc.isSuccess()){return 1;}
-    break;
-  }
+    case xAOD::Type::ParticleFlow: {
+      if (m_trigger) return 2;
+      auto sc = copyModRecordPFO();
+      if(!sc.isSuccess()) return 1;
+      break;
+    }
 
-  case xAOD::Type::TrackCaloCluster : {
-    auto sc  = copyModRecord(m_inTCCKey, 
-			       m_outTCCKey);
-    if(!sc.isSuccess()){return 1;}
-    break;}
+    case xAOD::Type::FlowElement: {
+      if (m_trigger) return 2;
+      auto sc = copyModRecordFE();
+      if(!sc.isSuccess()) return 1;
+      break;
+    }
 
-  default: {
-    ATH_MSG_WARNING( "Unsupported input type " << m_inputType );
-    return 1;
-  }
+    case xAOD::Type::TrackCaloCluster : {
+      auto sc  = copyModRecord(m_inTCCKey, 
+                                m_outTCCKey);
+      if(!sc.isSuccess()){return 1;}
+      break;
+    }
+
+    default: {
+      ATH_MSG_WARNING( "Unsupported input type " << m_inputType );
+      return 1;
+    }
     
   }
 
@@ -163,8 +194,8 @@ JetConstituentModSequence::copyModRecordPFO() const {
   auto inChargedPFOHandle = makeHandle(m_inChargedPFOKey);
   if(!inNeutralPFOHandle.isValid()){
     ATH_MSG_WARNING("Unable to retrieve input PFO containers \""
-		    << m_inNeutralPFOKey.key() << "\" and \""
-		    << m_inChargedPFOKey.key() << "\"");
+                    << m_inNeutralPFOKey.key() << "\" and \""
+                    << m_inChargedPFOKey.key() << "\"");
     return StatusCode::FAILURE;
   }
 
@@ -176,7 +207,7 @@ JetConstituentModSequence::copyModRecordPFO() const {
   //
   auto outNeutralPFOHandle = makeHandle(m_outNeutralPFOKey);
   ATH_CHECK(outNeutralPFOHandle.record(std::unique_ptr<xAOD::PFOContainer>(neutralCopy.first),
-				       std::unique_ptr<xAOD::ShallowAuxContainer>(neutralCopy.second)));
+                                       std::unique_ptr<xAOD::ShallowAuxContainer>(neutralCopy.second)));
   //    Charged PFOs
   std::pair<xAOD::PFOContainer*, xAOD::ShallowAuxContainer* > chargedCopy =
     xAOD::shallowCopyContainer(*inChargedPFOHandle);
@@ -184,7 +215,7 @@ JetConstituentModSequence::copyModRecordPFO() const {
   //
   auto outChargedPFOHandle = makeHandle(m_outChargedPFOKey);
   ATH_CHECK(outChargedPFOHandle.record(std::unique_ptr<xAOD::PFOContainer>(chargedCopy.first),
-				       std::unique_ptr<xAOD::ShallowAuxContainer>(chargedCopy.second)));
+                                       std::unique_ptr<xAOD::ShallowAuxContainer>(chargedCopy.second)));
 
   // 2. Set up output handle for merged (view) container and record
   auto outAllPFOHandle = makeHandle(m_outAllPFOKey);
@@ -192,8 +223,8 @@ JetConstituentModSequence::copyModRecordPFO() const {
   //    Merge charged & neutral PFOs into the viw container
   outAllPFOHandle->assign(outNeutralPFOHandle->begin(), outNeutralPFOHandle->end());
   outAllPFOHandle->insert(outAllPFOHandle->end(),
-			  outChargedPFOHandle->begin(), 
-			  outChargedPFOHandle->end());
+                          outChargedPFOHandle->begin(), 
+                          outChargedPFOHandle->end());
 
   // 3. Now process modifications on all PFOs
   for (auto t : m_modifiers) {ATH_CHECK(t->process(outAllPFOHandle.ptr()));}
@@ -201,6 +232,58 @@ JetConstituentModSequence::copyModRecordPFO() const {
   return StatusCode::SUCCESS;
 }
 
+StatusCode JetConstituentModSequence::copyModRecordFE() const {
+
+  // Cannot be handled the same way as other objects (e.g. clusters),
+  // because the data is split between two containers, but we need
+  // information from both to do the modifications.
+  //
+  // The logic is:
+  //   1. Copy the charged & neutral containers independently
+  //   2. Merge into a combined view container
+  //   3. Modify the combined container
+
+  // 1. Retrieve the input containers
+  SG::ReadHandle<xAOD::FlowElementContainer> inNeutralFEHandle = makeHandle(m_inNeutralFEKey);
+  SG::ReadHandle<xAOD::FlowElementContainer> inChargedFEHandle = makeHandle(m_inChargedFEKey);
+  if(!inNeutralFEHandle.isValid()){
+    ATH_MSG_WARNING("Unable to retrieve input FlowElement containers \""
+                    << m_inNeutralFEKey.key() << "\" and \""
+                    << m_inChargedFEKey.key() << "\"");
+    return StatusCode::FAILURE;
+  }
+
+  //    Copy the input containers individually, set I/O option and record
+  //    Neutral FEs
+  std::pair<xAOD::FlowElementContainer*, xAOD::ShallowAuxContainer* > neutralCopy = xAOD::shallowCopyContainer(*inNeutralFEHandle);
+  neutralCopy.second->setShallowIO(m_saveAsShallow);
+
+  SG::WriteHandle<xAOD::FlowElementContainer> outNeutralFEHandle = makeHandle(m_outNeutralFEKey);
+  ATH_CHECK(outNeutralFEHandle.record(std::unique_ptr<xAOD::FlowElementContainer>(neutralCopy.first),
+                                      std::unique_ptr<xAOD::ShallowAuxContainer>(neutralCopy.second)));
+  //    Charged FEs
+  std::pair<xAOD::FlowElementContainer*, xAOD::ShallowAuxContainer* > chargedCopy = xAOD::shallowCopyContainer(*inChargedFEHandle);
+  chargedCopy.second->setShallowIO(m_saveAsShallow);
+
+  SG::WriteHandle<xAOD::FlowElementContainer> outChargedFEHandle = makeHandle(m_outChargedFEKey);
+  ATH_CHECK(outChargedFEHandle.record(std::unique_ptr<xAOD::FlowElementContainer>(chargedCopy.first),
+                                      std::unique_ptr<xAOD::ShallowAuxContainer>(chargedCopy.second)));
+
+  // 2. Set up output handle for merged (view) container and record
+  SG::WriteHandle<xAOD::FlowElementContainer> outAllFEHandle = makeHandle(m_outAllFEKey);
+  ATH_CHECK(outAllFEHandle.record(std::make_unique<xAOD::FlowElementContainer>(SG::VIEW_ELEMENTS)));
+  //    Merge charged & neutral FEs into the view container
+  outAllFEHandle->assign(outNeutralFEHandle->begin(), outNeutralFEHandle->end());
+  outAllFEHandle->insert(outAllFEHandle->end(),
+                          outChargedFEHandle->begin(), 
+                          outChargedFEHandle->end());
+
+  // 3. Now process modifications on all FEs
+  for (auto t : m_modifiers) {ATH_CHECK(t->process(outAllFEHandle.ptr()));}
+
+  return StatusCode::SUCCESS;
+}
+
 void JetConstituentModSequence::setInputClusterCollection(const xAOD::IParticleContainer *cont) {
   m_trigInputConstits = cont;
 }
diff --git a/Reconstruction/Jet/JetRecTools/Root/JetConstituentModifierBase.cxx b/Reconstruction/Jet/JetRecTools/Root/JetConstituentModifierBase.cxx
index b87018b865b8e94aab56830ee5a74ff20d441b5e..16afb87a695c1cd49c44c9ce6e334f90e8fa900d 100644
--- a/Reconstruction/Jet/JetRecTools/Root/JetConstituentModifierBase.cxx
+++ b/Reconstruction/Jet/JetRecTools/Root/JetConstituentModifierBase.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 // Source file for the JetConstituentModifierBase.h
@@ -8,6 +8,7 @@
 #include "JetRecTools/JetConstituentModifierBase.h"
 #include "xAODCaloEvent/CaloCluster.h"
 #include "xAODPFlow/PFO.h"
+#include "xAODPFlow/FlowElement.h"
 #include "xAODPFlow/TrackCaloCluster.h"
 
 JetConstituentModifierBase::JetConstituentModifierBase(const std::string & name): asg::AsgTool(name) {
@@ -27,7 +28,7 @@ StatusCode JetConstituentModifierBase::process(xAOD::IParticleContainer* cont) c
   // be the ones passed to this method...
   if(!cont->empty() && cont->front()->type() != m_inputType) {
     ATH_MSG_ERROR("Object type mismatch! This tool expects " << m_inputType
-		  << ", but received " << cont->front()->type());
+               << ", but received " << cont->front()->type());
     return StatusCode::FAILURE;
   }
 
@@ -58,7 +59,7 @@ StatusCode JetConstituentModifierBase::setEtaPhi(xAOD::IParticle* obj, float eta
 }
 
 StatusCode JetConstituentModifierBase::setEnergyPt(xAOD::IParticle* obj, float e, float pt,
-						   const SG::AuxElement::Accessor<float>* weightAcc) const
+                                                   const SG::AuxElement::Accessor<float>* weightAcc) const
 {
   switch(m_inputType) {
   case xAOD::Type::CaloCluster:
@@ -74,11 +75,21 @@ StatusCode JetConstituentModifierBase::setEnergyPt(xAOD::IParticle* obj, float e
     {
       xAOD::PFO* pfo = static_cast<xAOD::PFO*>(obj);
       if( (m_applyToChargedPFO && pfo->isCharged()) || 
-	  (m_applyToNeutralPFO && !pfo->isCharged()) ) {
-	if(weightAcc) (*weightAcc)(*pfo) = pfo->pt() > 0. ? pt / pfo->pt() : 0.;
-	// KTJ: Temporary fix
-	// Defeats the purpose, but we need to use this to reset the 4-vec cache
-	pfo->setP4(pt, pfo->eta(), pfo->phi());
+          (m_applyToNeutralPFO && !pfo->isCharged()) ) {
+        if(weightAcc) (*weightAcc)(*pfo) = pfo->pt() > 0. ? pt / pfo->pt() : 0.;
+        // KTJ: Temporary fix
+        // Defeats the purpose, but we need to use this to reset the 4-vec cache
+        pfo->setP4(pt, pfo->eta(), pfo->phi());
+      }
+    }
+    break;
+  case xAOD::Type::FlowElement:
+    {
+      xAOD::FlowElement* pfo = static_cast<xAOD::FlowElement*>(obj);
+      if( (m_applyToChargedPFO && pfo->isCharged()) || 
+          (m_applyToNeutralPFO && !pfo->isCharged()) ) {
+        if(weightAcc) (*weightAcc)(*pfo) = pfo->pt() > 0. ? pt / pfo->pt() : 0.;
+        pfo->setP4(pt, pfo->eta(), pfo->phi(), 0.);
       }
     }
     break;
@@ -100,7 +111,7 @@ StatusCode JetConstituentModifierBase::setEnergyPt(xAOD::IParticle* obj, float e
 }
 
 StatusCode JetConstituentModifierBase::setP4(xAOD::IParticle* obj, const xAOD::JetFourMom_t& p4,
-					     const SG::AuxElement::Accessor<float>* weightAcc) const {
+                                             const SG::AuxElement::Accessor<float>* weightAcc) const {
   switch(m_inputType) {
   case xAOD::Type::CaloCluster:
     {
@@ -117,18 +128,27 @@ StatusCode JetConstituentModifierBase::setP4(xAOD::IParticle* obj, const xAOD::J
       xAOD::PFO* pfo = static_cast<xAOD::PFO*>(obj);
       // The PFO setter defaults to m=0
       if( (m_applyToChargedPFO && pfo->isCharged()) || 
-	  (m_applyToNeutralPFO && !pfo->isCharged()) ) {
-	if(weightAcc) (*weightAcc)(*pfo) = pfo->pt() > 0. ? p4.pt() / pfo->pt() : 0.;
-	pfo->setP4(p4.pt(),p4.eta(),p4.phi(),p4.mass());
+          (m_applyToNeutralPFO && !pfo->isCharged()) ) {
+        if(weightAcc) (*weightAcc)(*pfo) = pfo->pt() > 0. ? p4.pt() / pfo->pt() : 0.;
+        pfo->setP4(p4.pt(),p4.eta(),p4.phi(),p4.mass());
       }
-      break;
     }
-
+    break;
+  case xAOD::Type::FlowElement:
+    {
+      xAOD::FlowElement* pfo = static_cast<xAOD::FlowElement*>(obj);
+      if( (m_applyToChargedPFO && pfo->isCharged()) || 
+          (m_applyToNeutralPFO && !pfo->isCharged()) ) {
+        if(weightAcc) (*weightAcc)(*pfo) = pfo->pt() > 0. ? p4.pt() / pfo->pt() : 0.;
+        pfo->setP4(p4.pt(),p4.eta(),p4.phi(),p4.mass());
+      }
+    }
+    break;
   case xAOD::Type::TrackCaloCluster:
     {
       xAOD::TrackCaloCluster* tcc = static_cast<xAOD::TrackCaloCluster*>(obj);
       if( tcc->taste() != 0) {
-	      if(weightAcc) (*weightAcc)(*tcc) = tcc->pt() > 0. ? p4.pt() / tcc->pt() : 0.;
+              if(weightAcc) (*weightAcc)(*tcc) = tcc->pt() > 0. ? p4.pt() / tcc->pt() : 0.;
         tcc->setParameters(p4.pt(), p4.eta(), p4.phi(), p4.mass(), xAOD::TrackCaloCluster::Taste(tcc->taste()), tcc->trackParticleLink(), tcc->caloClusterLinks());
       }
       break;
diff --git a/Reconstruction/Jet/JetRecTools/Root/PuppiWeightTool.cxx b/Reconstruction/Jet/JetRecTools/Root/PuppiWeightTool.cxx
index 6f3e37700426e7832ab6b371b891f4d62cab2ea9..62b1da8bc69700555ef83850e36540b47a8c376b 100644
--- a/Reconstruction/Jet/JetRecTools/Root/PuppiWeightTool.cxx
+++ b/Reconstruction/Jet/JetRecTools/Root/PuppiWeightTool.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 #include "JetRecTools/PuppiWeightTool.h"
 
@@ -41,7 +41,7 @@ StatusCode PuppiWeightTool::initialize() {
   
   ATH_CHECK(m_vertexContainer_key.initialize());
 
-  if(m_inputType==xAOD::Type::ParticleFlow) {
+  if(m_inputType==xAOD::Type::ParticleFlow || m_inputType==xAOD::Type::FlowElement) {
     if(!m_applyToNeutralPFO) {
       ATH_MSG_ERROR("Incompatible configuration: ApplyToNeutralPFO=False -- what kind of pileup do you wish to suppress?");
       return StatusCode::FAILURE;
@@ -67,6 +67,10 @@ StatusCode PuppiWeightTool::finalize() {
 //------------------------------------------------------------------------------
 
 StatusCode PuppiWeightTool::process_impl(xAOD::IParticleContainer* cont) const {
+  if(m_inputType == xAOD::Type::FlowElement){
+    xAOD::FlowElementContainer* feCont = static_cast<xAOD::FlowElementContainer*> (cont);
+    return applyPuppiWeights(feCont);
+  }
   xAOD::PFOContainer* pfoCont = static_cast<xAOD::PFOContainer*> (cont);
   return applyPuppiWeights(pfoCont);
 }
@@ -103,8 +107,8 @@ StatusCode PuppiWeightTool::applyPuppiWeights(xAOD::PFOContainer* cont) const{
     else{     
       if(isCharged){
         bool matchedToPrimaryVertex=PVMatchedAcc(*ppfo);
-	if(matchedToPrimaryVertex) chargedHSVector.push_back(pj);
-	else chargedPUVector.push_back(pj);
+        if(matchedToPrimaryVertex) chargedHSVector.push_back(pj);
+        else chargedPUVector.push_back(pj);
       }
       else neutralVector.push_back(pj);
     }
@@ -152,3 +156,77 @@ StatusCode PuppiWeightTool::applyPuppiWeights(xAOD::PFOContainer* cont) const{
 
   return StatusCode::SUCCESS;
 }
+
+StatusCode PuppiWeightTool::applyPuppiWeights(xAOD::FlowElementContainer* cont) const{
+
+  const static SG::AuxElement::Accessor<bool> PVMatchedAcc("matchedToPV");
+  const static SG::AuxElement::Accessor<double> alphaAcc("PUPPI_alpha");
+  const static SG::AuxElement::Accessor<double> weightAcc("PUPPI_weight");
+
+  std::vector<fastjet::PseudoJet> chargedHSVector;
+  std::vector<fastjet::PseudoJet> chargedPUVector;
+  std::vector<fastjet::PseudoJet> neutralVector;
+  std::vector<fastjet::PseudoJet> forwardVector;
+
+  // Fill input particle vectors for puppi
+  for ( xAOD::FlowElement* pfe : *cont ) {
+    if (!PVMatchedAcc.isAvailable(*pfe)){
+      ATH_MSG_ERROR("Not known if FlowElement is matched to primary vertex.  Run CorrectPFOTool before ChargedHadronSubtractionTool");
+      return StatusCode::FAILURE;
+    }
+
+    if (pfe->pt()<=FLT_MIN) continue;
+
+    fastjet::PseudoJet pj(pfe->p4());
+
+    if(fabs(pfe->eta()) > m_etaBoundary) forwardVector.push_back(pj);
+    else{     
+      if(pfe->isCharged()){
+        bool matchedToPrimaryVertex=PVMatchedAcc(*pfe);
+        if(matchedToPrimaryVertex) chargedHSVector.push_back(pj);
+        else chargedPUVector.push_back(pj);
+      }
+      else neutralVector.push_back(pj);
+    }
+  }
+
+  //Count the number of primary vertices
+  const xAOD::VertexContainer* pvtxs = nullptr;
+  auto handle = SG::makeHandle(m_vertexContainer_key);
+  if (!handle.isValid()){
+    ATH_MSG_WARNING(" This event has no primary vertices " );
+    return StatusCode::FAILURE;
+  }
+    
+  pvtxs = handle.cptr();
+  if(pvtxs->empty()){
+    ATH_MSG_WARNING(" This event has no primary vertices " );
+    return StatusCode::FAILURE;
+  }
+
+  int nPV=0;
+  for( auto vtx_itr : *pvtxs ){
+    if((int)vtx_itr->nTrackParticles() < 2 ) continue;
+    ++nPV;
+  }
+
+  m_puppi->setParticles(chargedHSVector, chargedPUVector, neutralVector, forwardVector, nPV);
+
+  for ( xAOD::FlowElement* pfe : *cont ) {
+    bool isForward = (fabs(pfe->eta()) > m_etaBoundary);
+
+    fastjet::PseudoJet pj(pfe->p4());
+
+    double weight = m_puppi->getWeight(pj);
+    double alpha = m_puppi->getAlpha(pj);
+
+    if ((!(pfe->isCharged()) || isForward) && m_applyWeight) pfe->setP4(weight*pfe->p4());
+    alphaAcc(*pfe) = alpha;
+    weightAcc(*pfe) = weight;
+  }
+
+  ATH_MSG_DEBUG("Median: "<<m_puppi->getMedian());
+  ATH_MSG_DEBUG("RMS: "<<m_puppi->getRMS());
+
+  return StatusCode::SUCCESS;
+}
diff --git a/Reconstruction/Jet/JetRecTools/Root/SoftKillerWeightTool.cxx b/Reconstruction/Jet/JetRecTools/Root/SoftKillerWeightTool.cxx
index 60706f80a54aec8ce91cda6634a1e1773199060c..723179bdcdb96b7dcf564ed59deaea48afe628f3 100644
--- a/Reconstruction/Jet/JetRecTools/Root/SoftKillerWeightTool.cxx
+++ b/Reconstruction/Jet/JetRecTools/Root/SoftKillerWeightTool.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 #include "JetRecTools/SoftKillerWeightTool.h"
@@ -10,6 +10,7 @@
 #include "fastjet/contrib/SoftKiller.hh"
 
 #include "xAODPFlow/PFO.h"
+#include "xAODPFlow/FlowElement.h"
 #include "xAODPFlow/TrackCaloCluster.h"
 
 
@@ -48,7 +49,7 @@ SoftKillerWeightTool::SoftKillerWeightTool(const std::string& name) : JetConstit
 
 
 StatusCode SoftKillerWeightTool::initialize() {
-  if(m_isCaloSplit && (m_inputType!=xAOD::Type::ParticleFlow && m_inputType!=xAOD::Type::CaloCluster)) {
+  if(m_isCaloSplit && (m_inputType!=xAOD::Type::ParticleFlow && m_inputType!=xAOD::Type::CaloCluster && m_inputType!=xAOD::Type::FlowElement)) {
     ATH_MSG_ERROR("SoftKillerWeightTool requires CaloCluster or PFO inputs when isCaloSplit is true."
 		  << " It cannot apply split calorimeters on objects of type "
 		  << m_inputType);
@@ -131,6 +132,14 @@ double SoftKillerWeightTool::getSoftKillerMinPt(xAOD::IParticleContainer& cont)
       xAOD::TrackCaloCluster* tcc = static_cast<xAOD::TrackCaloCluster*>(part);
       accept = (tcc->taste()!= 0);
     }
+    if(m_inputType==xAOD::Type::FlowElement){
+      xAOD::FlowElement* fe = static_cast<xAOD::FlowElement*>(part);
+      if(fe->signalType() & xAOD::FlowElement::PFlow){
+        if(m_ignoreChargedPFOs) accept = !fe->isCharged();
+      }
+      else
+        accept = !fe->isCharged();
+    }
     if(accept) {
       partPJ.push_back( fastjet::PseudoJet( part->p4() ));
     }
@@ -164,6 +173,14 @@ std::pair<double,double> SoftKillerWeightTool::getSoftKillerMinPtSplit(xAOD::IPa
       xAOD::TrackCaloCluster* tcc = static_cast<xAOD::TrackCaloCluster*>(part);
       accept = (tcc->taste()!= 0);
     }
+    if(m_inputType==xAOD::Type::FlowElement){
+      xAOD::FlowElement* fe = static_cast<xAOD::FlowElement*>(part);
+      if(fe->signalType() & xAOD::FlowElement::PFlow){
+        if(m_ignoreChargedPFOs) accept = !fe->isCharged();
+      }
+      else
+        accept = !fe->isCharged();
+    }
     if(accept) {
       double center_lambda = acc_clambda.isAvailable(*part) ? acc_clambda(*part) : 0.;
       if( center_lambda < m_lambdaCalDivide) partPJ_ECal.push_back( fastjet::PseudoJet( part->p4() ));
diff --git a/Reconstruction/Jet/JetRecTools/Root/VoronoiWeightTool.cxx b/Reconstruction/Jet/JetRecTools/Root/VoronoiWeightTool.cxx
index b8bd6b3a16661120cb18e7adb2676093e55ccae8..551f8f225062202d96801ad8fe7ebde220b0f823 100644
--- a/Reconstruction/Jet/JetRecTools/Root/VoronoiWeightTool.cxx
+++ b/Reconstruction/Jet/JetRecTools/Root/VoronoiWeightTool.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 #include "JetRecTools/VoronoiWeightTool.h"
@@ -11,6 +11,7 @@
 
 #include "xAODPFlow/PFO.h"
 #include "xAODPFlow/TrackCaloCluster.h"
+#include "xAODPFlow/FlowElement.h"
 
 namespace SortHelper {
   //
@@ -123,6 +124,14 @@ StatusCode VoronoiWeightTool::process_impl(xAOD::IParticleContainer* particlesin
       xAOD::TrackCaloCluster* tcc = static_cast<xAOD::TrackCaloCluster*>(part);
       accept &= (tcc->taste()!= 0);
     }
+    if(m_inputType==xAOD::Type::FlowElement){
+      xAOD::FlowElement* fe = static_cast<xAOD::FlowElement*>(part);
+      if(fe->signalType() & xAOD::FlowElement::PFlow){
+        if(m_ignoreChargedPFOs) accept = !fe->isCharged();
+      }
+      else
+        accept = !fe->isCharged();
+    }
     if(accept) {
       particles.push_back( fastjet::PseudoJet(part->p4()) );
       // ATH_MSG_VERBOSE( "Accepted particle with pt " << part->pt() );
diff --git a/Reconstruction/Jet/JetUtils/CMakeLists.txt b/Reconstruction/Jet/JetUtils/CMakeLists.txt
index 12ca8d0e05f5aadde1cf023ea79ddbead1d242e0..4cd94d9879ec23e2fac3c63b7548fbaf30f5262d 100644
--- a/Reconstruction/Jet/JetUtils/CMakeLists.txt
+++ b/Reconstruction/Jet/JetUtils/CMakeLists.txt
@@ -8,13 +8,13 @@ if( XAOD_STANDALONE OR XAOD_ANALYSIS )
    atlas_add_library( JetUtils
       JetUtils/*.h JetUtils/*.icc Root/*.cxx
       PUBLIC_HEADERS JetUtils
-      LINK_LIBRARIES xAODCaloEvent xAODJet xAODTracking
+      LINK_LIBRARIES xAODCaloEvent xAODJet xAODTracking PFlowUtilsLib
       PRIVATE_LINK_LIBRARIES CaloGeoHelpers xAODPFlow )
 else()
    atlas_add_library( JetUtils
       JetUtils/*.h JetUtils/*.icc Root/*.cxx src/*.cxx
       PUBLIC_HEADERS JetUtils
-      LINK_LIBRARIES CaloEvent xAODCaloEvent xAODJet xAODTracking
+      LINK_LIBRARIES CaloEvent xAODCaloEvent xAODJet xAODTracking PFlowUtilsLib
       PRIVATE_LINK_LIBRARIES CaloGeoHelpers xAODPFlow TileEvent )
 endif()
 
diff --git a/Reconstruction/Jet/JetUtils/Root/JetCaloCalculations.cxx b/Reconstruction/Jet/JetUtils/Root/JetCaloCalculations.cxx
index ba027622237fa3adf51f169b26eac3e79f285d4b..e7767d649e3d2009492eb4c03f09bc4fef16d044 100644
--- a/Reconstruction/Jet/JetUtils/Root/JetCaloCalculations.cxx
+++ b/Reconstruction/Jet/JetUtils/Root/JetCaloCalculations.cxx
@@ -8,8 +8,10 @@
 #include "xAODCaloEvent/CaloCluster.h"
 
 #include "xAODPFlow/PFO.h"
+#include "xAODPFlow/FlowElement.h"
 
 #include "JetUtils/JetCaloCalculations.h"
+#include "PFlowUtils/FEHelpers.h"
 #include "xAODJet/JetAttributes.h"
 ////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////
@@ -77,6 +79,43 @@ namespace CaloConstitHelpers {
 
   };
 
+  ///****************************************************
+  /// \class FlowElementExtractor
+  /// 
+  /// Extract cluster moments when JetConstituents are FlowElements.
+  struct FlowElementExtractor : public CaloConstitExtractor {
+    virtual ~FlowElementExtractor(){}
+    virtual bool valid(JetConstitIterator & it ) {
+      const xAOD::FlowElement* fe = dynamic_cast<const xAOD::FlowElement*>(it->rawConstituent());
+      if (fe != nullptr) return (!fe->isCharged());
+      return false;
+    }
+
+    virtual double moment(JetConstitIterator & it , xAOD::CaloCluster::MomentType momentType){
+      float m;
+      const xAOD::FlowElement* fe = static_cast<const xAOD::FlowElement*>(it->rawConstituent());
+      FEHelpers::getClusterMoment(*fe, momentType, m);
+      return m;      
+    }
+
+    virtual double time(JetConstitIterator & it){
+      const xAOD::FlowElement* fe = static_cast<const xAOD::FlowElement*>(it->rawConstituent());
+      const static SG::AuxElement::ConstAccessor<float> accTiming("TIMING");
+      return accTiming(*fe);
+    }        
+
+    virtual double energyHEC(JetConstitIterator & it ){
+      const xAOD::FlowElement* fe = static_cast<const xAOD::FlowElement*>(it->rawConstituent());
+      // Add up the four individual HEC layers
+      const static SG::AuxElement::ConstAccessor<float> accHEC0("LAYERENERGY_HEC0");
+      const static SG::AuxElement::ConstAccessor<float> accHEC1("LAYERENERGY_HEC1");
+      const static SG::AuxElement::ConstAccessor<float> accHEC2("LAYERENERGY_HEC2");
+      const static SG::AuxElement::ConstAccessor<float> accHEC3("LAYERENERGY_HEC3");
+      return accHEC0(*fe) + accHEC1(*fe) + accHEC2(*fe) + accHEC3(*fe);
+    }
+
+  };
+
 
   /// returns a pointer to a CaloConstitExtractor for a given jet. Do not delete the pointer !
   /// WARNING : not entirely safe. Assumes that all jet constituents have the same type as 1st constit.
@@ -84,6 +123,7 @@ namespace CaloConstitHelpers {
     static CaloConstitExtractor nullEx;
     static CaloClusterExtractor clusteEx;
     static PFOExtractor pfoEx;
+    static FlowElementExtractor feEx;
 
     if(jet->numConstituents() == 0 ) return &nullEx;    
 
@@ -94,6 +134,8 @@ namespace CaloConstitHelpers {
       return &clusteEx;
     case xAOD::Type::ParticleFlow:
       return &pfoEx;
+    case xAOD::Type::FlowElement:
+      return &feEx;
     default:
       break;
     }
diff --git a/Reconstruction/PFlow/PFlowUtils/PFlowUtils/FEHelpers.h b/Reconstruction/PFlow/PFlowUtils/PFlowUtils/FEHelpers.h
new file mode 100644
index 0000000000000000000000000000000000000000..e9fff1b1f06e71d30405c546b7786035f251ec21
--- /dev/null
+++ b/Reconstruction/PFlow/PFlowUtils/PFlowUtils/FEHelpers.h
@@ -0,0 +1,27 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+// Author: Bill Balunas <bill.balunas@cern.ch>
+
+#ifndef PFLOWUTILS_FEHELPERS_H
+#define PFLOWUTILS_FEHELPERS_H 1
+
+#include "xAODPFlow/FlowElement.h"
+#include "xAODTracking/Vertex.h"
+#include "xAODCaloEvent/CaloCluster.h"
+
+namespace FEHelpers {
+
+  TLorentzVector getVertexCorrectedFourVec(const xAOD::FlowElement& fe, const xAOD::Vertex& vertexToCorrectTo);
+  TLorentzVector getVertexCorrectedFourVec(const xAOD::FlowElement& fe, const TVector3& vertexToCorrectTo);
+  void vertexCorrectTheFourVector(const xAOD::FlowElement& fe, const TVector3& vertexToCorrectTo, TLorentzVector& theFourVector);
+
+  bool getClusterMoment(const xAOD::FlowElement& fe, xAOD::CaloCluster::MomentType momentType, float& value);
+  std::vector<float> getEnergiesPerSampling(const xAOD::FlowElement& fe);
+
+  std::string getClusterMomentName(xAOD::CaloCluster::MomentType momentType);
+
+}
+
+#endif
\ No newline at end of file
diff --git a/Reconstruction/PFlow/PFlowUtils/PFlowUtils/IWeightPFOTool.h b/Reconstruction/PFlow/PFlowUtils/PFlowUtils/IWeightPFOTool.h
index c8c0056aefdf86a04ffc0255820aad24a3018465..b4a6b8b34c1b8614054abfc6e9d2d37291939a86 100644
--- a/Reconstruction/PFlow/PFlowUtils/PFlowUtils/IWeightPFOTool.h
+++ b/Reconstruction/PFlow/PFlowUtils/PFlowUtils/IWeightPFOTool.h
@@ -10,6 +10,7 @@
 #include "AsgTools/IAsgTool.h"
 
 #include "xAODPFlow/PFOContainer.h"
+#include "xAODPFlow/FlowElementContainer.h"
 #include "PFlowUtils/PFODefs.h"
 
 namespace CP {
@@ -23,6 +24,7 @@ namespace CP {
 
     /** given a PFO, extract weight */
     virtual StatusCode fillWeight( const xAOD::PFO& cpfo, float& weight ) const = 0;
+    virtual StatusCode fillWeight( const xAOD::FlowElement& cpfo, float& weight ) const = 0;
 
   };
 
diff --git a/Reconstruction/PFlow/PFlowUtils/PFlowUtils/WeightPFOTool.h b/Reconstruction/PFlow/PFlowUtils/PFlowUtils/WeightPFOTool.h
index 1e30a994e224158ba07b2742f65df20ce631fbac..44f3b14c8fa90b9300a5b34ec0622645eb4d18c9 100644
--- a/Reconstruction/PFlow/PFlowUtils/PFlowUtils/WeightPFOTool.h
+++ b/Reconstruction/PFlow/PFlowUtils/PFlowUtils/WeightPFOTool.h
@@ -26,6 +26,7 @@ namespace CP {
 
     // given a PFO, extract weight
     StatusCode fillWeight( const xAOD::PFO& cpfo, float& weight) const;
+    StatusCode fillWeight( const xAOD::FlowElement& cpfo, float& weight) const;
 
   private:
 
diff --git a/Reconstruction/PFlow/PFlowUtils/Root/FEHelpers.cxx b/Reconstruction/PFlow/PFlowUtils/Root/FEHelpers.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..14c3286782b43138dc8e0d991d67e2c1dd936c3f
--- /dev/null
+++ b/Reconstruction/PFlow/PFlowUtils/Root/FEHelpers.cxx
@@ -0,0 +1,174 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+// Author: Bill Balunas <bill.balunas@cern.ch>
+
+#include "PFlowUtils/FEHelpers.h"
+
+namespace FEHelpers {
+
+  TLorentzVector getVertexCorrectedFourVec(const xAOD::FlowElement& fe, const xAOD::Vertex& vertexToCorrectTo){
+    TVector3 theVertexVector(vertexToCorrectTo.x(), vertexToCorrectTo.y(), vertexToCorrectTo.z());
+    return getVertexCorrectedFourVec(fe, theVertexVector);
+  }
+
+  TLorentzVector getVertexCorrectedFourVec(const xAOD::FlowElement& fe, const TVector3& vertexToCorrectTo){
+
+    TLorentzVector theNewVector(0.0,0.0,0.0,0.0);
+    theNewVector.SetPtEtaPhiM(fe.pt(), fe.eta(), fe.phi(), fe.m());
+
+    vertexCorrectTheFourVector(fe, vertexToCorrectTo, theNewVector);
+    return theNewVector;
+  }
+
+  void vertexCorrectTheFourVector(const xAOD::FlowElement& fe, const TVector3& vertexToCorrectTo, TLorentzVector& theFourVector){
+
+    float clusterEta = theFourVector.Eta();
+
+    const static SG::AuxElement::ConstAccessor<float> accCenterMag("CENTER_MAG");
+    float centerMag = accCenterMag(fe);
+
+    float radius = centerMag/cosh(clusterEta);
+
+    float EtaVertexCorrection = 0.0, PhiVertexCorrection = 0.0;
+    float clusterPhi = theFourVector.Phi();
+
+    if (radius > 1.0 && centerMag > 1e-3){
+      EtaVertexCorrection = (-vertexToCorrectTo.Z()/cosh(clusterEta) + (vertexToCorrectTo.X()*cos(clusterPhi) + vertexToCorrectTo.Y()*sin(clusterPhi))*tanh(clusterEta))/radius;
+      PhiVertexCorrection = (vertexToCorrectTo.X()*sin(clusterPhi) - vertexToCorrectTo.Y()*cos(clusterPhi))/radius;
+    }
+
+    float etaVertexCorrected = clusterEta + EtaVertexCorrection;
+    float p = std::sqrt(theFourVector.E()*theFourVector.E()-theFourVector.M()*theFourVector.M());
+    float ptVertexCorrected = p/cosh(etaVertexCorrected); 
+    theFourVector.SetPtEtaPhiM(ptVertexCorrected, etaVertexCorrected, clusterPhi + PhiVertexCorrection, theFourVector.M());
+
+  }
+
+  bool getClusterMoment(const xAOD::FlowElement& fe, xAOD::CaloCluster::MomentType momentType, float& value){
+
+    // First try to retrieve directly from the cluster
+    const xAOD::CaloCluster* cluster = dynamic_cast<const xAOD::CaloCluster*>(fe.otherObject(0));
+    if(cluster != nullptr){
+      double tmpValue(-99.);
+      if(cluster->retrieveMoment(momentType, tmpValue)){
+        value = static_cast<float>(tmpValue);
+        return true;
+      }
+    }
+
+    // Cluster or moment from it was unavailable, try to retrieve from the FE itself instead
+    std::string momentName = getClusterMomentName(momentType);
+    if(!momentName.empty()){
+      const SG::AuxElement::ConstAccessor<float> acc(momentName);
+      if(acc.isAvailable(fe)){
+        value = acc(fe);
+        return true;
+      }
+    }
+
+    return false;
+  }
+
+  std::vector<float> getEnergiesPerSampling(const xAOD::FlowElement& fe){
+
+    const static SG::AuxElement::ConstAccessor< std::vector<float> > accEPerSampling("ePerSampling");
+    if(accEPerSampling.isAvailable(fe)) return accEPerSampling(fe);
+
+    // This wasn't stored as a vector, try to retrieve the individual elements
+    const static SG::AuxElement::ConstAccessor<float> accPreSamplerB("LAYERENERGY_PreSamplerB");
+    const static SG::AuxElement::ConstAccessor<float> accEMB1("LAYERENERGY_EMB1");
+    const static SG::AuxElement::ConstAccessor<float> accEMB2("LAYERENERGY_EMB2");
+    const static SG::AuxElement::ConstAccessor<float> accEMB3("LAYERENERGY_EMB3");
+    const static SG::AuxElement::ConstAccessor<float> accPreSamplerE("LAYERENERGY_PreSamplerE");
+    const static SG::AuxElement::ConstAccessor<float> accEME1("LAYERENERGY_EME1");
+    const static SG::AuxElement::ConstAccessor<float> accEME2("LAYERENERGY_EME2");
+    const static SG::AuxElement::ConstAccessor<float> accEME3("LAYERENERGY_EME3");
+    const static SG::AuxElement::ConstAccessor<float> accHEC0("LAYERENERGY_HEC0");
+    const static SG::AuxElement::ConstAccessor<float> accHEC1("LAYERENERGY_HEC1");
+    const static SG::AuxElement::ConstAccessor<float> accHEC2("LAYERENERGY_HEC2");
+    const static SG::AuxElement::ConstAccessor<float> accHEC3("LAYERENERGY_HEC3");
+    const static SG::AuxElement::ConstAccessor<float> accTileBar0("LAYERENERGY_TileBar0");
+    const static SG::AuxElement::ConstAccessor<float> accTileBar1("LAYERENERGY_TileBar1");
+    const static SG::AuxElement::ConstAccessor<float> accTileBar2("LAYERENERGY_TileBar2");
+    const static SG::AuxElement::ConstAccessor<float> accTileGap1("LAYERENERGY_TileGap1");
+    const static SG::AuxElement::ConstAccessor<float> accTileGap2("LAYERENERGY_TileGap2");
+    const static SG::AuxElement::ConstAccessor<float> accTileGap3("LAYERENERGY_TileGap3");
+    const static SG::AuxElement::ConstAccessor<float> accTileExt0("LAYERENERGY_TileExt0");
+    const static SG::AuxElement::ConstAccessor<float> accTileExt1("LAYERENERGY_TileExt1");
+    const static SG::AuxElement::ConstAccessor<float> accTileExt2("LAYERENERGY_TileExt2");
+    const static SG::AuxElement::ConstAccessor<float> accFCAL0("LAYERENERGY_FCAL0");
+    const static SG::AuxElement::ConstAccessor<float> accFCAL1("LAYERENERGY_FCAL1");
+    const static SG::AuxElement::ConstAccessor<float> accFCAL2("LAYERENERGY_FCAL2");
+    const static SG::AuxElement::ConstAccessor<float> accMINIFCAL0("LAYERENERGY_MINIFCAL0");
+    const static SG::AuxElement::ConstAccessor<float> accMINIFCAL1("LAYERENERGY_MINIFCAL1");
+    const static SG::AuxElement::ConstAccessor<float> accMINIFCAL2("LAYERENERGY_MINIFCAL2");
+    const static SG::AuxElement::ConstAccessor<float> accMINIFCAL3("LAYERENERGY_MINIFCAL3");
+
+    // Build the vector up from the individual samplings
+    std::vector<float> result(28, 0); 
+    result[0]  = accPreSamplerB(fe);
+    result[1]  = accEMB1(fe);
+    result[2]  = accEMB2(fe);
+    result[3]  = accEMB3(fe);
+    result[4]  = accPreSamplerE(fe);
+    result[5]  = accEME1(fe);
+    result[6]  = accEME2(fe);
+    result[7]  = accEME3(fe);
+    result[8]  = accHEC0(fe);
+    result[9]  = accHEC1(fe);
+    result[10] = accHEC2(fe);
+    result[11] = accHEC3(fe);
+    result[12] = accTileBar0(fe);
+    result[13] = accTileBar1(fe);
+    result[14] = accTileBar2(fe);
+    result[15] = accTileGap1(fe);
+    result[16] = accTileGap2(fe);
+    result[17] = accTileGap3(fe);
+    result[18] = accTileExt0(fe);
+    result[19] = accTileExt1(fe);
+    result[20] = accTileExt2(fe);
+    result[21] = accFCAL0(fe);
+    result[22] = accFCAL1(fe);
+    result[23] = accFCAL2(fe);
+    result[24] = accMINIFCAL0(fe);
+    result[25] = accMINIFCAL1(fe);
+    result[26] = accMINIFCAL2(fe);
+    result[27] = accMINIFCAL3(fe);
+
+    return result;
+  }
+
+  std::string getClusterMomentName(xAOD::CaloCluster::MomentType momentType){
+
+    switch(momentType){
+      case xAOD::CaloCluster::ENG_FRAC_CORE:       return "ENG_FRAC_CORE";
+      case xAOD::CaloCluster::FIRST_ENG_DENS:      return "FIRST_ENG_DENS";
+      case xAOD::CaloCluster::CENTER_LAMBDA:       return "CENTER_LAMBDA";
+      case xAOD::CaloCluster::SECOND_R:            return "SECOND_R";
+      case xAOD::CaloCluster::DELTA_ALPHA:         return "DELTA_ALPHA";
+      case xAOD::CaloCluster::LATERAL:             return "LATERAL";
+      case xAOD::CaloCluster::LONGITUDINAL:        return "LONGITUDINAL";
+      case xAOD::CaloCluster::SECOND_LAMBDA:       return "SECOND_LAMBDA";
+      case xAOD::CaloCluster::ISOLATION:           return "ISOLATION";
+      case xAOD::CaloCluster::ENG_FRAC_MAX:        return "ENG_FRAC_MAX";
+      case xAOD::CaloCluster::ENG_BAD_CELLS:       return "ENG_BAD_CELLS";
+      case xAOD::CaloCluster::N_BAD_CELLS:         return "N_BAD_CELLS";
+      case xAOD::CaloCluster::BADLARQ_FRAC:        return "BADLARQ_FRAC";
+      case xAOD::CaloCluster::ENG_POS:             return "ENG_POS";
+      case xAOD::CaloCluster::SIGNIFICANCE:        return "SIGNIFICANCE";
+      case xAOD::CaloCluster::CELL_SIGNIFICANCE:   return "CELL_SIGNIFICANCE";
+      case xAOD::CaloCluster::CELL_SIG_SAMPLING:   return "CELL_SIG_SAMPLING";
+      case xAOD::CaloCluster::AVG_LAR_Q:           return "AVG_LAR_Q";
+      case xAOD::CaloCluster::AVG_TILE_Q:          return "AVG_TILE_Q";
+      case xAOD::CaloCluster::EM_PROBABILITY:      return "EM_PROBABILITY";
+      case xAOD::CaloCluster::ENG_CALIB_TOT:       return "ENG_CALIB_TOT";
+      case xAOD::CaloCluster::ENG_CALIB_FRAC_EM:   return "ENG_CALIB_FRAC_EM";
+      case xAOD::CaloCluster::ENG_CALIB_FRAC_HAD:  return "ENG_CALIB_FRAC_HAD";
+      case xAOD::CaloCluster::ENG_CALIB_FRAC_REST: return "ENG_CALIB_FRAC_REST";
+      default: return "";
+    }
+  }
+
+}
\ No newline at end of file
diff --git a/Reconstruction/PFlow/PFlowUtils/Root/WeightPFOTool.cxx b/Reconstruction/PFlow/PFlowUtils/Root/WeightPFOTool.cxx
index 4621d5653df81d4c14a74073ec96cbf3d9cc449a..de19ee5cccb0e12d8f925d3bad4e8fb0cf5dad98 100644
--- a/Reconstruction/PFlow/PFlowUtils/Root/WeightPFOTool.cxx
+++ b/Reconstruction/PFlow/PFlowUtils/Root/WeightPFOTool.cxx
@@ -3,6 +3,7 @@
 */
 
 #include "PFlowUtils/WeightPFOTool.h"
+#include "PFlowUtils/FEHelpers.h"
 
 namespace CP {
 
@@ -48,7 +49,7 @@ namespace CP {
       ATH_MSG_WARNING("PFO with invalid pt " << cpfo.pt() << ", quitting.");
       return StatusCode::FAILURE;
     }
-	
+        
     int isInDenseEnvironment = false;
     float expectedEnergy = 0.0;
     bool gotVariable = cpfo.attribute(xAOD::PFODetails::PFOAttributes::eflowRec_isInDenseEnvironment,isInDenseEnvironment);
@@ -59,35 +60,99 @@ namespace CP {
     } else {
       //EM case first
       if (CP::EM == theNeutralPFOScale){
-	// Start by computing the correction as though we subtracted the calo energy
-	// This interpolates between the full track P and the expected calo E
-	float EoverP = expectedEnergy/cpfo.e(); // divide once only
-	if(m_doEoverPweight) {
-	  if(cpfo.pt()<30e3) {        // take full track
-	    weight = 1.;
-	  } else if(cpfo.pt()<60e3) { // linearly interpolate between 1 and E/P
-	    float interpolf = (1.0 - (cpfo.pt()-30000)/30000);
-	    weight = EoverP + interpolf * (1-EoverP);
-	  } else {                    // take the expected energy
-	    weight = EoverP;
-	  }
-	}
-
-	ATH_MSG_VERBOSE("cpfo in dense environment? " << isInDenseEnvironment);
-	ATH_MSG_VERBOSE("cpfo pt: " << cpfo.pt() << ", E/P: " << EoverP << ", weight: " << weight);
-
-	if(isInDenseEnvironment) {
-	  // In this case we further remove the expected deposited energy from the track
-	  weight -= EoverP;
-	}
+        // Start by computing the correction as though we subtracted the calo energy
+        // This interpolates between the full track P and the expected calo E
+        float EoverP = expectedEnergy/cpfo.e(); // divide once only
+        if(m_doEoverPweight) {
+          if(cpfo.pt()<30e3) {        // take full track
+            weight = 1.;
+          } else if(cpfo.pt()<60e3) { // linearly interpolate between 1 and E/P
+            float interpolf = (1.0 - (cpfo.pt()-30000)/30000);
+            weight = EoverP + interpolf * (1-EoverP);
+          } else {                    // take the expected energy
+            weight = EoverP;
+          }
+        }
+
+        ATH_MSG_VERBOSE("cpfo in dense environment? " << isInDenseEnvironment);
+        ATH_MSG_VERBOSE("cpfo pt: " << cpfo.pt() << ", E/P: " << EoverP << ", weight: " << weight);
+
+        if(isInDenseEnvironment) {
+          // In this case we further remove the expected deposited energy from the track
+          weight -= EoverP;
+        }
       }//EM Scale
       else if (CP::LC == theNeutralPFOScale){
-	if(!isInDenseEnvironment){
-	  if(cpfo.pt()<30e3 || cpfo.pt() >= 60e03) weight = 1;
-	}
+        if(!isInDenseEnvironment){
+          if(cpfo.pt()<30e3 || cpfo.pt() >= 60e03) weight = 1;
+        }
+      }
+    }
+        
+    ATH_MSG_VERBOSE("Weight before zero check: " << weight);
+    // If the weight went to 0, set it to the ghost scale, so that the cPFOs
+    // are always added to the track.
+    if (weight<1e-9) {weight = 1e-20;}
+    ATH_MSG_VERBOSE("Final weight: " << weight);
+    return StatusCode::SUCCESS;
+  }
+
+  StatusCode WeightPFOTool::fillWeight( const xAOD::FlowElement& cpfo, float& weight ) const {
+
+    if(!(cpfo.signalType() & xAOD::FlowElement::PFlow)){
+      ATH_MSG_WARNING("FlowElement was not a PFO. Signal type was: " << cpfo.signalType());
+      return StatusCode::FAILURE;
+    }
+
+    //we need to convert the string scale back to the enum
+    PFO_JetMETConfig_inputScale theNeutralPFOScale = CP::EM;
+    CP::inputScaleMapper inputScaleMapper;
+    bool answer = inputScaleMapper.getValue(m_theNeutralPFOScaleString,theNeutralPFOScale);
+    if (false == answer) ATH_MSG_FATAL("Invalid neutral PFO Scale has been specified in PFlowUtils::PFOWeightTool");
+
+    // Compute the weights internally
+    weight = 0.;
+    if(cpfo.pt()>100e3) {
+      ATH_MSG_WARNING("PFO with invalid pt " << cpfo.pt() << ", quitting.");
+      return StatusCode::FAILURE;
+    }
+
+    const static SG::AuxElement::ConstAccessor<int> accDenseEnv("IsInDenseEnvironment");
+    const static SG::AuxElement::ConstAccessor<float> accExpE("TracksExpectedEnergyDeposit");
+
+    int isInDenseEnvironment = accDenseEnv(cpfo);
+    float expectedEnergy = accExpE(cpfo);
+
+    //EM case first
+    if (CP::EM == theNeutralPFOScale){
+      // Start by computing the correction as though we subtracted the calo energy
+      // This interpolates between the full track P and the expected calo E
+      float EoverP = expectedEnergy/cpfo.e(); // divide once only
+      if(m_doEoverPweight) {
+        if(cpfo.pt()<30e3) {        // take full track
+          weight = 1.;
+        } else if(cpfo.pt()<60e3) { // linearly interpolate between 1 and E/P
+          float interpolf = (1.0 - (cpfo.pt()-30000)/30000);
+          weight = EoverP + interpolf * (1-EoverP);
+        } else {                    // take the expected energy
+          weight = EoverP;
+        }
+      }
+
+      ATH_MSG_VERBOSE("cpfo in dense environment? " << isInDenseEnvironment);
+      ATH_MSG_VERBOSE("cpfo pt: " << cpfo.pt() << ", E/P: " << EoverP << ", weight: " << weight);
+
+      if(isInDenseEnvironment) {
+        // In this case we further remove the expected deposited energy from the track
+        weight -= EoverP;
+      }
+    }//EM Scale
+    else if (CP::LC == theNeutralPFOScale){
+      if(!isInDenseEnvironment){
+        if(cpfo.pt()<30e3 || cpfo.pt() >= 60e03) weight = 1;
       }
     }
-	
+        
     ATH_MSG_VERBOSE("Weight before zero check: " << weight);
     // If the weight went to 0, set it to the ghost scale, so that the cPFOs
     // are always added to the track.
diff --git a/Reconstruction/eflowRec/src/PFChargedFlowElementCreatorAlgorithm.cxx b/Reconstruction/eflowRec/src/PFChargedFlowElementCreatorAlgorithm.cxx
index 3ec0847b065e2e58864e20aae14da605c10006a4..899a06a9256e868abf1c518262f7610568a66432 100644
--- a/Reconstruction/eflowRec/src/PFChargedFlowElementCreatorAlgorithm.cxx
+++ b/Reconstruction/eflowRec/src/PFChargedFlowElementCreatorAlgorithm.cxx
@@ -63,6 +63,7 @@ void PFChargedFlowElementCreatorAlgorithm::createChargedFlowElements(const eflow
 
     //Now set the charge
     thisFE->setCharge(efRecTrack->getTrack()->charge());
+    thisFE->setSignalType(xAOD::FlowElement::ChargedPFlow);
 
     std::pair<double,double> etaPhi(0.0,0.0);
 
diff --git a/Reconstruction/eflowRec/src/PFNeutralFlowElementCreatorAlgorithm.cxx b/Reconstruction/eflowRec/src/PFNeutralFlowElementCreatorAlgorithm.cxx
index 731f73c2ff3f0d87d2796c5b987d62d3c4db5fee..c07823b560d487a472fa5c7ad641b889a23deb21 100644
--- a/Reconstruction/eflowRec/src/PFNeutralFlowElementCreatorAlgorithm.cxx
+++ b/Reconstruction/eflowRec/src/PFNeutralFlowElementCreatorAlgorithm.cxx
@@ -89,69 +89,70 @@ StatusCode PFNeutralFlowElementCreatorAlgorithm::createNeutralFlowElement(const
     ATH_MSG_DEBUG("Created neutral FlowElement with E, pt, eta and phi of " << thisFE->e() << ", " << thisFE->pt() << ", " << thisFE->eta() << " and " << thisFE->phi());
     
     thisFE->setCharge(0);
-
-    this->addMoment(xAOD::CaloCluster::CENTER_MAG,"PF_CENTER_MAG",*cluster,*thisFE);
-    this->addMoment(xAOD::CaloCluster::SECOND_R,"PF_SECOND_R",*cluster,*thisFE);
-    this->addMoment(xAOD::CaloCluster::CENTER_LAMBDA,"PF_CENTER_LAMBDA",*cluster,*thisFE);
-    this->addMoment(xAOD::CaloCluster::ENG_BAD_CELLS,"PF_ENG_BAD_CELLS",*cluster,*thisFE);
-    this->addMoment(xAOD::CaloCluster::N_BAD_CELLS,"PF_N_BAD_CELLS",*cluster,*thisFE);
-    this->addMoment(xAOD::CaloCluster::BADLARQ_FRAC,"PF_BADLARQ_FRAC",*cluster,*thisFE);
-    this->addMoment(xAOD::CaloCluster::ENG_POS,"PF_ENG_POS",*cluster,*thisFE);
-    this->addMoment(xAOD::CaloCluster::AVG_LAR_Q,"PF_AVG_LAR_Q",*cluster,*thisFE);
-    this->addMoment(xAOD::CaloCluster::AVG_TILE_Q,"PF_AVG_TILE_Q",*cluster,*thisFE);
-    this->addMoment(xAOD::CaloCluster::ISOLATION,"PF_ISOLATION",*cluster,*thisFE);
-    this->addMoment(xAOD::CaloCluster::SECOND_LAMBDA,"PF_SECOND_LAMBDA",*cluster,*thisFE);
-    this->addMoment(xAOD::CaloCluster::EM_PROBABILITY,"PF_EM_PROBABILITY",*cluster,*thisFE);
+    thisFE->setSignalType(xAOD::FlowElement::NeutralPFlow);
+
+    this->addMoment(xAOD::CaloCluster::CENTER_MAG,"CENTER_MAG",*cluster,*thisFE);
+    this->addMoment(xAOD::CaloCluster::SECOND_R,"SECOND_R",*cluster,*thisFE);
+    this->addMoment(xAOD::CaloCluster::CENTER_LAMBDA,"CENTER_LAMBDA",*cluster,*thisFE);
+    this->addMoment(xAOD::CaloCluster::ENG_BAD_CELLS,"ENG_BAD_CELLS",*cluster,*thisFE);
+    this->addMoment(xAOD::CaloCluster::N_BAD_CELLS,"N_BAD_CELLS",*cluster,*thisFE);
+    this->addMoment(xAOD::CaloCluster::BADLARQ_FRAC,"BADLARQ_FRAC",*cluster,*thisFE);
+    this->addMoment(xAOD::CaloCluster::ENG_POS,"ENG_POS",*cluster,*thisFE);
+    this->addMoment(xAOD::CaloCluster::AVG_LAR_Q,"AVG_LAR_Q",*cluster,*thisFE);
+    this->addMoment(xAOD::CaloCluster::AVG_TILE_Q,"AVG_TILE_Q",*cluster,*thisFE);
+    this->addMoment(xAOD::CaloCluster::ISOLATION,"ISOLATION",*cluster,*thisFE);
+    this->addMoment(xAOD::CaloCluster::SECOND_LAMBDA,"SECOND_LAMBDA",*cluster,*thisFE);
+    this->addMoment(xAOD::CaloCluster::EM_PROBABILITY,"EM_PROBABILITY",*cluster,*thisFE);
 
     if (m_useCalibHitTruth){
-      this->addMoment(xAOD::CaloCluster::ENG_CALIB_TOT,"PF_ENG_CALIB_TOT",*cluster,*thisFE);
-      this->addMoment(xAOD::CaloCluster::ENG_CALIB_FRAC_EM,"PF_ENG_CALIB_FRAC_EM",*cluster,*thisFE);
-      this->addMoment(xAOD::CaloCluster::ENG_CALIB_FRAC_HAD,"PF_ENG_CALIB_FRAC_HAD",*cluster,*thisFE);
-      this->addMoment(xAOD::CaloCluster::ENG_CALIB_FRAC_REST,"PF_ENG_CALIB_FRAC_REST",*cluster,*thisFE);
+      this->addMoment(xAOD::CaloCluster::ENG_CALIB_TOT,"ENG_CALIB_TOT",*cluster,*thisFE);
+      this->addMoment(xAOD::CaloCluster::ENG_CALIB_FRAC_EM,"ENG_CALIB_FRAC_EM",*cluster,*thisFE);
+      this->addMoment(xAOD::CaloCluster::ENG_CALIB_FRAC_HAD,"ENG_CALIB_FRAC_HAD",*cluster,*thisFE);
+      this->addMoment(xAOD::CaloCluster::ENG_CALIB_FRAC_REST,"ENG_CALIB_FRAC_REST",*cluster,*thisFE);
     }
 
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::PreSamplerB,"PF_LAYERENERGY_PreSamplerB",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::EMB1,"PF_LAYERENERGY_EMB1",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::EMB2,"PF_LAYERENERGY_EMB2",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::EMB3,"PF_LAYERENERGY_EMB3",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::PreSamplerB,"LAYERENERGY_PreSamplerB",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::EMB1,"LAYERENERGY_EMB1",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::EMB2,"LAYERENERGY_EMB2",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::EMB3,"LAYERENERGY_EMB3",*cluster,*thisFE);
     
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::PreSamplerE,"PF_LAYERENERGY_PreSamplerE",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::EME1,"PF_LAYERENERGY_EME1",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::EME2,"PF_LAYERENERGY_EME2",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::EME3,"PF_LAYERENERGY_EME3",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::PreSamplerE,"LAYERENERGY_PreSamplerE",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::EME1,"LAYERENERGY_EME1",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::EME2,"LAYERENERGY_EME2",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::EME3,"LAYERENERGY_EME3",*cluster,*thisFE);
 
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::HEC0,"PF_LAYERENERGY_HEC0",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::HEC1,"PF_LAYERENERGY_HEC1",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::HEC2,"PF_LAYERENERGY_HEC2",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::HEC3,"PF_LAYERENERGY_HEC3",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::HEC0,"LAYERENERGY_HEC0",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::HEC1,"LAYERENERGY_HEC1",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::HEC2,"LAYERENERGY_HEC2",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::HEC3,"LAYERENERGY_HEC3",*cluster,*thisFE);
 
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileBar0,"PF_LAYERENERGY_TileBar0",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileBar1,"PF_LAYERENERGY_TileBar1",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileBar2,"PF_LAYERENERGY_TileBar2",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileBar0,"LAYERENERGY_TileBar0",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileBar1,"LAYERENERGY_TileBar1",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileBar2,"LAYERENERGY_TileBar2",*cluster,*thisFE);
 
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileGap1,"PF_LAYERENERGY_TileGap1",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileGap2,"PF_LAYERENERGY_TileGap2",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileGap3,"PF_LAYERENERGY_TileGap3",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileGap1,"LAYERENERGY_TileGap1",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileGap2,"LAYERENERGY_TileGap2",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileGap3,"LAYERENERGY_TileGap3",*cluster,*thisFE);
 
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileExt0,"PF_LAYERENERGY_TileExt0",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileExt1,"PF_LAYERENERGY_TileExt1",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileExt2,"PF_LAYERENERGY_TileExt2",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileExt0,"LAYERENERGY_TileExt0",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileExt1,"LAYERENERGY_TileExt1",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::TileExt2,"LAYERENERGY_TileExt2",*cluster,*thisFE);
 
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::FCAL0,"PF_LAYERENERGY_FCAL0",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::FCAL1,"PF_LAYERENERGY_FCAL1",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::FCAL2,"PF_LAYERENERGY_FCAL2",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::FCAL0,"LAYERENERGY_FCAL0",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::FCAL1,"LAYERENERGY_FCAL1",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::FCAL2,"LAYERENERGY_FCAL2",*cluster,*thisFE);
 
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::MINIFCAL0,"PF_LAYERENERGY_MINIFCAL0",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::MINIFCAL1,"PF_LAYERENERGY_MINIFCAL1",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::MINIFCAL2,"PF_LAYERENERGY_MINIFCAL2",*cluster,*thisFE);
-    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::MINIFCAL3,"PF_LAYERENERGY_MINIFCAL3",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::MINIFCAL0,"LAYERENERGY_MINIFCAL0",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::MINIFCAL1,"LAYERENERGY_MINIFCAL1",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::MINIFCAL2,"LAYERENERGY_MINIFCAL2",*cluster,*thisFE);
+    this->addSamplingEnergy(xAOD::CaloCluster::CaloSample::MINIFCAL3,"LAYERENERGY_MINIFCAL3",*cluster,*thisFE);
 
     float layerEnergy_TileBar0 = cluster->eSample(xAOD::CaloCluster::CaloSample::TileBar0);
     float layerEnergy_TileExt0 = cluster->eSample(xAOD::CaloCluster::CaloSample::TileExt0);
-    const static SG::AuxElement::Accessor<float> accFloatTIle0E("PF_LAYERENERGY_TILE0");
+    const static SG::AuxElement::Accessor<float> accFloatTIle0E("LAYERENERGY_TILE0");
     accFloatTIle0E(*thisFE) = layerEnergy_TileBar0 + layerEnergy_TileExt0;
 
-    const static SG::AuxElement::Accessor<float> accFloatTiming("PF_TIMING");
+    const static SG::AuxElement::Accessor<float> accFloatTiming("TIMING");
     accFloatTiming(*thisFE) = cluster->time();
  
   }//cluster loop
@@ -164,7 +165,7 @@ void PFNeutralFlowElementCreatorAlgorithm::addMoment(const xAOD::CaloCluster::Mo
   bool isRetrieved = theCluster.retrieveMoment(momentType, moment);
   if (true == isRetrieved) {    
     float float_moment = static_cast<float>(moment);
-    const static SG::AuxElement::Accessor<float> accFloat(feAttribute);
+    const SG::AuxElement::Accessor<float> accFloat(feAttribute);
     accFloat(theFE) = float_moment;
   }
   else ATH_MSG_WARNING(" Could not retrieve moment from the CaloCluster");
@@ -172,6 +173,6 @@ void PFNeutralFlowElementCreatorAlgorithm::addMoment(const xAOD::CaloCluster::Mo
 }
 
 void PFNeutralFlowElementCreatorAlgorithm::addSamplingEnergy(const xAOD::CaloCluster::CaloSample& sampling, const std::string& feAttribute, const xAOD::CaloCluster& theCluster, xAOD::FlowElement& theFE) const {
-  const static SG::AuxElement::Accessor<float> accFloat(feAttribute);
+  const SG::AuxElement::Accessor<float> accFloat(feAttribute);
   accFloat(theFE) = theCluster.eSample(sampling);
 }