diff --git a/PhysicsAnalysis/JetTagging/JetHitAssociation/CMakeLists.txt b/PhysicsAnalysis/JetTagging/JetHitAssociation/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..15d0938c1e8ed604987cff48adb021d605636be9
--- /dev/null
+++ b/PhysicsAnalysis/JetTagging/JetHitAssociation/CMakeLists.txt
@@ -0,0 +1,22 @@
+# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+
+# Declare the package name:
+atlas_subdir( JetHitAssociation )
+
+# External dependencies:
+find_package( ROOT COMPONENTS Core TVector3 )
+
+atlas_add_library( JetHitAssociationLib
+ src/*.h
+ INTERFACE
+ PUBLIC_HEADERS JetHitAssociation
+ LINK_LIBRARIES GaudiKernel AthenaBaseComps TrkValInterfaces xAODEventInfo AthContainers xAODJet xAODTracking InDetEventTPCnv )
+
+# Component(s) in the package:
+atlas_add_component( JetHitAssociation
+ src/*.cxx
+ src/components/*.cxx
+ INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
+ LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps EventInfo EventPrimitives GaudiKernel JetHitAssociationLib )
+
+atlas_install_joboptions( share/*.py )
diff --git a/PhysicsAnalysis/JetTagging/JetHitAssociation/JetHitAssociation/JetHitAssociation.h b/PhysicsAnalysis/JetTagging/JetHitAssociation/JetHitAssociation/JetHitAssociation.h
new file mode 100644
index 0000000000000000000000000000000000000000..65549157081224694d7f29295d7156ded870c606
--- /dev/null
+++ b/PhysicsAnalysis/JetTagging/JetHitAssociation/JetHitAssociation/JetHitAssociation.h
@@ -0,0 +1,64 @@
+/*
+ Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef JETHITASSOCIATION_H
+#define JETHITASSOCIATION_H
+
+// STL includes
+#include <string>
+
+// FrameWork includes
+#include "AthenaBaseComps/AthAlgorithm.h"
+
+// Containers
+#include "xAODJet/JetContainer.h"
+#include "xAODTracking/VertexContainer.h"
+#include "xAODTracking/TrackMeasurementValidationContainer.h"
+#include "xAODTracking/TrackMeasurementValidationAuxContainer.h"
+
+// Read and write handles
+#include "StoreGate/ReadHandleKey.h"
+#include "StoreGate/WriteHandleKey.h"
+
+class JetHitAssociation : public AthAlgorithm
+{
+ public:
+ /// Constructor with parameters:
+ JetHitAssociation (const std::string& name,
+ ISvcLocator* svcloc);
+
+
+ virtual StatusCode initialize() override;
+ virtual StatusCode finalize() override;
+ virtual StatusCode execute() override;
+
+
+ private:
+ // functions
+ StatusCode saveHits(const std::vector<const xAOD::Jet*> &jets,
+ const xAOD::Vertex* const vertex,
+ SG::ReadHandle<xAOD::TrackMeasurementValidationContainer> &hits,
+ const SG::WriteHandle<xAOD::TrackMeasurementValidationContainer> &writeHandle,
+ unsigned long long int &nStoredHits,
+ unsigned long long int &nTotalHits);
+
+ // Read and write handles
+ SG::ReadHandleKey<xAOD::JetContainer> m_jetCollectionName;
+ SG::ReadHandleKey<xAOD::VertexContainer> m_vertexCollectionName;
+ SG::ReadHandleKey<xAOD::TrackMeasurementValidationContainer> m_inputPixHitCollectionName;
+ SG::ReadHandleKey<xAOD::TrackMeasurementValidationContainer> m_inputSCTHitCollectionName;
+ SG::WriteHandleKey<xAOD::TrackMeasurementValidationContainer> m_outputPixHitCollectionName;
+ SG::WriteHandleKey<xAOD::TrackMeasurementValidationContainer> m_outputSCTHitCollectionName;
+
+ // algorithm parameters, counters
+ float m_jetPtThreshold;
+ float m_dRmatchHitToJet;
+ unsigned long long int m_nStoredPixHits;
+ unsigned long long int m_nStoredSCTHits;
+ unsigned long long int m_nTotalPixHits;
+ unsigned long long int m_nTotalSCTHits;
+
+};
+
+#endif
diff --git a/PhysicsAnalysis/JetTagging/JetHitAssociation/share/jetHitAssociation_config.py b/PhysicsAnalysis/JetTagging/JetHitAssociation/share/jetHitAssociation_config.py
new file mode 100644
index 0000000000000000000000000000000000000000..b73935b8e49578ce4233489d06c7b6862ccd9d44
--- /dev/null
+++ b/PhysicsAnalysis/JetTagging/JetHitAssociation/share/jetHitAssociation_config.py
@@ -0,0 +1,29 @@
+# For using GeV units
+import AthenaCommon.SystemOfUnits as Units
+
+# Needed in order that PixelPrepDataToxAOD does not generate a potential memory leak error.
+import InDetRecExample.TrackingCommon as TrackingCommon
+
+# We need to add this algorithm to get the TrackMeasurementValidationContainer
+from InDetPrepRawDataToxAOD.InDetPrepRawDataToxAODConf import PixelPrepDataToxAOD, SCT_PrepDataToxAOD
+#
+# Getting the Pixel cluster container is more complicated than geting the SCT cluster container
+# You MUST minimally, have these items - a 'name' and this ClusterSplitProbabilityName or you get a memory leak ERROR
+xAOD_PixelPrepDataToxAOD = PixelPrepDataToxAOD( name = "xAOD_PixelPrepDataToxAOD",
+ ClusterSplitProbabilityName = TrackingCommon.pixelClusterSplitProbName() )
+
+from AthenaCommon.AlgSequence import AlgSequence
+algSeq = AlgSequence()
+
+algSeq += xAOD_PixelPrepDataToxAOD
+algSeq += SCT_PrepDataToxAOD('SCT_PrepDataToxAOD')
+
+from JetHitAssociation.JetHitAssociationConf import JetHitAssociation
+jetHitAssociation = JetHitAssociation('JetHitAssociation')
+jetHitAssociation.jetCollectionName = "AntiKt4EMPFlowJets"
+jetHitAssociation.jetPtThreshold = 300 * Units.GeV
+jetHitAssociation.dRmatchHitToJet = 0.4
+#jetHitAssociation.OutputLevel = DEBUG
+
+# Add algorithm to sequence
+algSeq += jetHitAssociation
\ No newline at end of file
diff --git a/PhysicsAnalysis/JetTagging/JetHitAssociation/src/JetHitAssociation.cxx b/PhysicsAnalysis/JetTagging/JetHitAssociation/src/JetHitAssociation.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..ed6a9a65eaf4289e80361aa2b05f661aafad0894
--- /dev/null
+++ b/PhysicsAnalysis/JetTagging/JetHitAssociation/src/JetHitAssociation.cxx
@@ -0,0 +1,198 @@
+/*
+Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "JetHitAssociation/JetHitAssociation.h"
+
+//Include some helpful ROOT objects here.
+#include <TVector3.h>
+
+// Constructor
+JetHitAssociation::JetHitAssociation(const std::string& name, ISvcLocator* svcloc)
+ : AthAlgorithm(name, svcloc),
+ m_nStoredPixHits(0),
+ m_nStoredSCTHits(0),
+ m_nTotalPixHits(0),
+ m_nTotalSCTHits(0) {
+
+ declareProperty("jetCollectionName", m_jetCollectionName = "AntiKt4EMPFlowJets",
+ "Jet collection that will be used for matching the Hits");
+ declareProperty("vertexCollectionName", m_vertexCollectionName = "PrimaryVertices",
+ "Name of PV collection. Default is PrimaryVertices");
+ declareProperty("outputPixHitCollectionName", m_outputPixHitCollectionName = "JetAssociatedPixelClusters",
+ "Name of output pixel hit collection. JetAssociatedPixelClusters is default.");
+ declareProperty("outputSCTHitCollectionName", m_outputSCTHitCollectionName = "JetAssociatedSCTClusters",
+ "Name of output SCT hit collection. JetAssociatedSCTClusters is default.");
+ declareProperty("inputPixHitCollectionName", m_inputPixHitCollectionName = "PixelClusters",
+ "Name of input pixel hit collection. PixelClusters is default");
+ declareProperty("inputSCTHitCollectionName", m_inputSCTHitCollectionName = "SCT_Clusters",
+ "Name of input SCT hit collection. SCT_Clusters is default");
+ declareProperty("jetPtThreshold", m_jetPtThreshold = 300000.0,
+ "Hits are saved only if they match to jets with pT > jetPtThreshold [in MeV]");
+ declareProperty("dRmatchHitToJet", m_dRmatchHitToJet = 0.4,
+ "The radius used for matching hits to jets");
+
+ declare(m_outputPixHitCollectionName);
+ declare(m_outputSCTHitCollectionName);
+}
+
+StatusCode JetHitAssociation::initialize() {
+
+ ATH_CHECK( AthAlgorithm::initialize() );
+
+ ATH_MSG_DEBUG("jetCollectionName = " << m_jetCollectionName);
+ ATH_MSG_DEBUG("vertexCollectionName = " << m_vertexCollectionName);
+ ATH_MSG_DEBUG("inputPixHitCollectionName = " << m_inputPixHitCollectionName);
+ ATH_MSG_DEBUG("inputSCTHitCollectionName = " << m_inputSCTHitCollectionName);
+ ATH_MSG_DEBUG("outputPixHitCollectionName = " << m_outputPixHitCollectionName);
+ ATH_MSG_DEBUG("outputSCTHitCollectionName = " << m_outputSCTHitCollectionName);
+ ATH_MSG_DEBUG("jetPtThreshold [MeV] = " << m_jetPtThreshold);
+ ATH_MSG_DEBUG("dRmatchHitToJet = " << m_dRmatchHitToJet);
+
+ ATH_CHECK(m_inputPixHitCollectionName.initialize());
+ ATH_CHECK(m_inputSCTHitCollectionName.initialize());
+ ATH_CHECK(m_jetCollectionName.initialize());
+ ATH_CHECK(m_vertexCollectionName.initialize());
+ ATH_CHECK(m_outputPixHitCollectionName.initialize());
+ ATH_CHECK(m_outputSCTHitCollectionName.initialize());
+
+ ATH_MSG_DEBUG("JetHitAssociation: Initialized");
+
+ return StatusCode::SUCCESS;
+}
+
+
+// Function for saving hits matched to jets
+StatusCode JetHitAssociation::execute() {
+
+ const EventContext& ctx = Gaudi::Hive::currentContext();
+
+ // All jets
+ SG::ReadHandle<xAOD::JetContainer> jetReadHandle(m_jetCollectionName, ctx);
+ if ( !jetReadHandle.isValid() ) {
+ ATH_MSG_ERROR("Failed to retrieve jet container with key " << m_jetCollectionName.key() );
+ return StatusCode::FAILURE;
+ }
+
+ // Keep only jets above threashold
+ std::vector<const xAOD::Jet*> jets;
+ jets.reserve( jetReadHandle->size() );
+ for (const xAOD::Jet *jet : *jetReadHandle) {
+ if (jet->pt() > m_jetPtThreshold) jets.push_back(jet);
+ }
+
+// Create write handles
+ SG::WriteHandle<xAOD::TrackMeasurementValidationContainer> outputPixHits(m_outputPixHitCollectionName, ctx);
+ SG::WriteHandle<xAOD::TrackMeasurementValidationContainer> outputSCTHits(m_outputSCTHitCollectionName, ctx);
+
+
+ // Vertex collection
+ SG::ReadHandle<xAOD::VertexContainer> vertexReadHandle(m_vertexCollectionName, ctx);
+ if ( !vertexReadHandle.isValid() ) {
+ ATH_MSG_ERROR("Failed to retrieve PrimaryVertices container" );
+ return StatusCode::FAILURE;
+ }
+
+ // Find primary vertex
+ const xAOD::Vertex *primVtx = nullptr;
+ for (const xAOD::Vertex *vertex : *vertexReadHandle) {
+ if (vertex->vertexType() == xAOD::VxType::PriVtx) {
+ primVtx = vertex;
+ break;
+ }
+ }
+
+ // Read pixel hits
+ SG::ReadHandle<xAOD::TrackMeasurementValidationContainer> pixHitReadHandle(m_inputPixHitCollectionName, ctx);
+ if ( !pixHitReadHandle.isValid() ) {
+ ATH_MSG_ERROR("Failed to retrieve pixel hit container with key " << m_inputPixHitCollectionName);
+ return StatusCode::FAILURE;
+ }
+
+ // Read SCT hits
+ SG::ReadHandle<xAOD::TrackMeasurementValidationContainer> sctHitReadHandle(m_inputSCTHitCollectionName, ctx);
+ if ( !sctHitReadHandle.isValid() ) {
+ ATH_MSG_ERROR("Failed to retrieve SCT hit container with key " << m_inputSCTHitCollectionName);
+ return StatusCode::FAILURE;
+ }
+
+ // These subroutines check for viable jets and write out hits if any exist.
+ // A JetAssociatedPixelCluster or JetAssociatedSCTCluster container is created for each event regardless
+ ATH_CHECK( saveHits(jets, primVtx, pixHitReadHandle, outputPixHits, m_nStoredPixHits, m_nTotalPixHits) );
+ ATH_CHECK( saveHits(jets, primVtx, sctHitReadHandle, outputSCTHits, m_nStoredSCTHits, m_nTotalSCTHits) );
+
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode JetHitAssociation::saveHits(const std::vector<const xAOD::Jet*> &jets,
+ const xAOD::Vertex* const vertex,
+ SG::ReadHandle<xAOD::TrackMeasurementValidationContainer> &hits,
+ const SG::WriteHandle<xAOD::TrackMeasurementValidationContainer> &outputHits,
+ unsigned long long int &nStoredHits,
+ unsigned long long int &nTotalHits) {
+
+ // Collection where output hits are stored
+ auto outputHitCollection = std::make_unique<xAOD::TrackMeasurementValidationContainer>();
+ auto outputHitCollectionAux = std::make_unique<xAOD::TrackMeasurementValidationAuxContainer>();
+ outputHitCollection->setStore(outputHitCollectionAux.get());
+
+
+ // if there are no jets you need to skip this loop but you still need to write out the
+ // if the vertex is a null pointer (no primary vertex in the event) do the same
+ // JetAssociatedSCTClusters and JetAssociatedPixelClusters
+ if (!jets.empty() && vertex!=nullptr) {
+
+ // Get the x,y,z of the primary vertex
+ TVector3 PVposition(vertex->x(), vertex->y(), vertex->z());
+ ATH_MSG_DEBUG("JetHitAssociation: PrimaryVertex Z Position = " << PVposition.Z());
+
+ // Loop over hits
+ for (const xAOD::TrackMeasurementValidation *hit : *hits) {
+ // Get pixel globalX,Y,Z
+ float x = hit->globalX();
+ float y = hit->globalY();
+ float z = hit->globalZ();
+
+ // Put the x,y,z in a TVector3
+ TVector3 hitPosition(x,y,z);
+
+ // Correct for the PV position
+ hitPosition = hitPosition - PVposition;
+
+ // Loop over jets
+ for (const xAOD::Jet *jet : jets) {
+ // Calculate dR(hit,jet)
+ float dR = hitPosition.DeltaR(jet->p4().Vect());
+
+ // if the dR is smaller than the dR association threshold save the hit
+ if (dR < m_dRmatchHitToJet) {
+ outputHitCollection->push_back(std::make_unique<xAOD::TrackMeasurementValidation>(*hit));
+ *outputHitCollection->back() = *hit; //Without this line. entries are filled with zero instead of the correct values
+ ++nStoredHits;
+ break;
+ }
+ }
+ }
+ }
+ // count the total number of hits run
+ nTotalHits += hits->size();
+
+ // Write output container
+ if ( ! outputHits.put( std::move(outputHitCollection), std::move(outputHitCollectionAux) ) ) {
+ ATH_MSG_ERROR("Could not write output hit containers");
+ return StatusCode::FAILURE;
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode JetHitAssociation::finalize() {
+
+ ATH_MSG_DEBUG("JetHitAssociation: Total number of Pix Hits tested = " << m_nTotalPixHits);
+ ATH_MSG_DEBUG("JetHitAssociation: Total number of Pix Hits Stored = " << m_nStoredPixHits);
+ ATH_MSG_DEBUG("JetHitAssociation: Total number of SCT Hits tested = " << m_nTotalSCTHits);
+ ATH_MSG_DEBUG("JetHitAssociation: Total number of SCT Hits Stored = " << m_nStoredSCTHits);
+
+ return StatusCode::SUCCESS;
+}
diff --git a/PhysicsAnalysis/JetTagging/JetHitAssociation/src/components/JetHitAssociation_entries.cxx b/PhysicsAnalysis/JetTagging/JetHitAssociation/src/components/JetHitAssociation_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..95dc9f6cf29b08f390c4ca63cca62a0ce64153e2
--- /dev/null
+++ b/PhysicsAnalysis/JetTagging/JetHitAssociation/src/components/JetHitAssociation_entries.cxx
@@ -0,0 +1,3 @@
+#include "JetHitAssociation/JetHitAssociation.h"
+
+DECLARE_COMPONENT( JetHitAssociation )
diff --git a/PhysicsAnalysis/JetTagging/JetTagAlgs/BTagging/python/BTaggingFlags.py b/PhysicsAnalysis/JetTagging/JetTagAlgs/BTagging/python/BTaggingFlags.py
index bdafa3dd308f5b27b7da12f64192222864b3f41b..38119c60bb648415edbb3bb491cf0b44f251e065 100755
--- a/PhysicsAnalysis/JetTagging/JetTagAlgs/BTagging/python/BTaggingFlags.py
+++ b/PhysicsAnalysis/JetTagging/JetTagAlgs/BTagging/python/BTaggingFlags.py
@@ -24,6 +24,8 @@ class _BTaggingFlags:
btaggingESDList = [ ]
+ DoJetHitAssociation = False
+
def __init__ (self):
self.btaggingAODList = list()
diff --git a/PhysicsAnalysis/JetTagging/JetTagAlgs/BTagging/share/BTaggingReconstructionOutputAODList_jobOptions.py b/PhysicsAnalysis/JetTagging/JetTagAlgs/BTagging/share/BTaggingReconstructionOutputAODList_jobOptions.py
index cf2e3604082cb83dabbcb497404b7acc6ed30c2a..5bad5737826fb6d0bd40671d65b0c704663ae30f 100644
--- a/PhysicsAnalysis/JetTagging/JetTagAlgs/BTagging/share/BTaggingReconstructionOutputAODList_jobOptions.py
+++ b/PhysicsAnalysis/JetTagging/JetTagAlgs/BTagging/share/BTaggingReconstructionOutputAODList_jobOptions.py
@@ -28,6 +28,12 @@ if len(BTaggingAODList) == 0:
# BTaggingAODList += ["xAOD::BTagVertexAuxContainer#*"]
# BTaggingAODList += ["xAOD::VertexContainer#BTagging*"]
# BTaggingAODList += ["xAOD::VertexAuxContainer#BTagging*"]
+
+if BTaggingFlags.DoJetHitAssociation:
+ BTaggingAODList += ['xAOD::TrackMeasurementValidationContainer#JetAssociatedPixelClusters',
+ 'xAOD::TrackMeasurementValidationAuxContainer#JetAssociatedPixelClustersAux.']
+ BTaggingAODList += ['xAOD::TrackMeasurementValidationContainer#JetAssociatedSCTClusters',
+ 'xAOD::TrackMeasurementValidationAuxContainer#JetAssociatedSCTClustersAux.']
printfunc ("#BTAG# ESD output container list: " + str(BTaggingESDList))
printfunc ("#BTAG# AOD output container list: " + str(BTaggingAODList))
diff --git a/Reconstruction/RecExample/RecExCommon/share/CombinedRec_config.py b/Reconstruction/RecExample/RecExCommon/share/CombinedRec_config.py
index 79f8043311d194094c2652fc26ee8910f9c87b03..121797c1d4cbaab849b280d9b332d810d929eef7 100755
--- a/Reconstruction/RecExample/RecExCommon/share/CombinedRec_config.py
+++ b/Reconstruction/RecExample/RecExCommon/share/CombinedRec_config.py
@@ -157,6 +157,14 @@ if (jetOK or rec.readESD()) and rec.doBTagging() and DetFlags.ID_on() and DetFl
Configurable.configurableRun3Behavior=0
pass
+# Hits associated with high-pt jets for trackless b-tagging
+from BTagging.BTaggingFlags import BTaggingFlags
+if (jetOK or rec.readESD()) and DetFlags.ID_on() and rec.doWriteAOD() and BTaggingFlags.DoJetHitAssociation:
+ try:
+ include("JetHitAssociation/jetHitAssociation_config.py")
+ except Exception:
+ treatException("Could not set up jet hit association")
+
#
# functionality : tau identification
#