Muon R4 - Add infrastructure to digitize BI-RPC hits

Muon R4 - Add infrastructure to digitize BI-RPC hits

MuonSpectrometer/MuonPhaseII/MuonCnv/xAODMuonSimHitCnv/src/xAODSimHitToRpcMeasCnvAlg.cxx

 /*
-   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;
 }

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