Merge CKF changes

Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_emi_100GeV-101104.json

 {
-    "file_length": 100,
+    "file_length": 200,
     "mass": 0.511,
     "maxE": 100.0,
     "minE": 100.0,

Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_epl_100GeV-101103.json

 {
-    "file_length": 100,
+    "file_length": 200,
     "mass": 0.511,
     "maxE": 100.0,
     "minE": 100.0,

Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_gam_100GeV-102201.json

 {
-    "file_length": 100,
+    "file_length": 200,
     "mass": 0.0,
     "maxE": 100.0,
     "minE": 100.0,

Generators/DIFGenerator/python/DIFSampler.py

@@ -5,7 +5,7 @@
 import ParticleGun as PG
 from math import sqrt, sin, cos, acos
 from random import uniform
-from AthenaCommon.SystemOfUnits import TeV, MeV
+from AthenaCommon.SystemOfUnits import TeV, GeV, MeV
 from AthenaCommon.PhysicalConstants import pi
 from ROOT import TLorentzVector
 
@@ -90,8 +90,12 @@ class DIFSampler(PG.ParticleSampler):
             else:
                 self.mom = None
 
-    def __init__(self, daughter1_pid = 11, daughter2_pid = -11, mother_pid = None, mother_mom = PG.EThetaMPhiSampler(sqrt((1*TeV)**2 + (10*MeV)**2),0,10*MeV,0),mother_pos = CylinderSampler([0, 100**2],[0, 2*pi],[-1500, 0],0)):
-        self._mother_sampler = PG.ParticleSampler(pid = mother_pid, mom = mother_mom, pos = mother_pos)
+    # def __init__(self, daughter1_pid = 11, daughter2_pid = -11, mother_pid = None, mother_mom = PG.EThetaMPhiSampler(sqrt((1*TeV)**2 + (10*MeV)**2),0,10*MeV,0),mother_pos = CylinderSampler([0, 100**2],[0, 2*pi],[-1500, 0],0)):
+    def __init__(self, daughter1_pid=13, daughter2_pid=-13, mother_pid=None,
+                 mother_mom=PG.EThetaMPhiSampler(sqrt((1*TeV)**2 + (500*MeV)**2), [0, 0.0002], 500*MeV, [0, 2*pi]),
+                 my_z_position=-1500):
+                 # mother_pos=CylinderSampler([0, 100**2], [0, 2*pi], [-1500, 0], 0)):
+        self._mother_sampler = PG.ParticleSampler(pid = mother_pid, mom = mother_mom, pos=CylinderSampler([0, 100**2], [0, 2*pi], my_z_position, 0))
 
         self.daughter1 = self.particle(daughter1_pid)
         self.daughter2 = self.particle(daughter2_pid)

Generators/FaserParticleGun/python/FaserParticleGunConfig.py

@@ -190,7 +190,8 @@ def FaserParticleGunForeseeCfg(ConfigFlags, **kwargs) :
     return cfg
 
 def FaserParticleGunCfg(ConfigFlags) :
-    generator = ConfigFlags.Sim.Gun.setdefault("Generator", "SingleParticle")
+    # generator = ConfigFlags.Sim.Gun.setdefault("Generator", "SingleParticle")
+    generator = ConfigFlags.Sim.Gun.setdefault("Generator", "DecayInFlight")
     kwargs = ConfigFlags.Sim.Gun
     if generator == "SingleEcalParticle" :
         return FaserParticleGunSingleEcalParticleCfg(ConfigFlags, **kwargs)

Tracker/TrackerRecAlgs/FaserSpacePoints/CMakeLists.txt

-################################################################################
-# Package: FaserSpacePoints
-################################################################################
+atlas_subdir(FaserSpacePoints)
 
-# Declare the package name:
-atlas_subdir( FaserSpacePoints )
+atlas_add_component(
+  FaserSpacePoints
+  src/components/*.cxx src/*.cxx src/*.h
+  LINK_LIBRARIES AthenaBaseComps TrackerSimData TrackerSimEvent TrackerSpacePoint TrackerPrepRawData)
 
-atlas_add_component( FaserSpacePoints
-                     src/components/*.cxx src/*.cxx src/*.h
-	                   LINK_LIBRARIES AthenaBaseComps TrackerSpacePoint TrackerPrepRawData)
-
-atlas_install_python_modules( python/*.py )
\ No newline at end of file
+atlas_install_python_modules( python/*.py )
+atlas_install_scripts( test/*.py )

Tracker/TrackerRecAlgs/FaserSpacePoints/python/FaserSpacePointsConfig.py

+"""
+Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+"""
+
+from AthenaConfiguration.ComponentFactory import CompFactory
+from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
+
+
+def FaserSpacePointsCfg(flags, **kwargs):
+    acc = FaserSCT_GeometryCfg(flags)
+    Tracker__FaserSpacePoints = CompFactory.Tracker.FaserSpacePoints
+    acc.addEventAlgo(Tracker__FaserSpacePoints(**kwargs))
+
+    thistSvc = CompFactory.THistSvc()
+    thistSvc.Output += ["HIST DATAFILE='SpacePoints.root' OPT='RECREATE'"]
+    acc.addService(thistSvc)
+    return acc

Tracker/TrackerRecAlgs/FaserSpacePoints/python/FaserSpacePointsCosmics.py

@@ -22,7 +22,6 @@ Tracker__TruthSeededTrackFinder, THistSvc=CompFactory.getComps("Tracker::FaserSp
 def TruthSeededTrackFinderBasicCfg(flags, **kwargs):
     """Return ComponentAccumulator for TruthSeededTrackFinder"""
     acc = FaserSCT_GeometryCfg(flags)
-    kwargs.setdefault("SpacePointsSCTName", "SCT_SpacePointContainer")
     acc.addEventAlgo(Tracker__TruthSeededTrackFinder(**kwargs))
    # attach ToolHandles
     return acc

Tracker/TrackerRecAlgs/FaserSpacePoints/src/FaserSpacePoints.cxx

 #include "FaserSpacePoints.h"
 #include "Identifier/Identifier.h"
+#include "TrackerSimEvent/FaserSiHit.h"
 #include "TrackerSpacePoint/FaserSCT_SpacePoint.h"
 #include "TrackerPrepRawData/FaserSCT_ClusterCollection.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrackerReadoutGeometry/SCT_DetectorManager.h"
+
 
 namespace Tracker {
-  FaserSpacePoints::FaserSpacePoints(const std::string& name, ISvcLocator* pSvcLocator)
-    : AthReentrantAlgorithm(name, pSvcLocator) { }
 
-  
-  StatusCode FaserSpacePoints::initialize() {
-    ATH_MSG_INFO(name() << "::" << __FUNCTION__);
+FaserSpacePoints::FaserSpacePoints(const std::string& name, ISvcLocator* pSvcLocator)
+  : AthReentrantAlgorithm(name, pSvcLocator), AthHistogramming(name),
+    m_histSvc("THistSvc/THistSvc", name) {}
+
+
+StatusCode FaserSpacePoints::initialize() {
+  ATH_CHECK(m_spacePointContainerKey.initialize());
+  ATH_CHECK(m_siHitCollectionKey.initialize());
+  ATH_CHECK(m_simDataCollectionKey.initialize());
+  ATH_CHECK(detStore()->retrieve(m_idHelper, "FaserSCT_ID"));
+  ATH_CHECK(detStore()->retrieve(m_detMgr, "SCT"));
+  m_tree = new TTree("tree", "tree");
+  m_tree->Branch("run", &m_run, "run/I");
+  m_tree->Branch("event", &m_event, "event/I");
+  m_tree->Branch("station", &m_station, "station/I");
+  m_tree->Branch("layer", &m_layer, "layer/I");
+  m_tree->Branch("phi", &m_phi, "phi/I");
+  m_tree->Branch("eta", &m_eta, "eta/I");
+  m_tree->Branch("side", &m_side, "side/I");
+  m_tree->Branch("x", &m_x, "x/D");
+  m_tree->Branch("y", &m_y, "y/D");
+  m_tree->Branch("z", &m_z, "z/D");
+  m_tree->Branch("muon_hit", &m_muon_hit, "muon_hit/B");
+  m_tree->Branch("tx", &m_tx, "tx/D");
+  m_tree->Branch("ty", &m_ty, "ty/D");
+  m_tree->Branch("tz", &m_tz, "tz/D");
+  ATH_CHECK(histSvc()->regTree("/HIST/spacePoints", m_tree));
+  return StatusCode::SUCCESS;
+}
 
-    if ( m_sct_spcontainer.key().empty()) {
-      ATH_MSG_FATAL( "SCTs selected and no name set for SCT clusters");
-      return StatusCode::FAILURE;
-    }
-    ATH_CHECK( m_sct_spcontainer.initialize() );
 
-    ATH_CHECK(detStore()->retrieve(m_idHelper, "FaserSCT_ID"));
+StatusCode FaserSpacePoints::finalize() {
+  return StatusCode::SUCCESS;
+}
 
-    std::string filePath = m_filePath;
-    m_outputFile = TFile::Open(filePath.c_str(), "RECREATE");
-    if (m_outputFile == nullptr) {
-      ATH_MSG_ERROR("Unable to open output file at " << m_filePath);
-      return StatusCode::FAILURE;
-    }
-    m_outputFile->cd();
+bool FaserSpacePoints::muon_hit(const TrackerSimDataCollection* simDataMap, Identifier id) const {
+  // check if Identifier in simDataMap
+  if (!simDataMap->count(id)) {
+    return false;
+  }
+  // check if hit originates from a muon (if the pdg code of the hit is -13)
+  TrackerSimData simData = simDataMap->at(id);
+  const std::vector<TrackerSimData::Deposit>& deposits = simData.getdeposits();
+  auto it = std::find_if(
+      deposits.begin(), deposits.end(), [](TrackerSimData::Deposit deposit) {
+        return deposit.first->pdg_id() == -13;
+      });
+  return it != deposits.end();
+}
 
-    std::string treeName = m_treeName;
-    m_outputTree = new TTree(treeName.c_str(), "tree");
-    if (m_outputTree == nullptr) {
-      ATH_MSG_ERROR("Unable to create TTree");
-      return StatusCode::FAILURE;
-    }
+StatusCode FaserSpacePoints::execute(const EventContext& ctx) const {
+  m_run = ctx.eventID().run_number();
+  m_event = ctx.eventID().event_number();
 
-    FaserSpacePoints::initializeTree();
-    return StatusCode::SUCCESS;
-  }
+  SG::ReadHandle<FaserSCT_SpacePointContainer> spacePointContainer {m_spacePointContainerKey, ctx};
+  ATH_CHECK(spacePointContainer.isValid());
 
-  
-  StatusCode FaserSpacePoints::finalize() {
-    ATH_MSG_INFO(name() << "::" << __FUNCTION__);
-    m_outputFile->cd();
-    m_outputTree->Write();
-    return StatusCode::SUCCESS;
-  }
+  SG::ReadHandle<FaserSiHitCollection> siHitCollection {m_siHitCollectionKey, ctx};
+  ATH_CHECK(siHitCollection.isValid());
+
+  SG::ReadHandle<TrackerSimDataCollection> simDataCollection {m_simDataCollectionKey, ctx};
+  ATH_CHECK(simDataCollection.isValid());
 
-  
-  StatusCode FaserSpacePoints::execute(const EventContext& ctx) const {
-    ATH_MSG_INFO(name() << "::" << __FUNCTION__);
-  
-    m_run = ctx.eventID().run_number();
-    m_event = ctx.eventID().event_number();
-    ATH_MSG_DEBUG("run = " << m_run);
-    ATH_MSG_DEBUG("event = " << m_event);
-
-    SG::ReadHandle<FaserSCT_SpacePointContainer> sct_spcontainer( m_sct_spcontainer, ctx );
-    if (!sct_spcontainer.isValid()){
-      ATH_MSG_FATAL("Could not find the data object " << sct_spcontainer.name());
-      return StatusCode::RECOVERABLE;
+
+  // create map of the IdentifierHash of a SCT sensor and the true position of the corresponding hit
+  std::map<IdentifierHash, HepGeom::Point3D<double>> truePositions;
+  for (const FaserSiHit& hit: *siHitCollection) {
+    if ((!hit.particleLink()) || (hit.particleLink()->pdg_id() != -13)) {
+      continue;
+    }
+    Identifier waferId = m_idHelper->wafer_id(hit.getStation(), hit.getPlane(), hit.getRow(), hit.getModule(), hit.getSensor());
+    IdentifierHash waferHash = m_idHelper->wafer_hash(waferId);
+    const HepGeom::Point3D<double> localPosition = hit.localStartPosition();
+    const TrackerDD::SiDetectorElement* element = m_detMgr->getDetectorElement(waferId);
+    if (element) {
+      const HepGeom::Point3D<double> globalPosition =
+          Amg::EigenTransformToCLHEP(element->transformHit()) * localPosition;
+      truePositions[waferHash] = globalPosition;
     }
+  }
+
 
-    m_nCollections = sct_spcontainer->size();
-    ATH_MSG_DEBUG("nCollections = " << m_nCollections);
-    FaserSCT_SpacePointContainer::const_iterator it = sct_spcontainer->begin();
-    FaserSCT_SpacePointContainer::const_iterator it_end = sct_spcontainer->end();
-
-    for(; it != it_end; ++it) {
-      const FaserSCT_SpacePointCollection* col = *it;
-      m_nSpacepoints = col->size();
-      ATH_MSG_DEBUG("nSpacepoints = " << m_nSpacepoints);
-
-      FaserSCT_SpacePointCollection::const_iterator it_sp = col->begin();
-      FaserSCT_SpacePointCollection::const_iterator it_sp_end = col->end();
-      for (; it_sp != it_sp_end; ++it_sp) {
-        const FaserSCT_SpacePoint* sp = *it_sp;
-        const auto id = sp->clusterList().first->identify();
-        int station = m_idHelper->station(id);
-        int plane = m_idHelper->layer(id);
-        m_layer = (station - 1) * 3 + plane;
-
-        Amg::Vector3D position = sp->globalPosition();
-        m_x = position.x();
-        m_y = position.y();
-        m_z = position.z();
-
-        m_outputTree->Fill();
+  for (const FaserSCT_SpacePointCollection* spacePointCollection : *spacePointContainer) {
+    for (const FaserSCT_SpacePoint* spacePoint : *spacePointCollection) {
+      Identifier id1 = spacePoint->cluster1()->identify();
+      Identifier id2 = spacePoint->cluster2()->identify();
+      m_station = m_idHelper->station(id1);
+      m_layer = m_idHelper->layer(id1);
+      m_phi = m_idHelper->phi_module(id1);
+      m_eta = m_idHelper->eta_module(id1);
+      m_side = m_idHelper->side(id1);
+      // check if both clusters originate from a muon
+      m_muon_hit = (muon_hit(simDataCollection.get(), id1) && muon_hit(simDataCollection.get(), id2));
+      // global reconstructed position of space point
+      Amg::Vector3D position = spacePoint->globalPosition();
+      m_x = position.x();
+      m_y = position.y();
+      m_z = position.z();
+      // global true position of simulated hit
+      IdentifierHash waferHash = spacePoint->elementIdList().first;
+      if (truePositions.count(waferHash)) {
+        auto truePosition = truePositions[waferHash];
+        m_tx = truePosition.x();
+        m_ty = truePosition.y();
+        m_tz = truePosition.z();
+      } else {
+        m_tx = -1;
+        m_ty = -1;
+        m_tz = -1;
       }
+      m_tree->Fill();
     }
-    return StatusCode::SUCCESS;
   }
+  return StatusCode::SUCCESS;
+}
 
-  void FaserSpacePoints::initializeTree() {
-    m_outputTree->Branch("run", &m_run);
-    m_outputTree->Branch("event", &m_event);
-    m_outputTree->Branch("layer", &m_layer);
-    m_outputTree->Branch("n_collections", &m_nCollections);
-    m_outputTree->Branch("n_spacepoints", &m_nSpacepoints);
-    m_outputTree->Branch("x", &m_x);
-    m_outputTree->Branch("y", &m_y);
-    m_outputTree->Branch("z", &m_z);
-  }
 }

Tracker/TrackerRecAlgs/FaserSpacePoints/src/FaserSpacePoints.h

@@ -2,45 +2,70 @@
 #define FASERSPACEPOINTS_H
 
 #include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "AthenaBaseComps/AthHistogramming.h"
 #include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrackerSimEvent/FaserSiHitCollection.h"
 #include "TrackerSpacePoint/FaserSCT_SpacePointContainer.h"
+#include "TrackerSimData/TrackerSimDataCollection.h"
 #include <string>
 #include <atomic>
 
 #include "TTree.h"
 #include "TFile.h"
 
+class FaserSCT_ID;
+namespace  TrackerDD {
+class SCT_DetectorManager;
+}
+
+
 namespace Tracker {
-  class FaserSpacePoints : public AthReentrantAlgorithm {
-   public:
-    FaserSpacePoints(const std::string& name, ISvcLocator* pSvcLocator);
-    virtual ~FaserSpacePoints() = default;
-
-    virtual StatusCode initialize() override;
-    virtual StatusCode execute(const EventContext& ctx) const override;
-    virtual StatusCode finalize() override;
-
-   private:
-    void initializeTree();
-
-    mutable int m_run{0};
-    mutable int m_event{0};
-    mutable int m_layer{0};
-    mutable int m_nCollections{0};
-    mutable int m_nSpacepoints{0};
-    mutable float m_x{0};
-    mutable float m_y{0};
-    mutable float m_z{0};
-
-    const FaserSCT_ID* m_idHelper{nullptr};
-    SG::ReadHandleKey<FaserSCT_SpacePointContainer> m_sct_spcontainer{this, "SpacePointsSCTName", "SCT spacepoint container"};
-
-    Gaudi::Property<std::string> m_filePath{this, "FilePath", "spacepoints.root", "Output root file for spacepoints"};
-    Gaudi::Property<std::string> m_treeName{this, "TreeName", "tree", ""};
-
-    TFile* m_outputFile;
-    TTree* m_outputTree;
-  };
+
+class FaserSpacePoints : public AthReentrantAlgorithm, AthHistogramming {
+ public:
+  FaserSpacePoints(const std::string& name, ISvcLocator* pSvcLocator);
+  virtual ~FaserSpacePoints() = default;
+
+  virtual StatusCode initialize() override;
+  virtual StatusCode execute(const EventContext& ctx) const override;
+  virtual StatusCode finalize() override;
+
+  const ServiceHandle<ITHistSvc>& histSvc() const;
+
+ private:
+  const FaserSCT_ID* m_idHelper{nullptr};
+  const TrackerDD::SCT_DetectorManager* m_detMgr {nullptr};
+  bool muon_hit(const TrackerSimDataCollection* simDataMap, Identifier id) const;
+  ServiceHandle<ITHistSvc> m_histSvc;
+
+  SG::ReadHandleKey<FaserSCT_SpacePointContainer> m_spacePointContainerKey {
+    this, "SpacePoints", "SCT_SpacePointContainer"};
+  SG::ReadHandleKey<FaserSiHitCollection> m_siHitCollectionKey {
+    this, "FaserSiHitCollection", "SCT_Hits"};
+  SG::ReadHandleKey<TrackerSimDataCollection> m_simDataCollectionKey {
+    this, "TrackerSimDataCollection", "SCT_SDO_Map"};
+
+  mutable TTree* m_tree;
+  mutable int m_run {0};
+  mutable int m_event {0};
+  mutable int m_station {0};
+  mutable int m_layer {0};
+  mutable int m_phi {0};
+  mutable int m_eta {0};
+  mutable int m_side {0};
+  mutable double m_x {0};
+  mutable double m_y {0};
+  mutable double m_z {0};
+  mutable bool m_muon_hit {false};
+  mutable double m_tx {0};
+  mutable double m_ty {0};
+  mutable double m_tz {0};
+};
+
+inline const ServiceHandle<ITHistSvc>& FaserSpacePoints::histSvc() const {
+  return m_histSvc;
+}
+
 }
 
 #endif // FASERSPACEPOINTS_H

Tracker/TrackerRecAlgs/FaserSpacePoints/test/FaserSpacePointsDbg.py

+#!/usr/bin/env python
+import sys
+from AthenaCommon.Logging import log
+from AthenaCommon.Constants import DEBUG
+from AthenaCommon.Configurable import Configurable
+from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
+from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg
+from FaserSpacePoints.FaserSpacePointsConfig import FaserSpacePointsCfg
+
+log.setLevel(DEBUG)
+Configurable.configurableRun3Behavior = True
+
+ConfigFlags.Input.Files = ['my.RDO.pool.root']
+ConfigFlags.Output.ESDFileName = "spacePoints.ESD.pool.root"
+ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01"
+ConfigFlags.GeoModel.Align.Dynamic = False
+ConfigFlags.Beam.NumberOfCollisions = 0.
+ConfigFlags.lock()
+
+# Core components
+acc = MainServicesCfg(ConfigFlags)
+acc.merge(PoolReadCfg(ConfigFlags))
+acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags))
+acc.merge(TrackerSpacePointFinderCfg(ConfigFlags))
+acc.merge(FaserSpacePointsCfg(ConfigFlags))
+acc.getEventAlgo("Tracker::FaserSpacePoints").OutputLevel = DEBUG
+
+sc = acc.run(maxEvents=1000)
+sys.exit(not sc.isSuccess())

Tracker/TrackerRecAlgs/MCEvents/CMakeLists.txt

+atlas_subdir(MCEvents)
+
+atlas_add_component(
+        MCEvents
+        src/MCEventsAlg.h
+        src/MCEventsAlg.cxx
+        src/components/MCEvents_entries.cxx
+        LINK_LIBRARIES AthenaBaseComps GeneratorObjects StoreGateLib TrackerIdentifier TrackerSimEvent GaudiKernel)
+atlas_install_python_modules(python/*.py)
+atlas_install_scripts(test/*.py)

Tracker/TrackerRecAlgs/MCEvents/python/MCEventsConfig.py

+from AthenaConfiguration.ComponentFactory import CompFactory
+from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
+
+
+def MCEventsCfg(flags, **kwargs):
+    acc = FaserSCT_GeometryCfg(flags)
+    mc_events_alg = CompFactory.MCEventsAlg
+    acc.addEventAlgo(mc_events_alg (**kwargs))
+    return acc

Tracker/TrackerRecAlgs/MCEvents/src/MCEventsAlg.cxx

+#include "MCEventsAlg.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+
+
+MCEventsAlg::MCEventsAlg(const std::string& name, ISvcLocator* pSvcLocator)
+    : AthReentrantAlgorithm(name, pSvcLocator) {}
+
+
+StatusCode MCEventsAlg::initialize() {
+  ATH_CHECK(detStore()->retrieve(m_idHelper,"FaserSCT_ID"));
+  ATH_CHECK(m_mcEventKey.initialize());
+  ATH_CHECK(m_faserSiHitKey.initialize());
+
+  return StatusCode::SUCCESS;
+}
+
+
+StatusCode MCEventsAlg::execute(const EventContext& ctx) const {
+  SG::ReadHandle<FaserSiHitCollection> siHitCollection (m_faserSiHitKey, ctx);
+  std::map<int, std::vector<const FaserSiHit*>> particle_map {};
+  ATH_CHECK(siHitCollection.isValid());
+  for (const FaserSiHit& hit : *siHitCollection) {
+    int barcode = hit.particleLink().barcode();
+    particle_map[barcode].push_back(&hit);
+  }
+
+  SG::ReadHandle<McEventCollection> mcEvents (m_mcEventKey, ctx);
+  ATH_CHECK(mcEvents.isValid());
+  for (const HepMC::GenEvent* evt : *mcEvents) {
+    for (const HepMC::GenParticle* p : evt->particle_range()) {
+      const HepMC::FourVector vertex = p->production_vertex()->position();
+      const HepMC::FourVector momentum = p->momentum();
+      ATH_MSG_DEBUG("pdg: " << p->pdg_id() << ", barcode: " << p->barcode());
+      ATH_MSG_DEBUG("vertex: " << vertex.x() << ", " << vertex.y() << ", " << vertex.z());
+      ATH_MSG_DEBUG("momentum: " << momentum.px() << ", " << momentum.py() << ", " << momentum.pz());
+      auto hits = particle_map[p->barcode()];
+      for (const FaserSiHit* hit : hits) {
+        int station = hit->getStation();
+        int layer = hit->getPlane();
+        int phi = hit->getRow();
+        int eta = hit->getModule();
+        int side = hit->getSensor();
+        ATH_MSG_DEBUG(station << "/" << layer << "/" << phi << "/" << eta << "/" << side);
+      }
+    }
+  }
+  return StatusCode::SUCCESS;
+}
+
+
+StatusCode MCEventsAlg::finalize() {
+  return StatusCode::SUCCESS;
+}

Tracker/TrackerRecAlgs/MCEvents/src/MCEventsAlg.h

+#ifndef MCEVENTS_MCEVENTSALG_H
+#define MCEVENTS_MCEVENTSALG_H
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "TrackerSimEvent/FaserSiHitCollection.h"
+#include "GeneratorObjects/McEventCollection.h"
+
+class FaserSCT_ID;
+
+class MCEventsAlg : public AthReentrantAlgorithm {
+public:
+  MCEventsAlg(const std::string& name, ISvcLocator* pSvcLocator);
+  virtual ~MCEventsAlg() = default;
+
+  virtual StatusCode initialize() override;
+  virtual StatusCode execute(const EventContext& ctx) const override;
+  virtual StatusCode finalize() override;
+
+private:
+  const FaserSCT_ID* m_idHelper {nullptr};
+  SG::ReadHandleKey<McEventCollection> m_mcEventKey { this, "McEventCollection", "BeamTruthEvent" };
+  SG::ReadHandleKey<FaserSiHitCollection> m_faserSiHitKey { this, "FaserSiHitCollection", "SCT_Hits" };
+};
+
+
+#endif // MCEVENTS_MCEVENTSALG_H

Tracker/TrackerRecAlgs/MCEvents/src/components/MCEvents_entries.cxx

+#include "../MCEventsAlg.h"
+
+DECLARE_COMPONENT(MCEventsAlg)
\ No newline at end of file

Tracker/TrackerRecAlgs/MCEvents/test/MCEventsDbg.py

+#!/usr/bin/env python
+
+import sys
+from AthenaCommon.Logging import log
+from AthenaCommon.Constants import DEBUG
+from AthenaCommon.Configurable import Configurable
+from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+from MCEvents.MCEventsConfig import MCEventsCfg
+
+
+log.setLevel(DEBUG)
+Configurable.configurableRun3Behavior = True
+
+# Configure
+ConfigFlags.Input.Files = ['my.RDO.pool.root']
+ConfigFlags.Output.ESDFileName = "MCEvents.ESD.root"
+ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01"
+ConfigFlags.GeoModel.Align.Dynamic = False
+ConfigFlags.Beam.NumberOfCollisions = 0.
+ConfigFlags.lock()
+
+acc = MainServicesCfg(ConfigFlags)
+acc.merge(PoolReadCfg(ConfigFlags))
+acc.merge(MCEventsCfg(ConfigFlags))
+acc.getEventAlgo("MCEventsAlg").OutputLevel = DEBUG
+
+sc = acc.run(maxEvents=10)
+sys.exit(not sc.isSuccess())

Tracker/TrackerRecAlgs/MyClusters/CMakeLists.txt

+atlas_subdir(MyClusters)
+
+atlas_add_component(
+  MyClusters
+  src/MyClustersAlg.h
+  src/MyClustersAlg.cxx
+  src/components/MyClusters_entries.cxx
+  LINK_LIBRARIES AthenaBaseComps GeneratorObjects TrackerPrepRawData StoreGateLib GaudiKernel)
+
+atlas_install_python_modules(python/*.py)
+atlas_install_scripts(test/*.py)

Tracker/TrackerRecAlgs/MyClusters/python/MyClustersConfig.py

+"""
+Copyright (C) 2022 CERN for the benefit of the FASER collaboration
+"""
+
+from AthenaConfiguration.ComponentFactory import CompFactory
+from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
+
+
+def MyClustersCfg(flags, **kwargs):
+    acc = FaserSCT_GeometryCfg(flags)
+    my_clusters_alg = CompFactory.MyClustersAlg
+    acc.addEventAlgo(my_clusters_alg (**kwargs))
+
+    thistSvc = CompFactory.THistSvc()
+    thistSvc.Output += ["HIST DATAFILE='MyClusters.root' OPT='RECREATE'"]
+    acc.addService(thistSvc)
+    return acc