diff --git a/Tracker/TrackerRecAlgs/TrackerTruth/CMakeLists.txt b/Tracker/TrackerRecAlgs/TrackerTruth/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..b52a694c458489bb341116be7568c18ffd9de8e6 --- /dev/null +++ b/Tracker/TrackerRecAlgs/TrackerTruth/CMakeLists.txt @@ -0,0 +1,12 @@ +atlas_subdir(TrackerTruth) + +atlas_add_component( + TrackerTruth + src/TrackerTruthAlg.h + src/TrackerTruthAlg.cxx + src/components/TrackerTruth_entries.cxx + LINK_LIBRARIES AthenaBaseComps GeneratorObjects TrackerSimEvent TrackerIdentifier TrackerReadoutGeometry +) + +atlas_install_python_modules( python/*.py ) +atlas_install_scripts( test/*.py ) \ No newline at end of file diff --git a/Tracker/TrackerRecAlgs/TrackerTruth/ReadMe b/Tracker/TrackerRecAlgs/TrackerTruth/ReadMe new file mode 100644 index 0000000000000000000000000000000000000000..a763932b380cc4982774d5e409328dadea867411 --- /dev/null +++ b/Tracker/TrackerRecAlgs/TrackerTruth/ReadMe @@ -0,0 +1,6 @@ +This package can be used to write out the truth information of particles and hits from MC generated events. +For particles the run number, event number, pdg code, vertex and momentum is written out and for hits the +run number, event number, station, plane, module, row, sensor, local position, global position and pdg code of the +corresponding particle is written out. + +To write out truth information generate MC data using the G4FaserAlgConfigNew_Test.py and run TrackerTruthDbg.py diff --git a/Tracker/TrackerRecAlgs/TrackerTruth/python/TrackerTruthConfig.py b/Tracker/TrackerRecAlgs/TrackerTruth/python/TrackerTruthConfig.py new file mode 100644 index 0000000000000000000000000000000000000000..475e7ba6f60844cdea8c4b5da0ef4f43da70d2f5 --- /dev/null +++ b/Tracker/TrackerRecAlgs/TrackerTruth/python/TrackerTruthConfig.py @@ -0,0 +1,19 @@ +""" +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 TrackerTruthCfg(flags, **kwargs): + acc = FaserSCT_GeometryCfg(flags) + kwargs.setdefault("FaserSiHitCollection", "SCT_Hits") + TrackerTruthAlg = CompFactory.Tracker.TrackerTruthAlg + acc.addEventAlgo(TrackerTruthAlg(**kwargs)) + + thistSvc = CompFactory.THistSvc() + thistSvc.Output += ["HIST1 DATAFILE='TrackerTruth.root' OPT='RECREATE'"] + acc.addService(thistSvc) + + return acc diff --git a/Tracker/TrackerRecAlgs/TrackerTruth/python/__init__.py b/Tracker/TrackerRecAlgs/TrackerTruth/python/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/Tracker/TrackerRecAlgs/TrackerTruth/src/TrackerTruthAlg.cxx b/Tracker/TrackerRecAlgs/TrackerTruth/src/TrackerTruthAlg.cxx new file mode 100644 index 0000000000000000000000000000000000000000..cc600bad74f6e16341757becc305ea0343416889 --- /dev/null +++ b/Tracker/TrackerRecAlgs/TrackerTruth/src/TrackerTruthAlg.cxx @@ -0,0 +1,112 @@ +#include "TrackerTruthAlg.h" +#include "TrackerSimEvent/FaserSiHitCollection.h" +#include "TrackerIdentifier/FaserSCT_ID.h" +#include "TrackerReadoutGeometry/SCT_DetectorManager.h" +#include "TrackerReadoutGeometry/SiDetectorElement.h" +#include <TTree.h> + + +namespace Tracker { + +TrackerTruthAlg::TrackerTruthAlg(const std::string &name, ISvcLocator *pSvcLocator) + : AthReentrantAlgorithm(name, pSvcLocator), AthHistogramming(name), + m_sct(nullptr), m_sID(nullptr), m_histSvc("THistSvc/THistSvc", name) {} + + +StatusCode TrackerTruthAlg::initialize() { + ATH_CHECK(detStore()->retrieve(m_sct, "SCT")); + ATH_CHECK(detStore()->retrieve(m_sID, "FaserSCT_ID")); + ATH_CHECK(m_mcEventCollectionKey.initialize()); + ATH_CHECK(m_faserSiHitKey.initialize()); + + m_hits = new TTree("hits", "hits"); + m_hits->Branch("run_number", &m_run_number, "run_number/I"); + m_hits->Branch("event_number", &m_event_number, "event_number/I"); + m_hits->Branch("station", &m_station, "station/I"); + m_hits->Branch("plane", &m_plane, "plane/I"); + m_hits->Branch("module", &m_module, "module/I"); + m_hits->Branch("row", &m_row, "row/I"); + m_hits->Branch("sensor", &m_sensor, "sensor/I"); + m_hits->Branch("pdg", &m_pdg, "pdg/I"); + m_hits->Branch("local_x", &m_local_x, "local_x/D"); + m_hits->Branch("local_y", &m_local_y, "local_y/D"); + m_hits->Branch("global_x", &m_global_x, "global_x/D"); + m_hits->Branch("global_y", &m_global_y, "global_y/D"); + m_hits->Branch("global_z", &m_global_z, "global_z/D"); + ATH_CHECK(histSvc()->regTree("/HIST1/hits", m_hits)); + + m_particles = new TTree("particles", "particles"); + m_particles->Branch("run_number", &m_run_number, "run_number/I"); + m_particles->Branch("event_number", &m_event_number, "event_number/I"); + m_particles->Branch("pdg", &m_pdg, "pdg/I"); + m_particles->Branch("x", &m_vertex_x, "x/D"); + m_particles->Branch("y", &m_vertex_y, "y/D"); + m_particles->Branch("z", &m_vertex_z, "z/D"); + m_particles->Branch("px", &m_px, "px/D"); + m_particles->Branch("py", &m_px, "py/D"); + m_particles->Branch("pz", &m_px, "pz/D"); + ATH_CHECK(histSvc()->regTree("/HIST1/particles", m_particles)); + + return StatusCode::SUCCESS; +} + + +StatusCode TrackerTruthAlg::execute(const EventContext &ctx) const { + + m_run_number = ctx.eventID().run_number(); + m_event_number = ctx.eventID().event_number(); + + // Fill tree with simulated particles + SG::ReadHandle<McEventCollection> mcEventCollection(m_mcEventCollectionKey, ctx); + ATH_CHECK(mcEventCollection.isValid()); + for (const HepMC::GenEvent* event : *mcEventCollection) { + for (const HepMC::GenParticle* particle : event->particle_range()) { + m_pdg = particle->pdg_id(); + const auto vertex = particle->production_vertex()->position(); + m_vertex_x = vertex.x(); + m_vertex_y = vertex.y(); + m_vertex_z = vertex.z(); + const HepMC::FourVector& momentum = particle->momentum(); + m_px = momentum.x(); + m_py = momentum.y(); + m_pz = momentum.z(); + m_particles->Fill(); + } + } + + // Fill tree with simulated hits + SG::ReadHandle<FaserSiHitCollection> siHitCollection(m_faserSiHitKey, ctx); + ATH_CHECK(siHitCollection.isValid()); + for (const FaserSiHit &hit: *siHitCollection) { + if (!hit.particleLink().isValid()) + continue; + const HepGeom::Point3D<double> localStartPos = hit.localStartPosition(); + Identifier id = m_sID->wafer_id(hit.getStation(), hit.getPlane(), hit.getRow(), hit.getModule(), hit.getSensor()); + const TrackerDD::SiDetectorElement *geoelement = m_sct->getDetectorElement(id); + if (!geoelement) + continue; + const HepGeom::Point3D<double> globalStartPos = + Amg::EigenTransformToCLHEP(geoelement->transformHit()) * HepGeom::Point3D<double>(localStartPos); + m_station = hit.getStation(); + m_plane = hit.getPlane(); + m_module = hit.getModule(); + m_row = hit.getRow(); + m_sensor = hit.getSensor(); + m_pdg = hit.particleLink()->pdg_id(); + m_local_x = localStartPos.x(); + m_local_y = localStartPos.y(); + m_global_x = globalStartPos.x(); + m_global_y = globalStartPos.y(); + m_global_z = globalStartPos.z(); + m_hits->Fill(); + } + + return StatusCode::SUCCESS; +} + + +StatusCode TrackerTruthAlg::finalize() { + return StatusCode::SUCCESS; +} + +} // namespace Tracker diff --git a/Tracker/TrackerRecAlgs/TrackerTruth/src/TrackerTruthAlg.h b/Tracker/TrackerRecAlgs/TrackerTruth/src/TrackerTruthAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..7f21acc718500b5881bdb4956980fc1893c80747 --- /dev/null +++ b/Tracker/TrackerRecAlgs/TrackerTruth/src/TrackerTruthAlg.h @@ -0,0 +1,72 @@ +/* + Copyright (C) 2022 CERN for the benefit of the FASER collaboration +*/ + +#ifndef TRACKERTRUTH_TRACKERTRUTHALG_H +#define TRACKERTRUTH_TRACKERTRUTHALG_H + + +#include "AthenaBaseComps/AthReentrantAlgorithm.h" +#include "AthenaBaseComps/AthHistogramming.h" +#include "TrackerSimEvent/FaserSiHitCollection.h" +#include "GeneratorObjects/McEventCollection.h" + +class TTree; +class FaserSCT_ID; +namespace TrackerDD { +class SCT_DetectorManager; +} + + +namespace Tracker { + +class TrackerTruthAlg : public AthReentrantAlgorithm, AthHistogramming { +public: + TrackerTruthAlg(const std::string &name, ISvcLocator *pSvcLocator); + virtual ~TrackerTruthAlg() = default; + virtual StatusCode initialize() override; + virtual StatusCode execute(const EventContext &ctx) const override; + virtual StatusCode finalize() override; + const ServiceHandle <ITHistSvc> &histSvc() const; + +private: + SG::ReadHandleKey<McEventCollection> m_mcEventCollectionKey {this, "McEventCollection", "BeamTruthEvent"}; + SG::ReadHandleKey <FaserSiHitCollection> m_faserSiHitKey{this, "FaserSiHitCollection", "SCT_Hits"}; + + ServiceHandle <ITHistSvc> m_histSvc; + + const TrackerDD::SCT_DetectorManager *m_sct; + const FaserSCT_ID *m_sID; + + mutable TTree *m_particles; + mutable TTree *m_hits; + + mutable unsigned int m_run_number; + mutable unsigned int m_event_number; + mutable unsigned int m_station; + mutable unsigned int m_plane; + mutable unsigned int m_module; + mutable unsigned int m_row; + mutable unsigned int m_sensor; + mutable int m_pdg; + mutable double m_local_x; + mutable double m_local_y; + mutable double m_global_x; + mutable double m_global_y; + mutable double m_global_z; + mutable double m_vertex_x; + mutable double m_vertex_y; + mutable double m_vertex_z; + mutable double m_px; + mutable double m_py; + mutable double m_pz; +}; + + +inline const ServiceHandle <ITHistSvc> &TrackerTruthAlg::histSvc() const { + return m_histSvc; +} + +} // namespace Tracker + +#endif diff --git a/Tracker/TrackerRecAlgs/TrackerTruth/src/components/TrackerTruth_entries.cxx b/Tracker/TrackerRecAlgs/TrackerTruth/src/components/TrackerTruth_entries.cxx new file mode 100644 index 0000000000000000000000000000000000000000..0f49f0e7edefeccc3908ebc6822928ad08abbe7c --- /dev/null +++ b/Tracker/TrackerRecAlgs/TrackerTruth/src/components/TrackerTruth_entries.cxx @@ -0,0 +1,3 @@ +#include "../TrackerTruthAlg.h" + +DECLARE_COMPONENT(Tracker::TrackerTruthAlg ) \ No newline at end of file diff --git a/Tracker/TrackerRecAlgs/TrackerTruth/test/TrackerTruthDbg.py b/Tracker/TrackerRecAlgs/TrackerTruth/test/TrackerTruthDbg.py new file mode 100644 index 0000000000000000000000000000000000000000..479edc81bb284720c94ef4b16196a5aa7ac0f5d3 --- /dev/null +++ b/Tracker/TrackerRecAlgs/TrackerTruth/test/TrackerTruthDbg.py @@ -0,0 +1,38 @@ +#!/usr/bin/env python +"""" +Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration +""" +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 TrackerTruth.TrackerTruthConfig import TrackerTruthCfg + +log.setLevel(DEBUG) +Configurable.configurableRun3Behavior = True + +ConfigFlags.Input.Files = ["my.HITS.pool.root"] +ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01" +ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" +ConfigFlags.Input.ProjectName = "data22" +ConfigFlags.Input.isMC = True +ConfigFlags.GeoModel.FaserVersion = "FASER-01" +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)) + +acc.merge(TrackerTruthCfg(ConfigFlags)) +acc.getEventAlgo("Tracker::TrackerTruthAlg").OutputLevel = DEBUG + +# Execute and finish +sc = acc.run(maxEvents=-1) +sys.exit(not sc.isSuccess())