diff --git a/Tracker/TrackerRecAlgs/FaserSpacePoints/CMakeLists.txt b/Tracker/TrackerRecAlgs/FaserSpacePoints/CMakeLists.txt index 4978495c385798463bd84f20276920cdde578765..e88d445536b1770103382549c17e502e4f5f8b09 100644 --- a/Tracker/TrackerRecAlgs/FaserSpacePoints/CMakeLists.txt +++ b/Tracker/TrackerRecAlgs/FaserSpacePoints/CMakeLists.txt @@ -1,12 +1,9 @@ -################################################################################ -# 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 ) diff --git a/Tracker/TrackerRecAlgs/FaserSpacePoints/python/FaserSpacePointsConfig.py b/Tracker/TrackerRecAlgs/FaserSpacePoints/python/FaserSpacePointsConfig.py new file mode 100644 index 0000000000000000000000000000000000000000..36a2008a833453937a2adf9c4c1c2f3d81525362 --- /dev/null +++ b/Tracker/TrackerRecAlgs/FaserSpacePoints/python/FaserSpacePointsConfig.py @@ -0,0 +1,17 @@ +""" +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 diff --git a/Tracker/TrackerRecAlgs/FaserSpacePoints/python/FaserSpacePointsCosmics.py b/Tracker/TrackerRecAlgs/FaserSpacePoints/python/FaserSpacePointsCosmics.py index 48197207d0c1b321e658b74ddb2d40d74ee244dd..78bae901d8c5319dd35ad4f101ff5cb291740395 100644 --- a/Tracker/TrackerRecAlgs/FaserSpacePoints/python/FaserSpacePointsCosmics.py +++ b/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 diff --git a/Tracker/TrackerRecAlgs/FaserSpacePoints/src/FaserSpacePoints.cxx b/Tracker/TrackerRecAlgs/FaserSpacePoints/src/FaserSpacePoints.cxx index b9722a7c2f4318fd6d8e4cf2c76437ced6078d8e..caa4af80a79286b6bdb76e18754541f7ee498568 100644 --- a/Tracker/TrackerRecAlgs/FaserSpacePoints/src/FaserSpacePoints.cxx +++ b/Tracker/TrackerRecAlgs/FaserSpacePoints/src/FaserSpacePoints.cxx @@ -1,104 +1,128 @@ #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); - } } diff --git a/Tracker/TrackerRecAlgs/FaserSpacePoints/src/FaserSpacePoints.h b/Tracker/TrackerRecAlgs/FaserSpacePoints/src/FaserSpacePoints.h index 28807d73e7ac7939c5b830a758846028feedb815..1cc43f65e2f1a674b692487a26421ad526b466f3 100644 --- a/Tracker/TrackerRecAlgs/FaserSpacePoints/src/FaserSpacePoints.h +++ b/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 diff --git a/Tracker/TrackerRecAlgs/FaserSpacePoints/test/FaserSpacePointsDbg.py b/Tracker/TrackerRecAlgs/FaserSpacePoints/test/FaserSpacePointsDbg.py new file mode 100644 index 0000000000000000000000000000000000000000..a3d50f97a2127b4b33fec1c9cf50ab85218e8621 --- /dev/null +++ b/Tracker/TrackerRecAlgs/FaserSpacePoints/test/FaserSpacePointsDbg.py @@ -0,0 +1,32 @@ +#!/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())