From 822479463119a07de02008ee88f45c0ef1b437fd Mon Sep 17 00:00:00 2001
From: Johannes Junggeburth <johannes.josef.junggeburth@cern.ch>
Date: Tue, 20 Feb 2024 16:40:48 +0100
Subject: [PATCH] RpcGeoModelR4 - Use the stripLayer frame as reference to
 express the RE

RpcGeoModelR4 - Use the stripLayer frame as reference to express the RE
---
 .../MuonGeoModelTest/src/GeoModelRpcTest.cxx  | 10 +--
 .../MuonGeoModelTest/src/GeoModelRpcTest.h    |  5 +-
 .../src/MuonGeoUtilitiyTool.cxx               |  4 +-
 .../MuonGeoModelR4/src/RpcReadoutGeomTool.cxx | 84 +++++++++++++------
 .../MuonGeoModelTestR4/python/testGeoModel.py |  2 +-
 .../src/GeoModelMdtTest.cxx                   |  4 +-
 .../MuonGeoModelTestR4/src/GeoModelMdtTest.h  |  3 +-
 .../src/GeoModelRpcTest.cxx                   | 10 ++-
 .../MuonGeoModelTestR4/src/GeoModelRpcTest.h  |  5 +-
 .../util/runRpcGeoComparison.cxx              | 63 +++++++++-----
 .../src/TgcSensitiveDetector.cxx              | 19 +++--
 11 files changed, 136 insertions(+), 73 deletions(-)

