diff --git a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
index 21e1d51a6af7d92ebac09d051e80c919d1c39687..a778d2eddf7a5759c578bbb1e74600e0aaab48e6 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
+++ b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
@@ -29,6 +29,8 @@ parser.add_argument("-r", "--reco", default="",
help="Specify reco tag (to append to output filename)")
parser.add_argument("-n", "--nevents", type=int, default=-1,
help="Specify number of events to process (default: all)")
+parser.add_argument("--skip", type=int, default=0,
+ help="Specify number of events to skip (default: none)")
parser.add_argument("-v", "--verbose", action='store_true',
help="Turn on DEBUG output")
parser.add_argument("--isMC", action='store_true',
@@ -37,6 +39,8 @@ parser.add_argument("--MC_calibTag", default="",
help="Specify tag used to reconstruct MC calo energy: (WAVE-Calibration-01-LG-nofilt, WAVE-Calibration-01-LG, WAVE-Calibration-01-HG-nofilt, or WAVE-Calibration-01-HG) ")
parser.add_argument("--testBeam", action='store_true',
help="Set geometry for 2021 test beam")
+parser.add_argument("--isOverlay", action='store_true',
+ help="Set overlaid data input")
args = parser.parse_args()
@@ -81,6 +85,8 @@ else:
print(f"Starting reconstruction of {filepath.name} with type {runtype}")
if args.nevents > 0:
print(f"Reconstructing {args.nevents} events by command-line option")
+if args.skip > 0:
+ print(f"Skipping {args.skip} events by command-line option")
# Start reconstruction
@@ -103,6 +109,8 @@ else:
ConfigFlags.Input.ProjectName = "data20"
ConfigFlags.GeoModel.Align.Dynamic = False
+ConfigFlags.Exec.SkipEvents = args.skip
+
# Flags for later
useCKF = True
useCal = False
@@ -173,7 +181,7 @@ acc.merge(PoolWriteCfg(ConfigFlags))
#
# Set up RAW data access
-if args.isMC:
+if args.isMC or args.isOverlay:
from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
acc.merge(PoolReadCfg(ConfigFlags))
else:
@@ -185,21 +193,23 @@ else:
from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
acc.merge(FaserGeometryCfg(ConfigFlags))
-if useLHC:
+if useLHC and not args.isOverlay:
from LHCDataAlgs.LHCDataAlgConfig import LHCDataAlgCfg
acc.merge(LHCDataAlgCfg(ConfigFlags))
# Set up algorithms
-from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionCfg
-acc.merge(WaveformReconstructionCfg(ConfigFlags))
+if not args.isOverlay:
+ from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionCfg
+ acc.merge(WaveformReconstructionCfg(ConfigFlags))
-# Calorimeter Energy reconstruction
-if useCal:
- from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionCfg
- acc.merge(CalorimeterReconstructionCfg(ConfigFlags, MC_calibTag=args.MC_calibTag))
+ # Calorimeter Energy reconstruction
+ if useCal:
+ from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionCfg
+ acc.merge(CalorimeterReconstructionCfg(ConfigFlags, MC_calibTag=args.MC_calibTag))
# Tracker clusters
from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
+# acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags, DataObjectName="Pos_SCT_RDOs"))
acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags, DataObjectName="SCT_RDOs"))
#
@@ -226,7 +236,8 @@ if useCKF:
#
# Kalman Filter for tracking
from FaserActsKalmanFilter.CKF2Config import CKF2Cfg
- acc.merge(CKF2Cfg(ConfigFlags, noDiagnostics=True))
+ if not args.isOverlay:
+ acc.merge(CKF2Cfg(ConfigFlags, noDiagnostics=True))
# Add tracking collection with no IFT
acc.merge(CKF2Cfg(ConfigFlags, maskedLayers=[0, 1, 2], name="CKF_woIFT",
@@ -246,10 +257,10 @@ itemList = [ "xAOD::EventInfo#*"
, "TrackCollection#*"
]
#
-if useLHC:
+if useLHC and not args.isOverlay:
itemList.extend( ["xAOD::FaserLHCData#*", "xAOD::FaserLHCDataAux#*"] )
-if args.isMC:
+if args.isMC and not args.isOverlay:
# Make xAOD versions of truth
from Reconstruction.xAODTruthCnvAlgConfig import xAODTruthCnvAlgCfg
acc.merge(xAODTruthCnvAlgCfg(ConfigFlags))
@@ -266,13 +277,14 @@ tagBuilder = CompFactory.EventInfoTagBuilder()
tagBuilder.PropagateInput=False
acc.addEventAlgo(tagBuilder)
-# Waveform reconstruction output
-from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionOutputCfg
-acc.merge(WaveformReconstructionOutputCfg(ConfigFlags))
+if not args.isOverlay:
+ # Waveform reconstruction output
+ from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionOutputCfg
+ acc.merge(WaveformReconstructionOutputCfg(ConfigFlags))
-# Calorimeter reconstruction output
-from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionOutputCfg
-acc.merge(CalorimeterReconstructionOutputCfg(ConfigFlags))
+ # Calorimeter reconstruction output
+ from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionOutputCfg
+ acc.merge(CalorimeterReconstructionOutputCfg(ConfigFlags))
# Check what we have
print( "Writing out xAOD objects:" )
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt b/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8d8665ba924345e2fdcdb409384cdb66291446c3
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
@@ -0,0 +1,12 @@
+atlas_subdir(OverlayRDO)
+
+atlas_add_component(
+ OverlayRDO
+ src/OverlayRDOAlg.h src/SelectRDOAlg.h src/BinRDOAlg.h src/HistoRDOAlg.h
+ src/OverlayRDOAlg.cxx src/SelectRDOAlg.cxx src/BinRDOAlg.cxx src/HistoRDOAlg.cxx
+ src/component/OverlayRDO_entries.cxx
+ LINK_LIBRARIES AthenaBaseComps StoreGateLib TrackerRawData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack xAODFaserWaveform
+)
+
+atlas_install_python_modules(python/*.py)
+# atlas_install_scripts(test/*.py)
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..cfb2dfefe48927be61f6416bd206a78e0b754eab
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py
@@ -0,0 +1,127 @@
+"""
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+# from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+def HistoRDOAlgCfg(flags, **kwargs):
+
+ from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+ acc = FaserGeometryCfg(flags)
+ HistoRDOAlg = CompFactory.HistoRDOAlg("HistoRDOAlg",**kwargs)
+ acc.addEventAlgo(HistoRDOAlg)
+ HistoRDOAlg.OutputLevel = VERBOSE
+ # ItemList = []
+ # # ItemList += ["xAOD::EventInfo#" + chargePrefix + "EventInfo"]
+ # # ItemList += ["xAOD::EventAuxInfo#" + chargePrefix + "EventInfoAux."]
+ # ItemList += ["xAOD::EventInfo#EventInfo"]
+ # ItemList += ["xAOD::EventAuxInfo#EventInfoAux."]
+ # ItemList += ["TrackCollection#Orig_CKFTrackCollectionWithoutIFT"]
+ # ItemList += ["FaserSCT_RDO_Container#SCT_RDOs"]
+ # ItemList += ["FaserSCT_RDO_Container#SCT_EDGEMODE_RDOs"]
+ # # ItemList += ["Tracker::FaserSCT_ClusterContainer#" + chargePrefix + "SCT_ClusterContainer"]
+
+ # acc.merge(OutputStreamCfg(ConfigFlags,"RDO", ItemList=ItemList, disableEventTag=True))
+
+ thistSvc = CompFactory.THistSvc()
+ thistSvc.Output += ["HIST2 DATAFILE='RDOtree.root' OPT='RECREATE'"]
+ acc.addService(thistSvc)
+
+ return acc
+
+if __name__ == "__main__":
+
+ import sys
+ from AthenaCommon.Logging import log, logging
+ from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+ from AthenaCommon.Configurable import Configurable
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ from AthenaConfiguration.TestDefaults import defaultTestFiles
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ # from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+ # Set up logging and new style config
+ log.setLevel(DEBUG)
+ Configurable.configurableRun3Behavior = True
+
+ # Configure
+ # ConfigFlags.Overlay.DataOverlay = False
+ ConfigFlags.Input.Files = [
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin0RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin1RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin2RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin3RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin4RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin5RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin6RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin7RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myPos_Bin8RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin0RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin1RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin2RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin3RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin4RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin5RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin6RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin7RDO.pool.root',
+ '/home/dcasper/Work/faser/darkphoton/overlay/myNeg_Bin8RDO.pool.root',
+ ]
+ # ConfigFlags.Input.SecondaryFiles = [ 'Neg_RDO.pool.root' ]
+ # ConfigFlags.Input.Files = [ '/eos/experiment/faser/rec/2022/r0013/009171/Faser-Physics-009171-00006-r0013-xAOD.root']
+ # ConfigFlags.Input.SecondaryFiles = [ '/eos/experiment/faser/rec/2022/r0013/009166/Faser-Physics-009166-00485-r0013-xAOD.root' ]
+ # ConfigFlags.Output.RDOFileName = "Overlay.RDO.pool.root"
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03" # Always needed; must match FaserVersionS
+ ConfigFlags.IOVDb.DatabaseInstance = "CONDBR3" # Use data conditions for now
+ ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
+ ConfigFlags.Input.isMC = False # Needed to bypass autoconfig
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
+ ConfigFlags.Common.isOnline = False
+ ConfigFlags.GeoModel.Align.Dynamic = False
+ ConfigFlags.Beam.NumberOfCollisions = 0.
+
+ ConfigFlags.Detector.GeometryFaserSCT = True
+
+ ConfigFlags.lock()
+
+ # Core components
+ acc = MainServicesCfg(ConfigFlags)
+ acc.merge(PoolReadCfg(ConfigFlags))
+
+ # algorithm
+ acc.merge(HistoRDOAlgCfg(ConfigFlags))
+
+ # from SGComps.AddressRemappingConfig import AddressRemappingCfg
+ # acc.merge(AddressRemappingCfg([
+ # "xAOD::EventInfo#EventInfo->" + ConfigFlags.Overlay.SigPrefix + "EventInfo",
+ # "xAOD::EventAuxInfo#EventInfoAux.->" + ConfigFlags.Overlay.SigPrefix + "EventInfoAux.",
+ # ]))
+
+ # Hack to avoid problem with our use of MC databases when isMC = False
+ replicaSvc = acc.getService("DBReplicaSvc")
+ replicaSvc.COOLSQLiteVetoPattern = ""
+ replicaSvc.UseCOOLSQLite = True
+ replicaSvc.UseCOOLFrontier = False
+ replicaSvc.UseGeomSQLite = True
+
+ # Timing
+ #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
+
+ # Dump config
+ # logging.getLogger('forcomps').setLevel(VERBOSE)
+ # acc.foreach_component("*").OutputLevel = VERBOSE
+ # acc.foreach_component("*ClassID*").OutputLevel = INFO
+ # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+ # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+ # acc.getService("StoreGateSvc").Dump = True
+ # acc.getService("ConditionStore").Dump = True
+ # acc.printConfig(withDetails=True)
+ # ConfigFlags.dump()
+
+ # Execute and finish
+ sc = acc.run(maxEvents=-1)
+
+ # Success should be 0
+ sys.exit(not sc.isSuccess())
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..bce33065b965395662707d54a82606d1e7038b9e
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
@@ -0,0 +1,106 @@
+"""
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+def OverlayRDOAlgCfg(flags, **kwargs):
+
+ from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+ acc = FaserGeometryCfg(flags)
+ OverlayRDOAlg = CompFactory.OverlayRDOAlg("OverlayRDOAlg",**kwargs)
+ acc.addEventAlgo(OverlayRDOAlg)
+
+ ItemList = []
+ # ItemList += ["xAOD::EventInfo#" + chargePrefix + "EventInfo"]
+ # ItemList += ["xAOD::EventAuxInfo#" + chargePrefix + "EventInfoAux."]
+ ItemList += ["xAOD::EventInfo#EventInfo"]
+ ItemList += ["xAOD::EventAuxInfo#EventInfoAux."]
+ ItemList += ["TrackCollection#Orig_CKFTrackCollectionWithoutIFT"]
+ ItemList += ["FaserSCT_RDO_Container#SCT_RDOs"]
+ ItemList += ["FaserSCT_RDO_Container#SCT_EDGEMODE_RDOs"]
+ # ItemList += ["Tracker::FaserSCT_ClusterContainer#" + chargePrefix + "SCT_ClusterContainer"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,"RDO", ItemList=ItemList, disableEventTag=True))
+
+
+
+ return acc
+
+if __name__ == "__main__":
+
+ import sys
+ from AthenaCommon.Logging import log, logging
+ from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+ from AthenaCommon.Configurable import Configurable
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ from AthenaConfiguration.TestDefaults import defaultTestFiles
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ # from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+ # Set up logging and new style config
+ log.setLevel(DEBUG)
+ Configurable.configurableRun3Behavior = True
+
+ # Configure
+ ConfigFlags.Overlay.DataOverlay = False
+ ConfigFlags.Input.Files = [ 'Pos_RDO.pool.root']
+ ConfigFlags.Input.SecondaryFiles = [ 'Neg_RDO.pool.root' ]
+ # ConfigFlags.Input.Files = [ '/eos/experiment/faser/rec/2022/r0013/009171/Faser-Physics-009171-00006-r0013-xAOD.root']
+ # ConfigFlags.Input.SecondaryFiles = [ '/eos/experiment/faser/rec/2022/r0013/009166/Faser-Physics-009166-00485-r0013-xAOD.root' ]
+ ConfigFlags.Output.RDOFileName = "Overlay.RDO.pool.root"
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03" # Always needed; must match FaserVersionS
+ ConfigFlags.IOVDb.DatabaseInstance = "CONDBR3" # Use data conditions for now
+ ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
+ ConfigFlags.Input.isMC = False # Needed to bypass autoconfig
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
+ ConfigFlags.Common.isOnline = False
+ ConfigFlags.GeoModel.Align.Dynamic = False
+ ConfigFlags.Beam.NumberOfCollisions = 0.
+
+ ConfigFlags.Detector.GeometryFaserSCT = True
+
+ ConfigFlags.lock()
+
+ # Core components
+ acc = MainServicesCfg(ConfigFlags)
+ acc.merge(PoolReadCfg(ConfigFlags))
+
+ # algorithm
+ acc.merge(OverlayRDOAlgCfg(ConfigFlags))
+
+ # from SGComps.AddressRemappingConfig import AddressRemappingCfg
+ # acc.merge(AddressRemappingCfg([
+ # "xAOD::EventInfo#EventInfo->" + ConfigFlags.Overlay.SigPrefix + "EventInfo",
+ # "xAOD::EventAuxInfo#EventInfoAux.->" + ConfigFlags.Overlay.SigPrefix + "EventInfoAux.",
+ # ]))
+
+ # Hack to avoid problem with our use of MC databases when isMC = False
+ replicaSvc = acc.getService("DBReplicaSvc")
+ replicaSvc.COOLSQLiteVetoPattern = ""
+ replicaSvc.UseCOOLSQLite = True
+ replicaSvc.UseCOOLFrontier = False
+ replicaSvc.UseGeomSQLite = True
+
+ # Timing
+ #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
+
+ # Dump config
+ # logging.getLogger('forcomps').setLevel(VERBOSE)
+ # acc.foreach_component("*").OutputLevel = VERBOSE
+ # acc.foreach_component("*ClassID*").OutputLevel = INFO
+ # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+ # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+ # acc.getService("StoreGateSvc").Dump = True
+ # acc.getService("ConditionStore").Dump = True
+ # acc.printConfig(withDetails=True)
+ # ConfigFlags.dump()
+
+ # Execute and finish
+ sc = acc.run(maxEvents=-1)
+
+ # Success should be 0
+ sys.exit(not sc.isSuccess())
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..29534a9b3293511b430305770da134c073b10213
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
@@ -0,0 +1,156 @@
+"""
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+chargePrefix = "Pos_"
+
+def SelectRDOAlgCfg(flags, **kwargs):
+
+ from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+ acc = FaserGeometryCfg(flags)
+ SelectRDOAlg = CompFactory.SelectRDOAlg("SelectRDOAlg",**kwargs)
+ SelectRDOAlg.AcceptPositive = ("Pos_" in chargePrefix)
+ acc.addEventAlgo(SelectRDOAlg)
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin0RDO",Xmin=-66, Xmax=0, Ymin=-66, Ymax=0))
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin1RDO",Xmin=0, Xmax=66, Ymin=-66, Ymax=0))
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin2RDO",Xmin=-66, Xmax=0, Ymin=0, Ymax=66))
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin3RDO",Xmin=0, Xmax=66, Ymin=0, Ymax=66))
+
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin4RDO",Xmin=-33, Xmax=33, Ymin=-33, Ymax=33))
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin5RDO",Xmin=-33, Xmax=33, Ymin=33, Ymax=99))
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin6RDO",Xmin=-33, Xmax=33, Ymin=-99, Ymax=-33))
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin7RDO",Xmin=-99, Xmax=-33, Ymin=-33, Ymax=33))
+ acc.addEventAlgo(CompFactory.BinRDOAlg("Bin8RDO",Xmin=33, Xmax=99, Ymin=-33, Ymax=33))
+
+ ItemList = []
+ # ItemList += ["xAOD::EventInfo#" + chargePrefix + "EventInfo"]
+ # ItemList += ["xAOD::EventAuxInfo#" + chargePrefix + "EventInfoAux."]
+ ItemList += ["xAOD::EventInfo#EventInfo"]
+ ItemList += ["xAOD::EventAuxInfo#EventInfoAux."]
+ ItemList += ["TrackCollection#" + chargePrefix + "CKFTrackCollectionWithoutIFT"]
+ ItemList += ["FaserSCT_RDO_Container#" + chargePrefix + "SCT_RDOs"]
+ ItemList += ["FaserSCT_RDO_Container#" + chargePrefix + "SCT_EDGEMODE_RDOs"]
+ # ItemList += ["Tracker::FaserSCT_ClusterContainer#" + chargePrefix + "SCT_ClusterContainer"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin0RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo0 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin0RDO")
+ osrdo0.RequireAlgs += ["SelectRDOAlg","Bin0RDO"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin1RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo1 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin1RDO")
+ osrdo1.RequireAlgs += ["SelectRDOAlg","Bin1RDO"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin2RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo2 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin2RDO")
+ osrdo2.RequireAlgs += ["SelectRDOAlg","Bin2RDO"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin3RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo3 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin3RDO")
+ osrdo3.RequireAlgs += ["SelectRDOAlg","Bin3RDO"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin4RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo4 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin4RDO")
+ osrdo4.RequireAlgs += ["SelectRDOAlg","Bin4RDO"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin5RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo5 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin5RDO")
+ osrdo5.RequireAlgs += ["SelectRDOAlg","Bin5RDO"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin6RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo6 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin6RDO")
+ osrdo6.RequireAlgs += ["SelectRDOAlg","Bin6RDO"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin7RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo7 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin7RDO")
+ osrdo7.RequireAlgs += ["SelectRDOAlg","Bin7RDO"]
+
+ acc.merge(OutputStreamCfg(ConfigFlags,chargePrefix+"Bin8RDO", ItemList=ItemList, disableEventTag=True))
+ osrdo8 = acc.getEventAlgo("OutputStream" + chargePrefix + "Bin8RDO")
+ osrdo8.RequireAlgs += ["SelectRDOAlg","Bin8RDO"]
+
+ return acc
+
+
+
+if __name__ == "__main__":
+
+ import sys
+ from AthenaCommon.Logging import log, logging
+ from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+ from AthenaCommon.Configurable import Configurable
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ from AthenaConfiguration.TestDefaults import defaultTestFiles
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+
+ # Set up logging and new style config
+ log.setLevel(DEBUG)
+ Configurable.configurableRun3Behavior = True
+
+ # Configure
+ ConfigFlags.Input.Files = [ #'/eos/experiment/faser/rec/2022/r0013/009171/Faser-Physics-009171-00006-r0013-xAOD.root'
+ '/eos/experiment/faser/rec/2022/r0013/009166/Faser-Physics-009166-00485-r0013-xAOD.root'
+ ]
+ # ConfigFlags.Input.Files = [ 'Faser-Physics-009171-00006-r0013-xAOD.root',
+ # 'Faser-Physics-009166-00485-r0013-xAOD.root']
+ # ConfigFlags.Output.RDOFileName = chargePrefix + "RDO.pool.root"
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03" # Always needed; must match FaserVersionS
+ ConfigFlags.IOVDb.DatabaseInstance = "CONDBR3" # Use data conditions for now
+ ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
+ ConfigFlags.Input.isMC = False # Needed to bypass autoconfig
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
+ ConfigFlags.Common.isOnline = False
+ ConfigFlags.GeoModel.Align.Dynamic = False
+ ConfigFlags.Beam.NumberOfCollisions = 0.
+
+ ConfigFlags.Detector.GeometryFaserSCT = True
+
+ ConfigFlags.lock()
+
+ # Core components
+ acc = MainServicesCfg(ConfigFlags)
+ acc.merge(PoolReadCfg(ConfigFlags))
+
+ # algorithm
+ acc.merge(SelectRDOAlgCfg(ConfigFlags, TrackCollection = "CKFTrackCollectionWithoutIFT", OutputTrackCollection = chargePrefix + "CKFTrackCollectionWithoutIFT"))
+
+ from SGComps.AddressRemappingConfig import AddressRemappingCfg
+ acc.merge(AddressRemappingCfg([
+ # "xAOD::EventInfo#EventInfo->" + chargePrefix + "EventInfo",
+ # "xAOD::EventAuxInfo#EventInfoAux.->" + chargePrefix + "EventInfoAux.",
+ # "TrackCollection#CKFTrackCollectionWithoutIFT->" + chargePrefix + "CKFTrackCollectionWithoutIFT",
+ "FaserSCT_RDO_Container#SCT_RDOs->" + chargePrefix + "SCT_RDOs",
+ "FaserSCT_RDO_Container#SCT_EDGEMODE_RDOs->" + chargePrefix + "SCT_EDGEMODE_RDOs",
+ # "Tracker::FaserSCT_ClusterContainer#SCT_ClusterContainer->" + chargePrefix + "SCT_ClusterContainer"
+ ]))
+
+ # Hack to avoid problem with our use of MC databases when isMC = False
+ replicaSvc = acc.getService("DBReplicaSvc")
+ replicaSvc.COOLSQLiteVetoPattern = ""
+ replicaSvc.UseCOOLSQLite = True
+ replicaSvc.UseCOOLFrontier = False
+ replicaSvc.UseGeomSQLite = True
+
+ # Timing
+ #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
+
+ # Dump config
+ # logging.getLogger('forcomps').setLevel(VERBOSE)
+ # acc.foreach_component("*").OutputLevel = VERBOSE
+ # acc.foreach_component("*ClassID*").OutputLevel = INFO
+ # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+ # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+ # acc.getService("StoreGateSvc").Dump = True
+ # acc.getService("ConditionStore").Dump = True
+ # acc.printConfig(withDetails=True)
+ # ConfigFlags.dump()
+
+ # Execute and finish
+ sc = acc.run(maxEvents=-1)
+
+ # Success should be 0
+ sys.exit(not sc.isSuccess())
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a917447dde480b9b82a316653fb8e76cf76a2156
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.cxx
@@ -0,0 +1,55 @@
+#include "BinRDOAlg.h"
+#include "TrkTrack/Track.h"
+
+BinRDOAlg::BinRDOAlg(const std::string &name,
+ ISvcLocator *pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator)
+{}
+
+
+StatusCode BinRDOAlg::initialize()
+{
+ ATH_CHECK(m_trackCollection.initialize());
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode BinRDOAlg::execute(const EventContext &ctx) const
+{
+
+ setFilterPassed(false, ctx);
+ SG::ReadHandle<TrackCollection> trackCollection {m_trackCollection, ctx};
+ if (!trackCollection.isValid() || trackCollection->size() == 0) return StatusCode::SUCCESS;
+
+ // track cuts
+
+ double xUpstream {0.0};
+ double yUpstream {0.0};
+
+ for (const Trk::Track* track : *trackCollection)
+ {
+ if (track == nullptr) continue;
+ const Trk::TrackParameters* upstreamParameters {nullptr};
+ for (auto params : *(track->trackParameters()))
+ {
+ if (params->position().z() < 0) continue; // Ignore IFT hits
+ if (upstreamParameters == nullptr || params->position().z() < upstreamParameters->position().z()) upstreamParameters = params;
+ }
+ xUpstream = upstreamParameters->position().x();
+ yUpstream = upstreamParameters->position().y();
+ break;
+ }
+
+
+ // set filter passed
+
+ setFilterPassed(((xUpstream >= m_xMin) && (xUpstream <= m_xMax) && (yUpstream >= m_yMin) && (yUpstream <= m_yMax)), ctx);
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode BinRDOAlg::finalize()
+{
+ return StatusCode::SUCCESS;
+}
+
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..f8e86e844adf11e7a1a4021b9e8d7a11757ce58c
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.h
@@ -0,0 +1,23 @@
+#pragma once
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "TrkTrack/TrackCollection.h"
+
+class BinRDOAlg : public AthReentrantAlgorithm {
+public:
+ BinRDOAlg(const std::string &name, ISvcLocator *pSvcLocator);
+ virtual ~BinRDOAlg() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute(const EventContext &ctx) const override;
+ virtual StatusCode finalize() override;
+
+private:
+
+ SG::ReadHandleKey<TrackCollection> m_trackCollection { this, "TrackCollection", "CKFTrackCollectionWithoutIFT", "Input track collection name" };
+ DoubleProperty m_xMin {this, "Xmin", -25.0, "Minimum x position to accept"};
+ DoubleProperty m_xMax {this, "Xmax", 25.0, "Maximum x position to accept"};
+ DoubleProperty m_yMin {this, "Ymin", -25.0, "Minimum y position to accept"};
+ DoubleProperty m_yMax {this, "Ymax", 25.0, "Maximum y position to accept"};
+
+};
+
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..78096c60004a9e183bf83aa235ef72501ddee4c4
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx
@@ -0,0 +1,247 @@
+#include "HistoRDOAlg.h"
+#include "TrkTrack/Track.h"
+#include <TTree.h>
+
+
+HistoRDOAlg::HistoRDOAlg(const std::string &name,
+ ISvcLocator *pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator),
+ AthHistogramming(name),
+ m_histSvc("THistSvc/THistSvc", name) {}
+
+StatusCode HistoRDOAlg::initialize()
+{
+ m_tree = new TTree("overlay", "Overlay tree");
+
+ m_tree->Branch("longTracks1", &m_longTracks1, "longTracks1/I");
+ m_tree->Branch("x1", &m_x1);
+ m_tree->Branch("y1", &m_y1);
+ m_tree->Branch("p1", &m_p1);
+ m_tree->Branch("theta1", &m_theta1);
+ m_tree->Branch("phi1", &m_phi1);
+ m_tree->Branch("charge1", &m_charge1);
+ m_tree->Branch("chiSquared1", &m_chiSquared1);
+ m_tree->Branch("nDoF1", &m_nDoF1);
+ m_tree->Branch("pPolar1", &m_pTotPolar1, "pPolar1/D");
+
+ m_tree->Branch("longTracks2", &m_longTracks2, "longTracks2/I");
+ m_tree->Branch("x2", &m_x2);
+ m_tree->Branch("y2", &m_y2);
+ m_tree->Branch("p2", &m_p2);
+ m_tree->Branch("theta2", &m_theta2);
+ m_tree->Branch("phi2", &m_phi2);
+ m_tree->Branch("charge2", &m_charge2);
+ m_tree->Branch("chiSquared2", &m_chiSquared2);
+ m_tree->Branch("nDoF2", &m_nDoF2);
+ m_tree->Branch("pPolar2", &m_pTotPolar2, "pPolar2/D");
+
+ m_tree->Branch("yAgreement", &m_yAgreement, "yAgreement/D");
+ m_tree->Branch("matchFraction", &m_matchFraction, "matchFraction/D");
+
+ ATH_CHECK(histSvc()->regTree("/HIST2/tree", m_tree));
+
+ ATH_CHECK(m_trackCollection1.initialize());
+ ATH_CHECK(m_trackCollection2.initialize());
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode HistoRDOAlg::execute(const EventContext& ctx) const
+{
+
+ m_longTracks1 = 0;
+ m_x1.clear();
+ m_y1.clear();
+ m_p1.clear();
+ m_theta1.clear();
+ m_phi1.clear();
+ m_charge1.clear();
+ m_chiSquared1.clear();
+ m_nDoF1.clear();
+ m_pTotPolar1 = 0;
+
+ m_longTracks2 = 0;
+ m_x2.clear();
+ m_y2.clear();
+ m_p2.clear();
+ m_theta2.clear();
+ m_phi2.clear();
+ m_charge2.clear();
+ m_chiSquared2.clear();
+ m_nDoF2.clear();
+ m_pTotPolar2 = 0;
+
+ m_yAgreement = 0;
+ m_matchFraction = 0;
+
+ SG::ReadHandle<TrackCollection> trackCollection1 {m_trackCollection1, ctx};
+ // ATH_CHECK(trackCollection1.isValid());
+ SG::ReadHandle<TrackCollection> trackCollection2 {m_trackCollection2, ctx};
+ // ATH_CHECK(trackCollection2.isValid());
+
+ std::vector<const Trk::TrackParameters*> upstreamSingle {};
+ if (trackCollection1.isValid())
+ {
+ Amg::Vector3D pTot {0, 0, 0};
+ for (const Trk::Track* track : *trackCollection1)
+ {
+ if (track == nullptr) continue;
+ m_longTracks1++;
+ const Trk::TrackParameters* upstreamParameters {nullptr};
+ for (auto params : *(track->trackParameters()))
+ {
+ if (params->position().z() < 0) continue; // Ignore IFT hits
+ if (upstreamParameters == nullptr || params->position().z() < upstreamParameters->position().z()) upstreamParameters = params;
+ ATH_MSG_VERBOSE("TrackCollection 1, Track " << m_longTracks1 << ", z = " << params->position().z() << ", y = " << params->position().y() << ", q/p(TeV) = " << params->charge()/(params->momentum().mag()/1e6));
+ }
+ upstreamSingle.push_back(upstreamParameters);
+ m_x1.push_back(upstreamParameters->position().x());
+ m_y1.push_back(upstreamParameters->position().y());
+ m_p1.push_back(upstreamParameters->momentum().mag());
+ m_charge1.push_back(upstreamParameters->charge());
+ m_theta1.push_back(asin(upstreamParameters->momentum().perp()/upstreamParameters->momentum().mag()));
+ m_phi1.push_back(atan2(upstreamParameters->momentum().y(), upstreamParameters->momentum().x()));
+ m_chiSquared1.push_back(track->fitQuality()->chiSquared());
+ m_nDoF1.push_back(track->fitQuality()->doubleNumberDoF());
+ pTot += upstreamParameters->momentum();
+ }
+ m_pTotPolar1 = asin(pTot.perp()/pTot.mag());
+ }
+
+ std::vector<const Trk::TrackParameters*> upstreamOverlay {};
+ if (trackCollection2.isValid())
+ {
+ Amg::Vector3D pTot {0, 0, 0};
+ for (const Trk::Track* track : *trackCollection2)
+ {
+ if (track == nullptr) continue;
+ m_longTracks2++;
+ const Trk::TrackParameters* upstreamParameters {nullptr};
+ for (auto params : *(track->trackParameters()))
+ {
+ if (params->position().z() < 0) continue; // Ignore IFT hits
+ if (upstreamParameters == nullptr || params->position().z() < upstreamParameters->position().z()) upstreamParameters = params;
+ ATH_MSG_VERBOSE("TrackCollection 2, Track " << m_longTracks2 << ", z = " << params->position().z()<< ", y = " << params->position().y()<< ", q/p(TeV) = " << params->charge()/(params->momentum().mag()/1e6));
+ }
+ upstreamOverlay.push_back(upstreamParameters);
+ m_x2.push_back(upstreamParameters->position().x());
+ m_y2.push_back(upstreamParameters->position().y());
+ m_p2.push_back(upstreamParameters->momentum().mag());
+ m_charge2.push_back(upstreamParameters->charge());
+ m_theta2.push_back(asin(upstreamParameters->momentum().perp()/upstreamParameters->momentum().mag()));
+ m_phi2.push_back(atan2(upstreamParameters->momentum().y(), upstreamParameters->momentum().x()));
+ m_chiSquared2.push_back(track->fitQuality()->chiSquared());
+ m_nDoF2.push_back(track->fitQuality()->doubleNumberDoF());
+ pTot += upstreamParameters->momentum();
+ }
+ m_pTotPolar2 = asin(pTot.perp()/pTot.mag());
+ }
+
+ if (trackCollection1.isValid() && trackCollection2.isValid() && m_longTracks2 ==2)
+ {
+ double qOverPSingle1 = upstreamSingle[0]->charge()/upstreamSingle[0]->momentum().mag();
+ double qOverPCovSingle1 = (*upstreamSingle[0]->covariance())(4,4);
+ double qOverPSingle2 = upstreamSingle[1]->charge()/upstreamSingle[1]->momentum().mag();
+ double qOverPCovSingle2 = (*upstreamSingle[1]->covariance())(4,4);
+ double qOverPOverlay1 = upstreamOverlay[0]->charge()/upstreamOverlay[0]->momentum().mag();
+ double qOverPCovOverlay1 = (*upstreamOverlay[0]->covariance())(4,4);
+ double qOverPOverlay2 = upstreamOverlay[1]->charge()/upstreamOverlay[1]->momentum().mag();
+ double qOverPCovOverlay2 = (*upstreamOverlay[1]->covariance())(4,4);
+
+ double agreementDirect = pow(qOverPSingle1 - qOverPOverlay1, 2) / (qOverPCovSingle1 + qOverPCovOverlay1) +
+ pow(qOverPSingle2 - qOverPOverlay2, 2) / (qOverPCovSingle2 + qOverPCovOverlay2);
+ double agreementSwap = pow(qOverPSingle1 - qOverPOverlay2, 2) / (qOverPCovSingle1 + qOverPCovOverlay2) +
+ pow(qOverPSingle2 - qOverPOverlay1, 2) / (qOverPCovSingle2 + qOverPCovOverlay1);
+
+ const Trk::Track* firstSingle {nullptr};
+ const Trk::Track* firstOverlay {nullptr};
+ const Trk::Track* secondSingle {nullptr};
+ const Trk::Track* secondOverlay {nullptr};
+ if (agreementDirect <= agreementSwap)
+ {
+ firstSingle = (*trackCollection1)[0];
+ firstOverlay = (*trackCollection2)[0];
+ secondSingle = (*trackCollection1)[1];
+ secondOverlay = (*trackCollection2)[1];
+ ATH_MSG_VERBOSE("Matched q/p = " << qOverPSingle1 << " with " << qOverPOverlay1 << ", and q/p = " << qOverPSingle2 << " with " << qOverPOverlay2);
+ }
+ else
+ {
+ firstSingle = (*trackCollection1)[0];
+ firstOverlay = (*trackCollection2)[1];
+ secondSingle = (*trackCollection1)[1];
+ secondOverlay = (*trackCollection2)[0];
+ ATH_MSG_VERBOSE("Matched q/p = " << qOverPSingle1 << " with " << qOverPOverlay2 << ", and q/p = " << qOverPSingle2 << " with " << qOverPOverlay1);
+ }
+ int nZMatch {0};
+ auto singleParameters = firstSingle->trackParameters();
+ auto overlayParameters = firstOverlay->trackParameters();
+ for (size_t iSingle = 0, iOverlay = 0; (iSingle < singleParameters->size()) && (iOverlay < overlayParameters->size());)
+ {
+ auto singleState = (*singleParameters)[iSingle];
+ auto overlayState = (*overlayParameters)[iOverlay];
+ auto singleCov = singleState->covariance();
+ auto overlayCov = overlayState->covariance();
+ if (abs(singleState->position().z() - overlayState->position().z()) < 0.03)
+ {
+ nZMatch++;
+ m_yAgreement += pow(singleState->position().y() - overlayState->position().y(), 2)/((*singleCov)(1,1) + (*overlayCov)(1,1));
+ iSingle++;
+ iOverlay++;
+ ATH_MSG_VERBOSE("Matched (1): " << singleState->position().z() << " and " << overlayState->position().z());
+ }
+ else if (singleState->position().z() < overlayState->position().z())
+ {
+ ATH_MSG_VERBOSE("z not matched (1): " << singleState->position().z() << " vs " << overlayState->position().z());
+ iSingle++;
+ }
+ else if (overlayState->position().z() < singleState->position().z())
+ {
+ ATH_MSG_VERBOSE("z not matched (1): " << singleState->position().z() << " vs " << overlayState->position().z());
+ iOverlay++;
+ }
+ }
+ double matchFract1 { ((double) nZMatch)/singleParameters->size() };
+ ATH_MSG_VERBOSE("MatchFract1: " << matchFract1);
+ nZMatch = 0;
+ singleParameters = secondSingle->trackParameters();
+ overlayParameters = secondOverlay->trackParameters();
+ for (size_t iSingle = 0, iOverlay = 0; (iSingle < singleParameters->size()) && (iOverlay < overlayParameters->size());)
+ {
+ auto singleState = (*singleParameters)[iSingle];
+ auto overlayState = (*overlayParameters)[iOverlay];
+ auto singleCov = singleState->covariance();
+ auto overlayCov = overlayState->covariance();
+ if (abs(singleState->position().z() - overlayState->position().z()) < 0.03)
+ {
+ nZMatch++;
+ m_yAgreement += pow(singleState->position().y() - overlayState->position().y(), 2)/((*singleCov)(1,1) + (*overlayCov)(1,1));
+ iSingle++;
+ iOverlay++;
+ ATH_MSG_VERBOSE("Matched (2): " << singleState->position().z() << " and " << overlayState->position().z());
+ }
+ else if (singleState->position().z() < overlayState->position().z())
+ {
+ ATH_MSG_VERBOSE("z not matched (2): " << singleState->position().z() << " vs " << overlayState->position().z());
+ iSingle++;
+ }
+ else if (overlayState->position().z() < singleState->position().z())
+ {
+ ATH_MSG_VERBOSE("z not matched (2): " << singleState->position().z() << " vs " << overlayState->position().z());
+ iOverlay++;
+ }
+ }
+ double matchFract2 { ((double) nZMatch)/singleParameters->size() };
+ ATH_MSG_VERBOSE("MatchFract2: " << matchFract2);
+ m_matchFraction = std::min(matchFract1, matchFract2);
+ ATH_MSG_VERBOSE("MatchFraction: " << m_matchFraction);
+ }
+
+ m_tree->Fill();
+ return StatusCode::SUCCESS;
+}
+
+StatusCode HistoRDOAlg::finalize()
+{
+ return StatusCode::SUCCESS;
+}
\ No newline at end of file
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..4fe18f7ad71e075b8b37b56518058da6b7be6f0e
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.h
@@ -0,0 +1,58 @@
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "AthenaBaseComps/AthHistogramming.h"
+#include "TrkTrack/TrackCollection.h"
+
+#include <vector>
+
+class TTree;
+
+class HistoRDOAlg : public AthReentrantAlgorithm, AthHistogramming {
+public:
+ HistoRDOAlg(const std::string &name, ISvcLocator *pSvcLocator);
+ virtual ~HistoRDOAlg() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute(const EventContext &ctx) const override;
+ virtual StatusCode finalize() override;
+ const ServiceHandle <ITHistSvc> &histSvc() const;
+
+private:
+
+ ServiceHandle <ITHistSvc> m_histSvc;
+
+ SG::ReadHandleKey<TrackCollection> m_trackCollection1 { this, "TrackCollection1", "Pos_CKFTrackCollectionWithoutIFT", "Input track collection 1 name" };
+ SG::ReadHandleKey<TrackCollection> m_trackCollection2 { this, "TrackCollection2", "Neg_CKFTrackCollectionWithoutIFT", "Input track collection 2 name" };
+
+// mutable TTree* m_tree;
+
+ mutable TTree* m_tree;
+
+ mutable int m_longTracks1;
+ mutable std::vector<double> m_p1;
+ mutable std::vector<double> m_theta1;
+ mutable std::vector<double> m_phi1;
+ mutable std::vector<double> m_x1;
+ mutable std::vector<double> m_y1;
+ mutable std::vector<double> m_charge1;
+ mutable std::vector<double> m_chiSquared1;
+ mutable std::vector<double> m_nDoF1;
+ mutable double m_pTotPolar1;
+
+ mutable int m_longTracks2;
+ mutable std::vector<double> m_p2;
+ mutable std::vector<double> m_theta2;
+ mutable std::vector<double> m_phi2;
+ mutable std::vector<double> m_x2;
+ mutable std::vector<double> m_y2;
+ mutable std::vector<double> m_charge2;
+ mutable std::vector<double> m_chiSquared2;
+ mutable std::vector<double> m_nDoF2;
+ mutable double m_pTotPolar2;
+
+ mutable double m_yAgreement;
+ mutable double m_matchFraction;
+
+};
+
+inline const ServiceHandle <ITHistSvc> &HistoRDOAlg::histSvc() const {
+ return m_histSvc;
+}
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..9b88a96f9590462d5c03d6ffba601598185ec109
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
@@ -0,0 +1,186 @@
+#include "OverlayRDOAlg.h"
+#include "TrkTrack/Track.h"
+#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include "Identifier/Identifier.h"
+#include "TrackerReadoutGeometry/SCT_DetectorManager.h"
+#include "TrackerReadoutGeometry/SiDetectorElement.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include <cmath>
+
+
+
+OverlayRDOAlg::OverlayRDOAlg(const std::string &name,
+ ISvcLocator *pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator)
+ {}
+
+StatusCode OverlayRDOAlg::initialize()
+{
+ ATH_CHECK(detStore()->retrieve(m_sctHelper, "FaserSCT_ID"));
+ ATH_CHECK(m_posTrackKey.initialize());
+ ATH_CHECK(m_negTrackKey.initialize());
+ ATH_CHECK(m_posRdoKey.initialize());
+ ATH_CHECK(m_negRdoKey.initialize());
+ ATH_CHECK(m_posEdgeModeRdoKey.initialize());
+ ATH_CHECK(m_negEdgeModeRdoKey.initialize());
+ ATH_CHECK(m_outRdoKey.initialize());
+ ATH_CHECK(m_outEdgeModeRdoKey.initialize());
+ ATH_CHECK(m_outputTrackCollection.initialize());
+
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode OverlayRDOAlg::execute(const EventContext &ctx) const
+{
+
+ SG::ReadHandle<TrackCollection> posTrackCollection {m_posTrackKey, ctx};
+ ATH_CHECK(posTrackCollection.isValid());
+ SG::ReadHandle<TrackCollection> negTrackCollection {m_negTrackKey, ctx};
+ ATH_CHECK(negTrackCollection.isValid());
+ SG::ReadHandle<FaserSCT_RDO_Container> posRdoContainer {m_posRdoKey, ctx};
+ ATH_CHECK(posRdoContainer.isValid());
+ SG::ReadHandle<FaserSCT_RDO_Container> negRdoContainer {m_negRdoKey, ctx};
+ ATH_CHECK(negRdoContainer.isValid());
+ SG::ReadHandle<FaserSCT_RDO_Container> posEdgeModeRdoContainer {m_posEdgeModeRdoKey, ctx};
+ ATH_CHECK(posEdgeModeRdoContainer.isValid());
+ SG::ReadHandle<FaserSCT_RDO_Container> negEdgeModeRdoContainer {m_negEdgeModeRdoKey, ctx};
+ ATH_CHECK(negEdgeModeRdoContainer.isValid());
+
+ // SCT_RDOs
+
+ std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> > rdoDB {};
+ processContainer(*posRdoContainer, rdoDB);
+ processContainer(*negRdoContainer, rdoDB);
+
+ SG::WriteHandle<FaserSCT_RDO_Container> outRdoContainer {m_outRdoKey, ctx};
+ ATH_CHECK(outRdoContainer.record(std::make_unique<FaserSCT_RDO_Container>(m_sctHelper->wafer_hash_max())));
+
+ processDB(rdoDB, outRdoContainer, [](int i) {return true;});
+
+ // SCT_EDGEMODE_RDOs
+
+ rdoDB.clear();
+ processContainer(*posEdgeModeRdoContainer, rdoDB);
+ processContainer(*negEdgeModeRdoContainer, rdoDB);
+
+ SG::WriteHandle<FaserSCT_RDO_Container> outEdgeModeRdoContainer {m_outEdgeModeRdoKey, ctx};
+ ATH_CHECK(outEdgeModeRdoContainer.record(std::make_unique<FaserSCT_RDO_Container>(m_sctHelper->wafer_hash_max())));
+
+ processDB(rdoDB, outEdgeModeRdoContainer, [](int hitPattern) { return (((hitPattern & 0x2) == 0 ) || ((hitPattern & 0x4) != 0) ) ? false : true;});
+
+ // Tracks
+
+ SG::WriteHandle<TrackCollection> outputTrackCollection {m_outputTrackCollection, ctx};
+ std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
+
+ for (auto theTrack : *posTrackCollection)
+ {
+ outputTracks->push_back(cloneTrack(theTrack) );
+ }
+
+ for (auto theTrack : *negTrackCollection)
+ {
+ outputTracks->push_back(cloneTrack(theTrack) );
+ }
+
+ ATH_CHECK(outputTrackCollection.record(std::move(outputTracks)));
+
+ return StatusCode::SUCCESS;
+}
+
+
+void
+OverlayRDOAlg::processDB(const std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> >& rdoDB, SG::WriteHandle<FaserSCT_RDO_Container>& outRdoContainer, std::function<bool(int)> lambda) const
+{
+ for (auto& entry : rdoDB)
+ {
+ IdentifierHash waferHash = entry.first;
+ Identifier id = m_sctHelper->wafer_id(waferHash);
+ auto vec = entry.second;
+ if (vec.size() == 0)
+ {
+ std::cout << "Unexpected zero-length collection vector for ID = " << id << " ..." << std::endl;
+ continue;
+ }
+ std::unique_ptr<FaserSCT_RDO_Collection> current_collection = std::make_unique<FaserSCT_RDO_Collection>(waferHash);
+ current_collection->setIdentifier(id);
+ std::map<int, unsigned int> stripMap;
+ for (auto& collection : vec)
+ {
+ for (auto rawData : *collection)
+ {
+ // const FaserSCT3_RawData* data = dynamic_cast<const FaserSCT3_RawData*>(rawData);
+ Identifier stripID = rawData->identify();
+ int stripNumber = m_sctHelper->strip(stripID);
+ int groupSize = rawData->getGroupSize();
+ unsigned int dataWord = rawData->getWord();
+ // int time = data->getTimeBin();
+ // std::cout << "Word: " << std::hex << dataWord << std::dec << " Time: " << time << std::endl;
+ for (int i = stripNumber; i < stripNumber + groupSize; i++)
+ stripMap[i] |= dataWord;
+ }
+ }
+
+ for (auto stripEntry : stripMap)
+ {
+ Identifier rdoID {m_sctHelper->strip_id(id, stripEntry.first)};
+ if (!lambda(stripEntry.second>>22)) continue;
+ current_collection->emplace_back(new FaserSCT3_RawData(rdoID, stripEntry.second, std::vector<int> {}));
+ }
+
+ outRdoContainer->getWriteHandle(waferHash).addOrDelete(std::move(current_collection));
+ }
+}
+
+
+Trk::Track* OverlayRDOAlg::cloneTrack(const Trk::Track* theTrack) const
+{
+ Trk::TrackInfo i { theTrack->info() };
+ Trk::FitQuality* q = new Trk::FitQuality { *(theTrack->fitQuality()) };
+ DataVector<const Trk::TrackStateOnSurface>* s = new DataVector<const Trk::TrackStateOnSurface> {};
+
+ for (auto oldState : (*theTrack->trackStateOnSurfaces()))
+ {
+ const Trk::CurvilinearParameters* oldParam = dynamic_cast<const Trk::CurvilinearParameters*>(oldState->trackParameters());
+ std::unique_ptr<Trk::TrackParameters> newParam { new Trk::CurvilinearParameters { *oldParam } };
+ Trk::TrackStateOnSurface* newState = new Trk::TrackStateOnSurface { nullptr, newParam.release() };
+ s->push_back(newState);
+ }
+ std::unique_ptr<DataVector<const Trk::TrackStateOnSurface>> sink {s};
+ return new Trk::Track {i, std::move(*sink) , q } ;
+}
+
+void OverlayRDOAlg::processContainer(const FaserSCT_RDO_Container& container, std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> >& rdoDB) const
+{
+ size_t nLong{0};
+ size_t nData{0};
+ for( const auto& collection : container)
+ {
+ if (collection->empty()) continue;
+ Identifier id = collection->identify();
+ IdentifierHash hash = m_sctHelper->wafer_hash(id);
+ if (rdoDB.count(hash))
+ {
+ rdoDB[hash].push_back(collection);
+ }
+ else
+ {
+ rdoDB[hash] = std::vector<const FaserSCT_RDO_Collection*>{collection};
+ }
+ for(const auto& rawdata : *collection)
+ {
+ nData++;
+ if (rawdata->getGroupSize() > 1) nLong++;
+ }
+ }
+ std::cout << nLong << " long data out of " << nData << " total RDOs" << std::endl;
+}
+
+
+StatusCode OverlayRDOAlg::finalize()
+{
+ return StatusCode::SUCCESS;
+}
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..eba2f90d7cde46873bbe9124e17876fef5059e5a
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
@@ -0,0 +1,38 @@
+#pragma once
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "TrkTrack/TrackCollection.h"
+#include "TrackerRawData/FaserSCT_RDO_Container.h"
+
+class TTree;
+class FaserSCT_ID;
+
+class OverlayRDOAlg : public AthReentrantAlgorithm {
+public:
+ OverlayRDOAlg(const std::string &name, ISvcLocator *pSvcLocator);
+ virtual ~OverlayRDOAlg() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute(const EventContext &ctx) const override;
+ virtual StatusCode finalize() override;
+
+private:
+
+ void processContainer(const FaserSCT_RDO_Container& container, std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> >& rdoDB) const;
+ void processDB(const std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> >& db, SG::WriteHandle<FaserSCT_RDO_Container>& outContainer, std::function<bool(int)> lambda) const;
+ Trk::Track* cloneTrack(const Trk::Track* originalTrack) const;
+
+ SG::ReadHandleKey<TrackCollection> m_posTrackKey { this, "PosTrackCollection", "Pos_CKFTrackCollectionWithoutIFT", "Input track collection name" };
+ SG::ReadHandleKey<TrackCollection> m_negTrackKey { this, "NegTrackCollection", "Neg_CKFTrackCollectionWithoutIFT", "Input track collection name" };
+ SG::ReadHandleKey<FaserSCT_RDO_Container> m_posRdoKey { this, "PosRdoContainer", "Pos_SCT_RDOs"};
+ SG::ReadHandleKey<FaserSCT_RDO_Container> m_negRdoKey { this, "NegRdoContainer", "Neg_SCT_RDOs"};
+ SG::ReadHandleKey<FaserSCT_RDO_Container> m_posEdgeModeRdoKey { this, "PosEdgeModeRdoContainer", "Pos_SCT_EDGEMODE_RDOs"};
+ SG::ReadHandleKey<FaserSCT_RDO_Container> m_negEdgeModeRdoKey { this, "NegEdgeModeRdoContainer", "Neg_SCT_EDGEMODE_RDOs"};
+
+ SG::WriteHandleKey<FaserSCT_RDO_Container> m_outRdoKey{this, "OutputRDOObjectName", "SCT_RDOs", "Output RDO Object name"};
+ SG::WriteHandleKey<FaserSCT_RDO_Container> m_outEdgeModeRdoKey{this, "OutputEDGEMODEObjectName", "SCT_EDGEMODE_RDOs", "Output EDGEMODE Object name"};
+ SG::WriteHandleKey<TrackCollection> m_outputTrackCollection { this, "OutputTrackCollection", "Orig_CKFTrackCollectionWithoutIFT", "Output track collection name"};
+
+ const FaserSCT_ID* m_sctHelper;
+
+};
+
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b33bed9d18ac4d51f2e2dc7ae9fc52217cbffc67
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
@@ -0,0 +1,187 @@
+#include "SelectRDOAlg.h"
+#include "TrkTrack/Track.h"
+#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterContainer.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterCollection.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include "Identifier/Identifier.h"
+#include <cmath>
+
+SelectRDOAlg::SelectRDOAlg(const std::string &name,
+ ISvcLocator *pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator)
+{}
+
+
+StatusCode SelectRDOAlg::initialize()
+{
+ ATH_CHECK(m_trackCollection.initialize());
+ ATH_CHECK(m_segmentCollection.initialize());
+ ATH_CHECK(m_clusterContainer.initialize());
+ ATH_CHECK(m_triggerContainer.initialize());
+ ATH_CHECK(m_outputTrackCollection.initialize());
+
+ ATH_CHECK(detStore()->retrieve(m_sctHelper, "FaserSCT_ID"));
+
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode SelectRDOAlg::execute(const EventContext &ctx) const
+{
+
+ setFilterPassed(false, ctx);
+ SG::ReadHandle<TrackCollection> trackCollection {m_trackCollection, ctx};
+ if (!trackCollection.isValid() || trackCollection->size() == 0) return StatusCode::SUCCESS;
+
+ SG::ReadHandle<TrackCollection> segmentCollection {m_segmentCollection, ctx};
+ if (!segmentCollection.isValid()) return StatusCode::SUCCESS;
+
+ SG::ReadHandle<Tracker::FaserSCT_ClusterContainer> clusterContainer { m_clusterContainer, ctx };
+ if (!clusterContainer.isValid()) return StatusCode::SUCCESS;
+
+ SG::ReadHandle<xAOD::WaveformHitContainer> triggerContainer { m_triggerContainer, ctx };
+ if (!triggerContainer.isValid()) return StatusCode::SUCCESS;
+
+ // segment cuts
+
+ std::vector<int> segmentCount {0, 0, 0, 0};
+ int nSegments {0};
+
+ for (const Trk::Track* segment : *segmentCollection)
+ {
+ if (segment == nullptr) continue;
+ nSegments++;
+ for (auto measurement : *(segment->measurementsOnTrack()))
+ {
+ const Tracker::FaserSCT_ClusterOnTrack* cluster = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(measurement);
+ if (cluster != nullptr)
+ {
+ Identifier id = cluster->identify();
+ int station = m_sctHelper->station(id);
+ if (station >= 0 && station <= 3)
+ {
+ segmentCount[station]++;
+ break;
+ }
+ }
+ }
+ }
+ if ((nSegments > m_maxSegmentsTotal) || (segmentCount[1] > m_maxSegmentsStation) || (segmentCount[2] > m_maxSegmentsStation) || (segmentCount[3] > m_maxSegmentsStation))
+ return StatusCode::SUCCESS;
+
+ // track cuts
+
+ int nTracks {0};
+ int nLongTracks {0};
+ double chi2PerDoF {0};
+ double charge {0};
+ double maxRadius {0};
+ double momentum {0};
+ const Trk::Track* theTrack {nullptr};
+
+ for (const Trk::Track* track : *trackCollection)
+ {
+ if (track == nullptr) continue;
+ nTracks++;
+ std::set<int> stationMap;
+ std::set<std::pair<int, int>> layerMap;
+
+ // Check for hit in the three downstream stations
+ for (auto measurement : *(track->measurementsOnTrack()))
+ {
+ const Tracker::FaserSCT_ClusterOnTrack* cluster = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(measurement);
+ if (cluster != nullptr)
+ {
+ Identifier id = cluster->identify();
+ int station = m_sctHelper->station(id);
+ int layer = m_sctHelper->layer(id);
+ stationMap.emplace(station);
+ layerMap.emplace(station, layer);
+ }
+ }
+ if (stationMap.count(1) == 0 || stationMap.count(2) == 0 || stationMap.count(3) == 0) continue;
+
+ int nLayers = std::count_if(layerMap.begin(), layerMap.end(), [](std::pair<int,int> p){return p.first != 0;});
+ if (nLayers < m_minLayers) continue;
+ nLongTracks++;
+ const Trk::TrackParameters* upstreamParameters {nullptr};
+ for (auto params : *(track->trackParameters()))
+ {
+ if (params->position().z() < 0) continue; // Ignore IFT hits
+ double radius = sqrt(pow(params->position().x(), 2) + pow(params->position().y(), 2));
+ if (radius > maxRadius) maxRadius = radius;
+ if (upstreamParameters == nullptr || params->position().z() < upstreamParameters->position().z()) upstreamParameters = params;
+ }
+
+ momentum = upstreamParameters->momentum().mag()/1000;
+ chi2PerDoF = track->fitQuality()->chiSquared() / track->fitQuality()->numberDoF();
+ charge = upstreamParameters->charge();
+ theTrack = track;
+ }
+
+ if ((nTracks != 1) || (nLongTracks != 1) || (momentum < m_minMomentum) || (momentum > m_maxMomentum) || (chi2PerDoF > m_maxChi2PerDoF) || ((charge>0) != m_acceptPositive) || (maxRadius > m_maxRadius))
+ return StatusCode::SUCCESS;
+
+ // cluster cuts
+
+ // std::vector<int> clusterCount {0, 0, 0, 0};
+ // for (auto collection : *clusterContainer)
+ // {
+ // Identifier id = collection->identify();
+ // int station = m_sctHelper->station(id);
+ // if (station >= 0 && station <= 3) clusterCount[station] += collection->size();
+ // }
+ // if ((clusterCount[1] > m_maxClustersStation) || (clusterCount[2] > m_maxClustersStation) || (clusterCount[3] > m_maxClustersStation))
+ // return StatusCode::SUCCESS;
+
+ // waveform charge cuts
+
+ std::vector<double> timingCharge {0, 0, 0, 0};
+ const int baseChannel = 8;
+
+ for (auto hit : *triggerContainer) {
+ if ((hit->hit_status()&2)==0) { // dont store secondary hits as they can overwrite the primary hit
+ int ch=hit->channel();
+ timingCharge[ch - baseChannel] += hit->integral()/50;
+ }
+ }
+ if ((timingCharge[0] + timingCharge[1] > m_maxTimingCharge) || (timingCharge[2] + timingCharge[3] > m_maxTimingCharge))
+ return StatusCode::SUCCESS;
+
+ // passes all cuts
+
+ // copy track without associated clusters
+
+ SG::WriteHandle<TrackCollection> outputTrackCollection {m_outputTrackCollection, ctx};
+ std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
+
+ Trk::TrackInfo i { theTrack->info() };
+ Trk::FitQuality* q = new Trk::FitQuality { *(theTrack->fitQuality()) };
+ DataVector<const Trk::TrackStateOnSurface>* s = new DataVector<const Trk::TrackStateOnSurface> {};
+
+ for (auto oldState : (*theTrack->trackStateOnSurfaces()))
+ {
+ const Trk::CurvilinearParameters* oldParam = dynamic_cast<const Trk::CurvilinearParameters*>(oldState->trackParameters());
+ std::unique_ptr<Trk::TrackParameters> newParam { new Trk::CurvilinearParameters { *oldParam } };
+ Trk::TrackStateOnSurface* newState = new Trk::TrackStateOnSurface { nullptr, newParam.release() };
+ s->push_back(newState);
+ }
+ std::unique_ptr<DataVector<const Trk::TrackStateOnSurface>> sink {s};
+ outputTracks->push_back(new Trk::Track {i, std::move(*sink) , q } );
+
+ ATH_CHECK(outputTrackCollection.record(std::move(outputTracks)));
+
+ // set filter passed
+
+ setFilterPassed(true, ctx);
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode SelectRDOAlg::finalize()
+{
+ return StatusCode::SUCCESS;
+}
+
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..842282d534a23159625a1dabf59c3fbc0494c343
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
@@ -0,0 +1,43 @@
+#pragma once
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "TrkTrack/TrackCollection.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterContainer.h"
+#include "xAODFaserWaveform/WaveformHitContainer.h"
+#include "xAODFaserWaveform/WaveformHit.h"
+#include <xAODEventInfo/EventInfo.h>
+
+
+class FaserSCT_ID;
+
+class SelectRDOAlg : public AthReentrantAlgorithm {
+public:
+ SelectRDOAlg(const std::string &name, ISvcLocator *pSvcLocator);
+ virtual ~SelectRDOAlg() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute(const EventContext &ctx) const override;
+ virtual StatusCode finalize() override;
+
+private:
+
+ const FaserSCT_ID* m_sctHelper;
+
+ SG::ReadHandleKey<TrackCollection> m_trackCollection { this, "TrackCollection", "CKFTrackCollectionWithoutIFT", "Input track collection name" };
+ SG::ReadHandleKey<TrackCollection> m_segmentCollection { this, "SegmentCollection", "SegmentFit", "Track segment collection name" };
+ SG::ReadHandleKey<Tracker::FaserSCT_ClusterContainer> m_clusterContainer { this, "ClusterContainer", "SCT_ClusterContainer" };
+ SG::ReadHandleKey<xAOD::WaveformHitContainer> m_triggerContainer { this, "TriggerContainer", "TriggerWaveformHits", "Trigger hit container name" };
+ SG::WriteHandleKey<TrackCollection> m_outputTrackCollection { this, "OutputTrackCollection", "Pos_CKFTrackCollectionWithoutIFT", "Output track collection name"};
+
+ IntegerProperty m_minLayers {this, "MinLayers", 7, "Minimum hit layers to accept track"};
+ DoubleProperty m_maxChi2PerDoF {this, "MaxChi2PerDoF", 25.0, "Maximum chi2 per degree of freedom to accept track"};
+ DoubleProperty m_minMomentum {this, "MinMomentumGeV", 50.0, "Minimum momentum in GeV to accept track"};
+ DoubleProperty m_maxMomentum {this, "MaxMomentumGeV", 5000.0, "Maximum momentum in GeV to accept track"};
+ DoubleProperty m_maxRadius {this, "MaxRadiusMm", 95.0, "Maximum radius at first measurement to accept track"};
+ IntegerProperty m_maxSegmentsTotal {this, "MaxSegmentsTotal", 4, "Maximum number of segments in three stations to accept track"};
+ IntegerProperty m_maxSegmentsStation {this, "MaxSegmentsStation", 2, "Maximum number of segments in any single station to accept track"};
+ // IntegerProperty m_maxClustersStation {this, "MaxClustersStation", 9, "Maximum number of clusters in any single station to accept track"};
+ DoubleProperty m_maxTimingCharge {this, "MaxTimingChargePc", 70.0, "Maximum charge in pC recorded by upper or lower timing scintillator"};
+ BooleanProperty m_acceptPositive {this, "AcceptPositive", true, "Accept positive (true) or negative (false) charged tracks"};
+
+};
+
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..13e61cb84654d0c3e4096ffa01b0f6efbeb6184c
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx
@@ -0,0 +1,9 @@
+#include "../OverlayRDOAlg.h"
+#include "../SelectRDOAlg.h"
+#include "../BinRDOAlg.h"
+#include "../HistoRDOAlg.h"
+
+DECLARE_COMPONENT(OverlayRDOAlg)
+DECLARE_COMPONENT(SelectRDOAlg)
+DECLARE_COMPONENT(BinRDOAlg)
+DECLARE_COMPONENT(HistoRDOAlg)
diff --git a/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py b/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py
index f0de1582f2ccf47270203d4d1d6eb9a6bca341e9..eeaadd061bbc5b1c3c14af0a1e128e529d687cae 100644
--- a/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py
+++ b/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py
@@ -15,9 +15,9 @@ def SegmentFitAlgBasicCfg(flags, **kwargs):
Tracker__SegmentFitAlg=CompFactory.Tracker.SegmentFitAlg
acc.addEventAlgo(Tracker__SegmentFitAlg(**kwargs))
- thistSvc = CompFactory.THistSvc()
- thistSvc.Output += ["HIST DATAFILE='SegmentFitHistograms.root' OPT='RECREATE'"]
- acc.addService(thistSvc)
+ # thistSvc = CompFactory.THistSvc()
+ # thistSvc.Output += ["HIST DATAFILE='SegmentFitHistograms.root' OPT='RECREATE'"]
+ # acc.addService(thistSvc)
return acc
# with output defaults
diff --git a/Tracker/TrackerRecTools/FaserSiSpacePointTool/src/TrackerSpacePointMakerTool.cxx b/Tracker/TrackerRecTools/FaserSiSpacePointTool/src/TrackerSpacePointMakerTool.cxx
index 814fc31a5b82c9705472a283ba8abfd074ed067b..c8a7df00173d6503683ac47148baefd8bf54da4c 100644
--- a/Tracker/TrackerRecTools/FaserSiSpacePointTool/src/TrackerSpacePointMakerTool.cxx
+++ b/Tracker/TrackerRecTools/FaserSiSpacePointTool/src/TrackerSpacePointMakerTool.cxx
@@ -122,7 +122,7 @@ FaserSCT_SpacePoint* TrackerSpacePointMakerTool::makeSCT_SpacePoint(const Tracke
{
if (fabs(lambda1) > 1 + m_stripLengthTolerance)
{
- ATH_MSG_WARNING("Intersection lies outside the bounds of both strips");
+ ATH_MSG_DEBUG("Intersection lies outside the bounds of both strips");
ok = false;
}
}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
index f990be9458e0ee2ac0a5bdef8b670703f1424c26..14616c6f42dd0e60df67c5b73358da599093130c 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
@@ -50,6 +50,7 @@ StatusCode CKF2::initialize() {
ATH_CHECK(m_kalmanFitterTool1.retrieve());
ATH_CHECK(m_createTrkTrackTool.retrieve());
ATH_CHECK(m_trackCollection.initialize());
+ // ATH_CHECK(m_allTrackCollection.initialize());
ATH_CHECK(m_eventInfoKey.initialize());
if (m_performanceWriter && !m_noDiagnostics) {
ATH_CHECK(m_performanceWriterTool.retrieve());
@@ -85,6 +86,9 @@ StatusCode CKF2::execute() {
SG::WriteHandle trackContainer{m_trackCollection, ctx};
std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
+ // SG::WriteHandle allTrackContainer{m_allTrackCollection, ctx};
+ // std::unique_ptr<TrackCollection> outputAllTracks = std::make_unique<TrackCollection>();
+
std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry
= m_trackingGeometryTool->trackingGeometry();
@@ -157,10 +161,19 @@ StatusCode CKF2::execute() {
// TODO use status bits for different errors
// result.error() == Acts::CombinatorialKalmanFilterError::NoTrackFound
if (result.error() == Acts::PropagatorError::StepCountLimitReached ||
- result.error() == Acts::CombinatorialKalmanFilterError::PropagationReachesMaxSteps) {
- if (!eventInfo->updateErrorState(xAOD::EventInfo::SCT, xAOD::EventInfo::Error)) {
- ATH_MSG_WARNING (" cannot set error state for SCT");
- }
+ result.error() == Acts::CombinatorialKalmanFilterError::PropagationReachesMaxSteps)
+ {
+ try
+ {
+ if (!eventInfo->updateErrorState(xAOD::EventInfo::SCT, xAOD::EventInfo::Error))
+ {
+ ATH_MSG_WARNING ("Cannot set error state for SCT.");
+ }
+ }
+ catch (...)
+ {
+ ATH_MSG_DEBUG ("SCT error state is locked.");
+ }
}
continue;
}
@@ -183,6 +196,33 @@ StatusCode CKF2::execute() {
// select all tracks with at least 13 heats and with 6 or less shared hits, starting from the best track
// TODO use Gaudi parameters for the number of hits and shared hits
// TODO allow shared hits only in the first station?
+ // std::vector<FaserActsRecMultiTrajectory> rawTrajectories {};
+ // for (auto raw : allTrajectories)
+ // {
+ // rawTrajectories.push_back(raw.trajectory);
+ // }
+ // for (const FaserActsRecMultiTrajectory &traj : rawTrajectories) {
+ // const auto params = traj.trackParameters(traj.tips().front());
+ // ATH_MSG_DEBUG("Fitted parameters (raw)");
+ // ATH_MSG_DEBUG(" params: " << params.parameters().transpose());
+ // ATH_MSG_DEBUG(" position: " << params.position(gctx).transpose());
+ // ATH_MSG_DEBUG(" momentum: " << params.momentum().transpose());
+ // ATH_MSG_DEBUG(" charge: " << params.charge());
+ // std::unique_ptr<Trk::Track> track = m_createTrkTrackTool->createTrack(gctx, traj);
+ // if (track != nullptr) {
+ // m_numberOfSelectedTracks++;
+ // std::unique_ptr<Trk::Track> track2 = m_kalmanFitterTool1->fit(ctx, gctx, track.get(), Acts::BoundVector::Zero(), m_isMC, origin);
+ // if (track2) {
+ // outputAllTracks->push_back(std::move(track2));
+ // } else {
+ // outputAllTracks->push_back(std::move(track));
+ // ATH_MSG_WARNING("Re-Fit failed.");
+ // }
+ // } else {
+ // ATH_MSG_WARNING("CKF failed.");
+ // }
+ // }
+
std::vector<FaserActsRecMultiTrajectory> selectedTrajectories {};
while (not allTrajectories.empty()) {
TrajectoryInfo selected = allTrajectories.front();
@@ -224,6 +264,7 @@ StatusCode CKF2::execute() {
if (m_performanceWriter && !m_noDiagnostics) {
ATH_CHECK(m_performanceWriterTool->write(gctx, selectedTrajectories));
}
+ // ATH_CHECK(allTrackContainer.record(std::move(outputAllTracks)));
ATH_CHECK(trackContainer.record(std::move(outputTracks)));
return StatusCode::SUCCESS;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
index 170f0353904bef82095146c50bf98d5721ef2138..cd78bf72531931fc3ef7ef3b1b67edc854e46596 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
@@ -149,6 +149,7 @@ private:
ToolHandle<CreateTrkTrackTool> m_createTrkTrackTool {this, "CreateTrkTrackTool", "CreateTrkTrackTool"};
Gaudi::Property<bool> m_isMC {this, "isMC", false};
SG::WriteHandleKey<TrackCollection> m_trackCollection { this, "OutputCollection", "CKFTrackCollection", "Output track collection name" };
+ // SG::WriteHandleKey<TrackCollection> m_allTrackCollection { this, "AllTrackCollection", "CKFAllTrackCollection", "Output all track collection name" };
SG::WriteDecorHandleKey<xAOD::EventInfo> m_eventInfoKey{this, "EventInfoKey", "EventInfo"};
};
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
index 952f75287f594e5df31feb2b6a493bfe235a4086..bbae70e3b10be38ff5de9675796bbe509c8f2483 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
@@ -130,10 +130,16 @@ StatusCode CircleFitTrackSeedTool::run(std::vector<int> maskedLayers) {
std::list<Seed> allSeeds;
for (const Seed &seed : seeds) allSeeds.push_back(seed);
+ // allSeeds.sort([](const Seed &left, const Seed &right) {
+ // if (left.size > right.size) return true;
+ // if (left.size < right.size) return false;
+ // if (left.chi2 < right.chi2) return true;
+ // else return false;
+ // });
allSeeds.sort([](const Seed &left, const Seed &right) {
- if (left.size > right.size) return true;
- if (left.size < right.size) return false;
- if (left.chi2 < right.chi2) return true;
+ if (left.stations > right.stations) return true;
+ if (left.stations < right.stations) return false;
+ if (left.chi2/std::max(left.positions.size() + left.fakePositions.size() - left.constraints, 1UL) < right.chi2/std::max(right.positions.size() + right.fakePositions.size() - right.constraints, 1UL)) return true;
else return false;
});
@@ -141,6 +147,9 @@ StatusCode CircleFitTrackSeedTool::run(std::vector<int> maskedLayers) {
while (not allSeeds.empty()) {
Seed selected = allSeeds.front();
selectedSeeds.push_back(selected);
+ // allSeeds.remove_if([&](const Seed &p) {
+ // return ((p.clusterSet & selected.clusterSet).count() == p.clusterSet.count());
+ // });
allSeeds.remove_if([&](const Seed &p) {
return (p.size < 10) || ((p.clusterSet & selected.clusterSet).count() > 6);
});
@@ -161,7 +170,9 @@ StatusCode CircleFitTrackSeedTool::run(std::vector<int> maskedLayers) {
double origin = !selectedSeeds.empty() ? minSeed->minZ - 10 : 0;
m_targetZPosition = origin;
std::vector<Acts::CurvilinearTrackParameters> initParams {};
+ ATH_MSG_DEBUG("Sorted seed properties:");
for (const Seed &seed : selectedSeeds) {
+ ATH_MSG_DEBUG("seed size: " << seed.size << ", chi2: " << seed.chi2);
initParams.push_back(seed.get_params(origin, cov));
}
@@ -264,13 +275,16 @@ CircleFitTrackSeedTool::Seed::Seed(const std::vector<Segment> &segments) :
if (segments.size() > 1) {
fakeFit();
+ constraints = 5;
} else {
momentum = 9999999.;
charge = 1;
+ constraints = 2;
}
getChi2();
size = clusters.size();
+ stations = segments.size();
}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.h b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.h
index 0467c92620eebd955f0a78f11cfdea3fd0cb0b1d..a4d53c629c24f53e1dfe4488ca7524755f117ad3 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.h
@@ -66,7 +66,7 @@ private:
double c0, c1, cx, cy, r, chi2, momentum, charge, minZ;
Acts::Vector3 direction;
- size_t size;
+ size_t size, stations, constraints;
Acts::CurvilinearTrackParameters get_params(double origin, Acts::BoundSymMatrix cov) const;
private:
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx
index a7150f62b5daa6d951b5b27c112c8d71797c90b9..994689010203a269946c3c88392004a61a42adf4 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx
@@ -62,15 +62,16 @@ CreateTrkTrackTool::createTrack(const Acts::GeometryContext &gctx, const FaserAc
}
}
double nDoF = state.calibratedSize();
- const Trk::FitQualityOnSurface* quality = new Trk::FitQualityOnSurface(state.chi2(), nDoF);
- const Trk::TrackStateOnSurface *perState = new Trk::TrackStateOnSurface(
- clusterOnTrack,
- parm,
- quality,
- nullptr,
- typePattern);
- if (perState) {
- finalTrajectory->insert(finalTrajectory->begin(), perState);
+ const Trk::FitQualityOnSurface* quality = new Trk::FitQualityOnSurface(state.chi2(), nDoF);
+ const Trk::TrackStateOnSurface* perState = new Trk::TrackStateOnSurface(
+ clusterOnTrack,
+ parm,
+ quality,
+ nullptr,
+ typePattern);
+ if (perState)
+ {
+ finalTrajectory->insert(finalTrajectory->begin(), perState);
}
}
return;