From 4b40ebeb4182bb10e9d966632eafbfdd270425e2 Mon Sep 17 00:00:00 2001
From: Johannes Junggeburth <johannes.josef.junggeburth@cern.ch>
Date: Wed, 8 May 2024 09:36:56 +0200
Subject: [PATCH] Muon R4 - Add infrastructure to digitize BI-RPC hits
Muon R4 - Add infrastructure to digitize BI-RPC hits
---
.../src/xAODSimHitToRpcMeasCnvAlg.cxx | 248 ++++++++++++------
.../src/xAODSimHitToRpcMeasCnvAlg.h | 32 +++
2 files changed, 206 insertions(+), 74 deletions(-)
diff --git a/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.cxx b/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.cxx
index 32611e78c423..ffd7ee8eebea 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.cxx
+++ b/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.cxx
@@ -1,10 +1,12 @@
/*
- Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
#include "xAODSimHitToRpcMeasCnvAlg.h"
#include <xAODMuonPrepData/RpcStripAuxContainer.h>
+#include <xAODMuonPrepData/RpcStrip2DAuxContainer.h>
+
#include <MuonReadoutGeometryR4/RpcReadoutElement.h>
#include <MuonReadoutGeometryR4/MuonChamber.h>
@@ -21,6 +23,9 @@ namespace {
constexpr double percentage( unsigned int numerator, unsigned int denom) {
return 100. * numerator / std::max(denom, 1u);
}
+
+ using CheckVector2D = MuonGMR4::CheckVector2D;
+
}
xAODSimHitToRpcMeasCnvAlg::xAODSimHitToRpcMeasCnvAlg(const std::string& name,
@@ -31,7 +36,10 @@ StatusCode xAODSimHitToRpcMeasCnvAlg::initialize(){
ATH_CHECK(m_geoCtxKey.initialize());
ATH_CHECK(m_readKey.initialize());
ATH_CHECK(m_writeKey.initialize());
+ ATH_CHECK(m_writeKeyBI.initialize(!m_writeKeyBI.empty()));
+
ATH_CHECK(m_idHelperSvc.retrieve());
+ m_stIdxBIL = m_idHelperSvc->rpcIdHelper().stationNameIndex("BIL");
ATH_CHECK(detStore()->retrieve(m_DetMgr));
return StatusCode::SUCCESS;
}
@@ -41,6 +49,156 @@ StatusCode xAODSimHitToRpcMeasCnvAlg::finalize() {
<<percentage(m_acceptedHits[1], m_allHits[1]) <<" cases, the conversion was successful");
return StatusCode::SUCCESS;
}
+void xAODSimHitToRpcMeasCnvAlg::digitizeHit(const double hitTime,
+ const double locPosOnStrip,
+ const MuonGMR4::StripDesignPtr& designPtr,
+ const Identifier& gasGapId,
+ const bool measuresPhi,
+ xAOD::RpcStripContainer& prdContainer,
+ CLHEP::HepRandomEngine* rndEngine) const {
+
+ /// There're Rpc chambers without phi strips (BI)
+ if (!designPtr){
+ return;
+ }
+ const MuonGMR4::StripDesign& design{*designPtr};
+
+ const double uncert = design.stripPitch() / std::sqrt(12.);
+ const double smearedX = CLHEP::RandGaussZiggurat::shoot(rndEngine, locPosOnStrip, uncert);
+ const Amg::Vector2D locHitPos{smearedX * Amg::Vector2D::UnitX()};
+ ++(m_allHits[measuresPhi]);
+ if (!design.insideTrapezoid(locHitPos)) {
+ ATH_MSG_VERBOSE("The hit "<<Amg::toString(locHitPos)<<" is outside of the trapezoid bounds for "
+ <<m_idHelperSvc->toStringGasGap(gasGapId)<<", measuresPhi: "<<(measuresPhi ? "yay" : "nay"));
+ return;
+ }
+ const int strip = stripNumber(design, locHitPos);
+ if (strip < 0) {
+ ATH_MSG_VERBOSE("Hit " << Amg::toString(locHitPos) << " cannot trigger any signal in a strip for "
+ << m_idHelperSvc->toStringGasGap(gasGapId) <<", measuresPhi: "<<(measuresPhi ? "yay" : "nay"));
+ return;
+ }
+
+ const RpcIdHelper& id_helper{m_idHelperSvc->rpcIdHelper()};
+
+ bool isValid{false};
+ const Identifier prdId{id_helper.channelID(gasGapId,
+ id_helper.doubletZ(gasGapId),
+ id_helper.doubletPhi(gasGapId),
+ id_helper.gasGap(gasGapId),
+ measuresPhi, strip, isValid)};
+
+ if (!isValid) {
+ ATH_MSG_WARNING("Invalid hit identifier obtained for "<<m_idHelperSvc->toStringGasGap(gasGapId)
+ <<", eta strip "<<strip<<" & hit "<<Amg::toString(locHitPos,2 )
+ <<" /// "<<design);
+ return;
+ }
+ ++(m_acceptedHits[measuresPhi]);
+
+ xAOD::RpcStrip* prd = new xAOD::RpcStrip();
+ prdContainer.push_back(prd);
+
+ prd->setIdentifier(prdId.get_compact());
+ xAOD::MeasVector<1> lPos{smearedX};
+ xAOD::MeasMatrix<1> cov{uncert * uncert};
+ prd->setMeasurement<1>(m_idHelperSvc->detElementHash(prdId),
+ std::move(lPos), std::move(cov));
+ const MuonGMR4::RpcReadoutElement* readOutEle = m_DetMgr->getRpcReadoutElement(prdId);
+ prd->setReadoutElement(readOutEle);
+ prd->setStripNumber(strip);
+ prd->setGasGap(id_helper.gasGap(prdId));
+ prd->setDoubletPhi(id_helper.doubletPhi(prdId));
+ prd->setMeasuresPhi(id_helper.measuresPhi(prdId));
+ prd->setTime(hitTime);
+ prd->setAmbiguityFlag(0);
+}
+
+void xAODSimHitToRpcMeasCnvAlg::digitizeHit(const double hitTime,
+ const Amg::Vector2D& locPosOnStrip,
+ const MuonGMR4::StripDesignPtr& designPtr,
+ const Identifier& gasGapId,
+ xAOD::RpcStrip2DContainer& prdContainer,
+ CLHEP::HepRandomEngine* rndEngine) const{
+
+ /// There're Rpc chambers without phi strips (BI)
+ if (!designPtr){
+ return;
+ }
+ const MuonGMR4::StripDesign& design{*designPtr};
+
+ const double uncert = design.stripPitch() / std::sqrt(12.);
+ const double smearedX = CLHEP::RandGaussZiggurat::shoot(rndEngine, locPosOnStrip.x(), uncert);
+ const Amg::Vector2D locHitPos{smearedX, locPosOnStrip.y()};
+ ++(m_allHits[false]);
+ if (!design.insideTrapezoid(locHitPos)) {
+ ATH_MSG_VERBOSE("The 2D hit "<<Amg::toString(locHitPos)<<" is outside of the trapezoid bounds for "
+ <<m_idHelperSvc->toStringGasGap(gasGapId));
+ return;
+ }
+ const int strip = stripNumber(design, locHitPos);
+ if (strip < 0) {
+ ATH_MSG_VERBOSE("Hit " << Amg::toString(locHitPos) << " cannot trigger any signal in a strip for "
+ << m_idHelperSvc->toStringGasGap(gasGapId));
+ return;
+ }
+
+ const RpcIdHelper& id_helper{m_idHelperSvc->rpcIdHelper()};
+
+ bool isValid{false};
+ const Identifier prdId{id_helper.channelID(gasGapId,
+ id_helper.doubletZ(gasGapId),
+ id_helper.doubletPhi(gasGapId),
+ id_helper.gasGap(gasGapId),
+ false, strip, isValid)};
+
+ if (!isValid) {
+ ATH_MSG_WARNING("Invalid hit identifier obtained for "<<m_idHelperSvc->toStringGasGap(gasGapId)
+ <<", eta strip "<<strip<<" & hit "<<Amg::toString(locHitPos,2 )
+ <<" /// "<<design);
+ return;
+ }
+ ++(m_acceptedHits[false]);
+
+ xAOD::RpcStrip2D* prd = new xAOD::RpcStrip2D();
+ prdContainer.push_back(prd);
+
+ prd->setIdentifier(prdId.get_compact());
+ xAOD::MeasVector<2> lPos{};
+ xAOD::MeasMatrix<2> cov{xAOD::MeasMatrix<2>::Identity()};
+ cov(0,0) = uncert * uncert;
+ /// Dummy value for the moment
+ cov(1,1) = 1.;
+ prd->setMeasurement<2>(m_idHelperSvc->detElementHash(prdId),
+ xAOD::toStorage(locHitPos), std::move(cov));
+ const MuonGMR4::RpcReadoutElement* readOutEle = m_DetMgr->getRpcReadoutElement(prdId);
+ prd->setReadoutElement(readOutEle);
+ prd->setStripNumber(strip);
+ prd->setGasGap(id_helper.gasGap(prdId));
+ prd->setDoubletPhi(id_helper.doubletPhi(prdId));
+ prd->setTime(hitTime);
+ prd->setAmbiguityFlag(0);
+}
+
+int xAODSimHitToRpcMeasCnvAlg::stripNumber(const MuonGMR4::StripDesign& design,
+ const Amg::Vector2D& locHitPos) const {
+
+ int strip = design.stripNumber(locHitPos);
+ if (strip > 0) return strip;
+
+ const CheckVector2D firstStrip = design.center(1);
+ const CheckVector2D lastStrip = design.center(design.numStrips());
+ if (!firstStrip || !lastStrip) {
+ return -1;
+ }
+ if ( (*firstStrip).x() - 0.5 *design.stripPitch() < locHitPos.x()) {
+ return 1;
+ } else if ( (*lastStrip).x() + 0.5 * design.stripPitch() > locHitPos.x()) {
+ return design.numStrips();
+ }
+ return -1;
+}
+
StatusCode xAODSimHitToRpcMeasCnvAlg::execute(const EventContext& ctx) const {
SG::ReadHandle<xAOD::MuonSimHitContainer> simHitContainer{m_readKey, ctx};
if (!simHitContainer.isPresent()){
@@ -56,77 +214,15 @@ StatusCode xAODSimHitToRpcMeasCnvAlg::execute(const EventContext& ctx) const {
ATH_CHECK(prdContainer.record(std::make_unique<xAOD::RpcStripContainer>(),
std::make_unique<xAOD::RpcStripAuxContainer>()));
- const RpcIdHelper& id_helper{m_idHelperSvc->rpcIdHelper()};
- CLHEP::HepRandomEngine* rndEngine = getRandomEngine(ctx);
+ SG::WriteHandle<xAOD::RpcStrip2DContainer> prd2DContainer{};
+ if (!m_writeKeyBI.empty()) {
+ prd2DContainer = SG::WriteHandle<xAOD::RpcStrip2DContainer>{m_writeKeyBI, ctx};
+ ATH_CHECK(prd2DContainer.record(std::make_unique<xAOD::RpcStrip2DContainer>(),
+ std::make_unique<xAOD::RpcStrip2DAuxContainer>()));
- using CheckVector2D = MuonGMR4::CheckVector2D;
- auto digitizeHit = [&] (const double hitTime,
- const double locX,
- const MuonGMR4::StripDesignPtr& designPtr,
- const Identifier& hitId,
- bool measPhi) {
- /// There're Rpc chambers without phi strips (BI)
- if (!designPtr){
- return;
- }
- const MuonGMR4::StripDesign& design{*designPtr};
- const double uncert = design.stripPitch() / std::sqrt(12.);
- const double smearedX = CLHEP::RandGaussZiggurat::shoot(rndEngine, locX, uncert);
- const Amg::Vector2D locHitPos{locX * Amg::Vector2D::UnitX()};
- ++(m_allHits[measPhi]);
- if (!design.insideTrapezoid(locHitPos)) {
- ATH_MSG_VERBOSE("The hit "<<Amg::toString(locHitPos)<<" is outside of the trapezoid bounds for "
- <<m_idHelperSvc->toStringGasGap(hitId)<<", measuresPhi: "<<(measPhi ? "yay" : "nay"));
- return;
- }
- int stripNumber = design.stripNumber(locHitPos);
- /// There're subtle cases where the hit is smeared outside the boundaries
- if (stripNumber < 0) {
- const CheckVector2D firstStrip = design.center(1);
- const CheckVector2D lastStrip = design.center(design.numStrips());
- if (!firstStrip || !lastStrip) {
- return;
- }
- if ( (*firstStrip).x() - 0.5 *design.stripPitch() < locHitPos.x()) {
- stripNumber = 1;
- } else if ( (*lastStrip).x() + 0.5 * design.stripPitch() > locHitPos.x()) {
- stripNumber = design.numStrips();
- } else {
- ATH_MSG_VERBOSE("Hit " << Amg::toString(locHitPos) << " cannot trigger any signal in a strip for "
- << m_idHelperSvc->toStringGasGap(hitId) <<", measuresPhi: "<<(measPhi ? "yay" : "nay"));
- return;
- }
- }
- bool isValid{false};
- const Identifier prdId{id_helper.channelID(hitId,
- id_helper.doubletZ(hitId),
- id_helper.doubletPhi(hitId),
- id_helper.gasGap(hitId),
- measPhi, stripNumber, isValid)};
-
- if (!isValid) {
- ATH_MSG_WARNING("Invalid hit identifier obtained for "<<m_idHelperSvc->toStringGasGap(hitId)
- <<", eta strip "<<stripNumber<<" & hit "<<Amg::toString(locHitPos,2 )
- <<" /// "<<design);
- return;
- }
- ++(m_acceptedHits[measPhi]);
- xAOD::RpcStrip* prd = new xAOD::RpcStrip();
- prdContainer->push_back(prd);
- prd->setIdentifier(prdId.get_compact());
- xAOD::MeasVector<1> lPos{smearedX};
- xAOD::MeasMatrix<1> cov{uncert * uncert};
- prd->setMeasurement<1>(m_idHelperSvc->detElementHash(prdId),
- std::move(lPos), std::move(cov));
- const MuonGMR4::RpcReadoutElement* readOutEle = m_DetMgr->getRpcReadoutElement(prdId);
- prd->setReadoutElement(readOutEle);
- prd->setStripNumber(stripNumber);
- prd->setGasGap(id_helper.gasGap(prdId));
- prd->setDoubletPhi(id_helper.doubletPhi(prdId));
- prd->setMeasuresPhi(id_helper.measuresPhi(prdId));
- prd->setTime(hitTime);
- prd->setAmbiguityFlag(0);
- };
+ }
+
+ CLHEP::HepRandomEngine* rndEngine = getRandomEngine(ctx);
for (const xAOD::MuonSimHit* simHit : *simHitContainer) {
const Identifier hitId = simHit->identify();
@@ -136,9 +232,13 @@ StatusCode xAODSimHitToRpcMeasCnvAlg::execute(const EventContext& ctx) const {
const Amg::Vector3D locSimHitPos{xAOD::toEigen(simHit->localPosition())};
const double hitTime = simHit->globalTime() - invC *(readOutEle->localToGlobalTrans(gctx, hitId) * locSimHitPos).mag();
- digitizeHit(hitTime, locSimHitPos.x(), readOutEle->getParameters().etaDesign, hitId, false);
- digitizeHit(hitTime, locSimHitPos.y(), readOutEle->getParameters().phiDesign, hitId, true);
-
+ if (m_idHelperSvc->stationName(hitId) != m_stIdxBIL) {
+ digitizeHit(hitTime, locSimHitPos.x(), readOutEle->getParameters().etaDesign, hitId, false, *prdContainer, rndEngine);
+ digitizeHit(hitTime, locSimHitPos.y(), readOutEle->getParameters().phiDesign, hitId, true, *prdContainer, rndEngine);
+ } else if (!m_writeKeyBI.empty()) {
+ digitizeHit(hitTime, locSimHitPos.block<2,1>(0,0), readOutEle->getParameters().etaDesign, hitId,
+ *prd2DContainer, rndEngine);
+ }
}
return StatusCode::SUCCESS;
}
diff --git a/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.h b/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.h
index ec18e075d3a8..4b6e1f9adf87 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.h
+++ b/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.h
@@ -11,7 +11,10 @@
#include <StoreGate/WriteHandleKey.h>
#include <xAODMuonSimHit/MuonSimHitContainer.h>
+
#include <xAODMuonPrepData/RpcStripContainer.h>
+#include <xAODMuonPrepData/RpcStrip2DContainer.h>
+
#include <MuonIdHelpers/IMuonIdHelperSvc.h>
#include <MuonReadoutGeometryR4/MuonDetectorManager.h>
@@ -39,6 +42,31 @@ class xAODSimHitToRpcMeasCnvAlg : public AthReentrantAlgorithm {
private:
CLHEP::HepRandomEngine* getRandomEngine(const EventContext& ctx) const;
+ /** @brief Smears the local simHit position orthogonal to the strip and
+ * writes a 1D rpc strip measurement
+ */
+ void digitizeHit(const double hitTime,
+ const double locPosOnStrip,
+ const MuonGMR4::StripDesignPtr& designPtr,
+ const Identifier& gasGapId,
+ const bool measuresPhi,
+ xAOD::RpcStripContainer& prdContainer,
+ CLHEP::HepRandomEngine* engine) const;
+
+ void digitizeHit(const double hitTime,
+ const Amg::Vector2D& locPosOnStrip,
+ const MuonGMR4::StripDesignPtr& designPtr,
+ const Identifier& gasGapId,
+ xAOD::RpcStrip2DContainer& prdContainer,
+ CLHEP::HepRandomEngine* engine) const;
+
+ /** @brief Returns the number of the strip that's closest to the hit position on the strip plane
+ * -1 is returned if the hit is outside the strip panel acceptance
+ * */
+ int stripNumber(const MuonGMR4::StripDesign& design,
+ const Amg::Vector2D& locHitPosOnPlane) const;
+
+
SG::ReadHandleKey<ActsGeometryContext> m_geoCtxKey{this, "AlignmentKey", "ActsAlignment", "cond handle key"};
SG::ReadHandleKey<xAOD::MuonSimHitContainer> m_readKey{this, "InputCollection", "xRpcSimHits",
@@ -47,6 +75,8 @@ class xAODSimHitToRpcMeasCnvAlg : public AthReentrantAlgorithm {
SG::WriteHandleKey<xAOD::RpcStripContainer> m_writeKey{this, "OutputContainer", "xRpcStrips",
"Output container"};
+ SG::WriteHandleKey<xAOD::RpcStrip2DContainer> m_writeKeyBI{this, "OutputContainerBI", "xRpcBILStrips",
+ "Output container of the 2D BIL rpc strips"};
/// Access to the new readout geometry
const MuonGMR4::MuonDetectorManager* m_DetMgr{nullptr};
@@ -58,6 +88,8 @@ class xAODSimHitToRpcMeasCnvAlg : public AthReentrantAlgorithm {
mutable std::array<std::atomic<unsigned>, 2> m_allHits ATLAS_THREAD_SAFE{};
mutable std::array<std::atomic<unsigned>, 2> m_acceptedHits ATLAS_THREAD_SAFE{};
+ int m_stIdxBIL{-1};
+
};
#endif
\ No newline at end of file
--
GitLab