From 8b9980d848af9146b6fbd737cfbdb570c8e8a05d Mon Sep 17 00:00:00 2001
From: Johannes Junggeburth <johannes.josef.junggeburth@cern.ch>
Date: Mon, 26 Feb 2024 11:38:29 +0100
Subject: [PATCH] MuonR4 - Fast digitization restructure Rpc & Tgc
MuonR4 - Fast digitization restructure Rpc & Tgc
---
.../python/testMdtSimHitCnv.py | 7 +-
.../src/xAODSimHitToRpcMeasCnvAlg.cxx | 160 +++++++-------
.../src/xAODSimHitToRpcMeasCnvAlg.h | 7 +-
.../src/xAODSimHitToTgcMeasCnvAlg.cxx | 201 ++++++++++++------
.../src/xAODSimHitToTgcMeasCnvAlg.h | 10 +-
.../MuonReadoutGeometryR4/RadialStripDesign.h | 4 +
.../RadialStripDesign.icc | 9 +-
7 files changed, 248 insertions(+), 150 deletions(-)
diff --git a/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/python/testMdtSimHitCnv.py b/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/python/testMdtSimHitCnv.py
index f4c390f4e19e..668b34321df7 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/python/testMdtSimHitCnv.py
+++ b/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/python/testMdtSimHitCnv.py
@@ -5,12 +5,13 @@
if __name__=="__main__":
from MuonGeoModelTestR4.testGeoModel import setupGeoR4TestCfg, SetupArgParser, executeTest
parser = SetupArgParser()
- parser.set_defaults(nEvents = -1)
+ parser.set_defaults(nEvents = 1000)
+ parser.set_defaults(inputFile=["/cvmfs/atlas-nightlies.cern.ch/repo/data/data-art/MuonRecRTT/R4SimHits.pool.root"])
args = parser.parse_args()
flags, cfg = setupGeoR4TestCfg(args)
- from xAODMuonSimHitCnv.MuonSimHitCnvCfg import xAODtoMdtCnvAlgCfg
- cfg.merge(xAODtoMdtCnvAlgCfg(flags))
+ from xAODMuonSimHitCnv.MuonSimHitCnvCfg import MuonSimHitToMeasurementCfg
+ cfg.merge(MuonSimHitToMeasurementCfg(flags))
executeTest(cfg, num_events = args.nEvents)
diff --git a/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.cxx b/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.cxx
index a51c80c0123b..37d360f11cc0 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.cxx
+++ b/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.cxx
@@ -16,6 +16,9 @@
namespace {
constexpr double invC = 1./ Gaudi::Units::c_light;
+ constexpr double percentage( unsigned int numerator, unsigned int denom) {
+ return 100. * numerator / std::max(denom, 1u);
+ }
}
xAODSimHitToRpcMeasCnvAlg::xAODSimHitToRpcMeasCnvAlg(const std::string& name,
@@ -30,6 +33,12 @@ StatusCode xAODSimHitToRpcMeasCnvAlg::initialize(){
ATH_CHECK(detStore()->retrieve(m_DetMgr));
return StatusCode::SUCCESS;
}
+StatusCode xAODSimHitToRpcMeasCnvAlg::finalize() {
+ ATH_MSG_INFO("Tried to convert "<<m_allHits[0]<<"/"<<m_allHits[1]<<" hits. In, "
+ <<percentage(m_acceptedHits[0], m_allHits[0]) <<"/"
+ <<percentage(m_acceptedHits[1], m_allHits[1]) <<" cases, the conversion was successful");
+ return StatusCode::SUCCESS;
+}
StatusCode xAODSimHitToRpcMeasCnvAlg::execute(const EventContext& ctx) const {
SG::ReadHandle<xAOD::MuonSimHitContainer> simHitContainer{m_readKey, ctx};
if (!simHitContainer.isPresent()){
@@ -44,97 +53,92 @@ StatusCode xAODSimHitToRpcMeasCnvAlg::execute(const EventContext& ctx) const {
const RpcIdHelper& id_helper{m_idHelperSvc->rpcIdHelper()};
CLHEP::HepRandomEngine* rndEngine = getRandomEngine(ctx);
- for (const xAOD::MuonSimHit* simHit : *simHitContainer) {
- const Identifier hitId = simHit->identify();
- // ignore radiation for now
- if (std::abs(simHit->pdgId()) != 13) continue;
-
-
- const MuonGMR4::RpcReadoutElement* readOutEle = m_DetMgr->getRpcReadoutElement(hitId);
- const MuonGMR4::StripDesign& etaDesign{*readOutEle->getParameters().etaDesign};
-
- const Amg::Vector3D locSimHitPos{xAOD::toEigen(simHit->localPosition())};
-
- const double etaUncert = etaDesign.stripPitch() / std::sqrt(12);
- const Amg::Vector2D smearedEtaPos{CLHEP::RandGaussZiggurat::shoot(rndEngine, locSimHitPos.x(), etaUncert), 0.};
-
- const double hitTime = simHit->globalTime() -
- invC *(readOutEle->localToGlobalTrans(gctx, hitId) * locSimHitPos).mag();
- const unsigned int etaStripNum = etaDesign.stripNumber(smearedEtaPos.block<2,1>(0,0));
-
- ATH_MSG_VERBOSE("Convert simulated hit "<<m_idHelperSvc->toStringGasGap(hitId)<<" located in gas gap at "
- <<Amg::toString(locSimHitPos, 2)<<" eta strip number: "<<etaStripNum
- <<" strip position "<<Amg::toString(etaDesign.center(etaStripNum).value_or(Amg::Vector2D::Zero()), 2));
-
+ double hitTime{0.};
+ const MuonGMR4::RpcReadoutElement* readOutEle{nullptr};
- bool isValid{false};
- const Identifier etaHitId{id_helper.channelID(hitId, readOutEle->doubletZ(),
- id_helper.doubletPhi(hitId),
- id_helper.gasGap(hitId),
- false, etaStripNum, isValid)};
-
- if (!isValid) {
- ATH_MSG_WARNING("Invalid hit identifier obtained for "<<m_idHelperSvc->toStringGasGap(hitId)
- <<", eta strip "<<etaStripNum<<" & hit "<<Amg::toString(locSimHitPos,2 ));
- } else {
- xAOD::RpcStrip* prd = new xAOD::RpcStrip();
- prdContainer->push_back(prd);
- prd->setIdentifier(etaHitId.get_compact());
- xAOD::MeasVector<1> lPos{smearedEtaPos.x()};
- xAOD::MeasMatrix<1> cov{etaUncert * etaUncert};
- prd->setMeasurement<1>(m_idHelperSvc->detElementHash(etaHitId),
- std::move(lPos), std::move(cov));
- prd->setReadoutElement(readOutEle);
- prd->setStripNumber(etaStripNum);
- prd->setGasGap(id_helper.gasGap(etaHitId));
- prd->setDoubletPhi(id_helper.doubletPhi(etaHitId));
- prd->setMeasuresPhi(id_helper.measuresPhi(etaHitId));
- prd->setReadoutElement(readOutEle);
- prd->setTime(hitTime);
- prd->setAmbiguityFlag(0);
- const Amg::Vector3D strip3D = lPos.x() * Amg::Vector3D::UnitX();
- const Amg::Transform3D& globToCenter{m_surfaceProvTool->globalToChambCenter(gctx,etaHitId)};
- prd->setStripPosInStation(xAOD::toStorage(globToCenter * readOutEle->localToGlobalTrans(gctx,prd->layerHash()) * strip3D));
+ using CheckVector2D = MuonGMR4::CheckVector2D;
+ auto digitizeHit = [&] (const double locX,
+ const MuonGMR4::StripDesignPtr& designPtr,
+ const Identifier& hitId,
+ bool measPhi) {
+ /// There're Rpc chambers without phi strips (BI)
+ if (!designPtr){
+ return;
}
- /// Check whether the read out element contains phi strips or not.
- if (!readOutEle->nPhiStrips()) {
- continue;
+ 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] = m_allHits[measPhi] + 1;
+ 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;
}
- const MuonGMR4::StripDesign& phiDesign{*readOutEle->getParameters().phiDesign};
- const double phiUncert = phiDesign.stripPitch() / std::sqrt(12.);
- const Amg::Vector2D smearedPhiPos{CLHEP::RandGaussZiggurat::shoot(rndEngine, locSimHitPos.y(), phiUncert), 0.};
-
- const unsigned int phiStripNum = phiDesign.stripNumber(smearedPhiPos.block<2,1>(0,0));
- const Identifier phiHitId{id_helper.channelID(hitId, readOutEle->doubletZ(),
- id_helper.doubletPhi(hitId),
- id_helper.gasGap(hitId),
- true, phiStripNum, isValid)};
-
+ 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)
- <<", phi strip "<<phiStripNum<<" & hit "<<Amg::toString(locSimHitPos,2 ));
- continue;
+ <<", eta strip "<<stripNumber<<" & hit "<<Amg::toString(locHitPos,2 )
+ <<" /// "<<design);
+ return;
}
-
+ m_acceptedHits[measPhi] = m_acceptedHits[measPhi] + 1;
xAOD::RpcStrip* prd = new xAOD::RpcStrip();
prdContainer->push_back(prd);
- prd->setIdentifier(phiHitId.get_compact());
- xAOD::MeasVector<1> lPos{smearedPhiPos.x()};
- xAOD::MeasMatrix<1> cov{phiUncert * phiUncert};
- prd->setMeasurement<1>(m_idHelperSvc->detElementHash(phiHitId),
- std::move(lPos), std::move(cov));
+ 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));
prd->setReadoutElement(readOutEle);
- prd->setStripNumber(phiStripNum);
- prd->setGasGap(id_helper.gasGap(phiHitId));
- prd->setDoubletPhi(id_helper.doubletPhi(phiHitId));
- prd->setMeasuresPhi(id_helper.measuresPhi(phiHitId));
+ prd->setStripNumber(stripNumber);
+ prd->setGasGap(id_helper.gasGap(prdId));
+ prd->setDoubletPhi(id_helper.doubletPhi(prdId));
+ prd->setMeasuresPhi(id_helper.measuresPhi(prdId));
prd->setReadoutElement(readOutEle);
prd->setTime(hitTime);
prd->setAmbiguityFlag(0);
- const Amg::Vector3D strip3D = lPos.x() * Amg::Vector3D::UnitX();
- const Amg::Transform3D& globToCenter{m_surfaceProvTool->globalToChambCenter(gctx, phiHitId)};
- prd->setStripPosInStation(xAOD::toStorage(globToCenter * readOutEle->localToGlobalTrans(gctx,prd->layerHash()) * strip3D));
+ const Amg::Vector3D strip3D = lPos.x() * Amg::Vector3D::UnitX();
+ const Amg::Transform3D& globToCenter{m_surfaceProvTool->globalToChambCenter(gctx,prdId)};
+ prd->setStripPosInStation(xAOD::toStorage(globToCenter * readOutEle->localToGlobalTrans(gctx,prd->layerHash()) * strip3D));
+ };
+
+ for (const xAOD::MuonSimHit* simHit : *simHitContainer) {
+ const Identifier hitId = simHit->identify();
+ // ignore radiation for now
+ if (std::abs(simHit->pdgId()) != 13) continue;
+ readOutEle = m_DetMgr->getRpcReadoutElement(hitId);
+ const Amg::Vector3D locSimHitPos{xAOD::toEigen(simHit->localPosition())};
+ hitTime = simHit->globalTime() - invC *(readOutEle->localToGlobalTrans(gctx, hitId) * locSimHitPos).mag();
+
+ digitizeHit(locSimHitPos.x(), readOutEle->getParameters().etaDesign, hitId, false);
+ digitizeHit(locSimHitPos.y(), readOutEle->getParameters().phiDesign, hitId, true);
+
}
return StatusCode::SUCCESS;
}
diff --git a/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.h b/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.h
index 56f311a2589f..ca2c312b2782 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.h
+++ b/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.h
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
#ifndef XAODMUONSIMHITCNV_xAODSimHitToRpcMeasurementCnvAlg_H
#define XAODMUONSIMHITCNV_xAODSimHitToRpcMeasurementCnvAlg_H
@@ -19,6 +19,7 @@
#include <AthenaKernel/IAthRNGSvc.h>
#include <CLHEP/Random/RandomEngine.h>
+#include "CxxUtils/checker_macros.h"
/**
* The xAODSimHitToRpcMasCnvAlg is a short cut towards the RpcStrip measurement
* The RpcSimHits are taken and expressed w.r.t. eta & phi gas gaps.
@@ -34,6 +35,7 @@ class xAODSimHitToRpcMeasCnvAlg : public AthReentrantAlgorithm {
StatusCode execute(const EventContext& ctx) const override;
StatusCode initialize() override;
+ StatusCode finalize() override;
private:
CLHEP::HepRandomEngine* getRandomEngine(const EventContext& ctx) const;
@@ -54,6 +56,9 @@ class xAODSimHitToRpcMeasCnvAlg : public AthReentrantAlgorithm {
PublicToolHandle<MuonGMR4::IMuonStationLayerSurfaceTool> m_surfaceProvTool{this, "LayerGeoTool", ""};
+ mutable std::array<std::atomic<unsigned>, 2> m_allHits ATLAS_THREAD_SAFE{};
+ mutable std::array<std::atomic<unsigned>, 2> m_acceptedHits ATLAS_THREAD_SAFE{};
+
};
#endif
\ No newline at end of file
diff --git a/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToTgcMeasCnvAlg.cxx b/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToTgcMeasCnvAlg.cxx
index 74c061e89ce6..048086563c5a 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToTgcMeasCnvAlg.cxx
+++ b/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToTgcMeasCnvAlg.cxx
@@ -14,6 +14,12 @@
// Random Numbers
#include <AthenaKernel/RNGWrapper.h>
+namespace {
+ constexpr double percentage( unsigned int numerator, unsigned int denom) {
+ return 100. * numerator / std::max(denom, 1u);
+ }
+
+}
xAODSimHitToTgcMeasCnvAlg::xAODSimHitToTgcMeasCnvAlg(const std::string& name,
ISvcLocator* pSvcLocator):
AthReentrantAlgorithm{name, pSvcLocator} {}
@@ -26,6 +32,12 @@ StatusCode xAODSimHitToTgcMeasCnvAlg::initialize(){
ATH_CHECK(detStore()->retrieve(m_DetMgr));
return StatusCode::SUCCESS;
}
+StatusCode xAODSimHitToTgcMeasCnvAlg::finalize() {
+ ATH_MSG_INFO("Tried to convert "<<m_allHits[0]<<"/"<<m_allHits[1]<<" hits. In, "
+ <<percentage(m_acceptedHits[0], m_allHits[0]) <<"/"
+ <<percentage(m_acceptedHits[1], m_allHits[1]) <<" cases, the conversion was successful");
+ return StatusCode::SUCCESS;
+}
StatusCode xAODSimHitToTgcMeasCnvAlg::execute(const EventContext& ctx) const {
SG::ReadHandle<xAOD::MuonSimHitContainer> simHitContainer{m_readKey, ctx};
if (!simHitContainer.isPresent()){
@@ -40,84 +52,141 @@ StatusCode xAODSimHitToTgcMeasCnvAlg::execute(const EventContext& ctx) const {
const TgcIdHelper& id_helper{m_idHelperSvc->tgcIdHelper()};
CLHEP::HepRandomEngine* rndEngine = getRandomEngine(ctx);
- for (const xAOD::MuonSimHit* simHit : *simHitContainer) {
- const Identifier hitId = simHit->identify();
- // ignore radiation for now
- if (std::abs(simHit->pdgId()) != 13) continue;
-
- const unsigned int gasGap = id_helper.gasGap(hitId);
+ auto dumpHit = [&](const Identifier& hitId,
+ const double xLoc,
+ const double xUncert) {
+
+ xAOD::TgcStrip* prd = new xAOD::TgcStrip();
+ prdContainer->push_back(prd);
+ prd->setIdentifier(hitId.get_compact());
+ xAOD::MeasVector<1> lPos{xLoc};
+ xAOD::MeasMatrix<1> cov{xUncert * xUncert};
+ prd->setMeasurement<1>(m_idHelperSvc->detElementHash(hitId),
+ std::move(lPos), std::move(cov));
+ prd->setChannelNumber(id_helper.channel(hitId));
+ prd->setGasGap(id_helper.gasGap(hitId));
+ const bool measPhi{id_helper.measuresPhi(hitId)};
+ m_acceptedHits[measPhi] = m_acceptedHits[measPhi] + 1;
+ prd->setMeasuresPhi(measPhi);
const MuonGMR4::TgcReadoutElement* readOutEle = m_DetMgr->getTgcReadoutElement(hitId);
- const MuonGMR4::StripDesign& etaDesign{readOutEle->wireGangLayout(gasGap)};
-
- const Amg::Vector3D locSimHitPos{xAOD::toEigen(simHit->localPosition())};
+ prd->setReadoutElement(readOutEle);
+ const Amg::Vector3D strip3D = lPos.x() * Amg::Vector3D::UnitX();
+ const Amg::Transform3D& globToCenter{m_surfaceProvTool->globalToChambCenter(gctx, hitId)};
+ prd->setStripPosInStation(xAOD::toStorage(globToCenter * readOutEle->localToGlobalTrans(gctx,prd->layerHash()) * strip3D));
+ };
+
+ auto processEtaHit = [&] (const Amg::Vector3D& locSimHitPos,
+ const Identifier& hitId) {
+ m_allHits[false] = m_allHits[false] + 1;
+ const MuonGMR4::TgcReadoutElement* readOutEle = m_DetMgr->getTgcReadoutElement(hitId);
+ const unsigned int gasGap = id_helper.gasGap(hitId);
+ const MuonGMR4::WireGroupDesign& design{readOutEle->wireGangLayout(gasGap)};
+ if (!design.insideTrapezoid(locSimHitPos.block<2,1>(0,0))) {
+ ATH_MSG_DEBUG("The hit "<<Amg::toString(locSimHitPos)<<" in "<<m_idHelperSvc->toStringGasGap(hitId)
+ <<" is outside of the trapezoid "<<design);
+ return;
+ }
+ const int wireGrpNum = design.stripNumber(locSimHitPos.block<2,1>(0,0));
- const double etaUncert = etaDesign.stripPitch() / std::sqrt(12);
- const Amg::Vector2D smearedEtaPos{CLHEP::RandGaussZiggurat::shoot(rndEngine, locSimHitPos.x(), etaUncert), 0.};
- const unsigned int etaStripNum = etaDesign.stripNumber(smearedEtaPos.block<2,1>(0,0));
-
- ATH_MSG_VERBOSE("Convert simulated hit "<<m_idHelperSvc->toStringGasGap(hitId)<<" located in gas gap at "
- <<Amg::toString(locSimHitPos, 2)<<" eta strip number: "<<etaStripNum
- <<" strip position "<<Amg::toString(etaDesign.center(etaStripNum).value_or(Amg::Vector2D::Zero()), 2));
-
+ if (wireGrpNum < 0) {
+ ATH_MSG_DEBUG("True hit is not "<<Amg::toString(locSimHitPos)<<" "<<m_idHelperSvc->toStringGasGap(hitId));
+ return;
+ }
+
+ const double uncert = design.stripPitch() * design.numWiresInGroup(wireGrpNum) / std::sqrt(12);
+ const double locX = CLHEP::RandGaussZiggurat::shoot(rndEngine, locSimHitPos.x(), uncert);
+ /// Recalculate the strip number with the smeared hit -> Use the real Y to ensure that the
+ /// hit remains within the active trapzoid
+ const Amg::Vector2D smearedPos{locX, locSimHitPos.y()};
+ const int prdWireNum = design.stripNumber(smearedPos);
+ if (prdWireNum < 0) {
+ if (design.insideTrapezoid(smearedPos)) {
+ ATH_MSG_WARNING("True hit "<<Amg::toString(locSimHitPos, 2)<<" corresponding to "<<wireGrpNum<<" --> "
+ <<Amg::toString(smearedPos)<<" "<<uncert<<" is outside of "<<design);
+ }
+ return;
+ }
bool isValid{false};
- const Identifier etaHitId{id_helper.channelID(hitId, gasGap, false,
- etaStripNum, isValid)};
-
+ const Identifier prdId{id_helper.channelID(hitId, gasGap, false, prdWireNum, isValid)};
if (!isValid) {
- ATH_MSG_WARNING("Invalid hit identifier obtained for "<<m_idHelperSvc->toStringGasGap(hitId)
- <<", eta strip "<<etaStripNum<<" & hit "<<Amg::toString(locSimHitPos,2 ));
- } else {
- xAOD::TgcStrip* prd = new xAOD::TgcStrip();
- prdContainer->push_back(prd);
- prd->setIdentifier(etaHitId.get_compact());
- xAOD::MeasVector<1> lPos{smearedEtaPos.x()};
- xAOD::MeasMatrix<1> cov{etaUncert * etaUncert};
- prd->setMeasurement<1>(m_idHelperSvc->detElementHash(etaHitId),
- std::move(lPos), std::move(cov));
- prd->setReadoutElement(readOutEle);
- prd->setChannelNumber(etaStripNum);
- prd->setGasGap(gasGap);
- prd->setMeasuresPhi(false);
- prd->setReadoutElement(readOutEle);
- const Amg::Vector3D strip3D = lPos.x() * Amg::Vector3D::UnitX();
- const Amg::Transform3D& globToCenter{m_surfaceProvTool->globalToChambCenter(gctx,etaHitId)};
- prd->setStripPosInStation(xAOD::toStorage(globToCenter * readOutEle->localToGlobalTrans(gctx,prd->layerHash()) * strip3D));
+ ATH_MSG_WARNING("Invalid channel "<< m_idHelperSvc->toStringGasGap(hitId)<<", channel: "<<prdWireNum);
+ }
+ ATH_MSG_VERBOSE("Convert simulated hit "<<m_idHelperSvc->toStringGasGap(hitId)<<" located in gas gap at "
+ <<Amg::toString(locSimHitPos, 2)<<" eta strip number: "<<prdWireNum
+ <<" strip position "<<Amg::toString(design.center(prdWireNum).value_or(Amg::Vector2D::Zero()), 2));
+
+ dumpHit(prdId, locX, uncert);
+ };
+ auto processStripHit = [&] (const Amg::Vector3D& locSimHitPos,
+ const Identifier& hitId) {
+ m_allHits[true] = m_allHits[true] + 1;
+ const MuonGMR4::TgcReadoutElement* readOutEle = m_DetMgr->getTgcReadoutElement(hitId);
+ const unsigned int gasGap = id_helper.gasGap(hitId);
+ const MuonGMR4::RadialStripDesign& design{readOutEle->stripLayout(gasGap)};
+ if (!design.insideTrapezoid(locSimHitPos.block<2,1>(0,0))) {
+ ATH_MSG_DEBUG("The eta hit "<<Amg::toString(locSimHitPos)<<" in "<<m_idHelperSvc->toStringGasGap(hitId)
+ <<" is outside of the trapezoid "<<design);
+ return;
}
- /// Check whether the read out element contains phi strips or not.
- if (!readOutEle->numStrips(gasGap)) {
- continue;
+ int stripNum = design.stripNumber(locSimHitPos.block<2,1>(0,0));
+ if (stripNum < 0) {
+ ATH_MSG_WARNING("Strip hit "<<Amg::toString(locSimHitPos)<<" cannot be assigned to any active strip for "
+ <<m_idHelperSvc->toStringGasGap(hitId)<<". "<<design);
+ return;
}
- const MuonGMR4::StripDesign& phiDesign{readOutEle->stripLayout(gasGap)};
- const double phiUncert = phiDesign.stripPitch() / std::sqrt(12.);
- const Amg::Vector2D smearedPhiPos{CLHEP::RandGaussZiggurat::shoot(rndEngine, locSimHitPos.y(), phiUncert), 0.};
-
- const unsigned int phiStripNum = phiDesign.stripNumber(smearedPhiPos.block<2,1>(0,0));
- const Identifier phiHitId{id_helper.channelID(hitId, gasGap, true,
- phiStripNum, isValid)};
+ const double stripPitch = std::abs(*Amg::intersect<2>(design.stripLeftBottom(stripNum), design.stripLeftEdge(stripNum),
+ locSimHitPos.block<2,1>(0,0), design.stripNormal(stripNum))) +
+ std::abs(*Amg::intersect<2>(design.stripRightBottom(stripNum), design.stripRightEdge(stripNum),
+ locSimHitPos.block<2,1>(0,0), design.stripNormal(stripNum)));
+ const double uncert = stripPitch / std::sqrt(12);
+ const double locX = CLHEP::RandGaussZiggurat::shoot(rndEngine, locSimHitPos.x(), uncert);
+ /// Recalculate the strip number with the smeared hit -> Use the real Y to ensure that the
+ /// hit remains within the active trapzoid
+ const Amg::Vector2D smearedPos{locX, locSimHitPos.y()};
+ const int prdStripNum = design.stripNumber(smearedPos);
+
+ if (prdStripNum < 0) {
+ if (design.insideTrapezoid(smearedPos)) {
+ ATH_MSG_WARNING("True phi hit "<<Amg::toString(locSimHitPos, 2)<<" corresponding to "<<stripNum<<" --> "
+ <<Amg::toString(smearedPos)<<" "<<uncert<<" is outside of "<<design);
+ }
+ return;
+ }
+ bool isValid{false};
+ const Identifier prdId{id_helper.channelID(hitId, gasGap, true, prdStripNum, isValid)};
if (!isValid) {
- ATH_MSG_WARNING("Invalid hit identifier obtained for "<<m_idHelperSvc->toStringGasGap(hitId)
- <<", phi strip "<<phiStripNum<<" & hit "<<Amg::toString(locSimHitPos,2 ));
- continue;
+ ATH_MSG_WARNING("Invalid channel "<< m_idHelperSvc->toStringGasGap(hitId)<<", channel: "<<prdStripNum);
+ }
+ ATH_MSG_VERBOSE("Convert simulated hit "<<m_idHelperSvc->toStringGasGap(hitId)<<" located in gas gap at "
+ <<Amg::toString(locSimHitPos, 2)<<" eta strip number: "<<prdStripNum
+ <<" strip position "<<Amg::toString(design.center(prdStripNum).value_or(Amg::Vector2D::Zero()), 2));
+
+ dumpHit(prdId, locX, uncert);
+ };
+ for (const xAOD::MuonSimHit* simHit : *simHitContainer) {
+ const Identifier hitId = simHit->identify();
+ // ignore radiation for now
+ if (std::abs(simHit->pdgId()) != 13) continue;
+
+ const MuonGMR4::TgcReadoutElement* readOutEle = m_DetMgr->getTgcReadoutElement(hitId);
+ const Amg::Vector3D locSimHitPos{xAOD::toEigen(simHit->localPosition())};
+ const int gasGap = id_helper.gasGap(hitId);
+
+ if (readOutEle->numWireGangs(gasGap)) {
+ processEtaHit(locSimHitPos, hitId);
}
+ if (readOutEle->numStrips(gasGap)) {
+ const IdentifierHash stripHash{MuonGMR4::TgcReadoutElement::constructHash(0, gasGap, true)};
+ const IdentifierHash wireHash{MuonGMR4::TgcReadoutElement::constructHash(0, gasGap, false)};
- xAOD::TgcStrip* prd = new xAOD::TgcStrip();
- prdContainer->push_back(prd);
- prd->setIdentifier(phiHitId.get_compact());
- xAOD::MeasVector<1> lPos{smearedPhiPos.x()};
- xAOD::MeasMatrix<1> cov{phiUncert * phiUncert};
- prd->setMeasurement<1>(m_idHelperSvc->detElementHash(phiHitId),
- std::move(lPos), std::move(cov));
- prd->setReadoutElement(readOutEle);
- prd->setChannelNumber(phiStripNum);
- prd->setGasGap(gasGap);
- prd->setMeasuresPhi(true);
- prd->setReadoutElement(readOutEle);
- const Amg::Vector3D strip3D = lPos.x() * Amg::Vector3D::UnitX();
- const Amg::Transform3D& globToCenter{m_surfaceProvTool->globalToChambCenter(gctx, phiHitId)};
- prd->setStripPosInStation(xAOD::toStorage(globToCenter * readOutEle->localToGlobalTrans(gctx,prd->layerHash()) * strip3D));
+ const Amg::Transform3D toPhiRot{readOutEle->globalToLocalTrans(gctx, stripHash) *
+ readOutEle->localToGlobalTrans(gctx, wireHash)};
+
+ processStripHit(toPhiRot * locSimHitPos, hitId);
+ }
}
return StatusCode::SUCCESS;
}
diff --git a/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToTgcMeasCnvAlg.h b/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToTgcMeasCnvAlg.h
index bca0d2c4ad94..1f12d5d0a4a2 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToTgcMeasCnvAlg.h
+++ b/MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToTgcMeasCnvAlg.h
@@ -19,6 +19,9 @@
#include <AthenaKernel/IAthRNGSvc.h>
#include <CLHEP/Random/RandomEngine.h>
+
+#include "CxxUtils/checker_macros.h"
+
/**
* The xAODSimHitToTgcMasCnvAlg is a short cut towards the TgcStrip measurement
* The TgcSimHits are taken and expressed w.r.t. eta & phi gas gaps.
@@ -34,7 +37,7 @@ class xAODSimHitToTgcMeasCnvAlg : public AthReentrantAlgorithm {
StatusCode execute(const EventContext& ctx) const override;
StatusCode initialize() override;
-
+ StatusCode finalize() override;
private:
CLHEP::HepRandomEngine* getRandomEngine(const EventContext& ctx) const;
@@ -54,6 +57,11 @@ class xAODSimHitToTgcMeasCnvAlg : public AthReentrantAlgorithm {
PublicToolHandle<MuonGMR4::IMuonStationLayerSurfaceTool> m_surfaceProvTool{this, "LayerGeoTool", ""};
+
+ mutable std::array<std::atomic<unsigned>, 2> m_allHits ATLAS_THREAD_SAFE{};
+ mutable std::array<std::atomic<unsigned>, 2> m_acceptedHits ATLAS_THREAD_SAFE{};
+
+
};
#endif
\ No newline at end of file
diff --git a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonReadoutGeometryR4/MuonReadoutGeometryR4/RadialStripDesign.h b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonReadoutGeometryR4/MuonReadoutGeometryR4/RadialStripDesign.h
index 7e4eb9d2d06e..f003d0ae014b 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonReadoutGeometryR4/MuonReadoutGeometryR4/RadialStripDesign.h
+++ b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonReadoutGeometryR4/MuonReadoutGeometryR4/RadialStripDesign.h
@@ -37,6 +37,10 @@ class RadialStripDesign: public StripDesign {
* @param: Strip number in the global scheme [1- nStrips()]
*/
Amg::Vector2D stripDir(int stripNumber) const;
+ /** @brief: Returns the direction of the left edge */
+ Amg::Vector2D stripLeftEdge(int stripNumber) const;
+ /** @brief: Returns the direction of the right edge */
+ Amg::Vector2D stripRightEdge(int stripNumber) const;
/** @bief: Returns the vector perpendicular to the stripDir and pointing to the next strip*/
Amg::Vector2D stripNormal(int stripNumber) const;
/** @brief: Returns the intersection of the left strip edge at the bottom panel's edge*/
diff --git a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonReadoutGeometryR4/MuonReadoutGeometryR4/RadialStripDesign.icc b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonReadoutGeometryR4/MuonReadoutGeometryR4/RadialStripDesign.icc
index db16d05658f2..b628a59308a4 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonReadoutGeometryR4/MuonReadoutGeometryR4/RadialStripDesign.icc
+++ b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonReadoutGeometryR4/MuonReadoutGeometryR4/RadialStripDesign.icc
@@ -32,7 +32,14 @@ namespace MuonGMR4{
m_strips[stripNum +1].bottomMounting()));
}
-
+ inline Amg::Vector2D RadialStripDesign::stripLeftEdge(int stripNumber) const{
+ CHECK_STRIPRANGE(stripNumber);
+ return m_strips[stripCh].fromBottomToTop().unit();
+ }
+ inline Amg::Vector2D RadialStripDesign::stripRightEdge(int stripNumber) const {
+ CHECK_STRIPRANGE(stripNumber);
+ return m_strips[stripCh + 1].fromBottomToTop().unit();
+ }
inline Amg::Vector2D RadialStripDesign::stripDir(int stripNumber) const {
CHECK_STRIPRANGE(stripNumber);
return (m_strips[stripCh].fromBottomToTop() + m_strips[stripCh+1].fromBottomToTop()).unit();
--
GitLab