MuonReadoutGeomCnvAlg.cxx

/*
   Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
#include "MuonReadoutGeomCnvAlg.h"

#include <GeoPrimitives/GeoPrimitivesHelpers.h>
#include <StoreGate/WriteCondHandle.h>
#include <StoreGate/ReadCondHandle.h>
#include <GeoModelKernel/GeoFullPhysVol.h>

#include <MuonReadoutGeometryR4/MdtReadoutElement.h>
#include <MuonReadoutGeometryR4/RpcReadoutElement.h>
#include <MuonReadoutGeometryR4/TgcReadoutElement.h>
#include <MuonReadoutGeometryR4/MmReadoutElement.h>
#include <MuonReadoutGeometryR4/sTgcReadoutElement.h>


#include <MuonReadoutGeometryR4/MuonChamber.h>

#include <MuonAlignmentDataR4/MdtAlignmentStore.h>
#include <MuonAlignmentDataR4/sTgcAlignmentStore.h>
#include <MuonAlignmentDataR4/MmAlignmentStore.h>


#include <MuonReadoutGeometry/MuonStation.h>
#include <MuonReadoutGeometry/MdtReadoutElement.h>
#include <MuonReadoutGeometry/RpcReadoutElement.h>
#include <MuonReadoutGeometry/sTgcReadoutElement.h>
#include <MuonReadoutGeometry/MMReadoutElement.h>

#include <MuonReadoutGeometry/MuonChannelDesign.h>

#include <AthenaKernel/IOVInfiniteRange.h>
#include <GaudiKernel/SystemOfUnits.h>

#include <GeoModelHelpers/defineWorld.h>
#include <GeoModelHelpers/cloneVolume.h>
#include <GeoModelHelpers/getChildNodesWithTrf.h>
#include <GeoModelHelpers/TransformToStringConverter.h>
#include <GeoModelHelpers/GeoShapeUtils.h>
#include <map>

#include <GaudiKernel/SystemOfUnits.h>

namespace {
    Amg::Transform3D readOutToStation(const GeoVFullPhysVol*  readOutVol) {
        return  readOutVol->getAbsoluteTransform().inverse() *
                readOutVol->getParent()->getX();
    }
    using SubDetAlignment = ActsGeometryContext::AlignmentStorePtr;
}

MuonReadoutGeomCnvAlg::MuonReadoutGeomCnvAlg(const std::string& name, ISvcLocator* pSvcLocator):
    AthReentrantAlgorithm{name, pSvcLocator} {}

StatusCode MuonReadoutGeomCnvAlg::initialize()  {
    ATH_CHECK(m_idHelperSvc.retrieve());
    ATH_CHECK(m_writeKey.initialize());
    ATH_CHECK(m_alignStoreKeys.initialize());
    ATH_CHECK(detStore()->retrieve(m_detMgr));
    return StatusCode::SUCCESS;
}


StatusCode MuonReadoutGeomCnvAlg::execute(const EventContext& ctx) const {
    SG::WriteCondHandle<MuonGM::MuonDetectorManager> writeHandle{m_writeKey, ctx};
    if (writeHandle.isValid()) {
        ATH_MSG_DEBUG("The current readout geometry is still valid.");
        return StatusCode::SUCCESS;
    }
    writeHandle.addDependency(IOVInfiniteRange::infiniteRunLB());
    /// Prepare the Geometry context
    ActsGeometryContext geoContext{};
    using TrackingAlignment = ActsTrk::DetectorAlignStore::TrackingAlignStore;
    for (const SG::ReadCondHandleKey<ActsTrk::DetectorAlignStore>& key : m_alignStoreKeys) {
        SG::ReadCondHandle<ActsTrk::DetectorAlignStore> readHandle{key, ctx};
        if (!readHandle.isValid()) {
            ATH_MSG_FATAL("Failed to retrieve alignment store "<<key.fullKey());
            return StatusCode::FAILURE;
        }
        writeHandle.addDependency(readHandle);
        auto alignStore = std::make_unique<ActsTrk::DetectorAlignStore>(**readHandle);
        /// Ensure that the position & tracking alignment caches are split from the conditions object
        if (alignStore->geoModelAlignment) {
            alignStore->geoModelAlignment->clearPosCache();
        }
        alignStore->trackingAlignment = std::make_unique<TrackingAlignment>(alignStore->detType);
        geoContext.setStore(std::move(alignStore));
    }
    /// Check that for every detector technology there's an DetectorAlignStore in the geometry context
    /// Otherwise create an empty one.
    std::vector<ActsTrk::DetectorType> presentTechs = m_detMgr->getDetectorTypes();
    for (const ActsTrk::DetectorType detType : presentTechs) {
        if (geoContext.getStore(detType)) {
            continue;
        }
        ATH_MSG_WARNING("No external detector alignment has been defined for technology "<<ActsTrk::to_string(detType));
        geoContext.setStore(std::make_unique<ActsTrk::DetectorAlignStore>(detType));
    }

    std::unique_ptr<MuonGM::MuonDetectorManager> detMgr = std::make_unique<MuonGM::MuonDetectorManager>();
    PVLink world{createGeoWorld()};
    detMgr->addTreeTop(world);
    ATH_CHECK(buildMdt(geoContext, detMgr.get(), world));
    ATH_CHECK(buildSTGC(geoContext, detMgr.get(), world));
    ATH_CHECK(buildMM(geoContext, detMgr.get(), world));

    ATH_CHECK(buildRpc(geoContext, detMgr.get(), world));

    ATH_CHECK(writeHandle.record(std::move(detMgr)));
    return StatusCode::SUCCESS;
}
StatusCode MuonReadoutGeomCnvAlg::buildStation(const ActsGeometryContext& gctx,
                                               MuonGM::MuonDetectorManager& mgr,
                                               const Identifier& stationId,
                                               PVLink world) const {
    const std::string stName{m_idHelperSvc->stationNameString(stationId)};
    const int stEta{m_idHelperSvc->stationEta(stationId)};
    const int stPhi{m_idHelperSvc->stationPhi(stationId)};
    MuonGM::MuonStation* station = mgr.getMuonStation(stName, stEta, stPhi);
    if (station) {
        ATH_MSG_DEBUG("Station "<<stName<<" "<<stEta<<" "<<stPhi<<" already exists.");
        return StatusCode::SUCCESS;
    }
    /// Fetch the readout element to get its parent volume
    const MuonGMR4::MuonReadoutElement* copyMe = m_detMgr->getReadoutElement(stationId);
    /// Retrieve the full phyiscal volume
    const GeoVFullPhysVol* readOutVol = copyMe->getMaterialGeom();
    PVConstLink parentVolume = readOutVol->getParent();
    /// Construct the aligned station transformation 
    const Amg::Transform3D stationTransform =  copyMe->localToGlobalTrans(gctx) *
                                                readOutToStation(readOutVol);
    
    /// Copy the full physical volume of the muon station
    PVLink parentPhysVol{make_intrusive<GeoFullPhysVol>(parentVolume->getLogVol())};    
    /// Make sure to copy all the children from the original tree that're not FullPhysVols -> represent
    /// They represent the passive material inside the station and are needed for the TrackinGeometry building
    const std::vector<GeoChildNodeWithTrf> children = getChildrenWithRef(parentVolume, false);
    double minX{1.e9}, maxX{-1.e9}, minY1{1.e9}, maxY1{-1.e9}, minY2{1.e9}, maxY2{-1.e9}, minZ{1.e9}, maxZ{-1.e9};
    for (const GeoChildNodeWithTrf& child : children) {
        GeoVPhysVol* childVol = const_pointer_cast<GeoVPhysVol>(child.volume);
        std::vector<Amg::Vector3D> edges = getPolyShapeEdges(childVol->getLogVol()->getShape(),
                                                             child.transform);
        for (const Amg::Vector3D& edge : edges) {
            minX = std::min(minX, edge.x());
            maxX = std::max(maxX, edge.x());
            minZ = std::min(minZ, edge.z());
            maxZ = std::max(maxZ, edge.z());
            if (edge.z() < 0) {
                minY1 = std::min(minY1, edge.y());
                maxY1 = std::max(maxY1, edge.y());
            } else {
                minY2 = std::min(minY2, edge.y());
                maxY2 = std::max(maxY2, edge.y());
            }
        }

        /// Skip the full physical volumes as they represent the readout elements
        if (typeid(*childVol) == typeid(GeoFullPhysVol)) {
            continue;
        }
        world->add(make_intrusive<GeoTransform>(child.transform));
        world->add(cloneVolume(childVol));
    }
    /// Add the physical volume to the world
    world->add(make_intrusive<GeoTransform>(stationTransform));
    world->add(parentPhysVol);
    /// To create the muon station, we need to extract the dimensions
    ///  --> Recieve the edge points from the shapes
    const double shortS = (maxY1 - minY1);
    const double longS  = (maxY2 - minY2);
    const double lengthR = (maxX - minX);
    const double lengthZ = (maxZ - minZ);
    auto newStation = std::make_unique<MuonGM::MuonStation>(stName,
                                                            shortS, lengthR, lengthZ, /// S / R / Z size
                                                            longS, lengthR, lengthZ,  /// S / R / Z size (long)
                                                            stEta, stPhi, false);
    newStation->setPhysVol(parentPhysVol);
    mgr.addMuonStation(std::move(newStation));

    return StatusCode::SUCCESS;
}

StatusCode MuonReadoutGeomCnvAlg::buildRpc(const ActsGeometryContext& gctx,
                                           MuonGM::MuonDetectorManager* mgr,
                                           PVLink world) const {
    
    const std::vector<const MuonGMR4::RpcReadoutElement*> readoutEles = m_detMgr->getAllRpcReadoutElements();
    ATH_MSG_DEBUG("Going to build "<<readoutEles.size()<<" Rpc readout elements.");
    const RpcIdHelper& idHelper{m_idHelperSvc->rpcIdHelper()};
    for (const MuonGMR4::RpcReadoutElement* copyMe : readoutEles) {
        const Identifier reId = copyMe->identify();
        const MuonGMR4::RpcReadoutElement::parameterBook& pars{copyMe->getParameters()};
        /// Build the mother station if it's not already existing
        ATH_CHECK(buildStation(gctx, *mgr, reId, world));

        const std::string stName{m_idHelperSvc->stationNameString(reId)};
        MuonGM::MuonStation* station = mgr->getMuonStation(stName, 
                                                           m_idHelperSvc->stationEta(reId), 
                                                           m_idHelperSvc->stationPhi(reId));
        
        PVLink parentPhysVol{station->getPhysVol()};
        GeoIntrusivePtr<const GeoVFullPhysVol> readOutVol{copyMe->getMaterialGeom()};
        parentPhysVol->add(make_intrusive<GeoTransform>(readOutToStation(readOutVol).inverse()));
        PVLink clonedVol{cloneVolume(const_pointer_cast<GeoVFullPhysVol>(readOutVol))};
        GeoIntrusivePtr<GeoVFullPhysVol> physVol{dynamic_pointer_cast<GeoVFullPhysVol>(clonedVol)};
        parentPhysVol->add(physVol);

        std::unique_ptr<MuonGM::RpcReadoutElement> newElement = std::make_unique<MuonGM::RpcReadoutElement>(physVol, stName, 1, 1, false, mgr);
        const bool aSide{copyMe->stationEta() > 0};
        newElement->setDoubletPhi(copyMe->doubletPhi());
        newElement->setDoubletR(copyMe->doubletR());
        newElement->setDoubletZ(copyMe->doubletZ());
        newElement->setIdentifier(reId);
        newElement->setParentMuonStation(station);

        /// Define the dimensions
        newElement->setLongRsize(pars.halfLength);
        newElement->setLongSsize(pars.halfWidth);
        newElement->setLongZsize(pars.halfThickness);
        newElement->setRsize(pars.halfLength);
        newElement->setSsize(pars.halfWidth);
        newElement->setZsize(pars.halfThickness);

        newElement->m_nlayers = copyMe->nGasGaps();
        newElement->m_phistripwidth = copyMe->stripPhiWidth();
        newElement->m_etastripwidth = copyMe->stripEtaWidth();
        newElement->m_phistrippitch = copyMe->stripPhiPitch();
        newElement->m_etastrippitch =  (aSide > 0 ? 1. : -1.) *copyMe->stripEtaPitch();
        newElement->m_phistriplength = copyMe->stripPhiLength();
        newElement->m_etastriplength = copyMe->stripEtaLength();

        newElement->m_nphistripsperpanel = copyMe->nPhiStrips();
        newElement->m_netastripsperpanel = copyMe->nEtaStrips();
        newElement->m_nphistrippanels = copyMe->nPhiPanels();
        newElement->m_hasDEDontop = true;
        newElement->m_descratzneg = false;

         for (unsigned int gasGap = 1; gasGap <= copyMe->nGasGaps(); ++gasGap) {
            for (int doubPhi = copyMe->doubletPhiMax(); doubPhi >= copyMe->doubletPhi(); --doubPhi) {
                for (bool measPhi : {false, true}) {
                    const int channel = 1;
                    const Identifier gapId = idHelper.channelID(copyMe->identify(), 
                                                                copyMe->doubletZ(), 
                                                                doubPhi, gasGap, measPhi, 
                                                                channel);

                    const Amg::Vector3D locStripPos = copyMe->globalToLocalTrans(gctx) * copyMe->stripPosition(gctx, gapId);
                    ATH_MSG_VERBOSE("GasGap "<<m_idHelperSvc->toString(gapId)<<", local strip position: "<<Amg::toString(locStripPos));
                    newElement->m_gasGap_xPos[gasGap -1] = locStripPos.x();
                    /// Hack to assign the proper strip positions  for REs having doubletPhi =2
                    /// in their Identifier
                    const int dbPIdx = copyMe->doubletPhi() == 2 ? 1 : doubPhi;
                    if (measPhi) {
                        newElement->m_first_phistrip_s[dbPIdx -1] = locStripPos.y();
                        newElement->m_phistrip_z = locStripPos.z();
                    } else{
                        newElement->m_first_etastrip_z = locStripPos.z();
                        newElement->m_etastrip_s[dbPIdx-1] = locStripPos.y();
                    }
                }
            }
        }
        newElement->m_mirrored = true;
        newElement->fillCache();
        newElement->m_mirrored = false;
        ATH_CHECK(dumpAndCompare(gctx, *copyMe, *newElement));
        mgr->addRpcReadoutElement(std::move(newElement));
    }
    return StatusCode::SUCCESS;
}

StatusCode MuonReadoutGeomCnvAlg::dumpAndCompare(const ActsGeometryContext& gctx,
                                                 const MuonGMR4::MmReadoutElement& refEle,
                                                 const MuonGM::MMReadoutElement& testEle) const {
    if (!m_checkGeo) {
        return StatusCode::SUCCESS;
    }
    ATH_MSG_VERBOSE("Compare basic readout transforms"<<std::endl
                <<GeoTrf::toString(testEle.absTransform(),true)<<std::endl
                <<GeoTrf::toString(refEle.localToGlobalTrans(gctx), true));
    const MmIdHelper& idHelper{m_idHelperSvc->mmIdHelper()};
    for (unsigned int gasGap = 1; gasGap <= refEle.nGasGaps(); ++ gasGap) {
        const Identifier gapId = idHelper.channelID(refEle.identify(), refEle.multilayer(),  gasGap, 1);
        
        const Amg::Transform3D& refTrf{refEle.localToGlobalTrans(gctx, gapId)};
        const Amg::Transform3D& testTrf{testEle.transform(gapId)};
        if (!Amg::doesNotDeform(refTrf.inverse()*testTrf)) {
            ATH_MSG_FATAL("The layer "<<m_idHelperSvc->toStringGasGap(gapId)<<" does not transform equally"
                         <<GeoTrf::toString(refTrf, true) <<" vs. "<<GeoTrf::toString(testTrf, true));
            return StatusCode::FAILURE;
        }
        const MuonGMR4::StripDesign& stripDesign{refEle.stripLayer(gapId).design()};
        
        for (int strip = stripDesign.firstStripNumber(); strip <= stripDesign.numStrips(); ++strip) {
            const Identifier stripId = idHelper.channelID(refEle.identify(), refEle.multilayer(), gasGap, strip);
            const Amg::Vector3D refStripPos{refEle.stripPosition(gctx, stripId)};
            const Amg::Vector3D refStripDir{refEle.localToGlobalTrans(gctx, refEle.layerHash(stripId)).linear() * Amg::Vector3D::UnitX()};

            Amg::Vector3D testStripPos{Amg::Vector3D::Zero()};
            if (!testEle.stripGlobalPosition(stripId, testStripPos)) {
                ATH_MSG_FATAL("Failed to retrieve strip position "<<m_idHelperSvc->toString(stripId));
                return StatusCode::FAILURE;
            }
            const double dist = refStripDir.dot(refStripPos - testStripPos);
            if (std::abs(dist) > 10. * Gaudi::Units::micrometer) {
                ATH_MSG_FATAL("The strip "<<Amg::toString(testStripPos)<<" is not describing the same strip as "
                            <<Amg::toString(refStripPos)<<". Channel "<<m_idHelperSvc->toString(stripId)
                            <<" distance: "<<dist<<" "<<(dist / testEle.m_etaDesign[gasGap -1].inputWidth));
                return StatusCode::FAILURE;
            }
            ATH_MSG_VERBOSE("Channel postion "<<m_idHelperSvc->toString(stripId)<<" match between legacy & new");
        }
    }
    return StatusCode::SUCCESS;
}

StatusCode MuonReadoutGeomCnvAlg::dumpAndCompare(const ActsGeometryContext& gctx,
                                                 const MuonGMR4::RpcReadoutElement& refEle,
                                                 const MuonGM::RpcReadoutElement& testEle) const {
    if (!m_checkGeo) {
        return StatusCode::SUCCESS;
    }
    ATH_MSG_VERBOSE("Compare basic readout transforms"<<std::endl
                <<GeoTrf::toString(testEle.absTransform(),true)<<std::endl
                <<GeoTrf::toString(refEle.localToGlobalTrans(gctx), true));
    const RpcIdHelper& idHelper{m_idHelperSvc->rpcIdHelper()};
    for (unsigned int gasGap = 1; gasGap <= refEle.nGasGaps(); ++gasGap) {
        for (int doubPhi = refEle.doubletPhi(); doubPhi <= refEle.doubletPhiMax(); ++doubPhi) {
            for (bool measPhi : {false, true}) {
                for (int strip = 1; strip <= testEle.Nstrips(measPhi); ++strip) {
                    const Identifier stripId = idHelper.channelID(refEle.identify(), 
                                                                refEle.doubletZ(), 
                                                                doubPhi, gasGap, measPhi, strip);
                    
                    const Amg::Transform3D& refTrans{refEle.localToGlobalTrans(gctx, stripId)};
                    const Amg::Transform3D& testTrans{testEle.transform(stripId)};
                    if (strip == 1 && Amg::doesNotDeform(refTrans.inverse()*testTrans)) {
                        ATH_MSG_ERROR("Transformation for "<<m_idHelperSvc->toString(stripId)<<std::endl
                            <<" *** ref:  "<<GeoTrf::toString(refTrans, true)<<std::endl
                            <<" *** test: "<<GeoTrf::toString(testTrans, true));
                            return StatusCode::FAILURE;
                    }

                    const Amg::Vector3D refStripPos = refEle.stripPosition(gctx, stripId);
                    const Amg::Vector3D testStripPos = testEle.stripPos(stripId);
                    if ((refStripPos - testStripPos).mag() > std::numeric_limits<float>::epsilon()){
                        ATH_MSG_ERROR("Mismatch in strip positions "<<m_idHelperSvc->toString(stripId)
                                <<" ref: "<<Amg::toString(refStripPos)<<" test: "<<Amg::toString(testStripPos)
                                <<" local coordinates -- ref: "<<Amg::toString(testEle.absTransform().inverse()*refStripPos)
                                <<" test: "<<Amg::toString(testEle.absTransform().inverse()*testStripPos));
                        return StatusCode::FAILURE;
                    }
                    ATH_MSG_VERBOSE("Agreement between new and old geometry for channel "<<m_idHelperSvc->toString(stripId)
                                    <<" strip position "<<Amg::toString(refStripPos));
                }
            }
        }
    }
    return StatusCode::SUCCESS;
}

StatusCode MuonReadoutGeomCnvAlg::buildMM(const ActsGeometryContext& gctx,
                                          MuonGM::MuonDetectorManager* mgr,
                                          PVLink world) const {

    SubDetAlignment alignItr = gctx.getStore(ActsTrk::DetectorType::Mm);
    const auto alignStore = alignItr ?
                            static_cast<const MmAlignmentStore*>(alignItr->internalAlignment.get()) : nullptr;

    const std::vector<const MuonGMR4::MmReadoutElement*> mmReadouts{m_detMgr->getAllMmReadoutElements()};
    for (const MuonGMR4::MmReadoutElement* copyMe : mmReadouts) {
        const Identifier reId = copyMe->identify();
        GeoIntrusivePtr<const GeoVFullPhysVol> readOutVol{copyMe->getMaterialGeom()};
        PVLink clonedVol{cloneVolume(const_pointer_cast<GeoVFullPhysVol>(readOutVol))};
        GeoIntrusivePtr<GeoFullPhysVol> physVol{dynamic_pointer_cast<GeoFullPhysVol>(clonedVol)};
        world->add(make_intrusive<GeoTransform>(copyMe->localToGlobalTrans(gctx)));
        world->add(physVol);

        auto newRE = std::make_unique<MuonGM::MMReadoutElement>(physVol, 
                                                                m_idHelperSvc->stationNameString(reId),
                                                                copyMe->stationEta(),
                                                                copyMe->stationPhi(),
                                                                copyMe->multilayer(), mgr,
                                                                alignStore ? alignStore->passivation : nullptr);
        /// Loop over the gas gaps & efine the 
        for (unsigned int gasGap = 0; gasGap < copyMe->nGasGaps(); ++gasGap) {
            const MuonGMR4::StripLayer& stripLayer{copyMe->stripLayer(MuonGMR4::MmReadoutElement::createHash(gasGap +1, 0))};
            const MuonGMR4::StripDesign& designFrom{stripLayer.design()};
            
            newRE->m_Xlg[gasGap] = stripLayer.toOrigin() * 
                                   Amg::getRotateZ3D(-designFrom.stereoAngle()) * 
                                   Amg::getRotateY3D(90. * Gaudi::Units::deg);
            ATH_MSG_VERBOSE("Layer transform "<<gasGap<<" "<<GeoTrf::toString(newRE->m_Xlg[gasGap], true));
            
            MuonGM::MuonChannelDesign& designTo{newRE->m_etaDesign[gasGap]};
            designTo.defineTrapezoid(designFrom.shortHalfHeight(),
                                     designFrom.longHalfHeight(),
                                     designFrom.halfWidth(), 
                                     designFrom.stereoAngle());
            designTo.type = MuonGM::MuonChannelDesign::ChannelType::etaStrip;
            designTo.detType = MuonGM::MuonChannelDesign::DetType::MM;
            designTo.inputPitch = designFrom.stripPitch();
            designTo.inputWidth = designTo.inputPitch * std::cos(designTo.stereoAngle());
            designTo.nMissedBottomEta  = designTo.nMissedBottomStereo = designFrom.firstStripNumber() - 1;
            designTo.totalStrips = designFrom.numStrips();
            designTo.nch = designFrom.numStrips();
            
            designTo.setFirstPos(designFrom.firstStripPos().x() + 0.5*designTo.inputPitch);
        }

        newRE->fillCache();
        if (alignStore && alignStore->getBLine(reId)) {
            newRE->setBLinePar(*alignStore->getBLine(reId));
        }
        ATH_CHECK(dumpAndCompare(gctx, *copyMe, *newRE));
        mgr->addMMReadoutElement(std::move(newRE));

    }
    return StatusCode::SUCCESS;
}

StatusCode  MuonReadoutGeomCnvAlg::buildSTGC(const ActsGeometryContext& gctx,
                                             MuonGM::MuonDetectorManager* mgr,
                                             PVLink world) const{
    SubDetAlignment alignItr = gctx.getStore(ActsTrk::DetectorType::sTgc);
    auto alignStore = alignItr ? static_cast<const sTgcAlignmentStore*>(alignItr->internalAlignment.get()) : nullptr;

    const std::vector<const MuonGMR4::sTgcReadoutElement*> sTgcReadOuts{m_detMgr->getAllsTgcReadoutElements()};
    for (const MuonGMR4::sTgcReadoutElement* copyMe : sTgcReadOuts) {
        const Identifier reId = copyMe->identify();
        GeoIntrusivePtr<const GeoVFullPhysVol> readOutVol{copyMe->getMaterialGeom()};
        PVLink clonedVol{cloneVolume(const_pointer_cast<GeoVFullPhysVol>(readOutVol))};
        GeoIntrusivePtr<GeoFullPhysVol> physVol{dynamic_pointer_cast<GeoFullPhysVol>(clonedVol)};
        world->add(make_intrusive<GeoTransform>(copyMe->localToGlobalTrans(gctx)));
        world->add(physVol);

        auto newRE = std::make_unique<MuonGM::sTgcReadoutElement>(physVol, 
                                                                  m_idHelperSvc->stationNameString(reId).substr(1),
                                                                  copyMe->stationEta(),
                                                                  copyMe->stationPhi(),
                                                                  copyMe->multilayer(), mgr);
        
        if (alignStore && alignStore->getBLine(reId)) {
            newRE->setBLinePar(*alignStore->getBLine(reId));
        }
        for (unsigned int layer = 1; layer <= copyMe->numLayers(); ++layer) {
            using channelType = MuonGMR4::sTgcReadoutElement::ReadoutChannelType;
            using ChannelDesign =  MuonGM::MuonChannelDesign;
            const IdentifierHash layerHash = MuonGMR4::sTgcReadoutElement::createHash(layer,channelType::Strip,0);
            newRE->m_Xlg[layer -1] =  Amg::getTranslate3D(copyMe->globalToLocalTrans(gctx) *
                                                          copyMe->center(gctx, layerHash));

            const MuonGMR4::StripDesign& copyEtaDesign{copyMe->stripDesign(layerHash)}; 
            /// Initialize the eta design
            ChannelDesign& etaDesign{newRE->m_etaDesign[layer-1]};
            etaDesign.type =  ChannelDesign::ChannelType::etaStrip;
            etaDesign.detType = ChannelDesign::DetType::STGC;
            if (copyEtaDesign.yCutout()) {
                etaDesign.defineDiamond(copyEtaDesign.shortHalfHeight(),
                                        copyEtaDesign.longHalfHeight(),
                                        copyEtaDesign.halfWidth(),
                                        copyEtaDesign.yCutout());
            } else {
                etaDesign.defineTrapezoid(copyEtaDesign.shortHalfHeight(),
                                          copyEtaDesign.longHalfHeight(),
                                          copyEtaDesign.halfWidth());
            }
            etaDesign.inputPitch  = copyEtaDesign.stripPitch();
            etaDesign.inputWidth  = copyEtaDesign.stripWidth();
            etaDesign.nch = copyEtaDesign.numStrips();
            etaDesign.setFirstPos(copyEtaDesign.firstStripPos().x());
            /// Initialize the phi design

            const MuonGMR4::WireGroupDesign& copyPhiDesign{copyMe->wireDesign(layerHash)};
            
            ChannelDesign& phiDesign{newRE->m_phiDesign[layer-1]};
            phiDesign.type = ChannelDesign::ChannelType::phiStrip;
            phiDesign.detType = ChannelDesign::DetType::STGC;
            if (copyPhiDesign.yCutout() == 0.) {
                phiDesign.defineTrapezoid(copyPhiDesign.shortHalfHeight(),
                                          copyPhiDesign.longHalfHeight(),
                                          copyPhiDesign.halfWidth());
              } else { 
                phiDesign.defineDiamond(copyPhiDesign.shortHalfHeight(),
                                        copyPhiDesign.longHalfHeight(),
                                        copyPhiDesign.halfWidth(), 
                                        copyPhiDesign.yCutout());
              }
              phiDesign.inputPitch  = copyPhiDesign.stripPitch();
              phiDesign.inputWidth  = 0.015;
              phiDesign.setFirstPos(copyPhiDesign.firstStripPos().x()); // Position of 1st wire, accounts for staggering
            //   phiDesign.firstPitch = firstWireGroup[il];             // Number of Wires in 1st group, group staggering
              //phiDesign.groupWidth  = wireGroupWidth;                // Number of Wires normal group
              phiDesign.nGroups = copyPhiDesign.numStrips();                           // Number of Wire Groups
              phiDesign.wireCutout = copyPhiDesign.wireCutout();                       // Size of "active" wire region for digits
              phiDesign.nch = copyPhiDesign.nAllWires();

              const MuonGMR4::PadDesign& copyPadDesign{copyMe->padDesign(layerHash)};
              MuonGM::MuonPadDesign& padDesign{newRE->m_padDesign[layer-1]};
              padDesign.sPadWidth = 2.*copyPadDesign.shortHalfHeight();
              padDesign.lPadWidth = 2.*copyPadDesign.longHalfHeight();
              padDesign.Size =  2.*copyPadDesign.halfWidth();
        }     
        newRE->fillCache();
        mgr->addsTgcReadoutElement(std::move(newRE));
   }
    
    
    return StatusCode::SUCCESS;
}
StatusCode MuonReadoutGeomCnvAlg::buildMdt(const ActsGeometryContext& gctx,
                                           MuonGM::MuonDetectorManager* mgr,
                                           PVLink world) const {    
    /// Access the B-Line and As-built parameters
    SubDetAlignment alignItr = gctx.getStore(ActsTrk::DetectorType::Mdt);
    const MdtAlignmentStore* alignStore = alignItr ?
                                 static_cast<const MdtAlignmentStore*>(alignItr->internalAlignment.get()) : nullptr;

    const std::vector<const MuonGMR4::MdtReadoutElement*> mdtReadOuts{m_detMgr->getAllMdtReadoutElements()};
    ATH_MSG_INFO("Copy "<<mdtReadOuts.size()<<" Mdt readout elements to the legacy system");
    for (const MuonGMR4::MdtReadoutElement* copyMe : mdtReadOuts) {
        const Identifier reId = copyMe->identify();
        /// Build the mother station
        ATH_CHECK(buildStation(gctx, *mgr, reId, world));

        const std::string stName{m_idHelperSvc->stationNameString(reId)};
        MuonGM::MuonStation* station = mgr->getMuonStation(stName, m_idHelperSvc->stationEta(reId), m_idHelperSvc->stationPhi(reId));
        
        // cppcheck-suppress invalidLifetime; ok: mgr took ownership.
        PVLink parentPhysVol{station->getPhysVol()};
        GeoIntrusivePtr<const GeoVFullPhysVol> readOutVol{copyMe->getMaterialGeom()};
        PVLink clonedVol{cloneVolume(const_pointer_cast<GeoVFullPhysVol>(readOutVol))};
        GeoIntrusivePtr<GeoFullPhysVol> physVol{dynamic_pointer_cast<GeoFullPhysVol>(clonedVol)};
        parentPhysVol->add(make_intrusive<GeoTransform>(readOutToStation(readOutVol).inverse()));
        parentPhysVol->add(physVol);

        const MuonGMR4::MdtReadoutElement::parameterBook& pars{copyMe->getParameters()};

        std::unique_ptr<MuonGM::MdtReadoutElement> newElement = std::make_unique<MuonGM::MdtReadoutElement>(physVol, stName, mgr);
        newElement->setIdentifier(reId);
        newElement->setMultilayer(copyMe->multilayer());
        // cppcheck-suppress invalidLifetime; ok: mgr took ownership.
        newElement->setParentMuonStation(station);
        /// Define the dimensions
        newElement->setLongRsize(2*pars.halfY);
        /// 1 cm is added as safety margin to the Mdt multilayer envelope
        newElement->setLongSsize(2*pars.longHalfX - 1.*Gaudi::Units::cm);
        newElement->setLongZsize(2*pars.halfHeight);
        newElement->setRsize(2*pars.halfY);
        newElement->setSsize(2*pars.shortHalfX - 1.*Gaudi::Units::cm);
        newElement->setZsize(2*pars.halfHeight);

        newElement->m_nlayers = copyMe->numLayers();
        newElement->m_ntubesperlayer = copyMe->numTubesInLay();
        newElement->m_deadlength = pars.deadLength;
        newElement->m_endpluglength = pars.endPlugLength;
        newElement->m_innerRadius = pars.tubeInnerRad;
        newElement->m_tubeWallThickness = pars.tubeWall;
        newElement->m_tubepitch = pars.tubePitch;
        /// Need to check how to obtain this parameter from the new geometry
        /// newElement->m_cutoutShift;

        /// Determine the tube length's 
        const MuonGMR4::MdtTubeLayer& tubeLay{*pars.tubeLayers[0]};
        unsigned int step{1};
        double lastLength{2.*tubeLay.tubeHalfLength(1)}; 
        for (unsigned tube = 0; tube < copyMe->numTubesInLay(); ++tube) {
            const double currLength = 2.*tubeLay.tubeHalfLength(tube);
            if (std::abs(lastLength - currLength) > std::numeric_limits<float>::epsilon() ||
                tube == copyMe->numTubesInLay() -1) {
                newElement->m_tubelength[step-1] = lastLength;
                newElement->m_tubelength[step] = currLength;                
                if (step == 1) {
                    newElement->m_ntubesinastep = tube;
                }
                lastLength = currLength;
                ++step;
            }
        }
        newElement->m_nsteps = step;
        
        /// Define the tube staggering
        const Amg::Transform3D globToLoc{copyMe->globalToLocalTrans(gctx)};
        double xOffSet{pars.halfY}, yOffSet{pars.halfHeight};
        if (newElement->barrel())  std::swap(xOffSet, yOffSet);
        for (unsigned lay = 1; lay <= copyMe->numLayers(); ++lay) {
            const IdentifierHash tubeHash{copyMe->measurementHash(lay, 1)};
            const Amg::Vector3D locTube = globToLoc * copyMe->globalTubePos(gctx, tubeHash);
            newElement->m_firstwire_x[lay-1] = locTube.z() + xOffSet;
            newElement->m_firstwire_y[lay-1] = locTube.x() + yOffSet;
        }
        MdtAlignmentStore::chamberDistortions distort = alignStore ? alignStore->getDistortion(reId) : 
                                                        MdtAlignmentStore::chamberDistortions{};
        
        newElement->setBLinePar(distort.bLine);
        station->setMdtAsBuiltParams(distort.asBuilt);
        newElement->geoInitDone();
        newElement->fillCache();
        /// Add the readout element to the manager
        ATH_CHECK(dumpAndCompare(gctx, *copyMe, *newElement));
        mgr->addMdtReadoutElement(std::move(newElement));
    }
    return StatusCode::SUCCESS;
}

StatusCode MuonReadoutGeomCnvAlg::dumpAndCompare(const ActsGeometryContext& gctx,
                                                 const MuonGMR4::MdtReadoutElement& refEle,
                                                 const MuonGM::MdtReadoutElement& testEle) const {
    if (!m_checkGeo) {
        return StatusCode::SUCCESS;
    }
    
    ATH_MSG_VERBOSE("Detector element "<<m_idHelperSvc->toString(refEle.identify())
                <<std::endl<<GeoTrf::toString(refEle.localToGlobalTrans(gctx))                        
                <<std::endl<<GeoTrf::toString(testEle.getMaterialGeom()->getAbsoluteTransform())
                <<std::endl<<"r-size: "<<testEle.getRsize()<<"/"<<testEle.getLongRsize()
                            <<" s-size: "<<testEle.getSsize()<<"/"<<testEle.getLongSsize()
                            <<" z-size: "<<testEle.getZsize()<<"/"<<testEle.getLongZsize());
    for (unsigned int lay = 1; lay <= refEle.numLayers(); ++lay){
        for (unsigned int tube = 1; tube <= refEle.numTubesInLay(); ++tube) {
            const IdentifierHash tubeHash {refEle.measurementHash(lay,tube)};
            if (!refEle.isValid(tubeHash)) {
                ATH_MSG_VERBOSE("SKip layer / tube "<<lay <<","<<tube);
                continue;
            }
            const Amg::Transform3D globToLocal = refEle.globalToLocalTrans(gctx, tubeHash);

            const Amg::Vector3D refPos = refEle.globalTubePos(gctx, tubeHash);
            const Amg::Vector3D tubePos = testEle.tubePos(lay, tube);

            
            if ( (refPos - tubePos).mag() > Gaudi::Units::micrometer &&
                 (globToLocal*refPos - globToLocal * tubePos).perp() > Gaudi::Units::micrometer) {
                ATH_MSG_ERROR("Tube positions differ for "<<m_idHelperSvc->toString(refEle.measurementId(tubeHash))
                            <<" reference: "<<GeoTrf::toString(refPos)<<" vs. test: "
                            <<GeoTrf::toString(tubePos) <<" delta: "<<(refPos - tubePos).mag());
                return StatusCode::FAILURE;
            }
            ATH_MSG_VERBOSE("Tube positions layer: "<<lay<<", tube: "<<tube
                <<std::endl<<"reference: "<<GeoTrf::toString(refPos)
                <<std::endl<<"test:      "<<GeoTrf::toString(tubePos)
                <<std::endl<<testEle.tubeLength(lay, tube)<<"/"
                            <<testEle.getActiveTubeLength(lay, tube)<<"/"
                            <<testEle.getWireLength(lay,tube)
                            <<" vs. "<<refEle.tubeLength(tubeHash)<<"/"<<refEle.activeTubeLength(tubeHash)
                            <<"/"<<refEle.wireLength(tubeHash));
            if (std::abs(testEle.tubeLength(lay,tube) - refEle.tubeLength(tubeHash)) >
                std::numeric_limits<float>::epsilon() ) {
                ATH_MSG_WARNING("Different tube length's detected for "<<m_idHelperSvc->toStringDetEl(refEle.identify())
                                << " layer: "<<lay<<", tube: "<<tube<<" "<<testEle.tubeLength(lay,tube)<<" (new) vs. "
                                <<refEle.tubeLength(tubeHash)<<" (ref)");
            }
        }
    }

    return StatusCode::SUCCESS;
}