diff --git a/MuonSpectrometer/MuonGeoModelTest/src/GeoModelRpcTest.cxx b/MuonSpectrometer/MuonGeoModelTest/src/GeoModelRpcTest.cxx
index d03089777b9f..f258609d7fc9 100644
--- a/MuonSpectrometer/MuonGeoModelTest/src/GeoModelRpcTest.cxx
+++ b/MuonSpectrometer/MuonGeoModelTest/src/GeoModelRpcTest.cxx
@@ -103,8 +103,9 @@ StatusCode GeoModelRpcTest::dumpToTree(const EventContext& ctx, const RpcReadout
 
    m_numStripsEta = readoutEle->Nstrips(false);
    m_numStripsPhi = readoutEle->Nstrips(true);
-   m_numGasGapsEta = readoutEle->NgasGaps(false);
-   m_numGasGapsPhi = readoutEle->NgasGaps(true);   
+   m_numRpcLayers = readoutEle->numberOfLayers();
+   m_numGasGapsPhi = readoutEle->NgasGaps(true);
+   m_numPhiPanels = readoutEle->NphiStripPanels();   
 
    m_stripEtaPitch = readoutEle->StripPitch(false);
    m_stripPhiPitch = readoutEle->StripPitch(true);
@@ -127,12 +128,9 @@ StatusCode GeoModelRpcTest::dumpToTree(const EventContext& ctx, const RpcReadout
         m_ALineRotT   = station->getALine_rotz();
         m_ALineRotZ   = station->getALine_rott();
     }
-
-    const int numGaps = std::max(readoutEle->NgasGaps(false), 
-                                 readoutEle->NgasGaps(true));
     const int maxDoubPhi = std::max(readoutEle->getDoubletPhi(), readoutEle->NphiStripPanels());
     for (int doubPhi = readoutEle->getDoubletPhi(); doubPhi <= maxDoubPhi; ++doubPhi) {
-        for (int gap = 1; gap <= numGaps; ++gap) {   
+        for (int gap = 1; gap <= readoutEle->numberOfLayers(); ++gap) {   
             for (bool measPhi : {false, true}) {
                 unsigned int numStrip = readoutEle->Nstrips(measPhi);
                 for (unsigned int strip = 1; strip <= numStrip ; ++strip) {
diff --git a/MuonSpectrometer/MuonGeoModelTest/src/GeoModelRpcTest.h b/MuonSpectrometer/MuonGeoModelTest/src/GeoModelRpcTest.h
index fc8445058bee..0725f34b4f55 100644
--- a/MuonSpectrometer/MuonGeoModelTest/src/GeoModelRpcTest.h
+++ b/MuonSpectrometer/MuonGeoModelTest/src/GeoModelRpcTest.h
@@ -57,6 +57,7 @@ class GeoModelRpcTest : public AthHistogramAlgorithm {
     MuonVal::ScalarBranch<uint8_t>& m_doubletZ{m_tree.newScalar<uint8_t>("stationDoubletZ")};
     MuonVal::ScalarBranch<uint8_t>& m_doubletPhi{m_tree.newScalar<uint8_t>("stationDoubletPhi")};
     MuonVal::ScalarBranch<std::string>& m_chamberDesign{m_tree.newScalar<std::string>("chamberDesign")};
+
  
     /// Number of strips, strip pitch in eta & phi direction
     MuonVal::ScalarBranch<uint8_t>& m_numStripsEta{m_tree.newScalar<uint8_t>("numEtaStrips")};
@@ -70,9 +71,9 @@ class GeoModelRpcTest : public AthHistogramAlgorithm {
     MuonVal::ScalarBranch<float>& m_stripEtaLength{m_tree.newScalar<float>("stripEtaLength")};
     MuonVal::ScalarBranch<float>& m_stripPhiLength{m_tree.newScalar<float>("stripPhiLength")};
     /// Number of eta & phi gas gaps
-    MuonVal::ScalarBranch<uint8_t>& m_numGasGapsEta{m_tree.newScalar<uint8_t>("numEtaGasGaps")};
+    MuonVal::ScalarBranch<uint8_t>& m_numRpcLayers{m_tree.newScalar<uint8_t>("numRpcLayers")}; 
     MuonVal::ScalarBranch<uint8_t>& m_numGasGapsPhi{m_tree.newScalar<uint8_t>("numPhiGasGaps")};
-    
+    MuonVal::ScalarBranch<uint8_t>& m_numPhiPanels{m_tree.newScalar<uint8_t>("numPhiPanels")};
     /// Transformation of the readout element (Translation, ColX, ColY, ColZ)
     MuonVal::CoordTransformBranch m_readoutTransform{m_tree, "GeoModelTransform"};
     
diff --git a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelR4/src/MuonGeoUtilitiyTool.cxx b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelR4/src/MuonGeoUtilitiyTool.cxx
index f92fc47535e2..6897b6dab857 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelR4/src/MuonGeoUtilitiyTool.cxx
+++ b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelR4/src/MuonGeoUtilitiyTool.cxx
@@ -189,7 +189,9 @@ std::vector<MuonGeoUtilityTool::physVolWithTrans> MuonGeoUtilityTool::findAllLea
     PVConstLink childVol = aV.getVolume();
     const Amg::Transform3D childTrans{aV.getTransform()};
     /// The logical volume has precisely the name for what we're searching for
-    if (childVol->getLogVol()->getName() == volumeName) {
+    
+    if (childVol->getLogVol()->getName() == volumeName ||
+        aV.getName() == volumeName) {
         physVolWithTrans foundNode{};
         foundNode.physVol = childVol;
         foundNode.transform = childTrans;
diff --git a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelR4/src/RpcReadoutGeomTool.cxx b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelR4/src/RpcReadoutGeomTool.cxx
index eeb5787d089e..56ebd1dbd5b2 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelR4/src/RpcReadoutGeomTool.cxx
+++ b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelR4/src/RpcReadoutGeomTool.cxx
@@ -46,6 +46,13 @@ struct gapVolume: public physVolWithTrans {
     
 };
 
+inline bool layerSorter(const physVolWithTrans&a, const physVolWithTrans & b){
+    const Amg::Vector3D cA = a.transform.translation();
+    const Amg::Vector3D cB = b.transform.translation();
+    if (std::abs(cA.x() - cB.x()) > tolerance) return (cA.x() < cB.x());
+    return (cA.y() < cB.y());
+}
+
 
 RpcReadoutGeomTool::RpcReadoutGeomTool(const std::string& type,
                                        const std::string& name,
@@ -83,29 +90,53 @@ StatusCode RpcReadoutGeomTool::loadDimensions(RpcReadoutElement::defineArgs& def
     define.halfLength = box->getZHalfLength() * Gaudi::Units::mm;
     define.halfWidth = box->getYHalfLength() * Gaudi::Units::mm;
    
-    /// Navigate through the GeoModel tree to find all gas volume leaves
-    std::vector<physVolWithTrans> allGasGaps = m_geoUtilTool->findAllLeafNodesByName(define.physVol, "StripLayer");
+    /** Rpc are made up out of 2 or 3 gasGap singlet. A singlet module is a RPC gas gap
+     *  sandwiched by two strip layers. In large sectors, the gas gap may be split
+     *  into two gasGaps.  
+     *
+     *          | Strip layer  |  Strip layer |    |  Strip layer  |  Strip layer |
+     *          |           gas gap           |    |    Gas gap    |    Gas gap   | 
+     *          | Strip layer  |  Strip layer |    |  Strip layer  |  Strip layer |
+     *   
+    */
+    std::vector<physVolWithTrans> stripLayers = m_geoUtilTool->findAllLeafNodesByName(define.physVol, "bottomStripLayer");
+    if (stripLayers.empty()) {
+        ATH_MSG_FATAL("The volume "<<m_idHelperSvc->toStringDetEl(define.detElId)<<" does not have any childern 'bottomStripLayer'"
+            <<std::endl<<m_geoUtilTool->dumpVolume(define.physVol));
+        return StatusCode::FAILURE;
+    }   
+    std::vector<physVolWithTrans> allGasGaps = m_geoUtilTool->findAllLeafNodesByName(define.physVol, "RpcGasGap");
     if (allGasGaps.empty()) {
-        ATH_MSG_FATAL("The volume "<<m_idHelperSvc->toStringDetEl(define.detElId)<<" does not have any childern StripLayer"
+        ATH_MSG_FATAL("The volume "<<m_idHelperSvc->toStringDetEl(define.detElId)<<" does not have any childern 'RpcGasGap'"
             <<std::endl<<m_geoUtilTool->dumpVolume(define.physVol));
         return StatusCode::FAILURE;
     }
-    /// For one reason or another the x-axis points along the gasgap 
-    /// and y along doublet phi
-    std::stable_sort(allGasGaps.begin(), allGasGaps.end(), [](const physVolWithTrans&a, const physVolWithTrans & b){
-         const Amg::Vector3D cA = a.transform.translation();
-         const Amg::Vector3D cB = b.transform.translation();
-         if (std::abs(cA.x() - cB.x()) > tolerance) return (cA.x() < cB.x());
-         return (cA.y() < cB.y());
-    });
+    /// In the GeoModel world, the x-axis points in radial direction & y axis along the phi direction
+    std::stable_sort(allGasGaps.begin(), allGasGaps.end(), layerSorter);
+    std::stable_sort(stripLayers.begin(),stripLayers.end(), layerSorter);
+    /// The strip layers are used to express the dimensions of the strip layer. However, that's projected into the 
+    /// Center of the gasgap which may or maybe not be split into two --> Find the closest gas gap in x for each
+    /// strip layer and overwrite the x coordinate of the strip layerof that one.
+    for (physVolWithTrans& stripLayer : stripLayers){
+        /// Find the closest gas Gap
+        const Amg::Vector3D stripTrans = stripLayer.transform.translation();
+        std::vector<physVolWithTrans>::iterator closestGap =  std::min_element(allGasGaps.begin(), allGasGaps.end(), 
+                             [&stripTrans](const physVolWithTrans& a, const physVolWithTrans& b){
+                                return std::abs(stripTrans.x() - a.transform.translation().x()) <
+                                       std::abs(stripTrans.x() - b.transform.translation().x());                                
+                             });
+        stripLayer.transform.translation().x() = closestGap->transform.translation().x();
+    }
+
     /// Now we need to associate the gasGap volumes with the gas gap number &
     /// the doublet Phi
-    Amg::Vector3D prevGap{allGasGaps[0].transform.translation()};
+    Amg::Vector3D prevGap{stripLayers[0].transform.translation()};
     unsigned int gasGap{1}, doubletPhi{0};
     
     unsigned int modulePhi = m_idHelperSvc->rpcIdHelper().doubletPhi(define.detElId);
     std::vector<gapVolume> allGapsWithIdx{};
-    for (physVolWithTrans& gapVol : allGasGaps) {
+    const bool isAside{m_idHelperSvc->stationEta(define.detElId) > 0};
+    for (physVolWithTrans& gapVol : stripLayers) {
         Amg::Vector3D gCen = gapVol.transform.translation();
         /// The volume points to a new gasgap
         if (std::abs(gCen.x() - prevGap.x()) > tolerance) {
@@ -146,7 +177,7 @@ StatusCode RpcReadoutGeomTool::loadDimensions(RpcReadoutElement::defineArgs& def
                                      paramBook.numEtaStrips);
         /// Define the box layout
         etaDesign->defineTrapezoid(gapBox->getYHalfLength(), gapBox->getYHalfLength(), gapBox->getZHalfLength());
-        gapVol.transform = gapVol.transform * Amg::getRotateY3D(-90. * Gaudi::Units::degree);
+        gapVol.transform = gapVol.transform * Amg::getRotateY3D( (isAside ? -90. :  90.)* Gaudi::Units::degree);
         
         etaDesign = (*factoryCache.stripDesigns.emplace(etaDesign).first);
         const IdentifierHash etaHash {RpcReadoutElement::createHash(0, gapVol.gasGap, gapVol.doubPhi, false)};
@@ -160,9 +191,9 @@ StatusCode RpcReadoutGeomTool::loadDimensions(RpcReadoutElement::defineArgs& def
         if (!define.etaDesign) define.etaDesign = etaDesign;
         StripDesignPtr phiDesign = std::make_unique<StripDesign>();
         phiDesign->defineStripLayout(Amg::Vector2D{-gapBox->getYHalfLength() + paramBook.firstOffSetPhi, 0.},
-                                     paramBook.stripPitchPhi,
-                                     paramBook.stripWidthPhi,
-                                     paramBook.numPhiStrips);
+                                                   paramBook.stripPitchPhi,
+                                                   paramBook.stripWidthPhi,
+                                                   paramBook.numPhiStrips);
         phiDesign->defineTrapezoid(gapBox->getZHalfLength(), gapBox->getZHalfLength(), gapBox->getYHalfLength());
         /// Next build the phi layer
         phiDesign = (*factoryCache.stripDesigns.emplace(phiDesign).first);
@@ -222,7 +253,6 @@ StatusCode RpcReadoutGeomTool::buildReadOutElements(MuonDetectorManager& mgr) {
             ATH_MSG_FATAL("Failed to construct the station Identifier from "<<key);
             return StatusCode::FAILURE;
         }
-        
         defineArgs define{};        
         define.physVol = pv;
         define.chambDesign = key_tokens[1];
@@ -253,16 +283,16 @@ StatusCode RpcReadoutGeomTool::readParameterBook(FactoryCache& cache) {
     for (const IRDBRecord* record : *paramTable) {
         const std::string chambType = record->getString("WRPC_TYPE");
         wRPCTable& parBook = cache.parameterBook[record->getString("WRPC_TYPE")];
-        parBook.stripPitchEta = record->getDouble("STRIPPITCH_Z") * Gaudi::Units::cm;
-        parBook.stripPitchPhi = record->getDouble("STRIPPITCH_S") * Gaudi::Units::cm;
-        const double stripDeadWidth = record->getDouble("STRIPDEADSEP") * Gaudi::Units::cm;
+        parBook.stripPitchEta = record->getDouble("etaStripPitch") * Gaudi::Units::cm;
+        parBook.stripPitchPhi = record->getDouble("phiStripPitch") * Gaudi::Units::cm;
+        const double stripDeadWidth = record->getDouble("stripDeadWidth") * Gaudi::Units::cm;
         parBook.stripWidthEta = parBook.stripPitchEta - stripDeadWidth;
         parBook.stripWidthPhi = parBook.stripPitchPhi - stripDeadWidth;
-        parBook.numEtaStrips = record->getInt("NSTRIPS_Z");
-        parBook.numPhiStrips = record->getInt("NSTRIPS_S");
-        parBook.firstOffSetPhi = record->getDouble("STRIPOFFSET_S") * Gaudi::Units::cm + 
+        parBook.numEtaStrips = record->getInt("nEtaStrips");
+        parBook.numPhiStrips = record->getInt("nPhiStrips");
+        parBook.firstOffSetPhi = record->getDouble("phiStripOffSet") * Gaudi::Units::cm + 
                                  0.5 * parBook.stripPitchPhi;
-        parBook.firstOffSetEta = record->getDouble("STRIPOFFSET_Z") * Gaudi::Units::cm +
+        parBook.firstOffSetEta = record->getDouble("etaStripOffSet") * Gaudi::Units::cm +
                                  record->getDouble("TCKSSU") * Gaudi::Units::cm +
                                  0.5 * parBook.stripPitchEta;
         
@@ -271,8 +301,8 @@ StatusCode RpcReadoutGeomTool::readParameterBook(FactoryCache& cache) {
                        <<", strip pitch (eta/phi) "<<parBook.stripPitchEta<<"/"<<parBook.stripPitchPhi
                        <<", strip width (eta/phi): "<<parBook.stripWidthEta<<"/"<<parBook.stripWidthPhi
                        <<", strip offset (eta/phi): "<<parBook.firstOffSetEta<<"/"<<parBook.firstOffSetPhi
-                       <<", STRIPOFFSET_Z: "<<(record->getDouble("STRIPOFFSET_Z") * Gaudi::Units::cm)
-                       <<", STRIPOFFSET_S: "<<(record->getDouble("STRIPOFFSET_S") * Gaudi::Units::cm)
+                       <<", etaStripOffSet: "<<(record->getDouble("etaStripOffSet") * Gaudi::Units::cm)
+                       <<", phiStripOffSet: "<<(record->getDouble("phiStripOffSet") * Gaudi::Units::cm)
                        <<", TCKSSU: "<<(record->getDouble("TCKSSU")* Gaudi::Units::cm));
     }
     return StatusCode::SUCCESS;
diff --git a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/python/testGeoModel.py b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/python/testGeoModel.py
index 4254b7dc9227..d152b1131fed 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/python/testGeoModel.py
+++ b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/python/testGeoModel.py
@@ -182,7 +182,7 @@ if __name__=="__main__":
     cfg.merge(setupHistSvcCfg(flags, out_file = args.outRootFile))
     chambToTest =  args.chambers if len([x for x in args.chambers if x =="all"]) ==0 else []
     
-    if flags.Detector.GeometryMDT:
+    if False and flags.Detector.GeometryMDT:
         cfg.merge(GeoModelMdtTestCfg(flags, 
                                      TestStations = [ch for ch in chambToTest if ch[0] == "B" or ch[0] == "E"],
                                      ReadoutSideXML="ReadoutSides.xml"))
diff --git a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelMdtTest.cxx b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelMdtTest.cxx
index 7b2d01c46b3a..ee55187500d8 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelMdtTest.cxx
+++ b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelMdtTest.cxx
@@ -166,7 +166,9 @@ StatusCode GeoModelMdtTest::dumpToTree(const EventContext& ctx,
          const Amg::Transform3D& tubeTransform{readoutEle->localToGlobalTrans(gctx,measHash)};
          m_tubeLay.push_back(lay);
          m_tubeNum.push_back(tube);         
-         m_tubeTransform.push_back(tubeTransform);        
+         m_tubeTransform.push_back(tubeTransform);
+         m_tubePosInCh.push_back(m_surfaceProvTool->globalToChambCenter(gctx, readoutEle->identify()) * 
+                                 readoutEle->center(gctx, measHash));
          m_roPos.push_back(readoutEle->readOutPos(gctx, measHash));
          m_tubeLength.push_back(readoutEle->tubeLength(measHash));
          m_activeTubeLength.push_back(readoutEle->activeTubeLength(measHash));
diff --git a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelMdtTest.h b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelMdtTest.h
index 6f0d72fdb1fb..c50054bab450 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelMdtTest.h
+++ b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelMdtTest.h
@@ -81,7 +81,8 @@ class GeoModelMdtTest : public AthHistogramAlgorithm{
 
       /// Position of the readout
       MuonVal::ThreeVectorBranch m_roPos{m_tree, "readOutPos"};
-
+      /// Position of the tube in the chamber frame
+      MuonVal::ThreeVectorBranch m_tubePosInCh{m_tree, "chamberTubePos"};
 
       /// Alignment parameters
       MuonVal::ScalarBranch<float>& m_ALineTransS{m_tree.newScalar<float>("ALineTransS", 0.)};
diff --git a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelRpcTest.cxx b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelRpcTest.cxx
index 9801c37b0244..4c4dab000efb 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelRpcTest.cxx
+++ b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelRpcTest.cxx
@@ -74,8 +74,8 @@ StatusCode GeoModelRpcTest::execute() {
       ATH_MSG_FATAL("Failed to retrieve "<<m_geoCtxKey.fullKey());
       return StatusCode::FAILURE;
     }
-    const ActsGeometryContext& gctx{**geoContextHandle};
-
+    // const ActsGeometryContext& gctx{**geoContextHandle};
+    ActsGeometryContext gctx{};
     for (const Identifier& test_me : m_testStations) {
       ATH_MSG_DEBUG("Test retrieval of Rpc detector element "<<m_idHelperSvc->toStringDetEl(test_me));
       const RpcReadoutElement* reElement = m_detMgr->getRpcReadoutElement(test_me);
@@ -150,8 +150,10 @@ StatusCode GeoModelRpcTest::dumpToTree(const EventContext& ctx,
    m_doubletPhi = reElement->doubletPhi();
    m_chamberDesign = reElement->chamberDesign();
    
-   m_numGasGapsEta = reElement->nGasGaps();
-   m_numGasGapsPhi = reElement->nGasGaps();
+   m_numRpcLayers = reElement->nGasGaps();
+
+   m_numGasGapsPhi = reElement->nPhiPanels();
+   m_numPhiPanels = reElement->nPhiPanels();
 
    ///
    m_numStripsEta = reElement->nEtaStrips();
diff --git a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelRpcTest.h b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelRpcTest.h
index e276fd8bd803..a0c23fbc4461 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelRpcTest.h
+++ b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelRpcTest.h
@@ -68,9 +68,10 @@ class GeoModelRpcTest : public AthHistogramAlgorithm{
       MuonVal::ScalarBranch<float>& m_stripEtaLength{m_tree.newScalar<float>("stripEtaLength")};
       MuonVal::ScalarBranch<float>& m_stripPhiLength{m_tree.newScalar<float>("stripPhiLength")};    
       /// Number of eta & phi gas gaps
-      MuonVal::ScalarBranch<uint8_t>& m_numGasGapsEta{m_tree.newScalar<uint8_t>("numEtaGasGaps")};
-      MuonVal::ScalarBranch<uint8_t>& m_numGasGapsPhi{m_tree.newScalar<uint8_t>("numPhiGasGaps")};
       
+      MuonVal::ScalarBranch<uint8_t>& m_numRpcLayers{m_tree.newScalar<uint8_t>("numRpcLayers")};      
+      MuonVal::ScalarBranch<uint8_t>& m_numGasGapsPhi{m_tree.newScalar<uint8_t>("numPhiGasGaps")};
+      MuonVal::ScalarBranch<uint8_t>& m_numPhiPanels{m_tree.newScalar<uint8_t>("numPhiPanels")};
       /// Transformation of the readout element (Translation, ColX, ColY, ColZ)
       MuonVal::CoordTransformBranch m_readoutTransform{m_tree, "GeoModelTransform"};
       
diff --git a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/util/runRpcGeoComparison.cxx b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/util/runRpcGeoComparison.cxx
index 5091eb707345..9ab9841d3418 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/util/runRpcGeoComparison.cxx
+++ b/MuonSpectrometer/MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/util/runRpcGeoComparison.cxx
@@ -8,6 +8,7 @@
  *        python -m MuonGeoModelTestR4.runGeoModelTest 
  *  
 */
+#include <algorithm>
 #include <GeoPrimitives/GeoPrimitives.h>
 #include <GeoPrimitives/GeoPrimitivesHelpers.h>
 #include <GeoPrimitives/GeoPrimitivesToStringConverter.h>
@@ -23,6 +24,7 @@ using namespace ActsTrk;
 #include <TFile.h>
 #include <TTreeReader.h>
 
+
 constexpr double tolerance = 1.*Gaudi::Units::millimeter;
 
 /// Helper struct to represent a full Rpc chamber
@@ -52,9 +54,12 @@ struct RpcChamber{
 
     unsigned int numStripsEta{0};
     unsigned int numStripsPhi{0};
-
-    unsigned int numGasGapsEta{0};
+    /// Number of rpc singlets along the radial direction
+    unsigned int numLayers{0};
+    /// Number of rpc gasGaps along the phi direction
     unsigned int numGasGapsPhi{0};
+    /// Number of rpc readout panels along the phi direction
+    unsigned int numPhiPanels{0};
 
     struct RpcStrip{
         Amg::Vector3D position{Amg::Vector3D::Zero()};
@@ -155,8 +160,9 @@ std::set<RpcChamber> readTreeDump(const std::string& inputFile) {
     TTreeReaderValue<uint8_t> numStripsEta{treeReader, "numEtaStrips"};
     TTreeReaderValue<uint8_t> numStripsPhi{treeReader, "numPhiStrips"};
     /// Number of eta & phi gas gaps
-    TTreeReaderValue<uint8_t> numGasGapsEta{treeReader, "numEtaGasGaps"};
     TTreeReaderValue<uint8_t> numGasGapsPhi{treeReader, "numPhiGasGaps"};
+    TTreeReaderValue<uint8_t> numPhiPanels{treeReader, "numPhiPanels"};
+    TTreeReaderValue<uint8_t> numLayers{treeReader, "numRpcLayers"};
        
     /// Strip dimensions 
     TTreeReaderValue<float> stripEtaPitch{treeReader, "stripEtaPitch"};
@@ -171,6 +177,12 @@ std::set<RpcChamber> readTreeDump(const std::string& inputFile) {
     TTreeReaderValue<std::vector<float>> geoModelTransformY{treeReader, "GeoModelTransformY"};
     TTreeReaderValue<std::vector<float>> geoModelTransformZ{treeReader, "GeoModelTransformZ"};
 
+
+    TTreeReaderValue<std::vector<float>> stripRotTranslationX{treeReader, "stripRotTranslationX"};
+    TTreeReaderValue<std::vector<float>> stripRotTranslationY{treeReader, "stripRotTranslationY"};
+    TTreeReaderValue<std::vector<float>> stripRotTranslationZ{treeReader, "stripRotTranslationZ"};
+
+
     TTreeReaderValue<std::vector<float>> stripRotCol1X{treeReader, "stripRotLinearCol1X"};
     TTreeReaderValue<std::vector<float>> stripRotCol1Y{treeReader, "stripRotLinearCol1Y"};
     TTreeReaderValue<std::vector<float>> stripRotCol1Z{treeReader, "stripRotLinearCol1Z"};
@@ -215,8 +227,9 @@ std::set<RpcChamber> readTreeDump(const std::string& inputFile) {
 
         newchamber.numStripsEta = (*numStripsEta);
         newchamber.numStripsPhi = (*numStripsPhi);
-        newchamber.numGasGapsEta = (*numGasGapsEta);
         newchamber.numGasGapsPhi = (*numGasGapsPhi);
+        newchamber.numPhiPanels = (*numPhiPanels);
+        newchamber.numLayers = (*numLayers);
         
 
         Amg::Vector3D geoTrans{(*geoModelTransformX)[0], (*geoModelTransformY)[0], (*geoModelTransformZ)[0]};
@@ -243,10 +256,11 @@ std::set<RpcChamber> readTreeDump(const std::string& inputFile) {
             newLayer.gasGap = (*stripRotGasGap)[l];
             newLayer.doubletPhi = (*stripRotDblPhi)[l];
             Amg::RotationMatrix3D stripRot{Amg::RotationMatrix3D::Identity()};
+            Amg::Vector3D translation{(*stripRotTranslationX)[l],(*stripRotTranslationY)[l],(*stripRotTranslationZ)[l]};
             stripRot.col(0) = Amg::Vector3D((*stripRotCol1X)[l],(*stripRotCol1Y)[l], (*stripRotCol1Z)[l]);
             stripRot.col(1) = Amg::Vector3D((*stripRotCol2X)[l],(*stripRotCol2Y)[l], (*stripRotCol2Z)[l]);
             stripRot.col(2) = Amg::Vector3D((*stripRotCol3X)[l],(*stripRotCol3Y)[l], (*stripRotCol3Z)[l]);
-            newLayer.transform = Amg::getTransformFromRotTransl(std::move(stripRot), Amg::Vector3D::Zero());
+            newLayer.transform = Amg::getTransformFromRotTransl(std::move(stripRot), std::move(translation));
             newchamber.layers.insert(std::move(newLayer));
         } 
         
@@ -305,25 +319,28 @@ int main( int argc, char** argv ) {
         std::cerr<<"The file "<<testFile<<" should contain at least one chamber "<<std::endl;
         return EXIT_FAILURE;
     }
+    std::cout<<"Read "<<refChambers.size()<<" chambers from reference: "<<refFile
+             <<" & "<<testChambers.size()<<" from "<<testFile<<std::endl;
     int return_code = EXIT_SUCCESS;
     /// Start to loop over the chambers
     for (const RpcChamber& reference : refChambers) {
         std::set<RpcChamber>::const_iterator test_itr = testChambers.find(reference);
         
         if (test_itr == testChambers.end()) {
-            std::cerr<<"The chamber "<<reference<<" is not part of the testing "<<std::endl;
+            std::cerr<<"runRpcGeoComparison() "<<__LINE__<<": The chamber "<<reference<<" is not part of the testing "<<std::endl;
             return_code = EXIT_FAILURE;
             continue;
         }
         bool chamberOkay = true;
         const RpcChamber& test = {*test_itr};
+      
+        TEST_BASICPROP(numLayers, "number of rpc singlets");
+        TEST_BASICPROP(numGasGapsPhi, "number of phi gas gaps");
+        TEST_BASICPROP(numPhiPanels, "number of phi readout panels");
         
-        // TEST_BASICPROP(numGasGapsEta, "numer of eta gas gaps");
-        // chamberOkay = true;
-        TEST_BASICPROP(numGasGapsPhi, "numer of phi gas gaps");
         
-        TEST_BASICPROP(numStripsEta, "numer of eta strips");
-        TEST_BASICPROP(numStripsPhi, "numer of phi strips");
+        TEST_BASICPROP(numStripsEta, "number of eta strips");
+        TEST_BASICPROP(numStripsPhi, "number of phi strips");
         
         TEST_BASICPROP(stripPitchEta, "eta strip pitch");
         TEST_BASICPROP(stripPitchPhi, "phi strip pitch");
@@ -333,22 +350,29 @@ int main( int argc, char** argv ) {
 
         TEST_BASICPROP(stripLengthEta, "eta strip length");
         TEST_BASICPROP(stripLengthPhi, "phi strip length");
-        chamberOkay = true;
         using RpcLayer = RpcChamber::RpcLayer;
         for (const RpcLayer& refLayer : reference.layers) {
-            break;
             std::set<RpcLayer>::const_iterator lay_itr = test.layers.find(refLayer);
             if (lay_itr == test.layers.end()) {
-                std::cerr<<"runRpcGeoComparison() "<<__LINE__<<": in chamber "<<test<<" "
+                std::cerr<<"runRpcGeoComparison() "<<__LINE__<<": "<<test<<" "
                          <<refLayer<<" is not found. "<<std::endl;
                 chamberOkay = false;
                 continue;
             }
+            break;
             const RpcLayer& testLayer{*lay_itr};
             const Amg::Transform3D layAlignment = testLayer.transform.inverse() *
                                                   refLayer.transform;
-            if (Amg::doesNotDeform(layAlignment)) {
-                std::cerr<<"runRpcGeoComparison() "<<__LINE__<<": in chamber "<<test<<" "
+            if (layAlignment.translation().mag() > tolerance) {
+                std::cerr<<"runRpcGeoComparison() "<<__LINE__<<": "<<test<<" "
+                         <<"the layer "<<testLayer<<" is misplaced w.r.t. reference by "
+                         <<Amg::toString(layAlignment)<<std::endl;
+                chamberOkay = false;
+                continue;
+            }
+            continue;            
+            if (!Amg::doesNotDeform(layAlignment)) {
+                std::cerr<<"runRpcGeoComparison() "<<__LINE__<<": "<<test<<" "
                          <<"the layer "<<testLayer<<" is misaligned w.r.t. reference by "
                          <<Amg::toString(layAlignment)<<std::endl;
                 continue;
@@ -359,7 +383,7 @@ int main( int argc, char** argv ) {
         for (const RpcStrip& refStrip : reference.strips) {
             std::set<RpcStrip>::const_iterator strip_itr = test.strips.find(refStrip);
             if (strip_itr == test.strips.end()) {
-                std::cerr<<"runRpcGeoComparison() "<<__LINE__<<": in chamber "<<test<<" "
+                std::cerr<<"runRpcGeoComparison() "<<__LINE__<<": "<<test<<" "
                          <<refStrip<<" is not found. "<<std::endl;
                 chamberOkay = false;
                 continue;
@@ -367,9 +391,10 @@ int main( int argc, char** argv ) {
             const RpcStrip& testStrip{*strip_itr};
             const Amg::Vector3D diffStrip{testStrip.position - refStrip.position};
             if (diffStrip.mag() > tolerance) {
-                std::cerr<<"runRpcGeoComparison() "<<__LINE__<<": in chamber "<<test<<" "
+                std::cerr<<"runRpcGeoComparison() "<<__LINE__<<": "<<test<<" "
                          <<testStrip<<" should be located at "<<Amg::toString(refStrip.position, 2)
-                         <<" displacement: "<<Amg::toString(diffStrip,2)<<std::endl;
+                         <<" displacement: "<<Amg::toString(diffStrip,2)<<"Ich depp "<<diffStrip.perp()<<" "<<diffStrip.z()
+                            <<"  "<<diffStrip.mag()<<std::endl;
                 chamberOkay = false;
             }
         }
diff --git a/MuonSpectrometer/MuonPhaseII/MuonG4/MuonSensitiveDetectorsR4/src/TgcSensitiveDetector.cxx b/MuonSpectrometer/MuonPhaseII/MuonG4/MuonSensitiveDetectorsR4/src/TgcSensitiveDetector.cxx
index a7fb88f06063..16bbd2168a41 100644
--- a/MuonSpectrometer/MuonPhaseII/MuonG4/MuonSensitiveDetectorsR4/src/TgcSensitiveDetector.cxx
+++ b/MuonSpectrometer/MuonPhaseII/MuonG4/MuonSensitiveDetectorsR4/src/TgcSensitiveDetector.cxx
@@ -123,17 +123,18 @@ const MuonGMR4::TgcReadoutElement* TgcSensitiveDetector::getReadoutElement(const
     const std::string stationVolume = touchHist->GetVolume(2)->GetName();
     const std::vector<std::string> volumeTokens = CxxUtils::tokenize(stationVolume, "_");
     /// We should have a string which kind of looks like
-    ///     av_7088_impr_1_MuonR4::T1E1_Station2MuonStation_pv_172_T1E_Station2_-2_4
+    ///     av_319_impr_1_MuonR4::LogVolMuonStation_pv_190_T3E_Station2_-2_46
     /// Of interest are only the T1E part and the last 2 numbers
-    if (volumeTokens.size() != 12) {
-        ATH_MSG_FATAL(__FILE__<<":"<<__LINE__<<" Cannot deduce the station name from "<<stationVolume<<" "<<volumeTokens.size());
-        throw std::runtime_error("Invalid station Identifier");
-    }
+    const size_t nTokens{volumeTokens.size()};
     const TgcIdHelper& idHelper{m_detMgr->idHelperSvc()->tgcIdHelper()};
-    const std::string stationName{volumeTokens[8].substr(0,3)};
-    const int stationEta = CxxUtils::atoi(volumeTokens[10]);
-    const int stationPhi = CxxUtils::atoi(volumeTokens[11]);
-    const Identifier stationId = idHelper.elementID(stationName, stationEta, stationPhi);
+    const std::string stationName{volumeTokens[nTokens-4].substr(0,3)};
+    const int stationEta = CxxUtils::atoi(volumeTokens[nTokens-2]);
+    const int stationPhi = CxxUtils::atoi(volumeTokens[nTokens-1]);
+    bool isValid{false};
+    const Identifier stationId = idHelper.elementID(stationName, stationEta, stationPhi,  isValid);
+    if (!isValid) {
+       throw std::runtime_error("Failed to deduce station name from "+stationVolume);
+    }
     return m_detMgr->getTgcReadoutElement(stationId);
 }
 
-- 
GitLab