Skip to content
Snippets Groups Projects
Commit b3a4dd90 authored by Fang-Ying Tsai's avatar Fang-Ying Tsai Committed by Frank Winklmeier
Browse files

implement filter logic and integrate ML decision in trackoverlay

implement filter logic and integrate ML decision in trackoverlay
parent bf73e42c
No related branches found
No related tags found
No related merge requests found
Showing
with 538 additions and 29 deletions
......@@ -385,6 +385,11 @@ def initConfigFlags():
return createTrackingConfigFlags()
_addFlagsCategory(acf, "Tracking", __tracking, 'TrkConfig')
def __trackoverlay():
from TrackOverlayConfig.TrackOverlayConfigFlags import createTrackOverlayConfigFlags
return createTrackOverlayConfigFlags()
_addFlagsCategory(acf, "TrackOverlay", __trackoverlay, 'TrackOverlayConfig')
def __acts():
from ActsConfig.ActsConfigFlags import createActsConfigFlags
return createActsConfigFlags()
......
......@@ -25,6 +25,8 @@ def ITkXAODToInDetClusterConversionCfg(flags, name="ITkXAODToInDetClusterConvers
def PixelClusterizationCfg(flags, name = "InDetPixelClusterization", **kwargs):
acc = ComponentAccumulator()
doTrackOverlay = getattr(flags.TrackOverlay, "ActiveConfig.doTrackOverlay", None) or flags.Overlay.doTrackOverlay
prefix = 'Sig_' if doTrackOverlay else ''
if "clusteringTool" not in kwargs:
from InDetConfig.SiClusterizationToolConfig import MergedPixelsToolCfg
......@@ -36,14 +38,14 @@ def PixelClusterizationCfg(flags, name = "InDetPixelClusterization", **kwargs):
kwargs.setdefault("gangedAmbiguitiesFinder", acc.popToolsAndMerge(
PixelGangedAmbiguitiesFinderCfg(flags)))
if not flags.Overlay.doTrackOverlay:
if not doTrackOverlay:
kwargs.setdefault("DataObjectName", "PixelRDOs")
else:
#for track overlay, only run tracking on the HS RDOs
kwargs.setdefault("DataObjectName", flags.Overlay.SigPrefix + "PixelRDOs")
kwargs.setdefault("ClustersName", "PixelClusters")
acc.addEventAlgo(CompFactory.InDet.PixelClusterization(name, **kwargs))
acc.addEventAlgo(CompFactory.InDet.PixelClusterization(prefix+name, **kwargs))
return acc
def PixelClusterizationPUCfg(flags, name="InDetPixelClusterizationPU", **kwargs):
......@@ -114,6 +116,8 @@ def ITkTrigPixelClusterizationCfg(flags, name = "ITkTrigPixelClusterization", ro
def SCTClusterizationCfg(flags, name="InDetSCT_Clusterization", **kwargs):
acc = ComponentAccumulator()
doTrackOverlay = getattr(flags.TrackOverlay, "ActiveConfig.doTrackOverlay", None) or flags.Overlay.doTrackOverlay
prefix = 'Sig_' if doTrackOverlay else ''
if "conditionsTool" not in kwargs:
from SCT_ConditionsTools.SCT_ConditionsToolsConfig import SCT_ConditionsSummaryToolCfg
......@@ -130,14 +134,14 @@ def SCTClusterizationCfg(flags, name="InDetSCT_Clusterization", **kwargs):
kwargs.setdefault("clusteringTool", acc.popToolsAndMerge(
SCT_ClusteringToolCfg(flags)))
if not flags.Overlay.doTrackOverlay:
if not doTrackOverlay:
kwargs.setdefault("DataObjectName", 'SCT_RDOs')
else:
#for track overlay, only run tracking on the HS RDOs
kwargs.setdefault("DataObjectName", flags.Overlay.SigPrefix + "SCT_RDOs")
kwargs.setdefault("ClustersName", 'SCT_Clusters')
acc.addEventAlgo(CompFactory.InDet.SCT_Clusterization(name, **kwargs))
acc.addEventAlgo(CompFactory.InDet.SCT_Clusterization(prefix+name, **kwargs))
return acc
def SCTClusterizationPUCfg(flags, name="InDetSCT_ClusterizationPU", **kwargs):
......@@ -211,20 +215,20 @@ def ITkTrigStripClusterizationCfg(flags, name="ITkTrigStripClusterization", rois
def InDetTRT_RIO_MakerCfg(flags, name = "InDetTRT_RIO_Maker", **kwargs):
acc = ComponentAccumulator()
doTrackOverlay = getattr(flags.TrackOverlay, "ActiveConfig.doTrackOverlay", None) or flags.Overlay.doTrackOverlay
prefix = 'Sig_' if doTrackOverlay else ''
if "TRT_DriftCircleTool" not in kwargs:
from InDetConfig.TRT_DriftCircleToolConfig import TRT_DriftCircleToolCfg
kwargs.setdefault("TRT_DriftCircleTool", acc.popToolsAndMerge(
TRT_DriftCircleToolCfg(flags)))
if not flags.Overlay.doTrackOverlay:
kwargs.setdefault("TRTRIOLocation", 'TRT_DriftCircles')
if not doTrackOverlay:
kwargs.setdefault("TRTRDOLocation", 'TRT_RDOs')
else:
#for track overlay, only run tracking on the HS RDOs
kwargs.setdefault("TRTRDOLocation", flags.Overlay.SigPrefix + 'TRT_RDOs')
kwargs.setdefault("TRTRIOLocation", 'TRT_DriftCircles')
acc.addEventAlgo(CompFactory.InDet.TRT_RIO_Maker(name, **kwargs))
acc.addEventAlgo(CompFactory.InDet.TRT_RIO_Maker(prefix+name, **kwargs))
return acc
def InDetTRT_NoTime_RIO_MakerCfg(flags, name = "InDetTRT_NoTime_RIO_Maker", **kwargs):
......@@ -295,3 +299,4 @@ def AthenaTrkClusterizationCfg(flags):
acc.merge(ITkStripClusterizationCfg(flags))
return acc
......@@ -287,6 +287,7 @@ def TRTStandalonePassRecoCfg(flags,
def StoreTrackSeparateContainerCfg(flags, TrackContainer="",
ClusterSplitProbContainer=""):
result = ComponentAccumulator()
doTrackOverlay = getattr(flags.TrackOverlay, "ActiveConfig.doTrackOverlay", None) or flags.Overlay.doTrackOverlay
# Dummy Merger to fill additional info
# for PRD-associated pixel tracklets
......@@ -295,17 +296,18 @@ def StoreTrackSeparateContainerCfg(flags, TrackContainer="",
AssociationMapName = ""
extension = flags.Tracking.ActiveConfig.extension
doTrackOverlay = getattr(flags.TrackOverlay, "ActiveConfig.doTrackOverlay", None) or flags.Overlay.doTrackOverlay
if extension == "Disappearing" or flags.Overlay.doTrackOverlay:
if extension == "Disappearing" or doTrackOverlay:
if extension == "Disappearing":
InputTracks = [TrackContainer]
if flags.Overlay.doTrackOverlay:
if doTrackOverlay:
InputTracks += [flags.Overlay.BkgPrefix +
extension + "Tracks"]
TrackContainer = extension+"Tracks"
AssociationMapName = "PRDtoTrackMap" + TrackContainer
MergerOutputTracks = TrackContainer
elif flags.Overlay.doTrackOverlay:
elif doTrackOverlay:
# schedule merger to combine signal and background tracks
InputTracks = [flags.Overlay.SigPrefix+TrackContainer,
flags.Overlay.BkgPrefix+TrackContainer]
......@@ -316,7 +318,8 @@ def StoreTrackSeparateContainerCfg(flags, TrackContainer="",
from TrkConfig.TrkTrackCollectionMergerConfig import (
TrackCollectionMergerAlgCfg)
result.merge(TrackCollectionMergerAlgCfg(
flags, "InDetTrackCollectionMerger"+extension,
flags,
name = "TrackCollectionMergerAlgCfg"+extension,
InputCombinedTracks=InputTracks,
OutputCombinedTracks=MergerOutputTracks,
AssociationMapName=AssociationMapName))
......@@ -338,7 +341,8 @@ def StoreTrackSeparateContainerCfg(flags, TrackContainer="",
from xAODTrackingCnv.xAODTrackingCnvConfig import (
TrackParticleCnvAlgPIDCheckCfg)
result.merge(TrackParticleCnvAlgPIDCheckCfg(
flags, extension + "TrackParticleCnvAlg",
flags,
name = extension + "TrackParticleCnvAlg",
TrackContainerName=TrackContainer,
xAODTrackParticlesFromTracksContainerName=xAODTrackParticlesName,
ClusterSplitProbabilityName=ClusterSplitProbContainer,
......@@ -395,7 +399,7 @@ def TrackRecoPassCfg(flags, extension="",
StatTrackTruthCollections=None,
ClusterSplitProbContainer=""):
result = ComponentAccumulator()
doTrackOverlay = getattr(flags.TrackOverlay, "ActiveConfig.doTrackOverlay", None) or flags.Overlay.doTrackOverlay
if InputCombinedInDetTracks is None:
InputCombinedInDetTracks = []
if InputExtendedInDetTracks is None:
......@@ -409,7 +413,7 @@ def TrackRecoPassCfg(flags, extension="",
# for track overlay, save resolved track name
# for final merged track collection
if (flags.Overlay.doTrackOverlay and
if (doTrackOverlay and
flags.Tracking.ActiveConfig.storeSeparateContainer and
not flags.Tracking.ActiveConfig.useTRTExtension):
ResolvedTracks = flags.Overlay.SigPrefix + ResolvedTracks
......@@ -436,7 +440,7 @@ def TrackRecoPassCfg(flags, extension="",
ResolvedTracks+"TruthCollection"]
TrackContainer = ResolvedTracks
if (flags.Overlay.doTrackOverlay and
if (doTrackOverlay and
flags.Tracking.ActiveConfig.storeSeparateContainer):
TrackContainer = "Resolved" + extension + "Tracks"
......@@ -451,7 +455,7 @@ def TrackRecoPassCfg(flags, extension="",
ExtendedTracks = "ExtendedTracksDisappearing"
elif "LargeD0" in extension:
ExtendedTracks = "ExtendedLargeD0Tracks"
if flags.Overlay.doTrackOverlay:
if doTrackOverlay:
ExtendedTracks = flags.Overlay.SigPrefix+"ExtendedLargeD0Tracks"
ExtendedTracksMap = "ExtendedTracksMap" + extension
......@@ -463,7 +467,7 @@ def TrackRecoPassCfg(flags, extension="",
ExtendedTracksMap=ExtendedTracksMap))
TrackContainer = ExtendedTracks
if flags.Overlay.doTrackOverlay and "LargeD0" in extension:
if doTrackOverlay and "LargeD0" in extension:
TrackContainer = "ExtendedLargeD0Tracks"
StatTrackCollections += [ExtendedTracks]
StatTrackTruthCollections += [ExtendedTracks+"TruthCollection"]
......@@ -495,6 +499,7 @@ def TrackFinalCfg(flags,
StatTrackCollections=None,
StatTrackTruthCollections=None):
result = ComponentAccumulator()
doTrackOverlay = getattr(flags.TrackOverlay, "ActiveConfig.doTrackOverlay", None) or flags.Overlay.doTrackOverlay
if InputCombinedInDetTracks is None:
InputCombinedInDetTracks = []
......@@ -503,7 +508,7 @@ def TrackFinalCfg(flags,
if StatTrackTruthCollections is None:
StatTrackTruthCollections = []
if flags.Overlay.doTrackOverlay:
if doTrackOverlay:
InputCombinedInDetTracks += [flags.Overlay.BkgPrefix +
"CombinedInDetTracks"]
......@@ -987,3 +992,4 @@ if __name__ == "__main__":
sc = top_acc.run(1)
if sc.isFailure():
sys.exit(-1)
......@@ -9,6 +9,7 @@ atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )
atlas_install_joboptions( share/skeleton.RDOtoRDOtrigger*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )
# Install other job opts without flake8 check
atlas_install_joboptions( share/*.py EXCLUDE share/*RDOtoRDOtrigger*.py )
# Install scripts
atlas_install_runtime( scripts/*.py )
......
......@@ -19,6 +19,7 @@ _all_domains = [
"AFP",
"HI",
"PostProcessing",
"TrackOverlay",
]
......@@ -56,6 +57,8 @@ def createRecoConfigFlags():
flags.addFlag("Reco.EnableCombinedMuon",
lambda prevFlags: prevFlags.Detector.EnableMuon and
prevFlags.Reco.EnableTracking)
# Enable TrackOverlay Reconstruction
flags.addFlag("Reco.EnableTrackOverlay", lambda prevFlags: False)
# Enable PFlow Reconstruction
flags.addFlag("Reco.EnablePFlow", lambda prevFlags: (
prevFlags.Reco.EnableTracking
......@@ -122,7 +125,6 @@ def createRecoConfigFlags():
# Enable ZDC reconstruction if ZDC data is in Bytestream or simulated
flags.addFlag("Reco.EnableZDC",_recoZDC)
# Enable common thinning and other post-processing
flags.addFlag("Reco.EnablePostProcessing", True)
flags.addFlag("Reco.PostProcessing.ThinNegativeClusters",
......
......@@ -3,7 +3,6 @@
from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
from AthenaConfiguration.Enums import Format, MetadataCategory, HIMode
def RecoSteering(flags):
"""
Generates configuration of the reconstructions
......@@ -66,10 +65,15 @@ def RecoSteering(flags):
# ID / ITk
acc.flagPerfmonDomain('ID')
if flags.Reco.EnableTracking:
from InDetConfig.TrackRecoConfig import InDetTrackRecoCfg
acc.merge(InDetTrackRecoCfg(flags))
log.info("---------- Configured tracking")
if flags.Reco.EnableTrackOverlay:
from TrackOverlayRec.TrackOverlayRecoConfig import TrackOverlayRecoCfg
acc.merge(TrackOverlayRecoCfg(flags))
else:
from InDetConfig.TrackRecoConfig import InDetTrackRecoCfg
acc.merge(InDetTrackRecoCfg(flags))
log.info("---------- Configured tracking")
# HI
acc.flagPerfmonDomain('HI')
if flags.Reco.EnableHI:
......@@ -281,7 +285,7 @@ def RecoSteering(flags):
from PerfMonComps.PerfMonCompsConfig import PerfMonMTSvcCfg
acc.merge(PerfMonMTSvcCfg(flags))
log.info("---------- Configured PerfMon")
return acc
......
################################################################################
# Package: TrackOverlayConfig
################################################################################
# Declare the package name:
atlas_subdir( TrackOverlayConfig )
#install files from the package:
atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )
# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
from AthenaConfiguration.AthConfigFlags import AthConfigFlags
def createTrackOverlayConfigFlags():
mlcf = AthConfigFlags()
mlcf.addFlag("TrackOverlay.TrackOverlayConfig.doTrackOverlay", True)
mlcf.addFlag("TrackOverlay.MCOverlayConfig.doTrackOverlay", False)
mlcf.addFlag("TrackOverlay.MLThreshold", 0.74201)
return mlcf
if __name__ == "__main__":
flags = createTrackOverlayConfigFlags()
flags.dump()
# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
# Declare the package name:
atlas_subdir( TrackOverlayRec )
# External dependencies:
find_package( onnxruntime )
set(ALL_ONNX_LIBS ${onnxruntime_LIBRARY} ${onnxruntime_LIBRARIES} ${ONNXRUNTIME_LIBRARIES})
# Component(s) in the package:
atlas_add_component( TrackOverlayRec
src/*.cxx src/*.h
src/components/*.cxx
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${OOT_INCLUDE_DIRSRS} ${ONNXRUNTIME_INCLUDE_DIRS}
LINK_LIBRARIES ${ROOT_LIBRARIES} ${ONNXRUNTIME_LIBRARIES} EventBookkeeperToolsLib InDetPhysValMonitoringLib PathResolver AthOnnxruntimeServiceLib)
# Install python modules
atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )
# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
from AthenaPython.PyAthenaComps import Alg, StatusCode
class TestAlg (Alg):
def __init__(self, name):
super(TestAlg, self).__init__(name)
def execute(self):
self.msg.info("======Running "+self.name+"=========")
return StatusCode.Success
def DummyCfg(flags, **kwargs):
acc = ComponentAccumulator()
doTrackOverlay = getattr(flags.TrackOverlay, "ActiveConfig.doTrackOverlay", None)
from pprint import pprint
pprint(vars(flags))
acc.addEventAlgo(TestAlg(name=('Sig_' if doTrackOverlay else '')+'Dummy'))
return acc
# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
from AthenaConfiguration.ComponentFactory import CompFactory
def TrackOverlayDecisionAlgCfg(flags, name="TrackOverlayDecisionAlg", **kwargs):
acc = ComponentAccumulator()
kwargs.setdefault("MLThreshold", flags.TrackOverlay.MLThreshold)
acc.addEventAlgo(CompFactory.TrackOverlayDecisionAlg.TrackOverlayDecisionAlg(name, **kwargs))
return acc
def InvertedTrackOverlayDecisionAlgCfg(flags, name="InvertedTrackOverlayDecisionAlg", **kwargs):
kwargs.setdefault("InvertFilter", True)
return TrackOverlayDecisionAlgCfg(flags, name,**kwargs)
# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
from AthenaCommon.CFElements import seqAND, parOR
from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
from TrackOverlayRec.TrackOverlayEventFilterConfig import TrackOverlayDecisionAlgCfg, InvertedTrackOverlayDecisionAlgCfg
from InDetConfig.TrackRecoConfig import InDetTrackRecoCfg
def TrackOverlayRecoCfg(flags):
acc = ComponentAccumulator()
acc.addSequence(seqAND("MOSequence"), parentName='AthAlgSeq')
acc.merge(TrackOverlayDecisionAlgCfg(flags), sequenceName='MOSequence')
acc.addSequence(parOR('WorkMOSequence'), parentName='MOSequence')
flagsMO = flags.cloneAndReplace("TrackOverlay.ActiveConfig","TrackOverlay.MCOverlayConfig")
acc.merge(InDetTrackRecoCfg(flagsMO), sequenceName='WorkMOSequence')
acc.addSequence(seqAND("TOSequence"), parentName='AthAlgSeq')
acc.merge(InvertedTrackOverlayDecisionAlgCfg(flags), sequenceName='TOSequence')
acc.addSequence(parOR('WorkTOSequence'), parentName='TOSequence')
flagsTO = flags.cloneAndReplace("TrackOverlay.ActiveConfig","TrackOverlay.TrackOverlayConfig")
acc.merge(InDetTrackRecoCfg(flagsTO), sequenceName='WorkTOSequence')
return acc
/*
* * Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
* * */
//
#include <fstream>
#include <string>
#include <iostream>
#include <vector>
#include <sstream>
#include <unistd.h>
#include <Eigen/Core>
//
#include "GaudiKernel/SystemOfUnits.h"
#include "Gaudi/Property.h"
#include "EventInfo/EventInfo.h"
#include "EventInfo/EventID.h"
#include "xAODTruth/TruthParticle.h"
#include "xAODTruth/TruthPileupEvent.h"
#include "xAODTruth/TruthPileupEventAuxContainer.h"
#include "PathResolver/PathResolver.h"
// ONNX Runtime include(s).
#include <core/session/onnxruntime_cxx_api.h>
//
#include "TrackOverlayDecisionAlg.h"
#include "EventBookkeeperTools/FilterReporter.h"
namespace TrackOverlayDecisionAlg {
TrackOverlayDecisionAlg::TrackOverlayDecisionAlg( const std::string& name, ISvcLocator* pSvcLocator ) :
::AthReentrantAlgorithm( name, pSvcLocator )
{
}
StatusCode TrackOverlayDecisionAlg::initialize()
{
ATH_CHECK(m_filterParams.initialize(false));
ATH_CHECK( m_eventInfoContainerName.initialize() );
ATH_CHECK( m_truthParticleName.initialize( (m_pileupSwitch == "HardScatter" or m_pileupSwitch == "All") and not m_truthParticleName.key().empty() ) );
ATH_CHECK(m_truthSelectionTool.retrieve(EnableTool {not m_truthParticleName.key().empty()} ));
ATH_CHECK( m_truthEventName.initialize( (m_pileupSwitch == "HardScatter" or m_pileupSwitch == "All") and not m_truthEventName.key().empty() ) );
ATH_CHECK( m_truthPileUpEventName.initialize( (m_pileupSwitch == "PileUp" or m_pileupSwitch == "All") and not m_truthPileUpEventName.key().empty() ) );
ATH_CHECK(m_svc.retrieve());
std::string this_file = __FILE__;
const std::string model_path = PathResolverFindCalibFile("TrackOverlay/TrackOverlay_J7_model.onnx");
Ort::SessionOptions session_options;
m_session = std::make_unique<Ort::Session>(m_svc->env(), model_path.c_str(), session_options);
m_inputInfo = TrackOverlayDecisionAlg::GetInputNodeInfo(m_session);
m_outputInfo = TrackOverlayDecisionAlg::GetOutputNodeInfo(m_session);
return StatusCode::SUCCESS;
}
StatusCode TrackOverlayDecisionAlg::finalize()
{
ATH_MSG_VERBOSE ( "Finalizing ..." );
ATH_MSG_VERBOSE("-----------------------------------------------------------------");
ATH_MSG_VERBOSE("m_filterParams.summary()" << m_filterParams.summary());
ATH_MSG_VERBOSE("-----------------------------------------------------------------");
ATH_MSG_INFO(m_filterParams.summary());
ATH_MSG_VERBOSE(" =====================================================================");
return StatusCode::SUCCESS;
}
const std::vector<const xAOD::TruthParticle*> TrackOverlayDecisionAlg::getTruthParticles() const {
std::vector<const xAOD::TruthParticle*> tempVec {};
if (m_pileupSwitch == "All") {
if (m_truthParticleName.key().empty()) {
return tempVec;
}
SG::ReadHandle<xAOD::TruthParticleContainer> truthParticleContainer( m_truthParticleName);
if (not truthParticleContainer.isValid()) {
return tempVec;
}
tempVec.insert(tempVec.begin(), truthParticleContainer->begin(), truthParticleContainer->end());
} else {
if (m_pileupSwitch == "HardScatter") {
if (not m_truthEventName.key().empty()) {
ATH_MSG_VERBOSE("Getting TruthEvents container.");
SG::ReadHandle<xAOD::TruthEventContainer> truthEventContainer( m_truthEventName);
const xAOD::TruthEvent* event = (truthEventContainer.isValid()) ? truthEventContainer->at(0) : nullptr;
if (not event) {
return tempVec;
}
const auto& links = event->truthParticleLinks();
tempVec.reserve(event->nTruthParticles());
for (const auto& link : links) {
if (link.isValid()){
tempVec.push_back(*link);
}
}
}
}else if (m_pileupSwitch == "PileUp") {
if (not m_truthPileUpEventName.key().empty()) {
ATH_MSG_VERBOSE("getting TruthPileupEvents container");
// get truth particles from all pileup events
SG::ReadHandle<xAOD::TruthPileupEventContainer> truthPileupEventContainer(m_truthPileUpEventName);
if (truthPileupEventContainer.isValid()) {
const unsigned int nPileup = truthPileupEventContainer->size();
tempVec.reserve(nPileup * 200); // quick initial guess, will still save some time
for (unsigned int i(0); i != nPileup; ++i) {
const auto *eventPileup = truthPileupEventContainer->at(i);
// get truth particles from each pileup event
int ntruth = eventPileup->nTruthParticles();
ATH_MSG_VERBOSE("Adding " << ntruth << " truth particles from TruthPileupEvents container");
const auto& links = eventPileup->truthParticleLinks();
for (const auto& link : links) {
if (link.isValid()){
tempVec.push_back(*link);
}
}
}
} else {
ATH_MSG_ERROR("no entries in TruthPileupEvents container!");
}
}
} else {
ATH_MSG_ERROR("bad value for PileUpSwitch");
}
}
return tempVec;
}
StatusCode TrackOverlayDecisionAlg::execute(const EventContext &ctx) const
{
ATH_MSG_DEBUG ("Executing ...");
std::vector<const xAOD::TruthParticle*> truthParticlesVec = TrackOverlayDecisionAlg::getTruthParticles();
//Access truth info for the NN input
float eventPxSum = 0.0;
float eventPySum = 0.0;
float eventPt = 0.0;
float puEvents = 0.0;
std::vector<float> pxValues, pyValues, pzValues, eValues, etaValues, phiValues, ptValues;
float truthMultiplicity = 0.0;
const int truthParticles = truthParticlesVec.size();
for (int itruth = 0; itruth < truthParticles; itruth++) {
const xAOD::TruthParticle* thisTruth = truthParticlesVec[itruth];
const IAthSelectionTool::CutResult accept = m_truthSelectionTool->accept(thisTruth);
if(accept){
pxValues.push_back((thisTruth->px()*0.001-1.46988000e+03)* px_diff); //as MinMaxScaler: 1.46988000e+03 is the lowest value of px from a J7 sample; *(0.001) is used to convert unit rather than *(1/1000) to speed up.
pyValues.push_back((thisTruth->py()*0.001-1.35142000e+03)* py_diff); //the lowest value of py: 1.35142000e+03
pzValues.push_back((thisTruth->pz()*0.001-1.50464000e+03)* pz_diff); //the lowest value of pz: 1.50464000e+03
ptValues.push_back((thisTruth->pt()*0.001-5.00006000e-01)* pt_diff); //the lowest value of pt: 5.00006000e-01
etaValues.push_back(thisTruth->eta());
phiValues.push_back(thisTruth->phi());
eValues.push_back((thisTruth->e()*0.001-5.08307000e-01)*e_diff); //the lowest value of energy: 5.08307000e-01
eventPxSum += thisTruth->px();
eventPySum += thisTruth->py();
truthMultiplicity++;
}//accept
}//for itruth
SG::ReadHandle<xAOD::TruthPileupEventContainer> truthPileupEventContainer;
SG::ReadHandle<xAOD::EventInfo> pie = SG::ReadHandle<xAOD::EventInfo>(m_eventInfoContainerName, ctx);
if (!m_truthPileUpEventName.key().empty()) {
truthPileupEventContainer = SG::ReadHandle<xAOD::TruthPileupEventContainer>(m_truthPileUpEventName, ctx);
}
puEvents = !m_truthPileUpEventName.key().empty() and truthPileupEventContainer.isValid() ? static_cast<int>( truthPileupEventContainer->size() ) : pie.isValid() ? pie->actualInteractionsPerCrossing() : 0;
eventPt = std::sqrt(eventPxSum*eventPxSum + eventPySum*eventPySum)*0.001;
std::vector<float> puEventsVec(pxValues.size(), (puEvents-1.55000000e+01)*pu_diff); //min of puEvents= 15.5, max of puEvents=84.5
std::vector<float> truthMultiplicityVec(pxValues.size(), (truthMultiplicity-1.80000000e+01)*multi_diff);
std::vector<float> eventPtVec(pxValues.size(), (eventPt-3.42359395e-01)*eventPt_diff);
std::vector<float> predictions;
//Compute the distances using Eigen for Eigen's optimized operations. Initialize matirces. Observed a significant improvement on computing calculation.
Eigen::VectorXf ptEigen = Eigen::VectorXf::Map(ptValues.data(), ptValues.size());
Eigen::VectorXf phiEigen = Eigen::VectorXf::Map(phiValues.data(), phiValues.size());
Eigen::VectorXf etaEigen = Eigen::VectorXf::Map(etaValues.data(), etaValues.size());
for (std::size_t i = 0; i < truthMultiplicity; ++i) {
float multiplicity_0p05 = 0.0, multiplicity_0p2 = 0.0;
float sum_0p05 = 0.0, sum_0p2 = 0.0;
float pt_0p05 = 0.0, pt_0p2 = 0.0;
float deltaEtaI = etaEigen[i];
float phiI = phiEigen[i];
for (std::size_t j = 0; j < truthMultiplicity; ++j) {
if (i == j) continue; // Skip the particle itself
float deltaEta = deltaEtaI - etaEigen[j];
float deltaPhi = phiI - phiEigen[j];
if (deltaPhi > M_PI) {
deltaPhi -= 2.0 * M_PI;
}
float distances = std::sqrt(deltaEta * deltaEta + deltaPhi * deltaPhi);
if (distances < 0.05){
multiplicity_0p05++;
sum_0p05 += distances;
pt_0p05 += ptEigen[j];
}
if (distances < 0.2){
multiplicity_0p2++;
sum_0p2 += distances;
pt_0p2 += ptEigen[j];
}
}// for j
std::vector<float> featData;
featData.push_back(pxValues[i]);
featData.push_back(pyValues[i]);
featData.push_back(pzValues[i]);
featData.push_back(eValues[i]);
featData.push_back(ptValues[i]);
featData.push_back((multiplicity_0p2 * area0p2) * constant1);
featData.push_back((multiplicity_0p05 * area0p05) * constant2);
featData.push_back((sum_0p2 * area0p2) * constant3);
featData.push_back((sum_0p05 * area0p05) * constant4);
featData.push_back(pt_0p2 * constant5);
featData.push_back(pt_0p05 * constant6);
featData.push_back(puEventsVec[i]);
featData.push_back(truthMultiplicityVec[i]);
featData.push_back(eventPtVec[i]);
std::vector<int64_t> input_node_dims;
std::vector<char*> input_node_names;
input_node_dims = std::get<0>(m_inputInfo);
input_node_names = std::get<1>(m_inputInfo);
std::vector<int64_t> output_node_dims;
std::vector<char*> output_node_names;
output_node_dims = std::get<0>(m_outputInfo);
output_node_names = std::get<1>(m_outputInfo);
Ort::MemoryInfo memoryInfo = Ort::MemoryInfo::CreateCpu(OrtArenaAllocator, OrtMemTypeCPU);
input_node_dims[0]=1;
Ort::Value input_data = Ort::Value::CreateTensor(memoryInfo, featData.data(), featData.size(), input_node_dims.data(), input_node_dims.size());
Ort::RunOptions run_options(nullptr);
//Run the inference
Ort::Session& mysession ATLAS_THREAD_SAFE = *m_session;
auto output_values = mysession.Run(run_options, input_node_names.data(), &input_data, input_node_names.size(), output_node_names.data(), output_node_names.size());
float* predictionData = output_values[0].GetTensorMutableData<float>();
float prediction = predictionData[0];
predictions.push_back(prediction);
}//for i
float threshold = m_MLthreshold;
ATH_MSG_ALWAYS("ML threshold:" << threshold);
int badTracks = 0;
for (float prediction : predictions) {
if (prediction > threshold) {
badTracks++;
}
}
float rouletteScore = static_cast<float>(badTracks) / static_cast<float>(truthMultiplicity);
FilterReporter filter(m_filterParams, false, ctx);
bool pass = false;
int decision = rouletteScore == 0;
if (decision==0){ //if ML decision is False, it goes to the MC-overlay workflow
pass = true;
}
else{
pass = false;
}
if (m_invertfilter) {
pass =! pass;
}
filter.setPassed(pass);
ATH_MSG_ALWAYS("End TrackOverlayDecisionAlg, difference in filters: "<<(pass ? "found" : "not found")<<"="<<pass<<", invert="<<m_invertfilter);
return StatusCode::SUCCESS;
}
}// end namespace TrackOverlayDecisionAlg
/*
* Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
* */
#ifndef TRACKOVERLAYREC_TRACKOVERLAYDECISIONALG_H
#define TRACKOVERLAYREC_TRACKOVERLAYDECISIONALG_H
// STL includes
#include <string>
// FrameWork includes
#include "StoreGate/ReadHandle.h"
#include "GaudiKernel/ToolHandle.h"
#include "AthenaBaseComps/AthReentrantAlgorithm.h"
#include <EventBookkeeperTools/FilterReporterParams.h>
// local includes
#include "InDetPhysValMonitoring/IAthSelectionTool.h"
//#gaudi includes
#include "GaudiKernel/ToolHandle.h"
//
#include "xAODTruth/TruthParticleContainer.h"
#include "xAODTruth/TruthEventContainer.h"
#include "xAODTruth/TruthPileupEventContainer.h"
//ONNX Runtime include(s)
#include <core/session/onnxruntime_cxx_api.h>
#include "AthOnnxruntimeService/IONNXRuntimeSvc.h"
namespace TrackOverlayDecisionAlg{
const double M_TWOPI = 2.0 * M_PI;
// Calculate constants only once for feature scaling
const float px_diff = 0.00029734736; // (1/3363.07) px_max: 1893.19; px_min: -1.46988e+03
const float py_diff = 0.000335; //(1/2979.65) py_max: 1628.23; py_min: -1.35142e+03
const float pz_diff = 0.00025; //(1/3924.99) pz_max: 2420.35 pz_min: -1.50464e+03
const float pt_diff = 0.0005283; //(1/1892.74) pt_max: 1893.24 pt_min: 5.00006e-01
const float e_diff = 0.000372; //(1/2689.63) e_max: 2690.14 e_min: 5.08307e-01
const float pu_diff = 0.01449; //(1/69) pu_max: 84.5 pu_min: 15.5
const float multi_diff = 0.00342; //(1/292) multiplicity_max:310.0 min:1.80000000e+01
const float eventPt_diff = 0.0004798; //(1/2084.2) eventPt_max: 2084.61 min:3.42359395e-01
const float area0p2 = 127.388535032;//(1.0 / (M_PI * 0.05 * 0.05)) for R=0.05
const float area0p05 = 7.961783; // 1.0 / (M_PI * 0.2 * 0.2); for R=0.2
const float constant1 = 0.000417; //(1.0 / 2395.2); 1/max_of_multiplicity_0p2
const float constant2 = 0.00002626749; //(1.0 / 38069.86); 1/max_of_multiplicity_0p05
const float constant3 = 0.014797; //(1.0 / 67.58); 1/max_of_sum_0p2
const float constant4 = 0.006934; //(1.0 / 144.21); 1/max_of_sum_0p05
const float constant5 = 0.0004086; //(1.0 / 2447.24); 1/max_of_sumPtAroundTrack for R=0.2
const float constant6 = 0.0004186; //(1.0 / 2388.55); 1/max_of_sumPtAroundTrack for R=0.05
class TrackOverlayDecisionAlg : public AthReentrantAlgorithm{
public:
/** Constructor with parameters */
TrackOverlayDecisionAlg( const std::string& name, ISvcLocator* pSvcLocator );
/** Destructor */
virtual ~TrackOverlayDecisionAlg() = default;
/** Athena algorithm's interface method initialize() */
virtual StatusCode initialize() override final; /** Athena algorithm's interface method execute() */
virtual StatusCode execute(const EventContext& ctx) const override final;
/** Athena algorithm's interface method finalize() */
virtual StatusCode finalize() override final;
private:
ToolHandle<IAthSelectionTool> m_truthSelectionTool{this, "TruthSelectionTool","AthTruthSelectionTool", "Truth selection tool (for efficiencies and resolutions)"};
SG::AuxElement::Decorator<bool> m_dec_selectedByPileupSwitch{"selectedByPileupSwitch"};
bool m_usingSpecialPileupSwitch {false};
//set the "selectedByPileupSwitch" decoration for all particles in the passed vector
void markSelectedByPileupSwitch(const std::vector<const xAOD::TruthParticle*> & truthParticles) const;
BooleanProperty m_useTrackSelection {this, "useTrackSelection", false, "plot only tracks accepted by selection tool"};
StringProperty m_pileupSwitch {this, "PileupSwitch", "HardScatter", "Pileup truth strategy to use. May be \"All\", \"HardScatter\", or \"PileUp\""};
FloatProperty m_lowProb{this,"LowProb",0.5,"Truth match prob. cutoff for efficiency (lower bound) and fake (upper bound) classification."};
//EventInfo container name
SG::ReadHandleKey<xAOD::EventInfo> m_eventInfoContainerName{this,"EventInfoContainerName", "EventInfo", ""};
//TruthParticle container's name
SG::ReadHandleKey<xAOD::TruthPileupEventContainer> m_truthPileUpEventName{this, "TruthPileupEvents", "TruthPileupEvents","Name of the truth pileup events container probably TruthPileupEvent(s)"};
SG::ReadHandleKey<xAOD::TruthParticleContainer> m_truthParticleName{this, "TruthParticleContainerName", "TruthParticles", ""};
SG::ReadHandleKey<xAOD::TruthEventContainer> m_truthEventName{this, "TruthEvents", "TruthEvents","Name of the truth events container probably either TruthEvent or TruthEvents"};
// Get truth particles into a vector, possibly using the pileup from the event
const std::vector<const xAOD::TruthParticle *> getTruthParticles() const;
FilterReporterParams m_filterParams {this, "TrackOverlayDecisionAlg", "Decides whether events should be reconstructed in track-overlay workflow or MC-overlay."};
Gaudi::Property<bool> m_invertfilter{this, "InvertFilter", false, "Invert filter decision."}; //!< invert filter decision at the end
Gaudi::Property<float> m_MLthreshold{this, "MLThreshold", 0.74201, "ML threshold for bad/good tracks decision. ML scores larger than this threshold are considered as bad tracks."};
// Set up the ONNX Runtime session
ServiceHandle<AthONNX::IONNXRuntimeSvc> m_svc{this, "ONNXRuntimeSvc", "AthONNX::ONNXRuntimeSvc", "CaloMuonScoreTool ONNXRuntimeSvc"};
std::tuple<std::vector<int64_t>, std::vector<char*>> m_inputInfo;
std::tuple<std::vector<int64_t>, std::vector<char*>> m_outputInfo;
std::unique_ptr<Ort::Session> m_session;
inline std::tuple<std::vector<int64_t>, std::vector<char*> > GetInputNodeInfo(const std::unique_ptr< Ort::Session >& session) {
std::vector<int64_t> input_node_dims;
size_t num_input_nodes = session->GetInputCount();
std::vector<char*> input_node_names(num_input_nodes);
Ort::AllocatorWithDefaultOptions allocator;
for( std::size_t i = 0; i < num_input_nodes; i++ ) {
char* input_name = session->GetInputNameAllocated(i, allocator).release();
input_node_names[i] = input_name;
Ort::TypeInfo type_info = session->GetInputTypeInfo(i);
auto tensor_info = type_info.GetTensorTypeAndShapeInfo();
input_node_dims = tensor_info.GetShape();
}
return std::make_tuple(input_node_dims, input_node_names);
}
inline std::tuple<std::vector<int64_t>, std::vector<char*> > GetOutputNodeInfo(const std::unique_ptr< Ort::Session >& session){
std::vector<int64_t> output_node_dims;
size_t num_output_nodes = session->GetOutputCount();
std::vector<char*> output_node_names(num_output_nodes);
Ort::AllocatorWithDefaultOptions allocator;
for( std::size_t i = 0; i < num_output_nodes; i++ ) {
char* output_name = session->GetOutputNameAllocated(i, allocator).release();
output_node_names[i] = output_name;
Ort::TypeInfo type_info = session->GetOutputTypeInfo(i);
auto tensor_info = type_info.GetTensorTypeAndShapeInfo();
output_node_dims = tensor_info.GetShape();
}
return std::make_tuple(output_node_dims, output_node_names);
}
};
}
#endif
#include "../TrackOverlayDecisionAlg.h"
DECLARE_COMPONENT( TrackOverlayDecisionAlg::TrackOverlayDecisionAlg )
......@@ -9,6 +9,8 @@ def TrackCollectionMergerAlgCfg(flags, name="InDetTrackCollectionMerger",
OutputCombinedTracks="",
**kwargs):
result = ComponentAccumulator()
doTrackOverlay = getattr(flags.TrackOverlay, "ActiveConfig.doTrackOverlay", None) or flags.Overlay.doTrackOverlay
prefix = 'Sig_' if doTrackOverlay else ''
kwargs.setdefault("TracksLocation", InputCombinedTracks)
kwargs.setdefault("OutputTracksLocation", OutputCombinedTracks)
......@@ -17,9 +19,9 @@ def TrackCollectionMergerAlgCfg(flags, name="InDetTrackCollectionMerger",
from InDetConfig.InDetAssociationToolsConfig import InDetPRDtoTrackMapToolGangedPixelsCfg
kwargs.setdefault("AssociationTool", result.popToolsAndMerge(InDetPRDtoTrackMapToolGangedPixelsCfg(flags)))
kwargs.setdefault("DoTrackOverlay",flags.Overlay.doTrackOverlay)
kwargs.setdefault("DoTrackOverlay",doTrackOverlay)
result.addEventAlgo(CompFactory.Trk.TrackCollectionMerger(name, **kwargs))
result.addEventAlgo(CompFactory.Trk.TrackCollectionMerger(prefix+name, **kwargs))
return result
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment