From 5e3e7f86a2c340251ac45c2adbe34ae471348a22 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Mon, 27 Feb 2023 20:08:57 -0800
Subject: [PATCH 01/21] Work in progress on overlay

---
 .../TrackerRecAlgs/OverlayRDO/CMakeLists.txt  |  12 ++
 .../OverlayRDO/python/OverlayRDOConfig.py     |  91 +++++++++++++++
 .../OverlayRDO/python/SelectRDOConfig.py      | 108 ++++++++++++++++++
 .../OverlayRDO/src/OverlayRDOAlg.cxx          |  51 +++++++++
 .../OverlayRDO/src/OverlayRDOAlg.h            |  29 +++++
 .../OverlayRDO/src/SelectRDOAlg.cxx           |  78 +++++++++++++
 .../OverlayRDO/src/SelectRDOAlg.h             |  32 ++++++
 .../src/component/OverlayRDO_entries.cxx      |   5 +
 8 files changed, 406 insertions(+)
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx

diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt b/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
new file mode 100644
index 000000000..50aecc583
--- /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/OverlayRDOAlg.cxx src/SelectRDOAlg.cxx
+        src/component/OverlayRDO_entries.cxx
+        LINK_LIBRARIES AthenaBaseComps StoreGateLib TrackerRawData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack
+)
+
+atlas_install_python_modules(python/*.py)
+# atlas_install_scripts(test/*.py)
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
new file mode 100644
index 000000000..f31a7cbd5
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
@@ -0,0 +1,91 @@
+"""
+    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)
+
+    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.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=10)
+
+    # 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 000000000..83d90a541
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
@@ -0,0 +1,108 @@
+"""
+    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)
+
+    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_EDGEMODE_RDOs"]
+    # ItemList += ["Tracker::FaserSCT_ClusterContainer#" + chargePrefix + "SCT_ClusterContainer"]
+
+    acc.merge(OutputStreamCfg(ConfigFlags,"RDO", ItemList=ItemList, disableEventTag=True))
+    osrdo = acc.getEventAlgo("OutputStreamRDO")
+    osrdo.AcceptAlgs += ["SelectRDOAlg"]
+
+    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']
+    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"))
+
+    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_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/OverlayRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
new file mode 100644
index 000000000..afac3427b
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
@@ -0,0 +1,51 @@
+#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());
+
+  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());
+
+
+  return StatusCode::SUCCESS;
+}
+
+
+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 000000000..365ef69aa
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
@@ -0,0 +1,29 @@
+#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:
+
+  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_EDGEMODE_RDOs"};
+  SG::ReadHandleKey<FaserSCT_RDO_Container> m_negRdoKey { this, "NegRdoContainer", "Neg_SCT_EDGEMODE_RDOs"};
+
+
+  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 000000000..fc61c5b2a
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
@@ -0,0 +1,78 @@
+#include "SelectRDOAlg.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>
+
+SelectRDOAlg::SelectRDOAlg(const std::string &name, 
+                           ISvcLocator *pSvcLocator)
+    : AthReentrantAlgorithm(name, pSvcLocator)
+{}
+
+
+StatusCode SelectRDOAlg::initialize() 
+{
+  ATH_CHECK(m_trackCollection.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()) return StatusCode::SUCCESS;
+
+  // const Trk::TrackParameters* candidateParameters {nullptr};
+  // const Trk::Track* candidateTrack {nullptr};
+
+  int nTracks {0};
+  int nLongTracks {0};
+
+  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++;
+  }
+
+  setFilterPassed(nTracks == 1 && nLongTracks == 1, ctx);
+  // std::cout << "Tracks: " << nTracks << ", Long Tracks: " << nLongTracks << std::endl;
+  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 000000000..226d1b3a7
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
@@ -0,0 +1,32 @@
+#pragma once
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "TrkTrack/TrackCollection.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", "CKFTrackCollectionWithoutIDF", "Input 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", 1000.0, "Maximum momentum in GeV to accept track"};
+  DoubleProperty    m_maxRadius {this, "MaxRadiusMm", 85.0, "Maximum radius at first measurement to accept track"};
+  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 000000000..034c80451
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx
@@ -0,0 +1,5 @@
+#include "../OverlayRDOAlg.h"
+#include "../SelectRDOAlg.h"
+
+DECLARE_COMPONENT(OverlayRDOAlg)
+DECLARE_COMPONENT(SelectRDOAlg)
\ No newline at end of file
-- 
GitLab


From b044e939ce28f712101816c0544148c010d375fa Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Mon, 27 Feb 2023 21:29:36 -0800
Subject: [PATCH 02/21] Fixes and progress

---
 .../OverlayRDO/python/OverlayRDOConfig.py     | 16 ++++++++-
 .../OverlayRDO/src/OverlayRDOAlg.cxx          | 34 +++++++++++++++++++
 .../OverlayRDO/src/OverlayRDOAlg.h            |  3 ++
 .../OverlayRDO/src/SelectRDOAlg.cxx           | 17 ++++++++--
 4 files changed, 67 insertions(+), 3 deletions(-)

diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
index f31a7cbd5..3684f09ba 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
@@ -13,6 +13,19 @@ def OverlayRDOAlgCfg(flags, **kwargs):
     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#" + chargePrefix + "CKFTrackCollectionWithoutIFT"]
+    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__":
@@ -37,6 +50,7 @@ if __name__ == "__main__":
     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
@@ -85,7 +99,7 @@ if __name__ == "__main__":
     # ConfigFlags.dump()
 
     # Execute and finish
-    sc = acc.run(maxEvents=10)
+    sc = acc.run(maxEvents=-1)
 
     # Success should be 0
     sys.exit(not sc.isSuccess())    
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
index afac3427b..4fc9699d4 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
@@ -23,6 +23,7 @@ StatusCode OverlayRDOAlg::initialize()
   // ATH_CHECK(m_negTrackKey.initialize());
   ATH_CHECK(m_posRdoKey.initialize());
   ATH_CHECK(m_negRdoKey.initialize());
+  ATH_CHECK(m_outRdoKey.initialize());
 
   return StatusCode::SUCCESS;
 }
@@ -40,11 +41,44 @@ StatusCode OverlayRDOAlg::execute(const EventContext &ctx) const
   SG::ReadHandle<FaserSCT_RDO_Container> negRdoContainer {m_negRdoKey, ctx};
   ATH_CHECK(negRdoContainer.isValid());
 
+  std::map<Identifier, 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())));
+
 
   return StatusCode::SUCCESS;
 }
 
 
+void OverlayRDOAlg::processContainer(const FaserSCT_RDO_Container& container, std::map<Identifier, 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();
+      if (rdoDB.count(id))
+      {
+        rdoDB[id].push_back(collection);
+      }
+      else
+      {
+        rdoDB[id] = 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
index 365ef69aa..de9f5f1ff 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
@@ -17,11 +17,14 @@ public:
 
 private:
 
+  void processContainer(const FaserSCT_RDO_Container& container, std::map<Identifier, std::vector<const FaserSCT_RDO_Collection*> >& rdoDB) 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_EDGEMODE_RDOs"};
   SG::ReadHandleKey<FaserSCT_RDO_Container> m_negRdoKey { this, "NegRdoContainer", "Neg_SCT_EDGEMODE_RDOs"};
 
+  SG::WriteHandleKey<FaserSCT_RDO_Container> m_outRdoKey{this, "OutputObjectName", "SCT_EDGEMODE_RDOs", "Output Object name"};
 
   const FaserSCT_ID* m_sctHelper;
 
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
index fc61c5b2a..70714a842 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
@@ -37,6 +37,9 @@ StatusCode SelectRDOAlg::execute(const EventContext &ctx) const
 
   int nTracks {0};
   int nLongTracks {0};
+  double chi2PerDoF {0};
+  double charge {0};
+  double radius {0};
 
   for (const Trk::Track* track : *trackCollection)
   {
@@ -63,10 +66,20 @@ StatusCode SelectRDOAlg::execute(const EventContext &ctx) const
     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
+        if (upstreamParameters == nullptr || params->position().z() < upstreamParameters->position().z()) upstreamParameters = params;
+    }
+
+    chi2PerDoF = track->fitQuality()->chiSquared() / track->fitQuality()->numberDoF();
+    charge = upstreamParameters->charge();
+    radius = sqrt(pow(upstreamParameters->position().x(), 2) + pow(upstreamParameters->position().y(), 2));
   }
 
-  setFilterPassed(nTracks == 1 && nLongTracks == 1, ctx);
-  // std::cout << "Tracks: " << nTracks << ", Long Tracks: " << nLongTracks << std::endl;
+  setFilterPassed( (nTracks == 1) && (nLongTracks == 1)  && (chi2PerDoF <= m_maxChi2PerDoF) && ((charge > 0) == m_acceptPositive) && (radius <= m_maxRadius), ctx);
   return StatusCode::SUCCESS;
 }
 
-- 
GitLab


From 821b1f0b75abcb6ab58ddb3194c339c356c56d58 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Mon, 27 Feb 2023 23:07:04 -0800
Subject: [PATCH 03/21] Working overlay of RDO

---
 .../OverlayRDO/src/OverlayRDOAlg.cxx          | 48 +++++++++++++++++--
 .../OverlayRDO/src/OverlayRDOAlg.h            |  2 +-
 2 files changed, 44 insertions(+), 6 deletions(-)

diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
index 4fc9699d4..b5099c066 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
@@ -41,19 +41,56 @@ StatusCode OverlayRDOAlg::execute(const EventContext &ctx) const
   SG::ReadHandle<FaserSCT_RDO_Container> negRdoContainer {m_negRdoKey, ctx};
   ATH_CHECK(negRdoContainer.isValid());
 
-  std::map<Identifier, std::vector<const FaserSCT_RDO_Collection*> > rdoDB {};
+  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())));
 
+  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;
+        stripMap[stripNumber] |= dataWord;
+      }
+    }
+
+    for (auto stripEntry : stripMap)
+    {
+      Identifier rdoID {m_sctHelper->strip_id(id, stripEntry.first)};
+      current_collection->emplace_back(new FaserSCT3_RawData(rdoID, stripEntry.second, std::vector<int> {}));
+    }
+
+    outRdoContainer->getWriteHandle(waferHash).addOrDelete(std::move(current_collection));
+  }
+
 
   return StatusCode::SUCCESS;
 }
 
 
-void OverlayRDOAlg::processContainer(const FaserSCT_RDO_Container& container, std::map<Identifier, std::vector<const FaserSCT_RDO_Collection*> >& rdoDB) const
+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};
@@ -61,13 +98,14 @@ void OverlayRDOAlg::processContainer(const FaserSCT_RDO_Container& container, st
   {
       if (collection->empty()) continue;
       Identifier id = collection->identify();
-      if (rdoDB.count(id))
+      IdentifierHash hash = m_sctHelper->wafer_hash(id);
+      if (rdoDB.count(hash))
       {
-        rdoDB[id].push_back(collection);
+        rdoDB[hash].push_back(collection);
       }
       else
       {
-        rdoDB[id] = std::vector<const FaserSCT_RDO_Collection*>{collection};
+        rdoDB[hash] = std::vector<const FaserSCT_RDO_Collection*>{collection};
       }
       for(const auto& rawdata : *collection)
       {
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
index de9f5f1ff..ab7243c4b 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
@@ -17,7 +17,7 @@ public:
 
 private:
 
-  void processContainer(const FaserSCT_RDO_Container& container, std::map<Identifier, std::vector<const FaserSCT_RDO_Collection*> >& rdoDB) const;
+  void processContainer(const FaserSCT_RDO_Container& container, std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> >& rdoDB) 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" };
-- 
GitLab


From 412c7343efaf1ca9884688cb70f63e0cf318054f Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Tue, 28 Feb 2023 11:02:34 -0800
Subject: [PATCH 04/21] Add segment, cluster and waveform cuts

---
 .../TrackerRecAlgs/OverlayRDO/CMakeLists.txt  |  2 +-
 .../OverlayRDO/python/SelectRDOConfig.py      |  5 +-
 .../OverlayRDO/src/SelectRDOAlg.cxx           | 87 +++++++++++++++++--
 .../OverlayRDO/src/SelectRDOAlg.h             | 12 ++-
 4 files changed, 94 insertions(+), 12 deletions(-)

diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt b/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
index 50aecc583..d1213ed93 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
@@ -5,7 +5,7 @@ atlas_add_component(
         src/OverlayRDOAlg.h src/SelectRDOAlg.h
         src/OverlayRDOAlg.cxx src/SelectRDOAlg.cxx
         src/component/OverlayRDO_entries.cxx
-        LINK_LIBRARIES AthenaBaseComps StoreGateLib TrackerRawData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack
+        LINK_LIBRARIES AthenaBaseComps StoreGateLib TrackerRawData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack xAODFaserWaveform
 )
 
 atlas_install_python_modules(python/*.py)
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
index 83d90a541..875529769 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
@@ -49,7 +49,10 @@ if __name__ == "__main__":
     Configurable.configurableRun3Behavior = True
 
     # Configure
-    ConfigFlags.Input.Files = [ '/eos/experiment/faser/rec/2022/r0013/009171/Faser-Physics-009171-00006-r0013-xAOD.root']
+    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
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
index 70714a842..c9b1db69a 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
@@ -2,11 +2,10 @@
 #include "TrkTrack/Track.h"
 #include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h"
 #include "TrackerIdentifier/FaserSCT_ID.h"
-// #include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterContainer.h"
+#include "TrackerPrepRawData/FaserSCT_ClusterCollection.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>
 
 SelectRDOAlg::SelectRDOAlg(const std::string &name, 
@@ -18,6 +17,9 @@ SelectRDOAlg::SelectRDOAlg(const std::string &name,
 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(detStore()->retrieve(m_sctHelper, "FaserSCT_ID"));
 
@@ -30,16 +32,52 @@ StatusCode SelectRDOAlg::execute(const EventContext &ctx) const
 
   setFilterPassed(false, ctx);
   SG::ReadHandle<TrackCollection> trackCollection {m_trackCollection, ctx};
-  if (!trackCollection.isValid()) return StatusCode::SUCCESS;
+  if (!trackCollection.isValid() || trackCollection->size() == 0) return StatusCode::SUCCESS;
 
-  // const Trk::TrackParameters* candidateParameters {nullptr};
-  // const Trk::Track* candidateTrack {nullptr};
+  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 radius {0};
+  double momentum {0};
 
   for (const Trk::Track* track : *trackCollection)
   {
@@ -66,7 +104,6 @@ StatusCode SelectRDOAlg::execute(const EventContext &ctx) const
     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()))
     {
@@ -74,12 +111,44 @@ StatusCode SelectRDOAlg::execute(const EventContext &ctx) const
         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();
     radius = sqrt(pow(upstreamParameters->position().x(), 2) + pow(upstreamParameters->position().y(), 2));
   }
 
-  setFilterPassed( (nTracks == 1) && (nLongTracks == 1)  && (chi2PerDoF <= m_maxChi2PerDoF) && ((charge > 0) == m_acceptPositive) && (radius <= m_maxRadius), ctx);
+  if ((nTracks != 1) || (nLongTracks != 1) || (momentum < m_minMomentum) || (momentum > m_maxMomentum) || (chi2PerDoF > m_maxChi2PerDoF) || ((charge>0) != m_acceptPositive) || (radius > 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
+
+  setFilterPassed(true, ctx);
   return StatusCode::SUCCESS;
 }
 
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
index 226d1b3a7..d1867a9e1 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
@@ -2,6 +2,9 @@
 
 #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>
 
 
@@ -19,13 +22,20 @@ private:
 
   const FaserSCT_ID* m_sctHelper;
 
-  SG::ReadHandleKey<TrackCollection> m_trackCollection { this, "TrackCollection", "CKFTrackCollectionWithoutIDF", "Input track collection name" };
+  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" };
 
   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", 1000.0, "Maximum momentum in GeV to accept track"};
   DoubleProperty    m_maxRadius {this, "MaxRadiusMm", 85.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", 52.5, "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"};
 
 };
-- 
GitLab


From 31b931680a4b68b5c9c238657c398d87b64c6ba9 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Tue, 28 Feb 2023 12:27:16 -0800
Subject: [PATCH 05/21] Mods to get reco working; locked store exception occurs

---
 .../Reconstruction/scripts/faser_reco.py      | 30 +++++++++++--------
 .../OverlayRDO/python/OverlayRDOConfig.py     |  2 +-
 .../OverlayRDO/python/SelectRDOConfig.py      |  9 +++---
 .../OverlayRDO/src/OverlayRDOAlg.h            |  6 ++--
 4 files changed, 26 insertions(+), 21 deletions(-)

diff --git a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
index 21e1d51a6..3a85e15b1 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
+++ b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
@@ -37,6 +37,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()
 
@@ -173,7 +175,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:    
@@ -190,13 +192,14 @@ if useLHC:
     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
@@ -266,13 +269,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/python/OverlayRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
index 3684f09ba..d8cf85688 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
@@ -19,7 +19,7 @@ def OverlayRDOAlgCfg(flags, **kwargs):
     ItemList += ["xAOD::EventInfo#EventInfo"]
     ItemList += ["xAOD::EventAuxInfo#EventInfoAux."]
     # ItemList += ["TrackCollection#" + chargePrefix + "CKFTrackCollectionWithoutIFT"]
-    ItemList += ["FaserSCT_RDO_Container#SCT_EDGEMODE_RDOs"]
+    ItemList += ["FaserSCT_RDO_Container#SCT_RDOs"]
     # ItemList += ["Tracker::FaserSCT_ClusterContainer#" + chargePrefix + "SCT_ClusterContainer"]
 
     acc.merge(OutputStreamCfg(ConfigFlags,"RDO", ItemList=ItemList, disableEventTag=True))
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
index 875529769..f2ea95a02 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
@@ -22,7 +22,7 @@ def SelectRDOAlgCfg(flags, **kwargs):
     ItemList += ["xAOD::EventInfo#EventInfo"]
     ItemList += ["xAOD::EventAuxInfo#EventInfoAux."]
     # ItemList += ["TrackCollection#" + chargePrefix + "CKFTrackCollectionWithoutIFT"]
-    ItemList += ["FaserSCT_RDO_Container#" + chargePrefix + "SCT_EDGEMODE_RDOs"]
+    ItemList += ["FaserSCT_RDO_Container#" + chargePrefix + "SCT_RDOs"]
     # ItemList += ["Tracker::FaserSCT_ClusterContainer#" + chargePrefix + "SCT_ClusterContainer"]
 
     acc.merge(OutputStreamCfg(ConfigFlags,"RDO", ItemList=ItemList, disableEventTag=True))
@@ -49,8 +49,9 @@ if __name__ == "__main__":
     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 = [ #'/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"
@@ -79,7 +80,7 @@ if __name__ == "__main__":
         # "xAOD::EventInfo#EventInfo->" + chargePrefix + "EventInfo",
         # "xAOD::EventAuxInfo#EventInfoAux.->" + chargePrefix + "EventInfoAux.",
         # "TrackCollection#CKFTrackCollectionWithoutIFT->" + chargePrefix + "CKFTrackCollectionWithoutIFT",
-        "FaserSCT_RDO_Container#SCT_EDGEMODE_RDOs->" + chargePrefix + "SCT_EDGEMODE_RDOs",
+        "FaserSCT_RDO_Container#SCT_RDOs->" + chargePrefix + "SCT_RDOs",
         # "Tracker::FaserSCT_ClusterContainer#SCT_ClusterContainer->" + chargePrefix + "SCT_ClusterContainer"
     ]))
 
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
index ab7243c4b..cf829e5a4 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
@@ -21,10 +21,10 @@ private:
 
   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_EDGEMODE_RDOs"};
-  SG::ReadHandleKey<FaserSCT_RDO_Container> m_negRdoKey { this, "NegRdoContainer", "Neg_SCT_EDGEMODE_RDOs"};
+  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::WriteHandleKey<FaserSCT_RDO_Container> m_outRdoKey{this, "OutputObjectName", "SCT_EDGEMODE_RDOs", "Output Object name"};
+  SG::WriteHandleKey<FaserSCT_RDO_Container> m_outRdoKey{this, "OutputObjectName", "SCT_RDOs", "Output Object name"};
 
   const FaserSCT_ID* m_sctHelper;
 
-- 
GitLab


From a97594af65f3b11f436880f8be18b1f63d4d7588 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Wed, 1 Mar 2023 08:51:57 -0800
Subject: [PATCH 06/21] Handle locked error word, don't fit IFT

---
 .../Reconstruction/scripts/faser_reco.py        |  3 ++-
 .../Acts/FaserActsKalmanFilter/src/CKF2.cxx     | 17 +++++++++++++----
 2 files changed, 15 insertions(+), 5 deletions(-)

diff --git a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
index 3a85e15b1..279127fb4 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
+++ b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
@@ -229,7 +229,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", 
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
index f990be945..6b02bc244 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
@@ -157,10 +157,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_WARNING ("SCT error state is locked.");
+          }
       }
       continue;
     }
-- 
GitLab


From 4e5558a7819ef85ccfa2ffc183f1748cd30f34f8 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Wed, 1 Mar 2023 08:55:16 -0800
Subject: [PATCH 07/21] Tweak reco for overlay

---
 Control/CalypsoExample/Reconstruction/scripts/faser_reco.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
index 279127fb4..cc608ff03 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
+++ b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
@@ -187,7 +187,7 @@ 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))
 
@@ -250,7 +250,7 @@ itemList = [ "xAOD::EventInfo#*"
              , "TrackCollection#*"
 ]
 #
-if useLHC:
+if useLHC and not args.isOverlay:
     itemList.extend( ["xAOD::FaserLHCData#*", "xAOD::FaserLHCDataAux#*"] )
 
 if args.isMC:
-- 
GitLab


From 99544bba2c9bbfebd1f3a151627ecd4c22586871 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Thu, 2 Mar 2023 08:28:51 -0800
Subject: [PATCH 08/21] Preserve single-track fit info during overlay

---
 .../OverlayRDO/python/OverlayRDOConfig.py     |  2 +-
 .../OverlayRDO/python/SelectRDOConfig.py      |  4 +-
 .../OverlayRDO/src/OverlayRDOAlg.cxx          | 43 ++++++++++++++++---
 .../OverlayRDO/src/OverlayRDOAlg.h            |  2 +
 .../OverlayRDO/src/SelectRDOAlg.cxx           | 26 +++++++++++
 .../OverlayRDO/src/SelectRDOAlg.h             |  1 +
 6 files changed, 69 insertions(+), 9 deletions(-)

diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
index d8cf85688..a9f643e3c 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
@@ -18,7 +18,7 @@ def OverlayRDOAlgCfg(flags, **kwargs):
     # ItemList += ["xAOD::EventAuxInfo#" + chargePrefix + "EventInfoAux."]
     ItemList += ["xAOD::EventInfo#EventInfo"]
     ItemList += ["xAOD::EventAuxInfo#EventInfoAux."]
-    # ItemList += ["TrackCollection#" + chargePrefix + "CKFTrackCollectionWithoutIFT"]
+    ItemList += ["TrackCollection#Orig_CKFTrackCollectionWithoutIFT"]
     ItemList += ["FaserSCT_RDO_Container#SCT_RDOs"]
     # ItemList += ["Tracker::FaserSCT_ClusterContainer#" + chargePrefix + "SCT_ClusterContainer"]
 
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
index f2ea95a02..95b7b692d 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
@@ -21,7 +21,7 @@ def SelectRDOAlgCfg(flags, **kwargs):
     # ItemList += ["xAOD::EventAuxInfo#" + chargePrefix + "EventInfoAux."]
     ItemList += ["xAOD::EventInfo#EventInfo"]
     ItemList += ["xAOD::EventAuxInfo#EventInfoAux."]
-    # ItemList += ["TrackCollection#" + chargePrefix + "CKFTrackCollectionWithoutIFT"]
+    ItemList += ["TrackCollection#" + chargePrefix + "CKFTrackCollectionWithoutIFT"]
     ItemList += ["FaserSCT_RDO_Container#" + chargePrefix + "SCT_RDOs"]
     # ItemList += ["Tracker::FaserSCT_ClusterContainer#" + chargePrefix + "SCT_ClusterContainer"]
 
@@ -73,7 +73,7 @@ if __name__ == "__main__":
     acc.merge(PoolReadCfg(ConfigFlags))
 
     # algorithm
-    acc.merge(SelectRDOAlgCfg(ConfigFlags, TrackCollection = "CKFTrackCollectionWithoutIFT"))
+    acc.merge(SelectRDOAlgCfg(ConfigFlags, TrackCollection = "CKFTrackCollectionWithoutIFT", OutputTrackCollection = chargePrefix + "CKFTrackCollectionWithoutIFT"))
 
     from SGComps.AddressRemappingConfig import AddressRemappingCfg
     acc.merge(AddressRemappingCfg([
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
index b5099c066..0384d5d03 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
@@ -19,11 +19,12 @@ OverlayRDOAlg::OverlayRDOAlg(const std::string &name,
 StatusCode OverlayRDOAlg::initialize() 
 {
   ATH_CHECK(detStore()->retrieve(m_sctHelper, "FaserSCT_ID"));
-  // ATH_CHECK(m_posTrackKey.initialize());
-  // ATH_CHECK(m_negTrackKey.initialize());
+  ATH_CHECK(m_posTrackKey.initialize());
+  ATH_CHECK(m_negTrackKey.initialize());
   ATH_CHECK(m_posRdoKey.initialize());
   ATH_CHECK(m_negRdoKey.initialize());
   ATH_CHECK(m_outRdoKey.initialize());
+  ATH_CHECK(m_outputTrackCollection.initialize());
 
   return StatusCode::SUCCESS;
 }
@@ -32,10 +33,10 @@ StatusCode OverlayRDOAlg::initialize()
 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<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};
@@ -85,10 +86,40 @@ StatusCode OverlayRDOAlg::execute(const EventContext &ctx) const
     outRdoContainer->getWriteHandle(waferHash).addOrDelete(std::move(current_collection));
   }
 
+  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;
 }
 
+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
 {
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
index cf829e5a4..57900ef49 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
@@ -18,6 +18,7 @@ public:
 private:
 
   void processContainer(const FaserSCT_RDO_Container& container, std::map<IdentifierHash, std::vector<const FaserSCT_RDO_Collection*> >& rdoDB) 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" };
@@ -25,6 +26,7 @@ private:
   SG::ReadHandleKey<FaserSCT_RDO_Container> m_negRdoKey { this, "NegRdoContainer", "Neg_SCT_RDOs"};
 
   SG::WriteHandleKey<FaserSCT_RDO_Container> m_outRdoKey{this, "OutputObjectName", "SCT_RDOs", "Output 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
index c9b1db69a..2bb95c96c 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
@@ -20,6 +20,7 @@ StatusCode SelectRDOAlg::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"));
 
@@ -78,6 +79,7 @@ StatusCode SelectRDOAlg::execute(const EventContext &ctx) const
   double charge {0};
   double radius {0};
   double momentum {0};
+  const Trk::Track* theTrack {nullptr};
 
   for (const Trk::Track* track : *trackCollection)
   {
@@ -115,6 +117,7 @@ StatusCode SelectRDOAlg::execute(const EventContext &ctx) const
     chi2PerDoF = track->fitQuality()->chiSquared() / track->fitQuality()->numberDoF();
     charge = upstreamParameters->charge();
     radius = sqrt(pow(upstreamParameters->position().x(), 2) + pow(upstreamParameters->position().y(), 2));
+    theTrack = track;
   }
 
   if ((nTracks != 1) || (nLongTracks != 1) || (momentum < m_minMomentum) || (momentum > m_maxMomentum) || (chi2PerDoF > m_maxChi2PerDoF) || ((charge>0) != m_acceptPositive) || (radius > m_maxRadius))
@@ -148,6 +151,29 @@ StatusCode SelectRDOAlg::execute(const EventContext &ctx) const
 
   // 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;
 }
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
index d1867a9e1..4c3bee533 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
@@ -26,6 +26,7 @@ private:
   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"};
-- 
GitLab


From 2d77fc67f078a1aa03c81700d25298067c4fbe3c Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Thu, 2 Mar 2023 18:39:56 +0100
Subject: [PATCH 09/21] Algorithms to overlay tracker RDO data

---
 .../Reconstruction/scripts/faser_reco.py      |  37 ++--
 .../TrackerRecAlgs/OverlayRDO/CMakeLists.txt  |  12 ++
 .../OverlayRDO/python/OverlayRDOConfig.py     | 105 ++++++++++
 .../OverlayRDO/python/SelectRDOConfig.py      | 112 +++++++++++
 .../OverlayRDO/src/OverlayRDOAlg.cxx          | 154 +++++++++++++++
 .../OverlayRDO/src/OverlayRDOAlg.h            |  34 ++++
 .../OverlayRDO/src/SelectRDOAlg.cxx           | 186 ++++++++++++++++++
 .../OverlayRDO/src/SelectRDOAlg.h             |  43 ++++
 .../src/component/OverlayRDO_entries.cxx      |   5 +
 .../Acts/FaserActsKalmanFilter/src/CKF2.cxx   |  17 +-
 10 files changed, 685 insertions(+), 20 deletions(-)
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx

diff --git a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
index 21e1d51a6..cc608ff03 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
+++ b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
@@ -37,6 +37,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()
 
@@ -173,7 +175,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,18 +187,19 @@ 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
@@ -226,7 +229,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,7 +250,7 @@ itemList = [ "xAOD::EventInfo#*"
              , "TrackCollection#*"
 ]
 #
-if useLHC:
+if useLHC and not args.isOverlay:
     itemList.extend( ["xAOD::FaserLHCData#*", "xAOD::FaserLHCDataAux#*"] )
 
 if args.isMC:
@@ -266,13 +270,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 000000000..d1213ed93
--- /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/OverlayRDOAlg.cxx src/SelectRDOAlg.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/OverlayRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
new file mode 100644
index 000000000..a9f643e3c
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
@@ -0,0 +1,105 @@
+"""
+    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 += ["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 000000000..95b7b692d
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
@@ -0,0 +1,112 @@
+"""
+    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)
+
+    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 += ["Tracker::FaserSCT_ClusterContainer#" + chargePrefix + "SCT_ClusterContainer"]
+
+    acc.merge(OutputStreamCfg(ConfigFlags,"RDO", ItemList=ItemList, disableEventTag=True))
+    osrdo = acc.getEventAlgo("OutputStreamRDO")
+    osrdo.AcceptAlgs += ["SelectRDOAlg"]
+
+    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",
+        # "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/OverlayRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
new file mode 100644
index 000000000..0384d5d03
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
@@ -0,0 +1,154 @@
+#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_outRdoKey.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());
+
+  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())));
+
+  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;
+        stripMap[stripNumber] |= dataWord;
+      }
+    }
+
+    for (auto stripEntry : stripMap)
+    {
+      Identifier rdoID {m_sctHelper->strip_id(id, stripEntry.first)};
+      current_collection->emplace_back(new FaserSCT3_RawData(rdoID, stripEntry.second, std::vector<int> {}));
+    }
+
+    outRdoContainer->getWriteHandle(waferHash).addOrDelete(std::move(current_collection));
+  }
+
+  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;
+}
+
+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 000000000..57900ef49
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
@@ -0,0 +1,34 @@
+#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;
+  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::WriteHandleKey<FaserSCT_RDO_Container> m_outRdoKey{this, "OutputObjectName", "SCT_RDOs", "Output 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 000000000..2bb95c96c
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
@@ -0,0 +1,186 @@
+#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 radius {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
+        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();
+    radius = sqrt(pow(upstreamParameters->position().x(), 2) + pow(upstreamParameters->position().y(), 2));
+    theTrack = track;
+  }
+
+  if ((nTracks != 1) || (nLongTracks != 1) || (momentum < m_minMomentum) || (momentum > m_maxMomentum) || (chi2PerDoF > m_maxChi2PerDoF) || ((charge>0) != m_acceptPositive) || (radius > 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 000000000..4c3bee533
--- /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", 1000.0, "Maximum momentum in GeV to accept track"};
+  DoubleProperty    m_maxRadius {this, "MaxRadiusMm", 85.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", 52.5, "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 000000000..034c80451
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx
@@ -0,0 +1,5 @@
+#include "../OverlayRDOAlg.h"
+#include "../SelectRDOAlg.h"
+
+DECLARE_COMPONENT(OverlayRDOAlg)
+DECLARE_COMPONENT(SelectRDOAlg)
\ No newline at end of file
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
index f990be945..6b02bc244 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
@@ -157,10 +157,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_WARNING ("SCT error state is locked.");
+          }
       }
       continue;
     }
-- 
GitLab


From 38a49da85b3aaf8b7454dad5bcab5764b2dbe72c Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Fri, 3 Mar 2023 08:43:02 -0800
Subject: [PATCH 10/21] Carry through EDGEMODE RDOs as well

---
 .../OverlayRDO/python/OverlayRDOConfig.py     |  1 +
 .../OverlayRDO/python/SelectRDOConfig.py      |  2 +
 .../OverlayRDO/src/OverlayRDOAlg.cxx          | 65 ++++++++++++++-----
 .../OverlayRDO/src/OverlayRDOAlg.h            |  6 +-
 4 files changed, 56 insertions(+), 18 deletions(-)

diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
index a9f643e3c..bce33065b 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
@@ -20,6 +20,7 @@ def OverlayRDOAlgCfg(flags, **kwargs):
     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))
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
index 95b7b692d..42b711789 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
@@ -23,6 +23,7 @@ def SelectRDOAlgCfg(flags, **kwargs):
     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,"RDO", ItemList=ItemList, disableEventTag=True))
@@ -81,6 +82,7 @@ if __name__ == "__main__":
         # "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"
     ]))
 
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
index 0384d5d03..f73f8fd3e 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
@@ -23,7 +23,10 @@ StatusCode OverlayRDOAlg::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;
@@ -41,6 +44,12 @@ StatusCode OverlayRDOAlg::execute(const EventContext &ctx) const
   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);
@@ -49,6 +58,43 @@ StatusCode OverlayRDOAlg::execute(const EventContext &ctx) const
   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;
@@ -80,30 +126,15 @@ StatusCode OverlayRDOAlg::execute(const EventContext &ctx) const
     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));
   }
-
-  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;
 }
 
+
 Trk::Track* OverlayRDOAlg::cloneTrack(const Trk::Track* theTrack) const
 {
   Trk::TrackInfo i { theTrack->info() };
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
index 57900ef49..eba2f90d7 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
@@ -18,14 +18,18 @@ public:
 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, "OutputObjectName", "SCT_RDOs", "Output Object name"};
+  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;
-- 
GitLab


From 0c18fcdb969987a72ba1ab601225e238aeaebc15 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Fri, 3 Mar 2023 10:59:58 -0800
Subject: [PATCH 11/21] Multiple output streams and binning in (x,y)

---
 .../TrackerRecAlgs/OverlayRDO/CMakeLists.txt  |  4 +-
 .../OverlayRDO/python/SelectRDOConfig.py      | 50 +++++++++++++++--
 .../OverlayRDO/src/BinRDOAlg.cxx              | 55 +++++++++++++++++++
 .../TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.h | 23 ++++++++
 .../OverlayRDO/src/SelectRDOAlg.cxx           | 25 +++++----
 .../OverlayRDO/src/SelectRDOAlg.h             |  8 +--
 .../src/component/OverlayRDO_entries.cxx      |  4 +-
 7 files changed, 146 insertions(+), 23 deletions(-)
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.cxx
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.h

diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt b/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
index d1213ed93..75394179e 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
@@ -2,8 +2,8 @@ atlas_subdir(OverlayRDO)
 
 atlas_add_component(
         OverlayRDO
-        src/OverlayRDOAlg.h src/SelectRDOAlg.h
-        src/OverlayRDOAlg.cxx src/SelectRDOAlg.cxx
+        src/OverlayRDOAlg.h src/SelectRDOAlg.h src/BinRDOAlg.h 
+        src/OverlayRDOAlg.cxx src/SelectRDOAlg.cxx src/BinRDOAlg.cxx
         src/component/OverlayRDO_entries.cxx
         LINK_LIBRARIES AthenaBaseComps StoreGateLib TrackerRawData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack xAODFaserWaveform
 )
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
index 42b711789..29534a9b3 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
@@ -15,6 +15,16 @@ def SelectRDOAlgCfg(flags, **kwargs):
     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"]
@@ -26,9 +36,41 @@ def SelectRDOAlgCfg(flags, **kwargs):
     ItemList += ["FaserSCT_RDO_Container#" + chargePrefix + "SCT_EDGEMODE_RDOs"]
     # ItemList += ["Tracker::FaserSCT_ClusterContainer#" + chargePrefix + "SCT_ClusterContainer"]
 
-    acc.merge(OutputStreamCfg(ConfigFlags,"RDO", ItemList=ItemList, disableEventTag=True))
-    osrdo = acc.getEventAlgo("OutputStreamRDO")
-    osrdo.AcceptAlgs += ["SelectRDOAlg"]
+    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
 
@@ -55,7 +97,7 @@ if __name__ == "__main__":
                               ]
     # 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.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
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.cxx
new file mode 100644
index 000000000..a917447dd
--- /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 000000000..f8e86e844
--- /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/SelectRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
index 2bb95c96c..b33bed9d1 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
@@ -77,7 +77,7 @@ StatusCode SelectRDOAlg::execute(const EventContext &ctx) const
   int nLongTracks {0};
   double chi2PerDoF {0};
   double charge {0};
-  double radius {0};
+  double maxRadius {0};
   double momentum {0};
   const Trk::Track* theTrack {nullptr};
 
@@ -110,30 +110,31 @@ StatusCode SelectRDOAlg::execute(const EventContext &ctx) const
     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();
-    radius = sqrt(pow(upstreamParameters->position().x(), 2) + pow(upstreamParameters->position().y(), 2));
     theTrack = track;
   }
 
-  if ((nTracks != 1) || (nLongTracks != 1) || (momentum < m_minMomentum) || (momentum > m_maxMomentum) || (chi2PerDoF > m_maxChi2PerDoF) || ((charge>0) != m_acceptPositive) || (radius > m_maxRadius))
+  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;
+  // 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
 
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
index 4c3bee533..842282d53 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
@@ -31,12 +31,12 @@ private:
   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", 1000.0, "Maximum momentum in GeV to accept track"};
-  DoubleProperty    m_maxRadius {this, "MaxRadiusMm", 85.0, "Maximum radius at first measurement 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", 52.5, "Maximum charge in pC recorded by upper or lower timing scintillator"};
+  // 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
index 034c80451..93a111b85 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx
@@ -1,5 +1,7 @@
 #include "../OverlayRDOAlg.h"
 #include "../SelectRDOAlg.h"
+#include "../BinRDOAlg.h"
 
 DECLARE_COMPONENT(OverlayRDOAlg)
-DECLARE_COMPONENT(SelectRDOAlg)
\ No newline at end of file
+DECLARE_COMPONENT(SelectRDOAlg)
+DECLARE_COMPONENT(BinRDOAlg)
\ No newline at end of file
-- 
GitLab


From 9b2a3c3a02a18e6dd72ebbd404f2fd743a05821f Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Fri, 3 Mar 2023 22:39:23 +0100
Subject: [PATCH 12/21] Add multiple output streams and (x,y) binning

---
 .../TrackerRecAlgs/OverlayRDO/CMakeLists.txt  |  4 +-
 .../OverlayRDO/python/OverlayRDOConfig.py     |  1 +
 .../OverlayRDO/python/SelectRDOConfig.py      | 52 +++++++++++++--
 .../OverlayRDO/src/BinRDOAlg.cxx              | 55 ++++++++++++++++
 .../TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.h | 23 +++++++
 .../OverlayRDO/src/OverlayRDOAlg.cxx          | 65 ++++++++++++++-----
 .../OverlayRDO/src/OverlayRDOAlg.h            |  6 +-
 .../OverlayRDO/src/SelectRDOAlg.cxx           | 25 +++----
 .../OverlayRDO/src/SelectRDOAlg.h             |  8 +--
 .../src/component/OverlayRDO_entries.cxx      |  4 +-
 10 files changed, 202 insertions(+), 41 deletions(-)
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.cxx
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.h

diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt b/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
index d1213ed93..75394179e 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
@@ -2,8 +2,8 @@ atlas_subdir(OverlayRDO)
 
 atlas_add_component(
         OverlayRDO
-        src/OverlayRDOAlg.h src/SelectRDOAlg.h
-        src/OverlayRDOAlg.cxx src/SelectRDOAlg.cxx
+        src/OverlayRDOAlg.h src/SelectRDOAlg.h src/BinRDOAlg.h 
+        src/OverlayRDOAlg.cxx src/SelectRDOAlg.cxx src/BinRDOAlg.cxx
         src/component/OverlayRDO_entries.cxx
         LINK_LIBRARIES AthenaBaseComps StoreGateLib TrackerRawData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack xAODFaserWaveform
 )
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
index a9f643e3c..bce33065b 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/OverlayRDOConfig.py
@@ -20,6 +20,7 @@ def OverlayRDOAlgCfg(flags, **kwargs):
     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))
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
index 95b7b692d..29534a9b3 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/SelectRDOConfig.py
@@ -15,6 +15,16 @@ def SelectRDOAlgCfg(flags, **kwargs):
     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"]
@@ -23,11 +33,44 @@ def SelectRDOAlgCfg(flags, **kwargs):
     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,"RDO", ItemList=ItemList, disableEventTag=True))
-    osrdo = acc.getEventAlgo("OutputStreamRDO")
-    osrdo.AcceptAlgs += ["SelectRDOAlg"]
+    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
 
@@ -54,7 +97,7 @@ if __name__ == "__main__":
                               ]
     # 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.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
@@ -81,6 +124,7 @@ if __name__ == "__main__":
         # "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"
     ]))
 
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/BinRDOAlg.cxx
new file mode 100644
index 000000000..a917447dd
--- /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 000000000..f8e86e844
--- /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/OverlayRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
index 0384d5d03..f73f8fd3e 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
@@ -23,7 +23,10 @@ StatusCode OverlayRDOAlg::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;
@@ -41,6 +44,12 @@ StatusCode OverlayRDOAlg::execute(const EventContext &ctx) const
   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);
@@ -49,6 +58,43 @@ StatusCode OverlayRDOAlg::execute(const EventContext &ctx) const
   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;
@@ -80,30 +126,15 @@ StatusCode OverlayRDOAlg::execute(const EventContext &ctx) const
     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));
   }
-
-  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;
 }
 
+
 Trk::Track* OverlayRDOAlg::cloneTrack(const Trk::Track* theTrack) const
 {
   Trk::TrackInfo i { theTrack->info() };
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
index 57900ef49..eba2f90d7 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.h
@@ -18,14 +18,18 @@ public:
 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, "OutputObjectName", "SCT_RDOs", "Output Object name"};
+  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
index 2bb95c96c..b33bed9d1 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.cxx
@@ -77,7 +77,7 @@ StatusCode SelectRDOAlg::execute(const EventContext &ctx) const
   int nLongTracks {0};
   double chi2PerDoF {0};
   double charge {0};
-  double radius {0};
+  double maxRadius {0};
   double momentum {0};
   const Trk::Track* theTrack {nullptr};
 
@@ -110,30 +110,31 @@ StatusCode SelectRDOAlg::execute(const EventContext &ctx) const
     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();
-    radius = sqrt(pow(upstreamParameters->position().x(), 2) + pow(upstreamParameters->position().y(), 2));
     theTrack = track;
   }
 
-  if ((nTracks != 1) || (nLongTracks != 1) || (momentum < m_minMomentum) || (momentum > m_maxMomentum) || (chi2PerDoF > m_maxChi2PerDoF) || ((charge>0) != m_acceptPositive) || (radius > m_maxRadius))
+  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;
+  // 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
 
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
index 4c3bee533..842282d53 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/SelectRDOAlg.h
@@ -31,12 +31,12 @@ private:
   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", 1000.0, "Maximum momentum in GeV to accept track"};
-  DoubleProperty    m_maxRadius {this, "MaxRadiusMm", 85.0, "Maximum radius at first measurement 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", 52.5, "Maximum charge in pC recorded by upper or lower timing scintillator"};
+  // 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
index 034c80451..330b7e130 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx
@@ -1,5 +1,7 @@
 #include "../OverlayRDOAlg.h"
 #include "../SelectRDOAlg.h"
+#include "../BinRDOAlg.h"
 
 DECLARE_COMPONENT(OverlayRDOAlg)
-DECLARE_COMPONENT(SelectRDOAlg)
\ No newline at end of file
+DECLARE_COMPONENT(SelectRDOAlg)
+DECLARE_COMPONENT(BinRDOAlg)
-- 
GitLab


From 306895e643bcc424df6668fc42cbe470c8bc9315 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Sun, 5 Mar 2023 09:36:23 -0800
Subject: [PATCH 13/21] Suppress unneeded histogram file

---
 .../TrackerSegmentFit/python/TrackerSegmentFitConfig.py     | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py b/Tracker/TrackerRecAlgs/TrackerSegmentFit/python/TrackerSegmentFitConfig.py
index f0de1582f..eeaadd061 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
-- 
GitLab


From e76390c6b063becd88834d21dc5b0de8ac6506f4 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Mon, 6 Mar 2023 10:29:03 -0800
Subject: [PATCH 14/21] Add tree-filling algorithm

---
 .../TrackerRecAlgs/OverlayRDO/CMakeLists.txt  |   4 +-
 .../OverlayRDO/python/HistoRDOConfig.py       | 128 ++++++++++++++++++
 .../OverlayRDO/src/HistoRDOAlg.cxx            | 108 +++++++++++++++
 .../OverlayRDO/src/HistoRDOAlg.h              |  48 +++++++
 .../src/component/OverlayRDO_entries.cxx      |   2 +
 5 files changed, 288 insertions(+), 2 deletions(-)
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx
 create mode 100644 Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.h

diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt b/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
index 75394179e..8d8665ba9 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/CMakeLists.txt
@@ -2,8 +2,8 @@ atlas_subdir(OverlayRDO)
 
 atlas_add_component(
         OverlayRDO
-        src/OverlayRDOAlg.h src/SelectRDOAlg.h src/BinRDOAlg.h 
-        src/OverlayRDOAlg.cxx src/SelectRDOAlg.cxx src/BinRDOAlg.cxx
+        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
 )
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py
new file mode 100644
index 000000000..8cd2109f2
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py
@@ -0,0 +1,128 @@
+"""
+    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)
+    acc = ComponentAccumulator()
+    HistoRDOAlg = CompFactory.HistoRDOAlg("HistoRDOAlg",**kwargs)
+    acc.addEventAlgo(HistoRDOAlg)
+
+    # 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/src/HistoRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx
new file mode 100644
index 000000000..d7c954aeb
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx
@@ -0,0 +1,108 @@
+#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("charge1", &m_charge1);
+
+    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("charge2", &m_charge2);
+
+    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_charge1.clear();
+
+    m_longTracks2 = 0;
+    m_x2.clear();
+    m_y2.clear();
+    m_p2.clear();
+    m_theta2.clear();
+    m_charge2.clear();
+
+    SG::ReadHandle<TrackCollection> trackCollection1 {m_trackCollection1, ctx};
+    // ATH_CHECK(trackCollection1.isValid());
+    SG::ReadHandle<TrackCollection> trackCollection2 {m_trackCollection2, ctx};
+    // ATH_CHECK(trackCollection2.isValid());
+
+    if (trackCollection1.isValid())
+    {
+        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;
+            }
+            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()));
+        }
+    }
+
+    if (trackCollection2.isValid())
+    {
+        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;
+            }
+            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_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 000000000..674b07031
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.h
@@ -0,0 +1,48 @@
+#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_x1;
+  mutable std::vector<double> m_y1;
+  mutable std::vector<double> m_charge1;
+
+  mutable int m_longTracks2;
+  mutable std::vector<double> m_p2;
+  mutable std::vector<double> m_theta2;
+  mutable std::vector<double> m_x2;
+  mutable std::vector<double> m_y2;
+  mutable std::vector<double> m_charge2;
+
+
+};
+
+inline const ServiceHandle <ITHistSvc> &HistoRDOAlg::histSvc() const {
+  return m_histSvc;
+}
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx
index 330b7e130..13e61cb84 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/component/OverlayRDO_entries.cxx
@@ -1,7 +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)
-- 
GitLab


From f6a7d3e4d81ec64b7a4ad9930b4cc5d5837640c0 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Tue, 7 Mar 2023 19:45:04 +0100
Subject: [PATCH 15/21] Handle group size > 1

---
 Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
index f73f8fd3e..9b88a96f9 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/OverlayRDOAlg.cxx
@@ -115,11 +115,12 @@ OverlayRDOAlg::processDB(const std::map<IdentifierHash, std::vector<const FaserS
         // const FaserSCT3_RawData* data = dynamic_cast<const FaserSCT3_RawData*>(rawData);
         Identifier stripID = rawData->identify();
         int stripNumber = m_sctHelper->strip(stripID);
-        // int groupSize = rawData->getGroupSize();
+        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;
-        stripMap[stripNumber] |= dataWord;
+	for (int i = stripNumber; i < stripNumber + groupSize; i++)
+	  stripMap[i] |= dataWord;
       }
     }
 
-- 
GitLab


From 2280ae4143bbb380ff2dfad307e48ddd84e00bc6 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Wed, 8 Mar 2023 11:12:44 -0800
Subject: [PATCH 16/21] Add matching quality information

---
 .../Reconstruction/scripts/faser_reco.py      |   2 +-
 .../OverlayRDO/python/HistoRDOConfig.py       |  19 ++-
 .../OverlayRDO/src/HistoRDOAlg.cxx            | 139 ++++++++++++++++++
 .../OverlayRDO/src/HistoRDOAlg.h              |  10 ++
 4 files changed, 159 insertions(+), 11 deletions(-)

diff --git a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
index cc608ff03..8cf58dab0 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
+++ b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
@@ -253,7 +253,7 @@ itemList = [ "xAOD::EventInfo#*"
 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))
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py b/Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py
index 8cd2109f2..cfb2dfefe 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/python/HistoRDOConfig.py
@@ -8,12 +8,11 @@ from AthenaConfiguration.ComponentFactory import CompFactory
 
 def HistoRDOAlgCfg(flags, **kwargs):
 
-    # from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
-    # acc = FaserGeometryCfg(flags)
-    acc = ComponentAccumulator()
+    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."]
@@ -83,7 +82,7 @@ if __name__ == "__main__":
     ConfigFlags.GeoModel.Align.Dynamic = False
     ConfigFlags.Beam.NumberOfCollisions = 0.
 
-    # ConfigFlags.Detector.GeometryFaserSCT = True
+    ConfigFlags.Detector.GeometryFaserSCT = True
 
     ConfigFlags.lock()
 
@@ -101,11 +100,11 @@ if __name__ == "__main__":
     # ]))
 
     # 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
+    replicaSvc = acc.getService("DBReplicaSvc")
+    replicaSvc.COOLSQLiteVetoPattern = ""
+    replicaSvc.UseCOOLSQLite = True
+    replicaSvc.UseCOOLFrontier = False
+    replicaSvc.UseGeomSQLite = True
 
     # Timing
     #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx
index d7c954aeb..041b4c17b 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx
@@ -18,14 +18,25 @@ StatusCode HistoRDOAlg::initialize()
     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));
 
@@ -43,22 +54,35 @@ StatusCode HistoRDOAlg::execute(const EventContext& ctx) const
     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;
@@ -68,17 +92,26 @@ StatusCode HistoRDOAlg::execute(const EventContext& ctx) const
             {
                 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;
@@ -88,15 +121,121 @@ StatusCode HistoRDOAlg::execute(const EventContext& ctx) const
             {
                 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_ALWAYS("Matched (1): " << singleState->position().z() << " and " << overlayState->position().z());
+            }
+            else if (singleState->position().z() < overlayState->position().z())
+            {
+                ATH_MSG_ALWAYS("z not matched (1): " << singleState->position().z() << " vs " << overlayState->position().z());
+                iSingle++;
+            }
+            else if (overlayState->position().z() < singleState->position().z())
+            {
+                ATH_MSG_ALWAYS("z not matched (1): " << singleState->position().z() << " vs " << overlayState->position().z());
+                iOverlay++;
+            }
+        }
+        double matchFract1 { ((double) nZMatch)/singleParameters->size() };
+        ATH_MSG_ALWAYS("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_ALWAYS("Matched (2): " << singleState->position().z() << " and " << overlayState->position().z());
+            }
+            else if (singleState->position().z() < overlayState->position().z())
+            {
+                ATH_MSG_ALWAYS("z not matched (2): " << singleState->position().z() << " vs " << overlayState->position().z());
+                iSingle++;
+            }
+            else if (overlayState->position().z() < singleState->position().z())
+            {
+                ATH_MSG_ALWAYS("z not matched (2): " << singleState->position().z() << " vs " << overlayState->position().z());
+                iOverlay++;
+            }
+        }
+        double matchFract2 { ((double) nZMatch)/singleParameters->size() };
+        ATH_MSG_ALWAYS("MatchFract2: " << matchFract2);
+        m_matchFraction = std::min(matchFract1, matchFract2);
+        ATH_MSG_ALWAYS("MatchFraction: " << m_matchFraction);
+    }
 
     m_tree->Fill();
     return StatusCode::SUCCESS;
diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.h b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.h
index 674b07031..4fe18f7ad 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.h
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.h
@@ -29,17 +29,27 @@ private:
   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;
 
 };
 
-- 
GitLab


From cde6702be6cd4ca063f0610cd9b75a56d89fdf10 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Wed, 8 Mar 2023 11:40:16 -0800
Subject: [PATCH 17/21] Suppress debugging printout

---
 .../OverlayRDO/src/HistoRDOAlg.cxx             | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx
index 041b4c17b..78096c600 100644
--- a/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx
+++ b/Tracker/TrackerRecAlgs/OverlayRDO/src/HistoRDOAlg.cxx
@@ -188,21 +188,21 @@ StatusCode HistoRDOAlg::execute(const EventContext& ctx) const
                 m_yAgreement += pow(singleState->position().y() - overlayState->position().y(), 2)/((*singleCov)(1,1) + (*overlayCov)(1,1));
                 iSingle++;
                 iOverlay++;
-                ATH_MSG_ALWAYS("Matched (1): " << singleState->position().z() << " and " << overlayState->position().z());
+                ATH_MSG_VERBOSE("Matched (1): " << singleState->position().z() << " and " << overlayState->position().z());
             }
             else if (singleState->position().z() < overlayState->position().z())
             {
-                ATH_MSG_ALWAYS("z not matched (1): " << singleState->position().z() << " vs " << 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_ALWAYS("z not matched (1): " << singleState->position().z() << " vs " << overlayState->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_ALWAYS("MatchFract1: " << matchFract1);
+        ATH_MSG_VERBOSE("MatchFract1: " << matchFract1);
         nZMatch = 0;
         singleParameters = secondSingle->trackParameters();
         overlayParameters = secondOverlay->trackParameters();
@@ -218,23 +218,23 @@ StatusCode HistoRDOAlg::execute(const EventContext& ctx) const
                 m_yAgreement += pow(singleState->position().y() - overlayState->position().y(), 2)/((*singleCov)(1,1) + (*overlayCov)(1,1));
                 iSingle++;
                 iOverlay++;
-                ATH_MSG_ALWAYS("Matched (2): " << singleState->position().z() << " and " << overlayState->position().z());
+                ATH_MSG_VERBOSE("Matched (2): " << singleState->position().z() << " and " << overlayState->position().z());
             }
             else if (singleState->position().z() < overlayState->position().z())
             {
-                ATH_MSG_ALWAYS("z not matched (2): " << singleState->position().z() << " vs " << 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_ALWAYS("z not matched (2): " << singleState->position().z() << " vs " << overlayState->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_ALWAYS("MatchFract2: " << matchFract2);
+        ATH_MSG_VERBOSE("MatchFract2: " << matchFract2);
         m_matchFraction = std::min(matchFract1, matchFract2);
-        ATH_MSG_ALWAYS("MatchFraction: " << m_matchFraction);
+        ATH_MSG_VERBOSE("MatchFraction: " << m_matchFraction);
     }
 
     m_tree->Fill();
-- 
GitLab


From c99e2abe0f28f7c015992f7a30c3e29308238b2d Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Sat, 11 Mar 2023 08:56:16 -0800
Subject: [PATCH 18/21] Accept more track seeds

---
 .../Reconstruction/scripts/faser_reco.py      |  7 ++++
 .../Acts/FaserActsKalmanFilter/src/CKF2.cxx   | 32 +++++++++++++++++++
 .../Acts/FaserActsKalmanFilter/src/CKF2.h     |  1 +
 .../src/CircleFitTrackSeedTool.cxx            |  7 +++-
 4 files changed, 46 insertions(+), 1 deletion(-)

diff --git a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
index 8cf58dab0..a778d2edd 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',
@@ -83,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
 
@@ -105,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
@@ -203,6 +209,7 @@ if not args.isOverlay:
 
 # 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"))
 
 #
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
index 6b02bc244..2b94744be 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();
 
@@ -192,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();
@@ -233,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 170f03539..cd78bf725 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 952f75287..56960d458 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
@@ -142,8 +142,11 @@ StatusCode CircleFitTrackSeedTool::run(std::vector<int> maskedLayers) {
     Seed selected = allSeeds.front();
     selectedSeeds.push_back(selected);
     allSeeds.remove_if([&](const Seed &p) {
-      return (p.size < 10) || ((p.clusterSet & selected.clusterSet).count() > 6);
+      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);
+    // });
   }
 
   Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero();
@@ -161,7 +164,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));
   }
 
-- 
GitLab


From cc6ec485b2d0f77810c40d3cb8ab06aa88f72c9e Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Fri, 17 Mar 2023 11:42:16 -0700
Subject: [PATCH 19/21] Suppress warning spam

---
 .../FaserSiSpacePointTool/src/TrackerSpacePointMakerTool.cxx    | 2 +-
 Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx                | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/Tracker/TrackerRecTools/FaserSiSpacePointTool/src/TrackerSpacePointMakerTool.cxx b/Tracker/TrackerRecTools/FaserSiSpacePointTool/src/TrackerSpacePointMakerTool.cxx
index 814fc31a5..c8a7df001 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 2b94744be..14616c6f4 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
@@ -172,7 +172,7 @@ StatusCode CKF2::execute() {
           }
           catch (...)
           {
-              ATH_MSG_WARNING ("SCT error state is locked.");
+              ATH_MSG_DEBUG ("SCT error state is locked.");
           }
       }
       continue;
-- 
GitLab


From 9c7ce3ee2e4857cc18cfdc5591f5d2e4d5c2bb52 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Fri, 17 Mar 2023 12:35:35 -0700
Subject: [PATCH 20/21] Attempt at improved seed selection

---
 .../src/CircleFitTrackSeedTool.cxx            | 23 +++++++++++++------
 .../src/CircleFitTrackSeedTool.h              |  2 +-
 2 files changed, 17 insertions(+), 8 deletions(-)

diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CircleFitTrackSeedTool.cxx
index 56960d458..bbae70e3b 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,12 +147,12 @@ 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);
+    //   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);
+    });
   }
 
   Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero();
@@ -269,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 0467c9262..a4d53c629 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:
-- 
GitLab


From 03acb2677185f980298fa06c69024ce15ada901f Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Wed, 19 Apr 2023 13:40:44 +0200
Subject: [PATCH 21/21] Incorporate Ke's fix to refit

---
 .../FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx     | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx
index 4928fc917..83be3bbac 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx
@@ -61,7 +61,14 @@ CreateTrkTrackTool::createTrack(const Acts::GeometryContext &gctx, const FaserAc
           };
         }
       }
-      const Trk::TrackStateOnSurface *perState = new Trk::TrackStateOnSurface(clusterOnTrack, parm);
+      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);
       }
-- 
GitLab