MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelMdtTest.cxx

Go to the documentation of this file.

 
/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
#include "GeoModelMdtTest.h"
#include <ActsGeometryInterfaces/GeometryContext.h>
#include <MuonReadoutGeometryR4/MdtReadoutElement.h>
#include <MuonReadoutGeometryR4/SpectrometerSector.h>
#include <GeoPrimitives/GeoPrimitivesToStringConverter.h>
#include <MuonVisualizationHelpersR4/VisualizationHelpers.h>
 
#include <fstream>
 
using namespace ActsTrk;
 
namespace MuonGMR4{
 
StatusCode GeoModelMdtTest::initialize() {
    ATH_CHECK(m_idHelperSvc.retrieve());
    ATH_CHECK(m_geoCtxKey.initialize()); 
    ATH_CHECK(m_cablingKey.initialize(!m_cablingKey.empty()));   
    ATH_CHECK(m_tree.init(this));
    if(m_visualTubes || m_visualStaggering) {
        ATH_CHECK(m_visualSvc.retrieve());
    }
    if (m_visualTubes) {
        m_clientTokenLayerVis.preFixName="GeoModelMdtTest_Layer";
        m_clientTokenLayerVis.subDirectory = "MdtLayerPlots";
        m_clientTokenLayerVis.canvasLimit = -1;
        ATH_CHECK(m_visualSvc->registerClient(m_clientTokenLayerVis));
    }
    if(m_visualStaggering) {
        m_clientTokenStaggeringVis.preFixName="GeoModelMdtTest_Staggering";
        m_clientTokenStaggeringVis.subDirectory = "MdtStaggeringPlots";
        m_clientTokenStaggeringVis.canvasLimit = -1;
        ATH_CHECK(m_visualSvc->registerClient(m_clientTokenStaggeringVis));
    }
    const MdtIdHelper& idHelper{m_idHelperSvc->mdtIdHelper()};
    auto translateTokenList = [this, &idHelper](const std::vector<std::string>& chNames){
        
        std::set<Identifier> transcriptedIds{};
        for (const std::string& token : chNames) { 
            if (token.size() != 6) {
                ATH_MSG_WARNING("Wrong format given for "<<token<<". Expecting 6 characters");
                continue;
            }
            const std::string statName = token.substr(0, 3);
            const unsigned statEta = std::atoi(token.substr(3, 1).c_str()) * (token[4] == 'A' ? 1 : -1);
            const unsigned statPhi = std::atoi(token.substr(5, 1).c_str());
            bool isValid{false};
            const Identifier eleId = idHelper.elementID(statName, statEta, statPhi, isValid);
            if (!isValid) {
                ATH_MSG_WARNING("Failed to deduce a station name for " << token);
                continue;
            }
            transcriptedIds.insert(eleId);
            const Identifier secMlId = idHelper.multilayerID(eleId, 2, isValid);
            if (isValid){
                transcriptedIds.insert(secMlId);
            }
        }
        return transcriptedIds;
    };
 
    std::vector <std::string>& selectedSt = m_selectStat.value();
    const std::vector <std::string>& excludedSt = m_excludeStat.value();
    selectedSt.erase(std::remove_if(selectedSt.begin(), selectedSt.end(),
                     [&excludedSt](const std::string& token){
                        return std::ranges::find(excludedSt, token) != excludedSt.end();
                     }), selectedSt.end());
 
    if (selectedSt.size()) {
        m_testStations = translateTokenList(selectedSt);
        std::stringstream sstr{};
        for (const Identifier& id : m_testStations) {
            sstr<<" *** "<<m_idHelperSvc->toString(id)<<std::endl;
        }
        ATH_MSG_INFO("Test only the following stations "<<std::endl<<sstr.str());
    } else {
        const std::set<Identifier> excluded = translateTokenList(excludedSt);
        for(auto itr = idHelper.detectorElement_begin();
                 itr!= idHelper.detectorElement_end();++itr){
            if (!excluded.count(*itr)) {
               m_testStations.insert(*itr);
            }
        }
        if (!excluded.empty()) {
            std::stringstream excluded_report{};
            for (const Identifier& id : excluded){
                excluded_report << " *** " << m_idHelperSvc->toStringDetEl(id) << std::endl;
            }
            ATH_MSG_INFO("Test all station except the following excluded ones " << std::endl << excluded_report.str());
        }
    }
    ATH_CHECK(detStore()->retrieve(m_detMgr));
    return StatusCode::SUCCESS;
}
void GeoModelMdtTest::visualizeTubeLayer(const EventContext& ctx,
                                         const MuonGMR4::MdtReadoutElement& reEle,
                                         const unsigned layer) const {
    if (!m_visualTubes) {
        return;
    }
     
    const double h = reEle.moduleHeight();
    const double wL = reEle.moduleWidthL();
    const double wS = reEle.moduleWidthS();
 
    const MdtIdHelper& idHelper{m_idHelperSvc->mdtIdHelper()};
    const Identifier detId = reEle.identify();
    const std::string chName = std::format("{:}{:}{:}{:}M{:}T{:}",
                                           m_idHelperSvc->stationNameString(detId),
                                           std::abs(m_idHelperSvc->stationEta(detId)),
                                           m_idHelperSvc->stationEta(detId) > 0 ? 'A' : 'C',
                                           m_idHelperSvc->stationPhi(detId),
                                           idHelper.multilayer(detId), layer);
    auto canvas = m_visualSvc->prepareCanvas(ctx, m_clientTokenLayerVis, chName);
    canvas->expandPad(-0.5*wL, -0.5*h);
    canvas->expandPad( 0.5*wL, 0.5*h);
    canvas->setAxisTitles("x [mm]", "y [mm]");
    canvas->setRangeScale(1.1);
    using namespace MuonValR4;
    canvas->add(drawLine(Amg::Vector3D{0., -0.5*wL, 0.5*h}, Amg::Vector3D{0., 0.5*wL, 0.5*h}, kBlack, kSolid));
    canvas->add(drawLine(Amg::Vector3D{0., -0.5*wS, -0.5*h}, Amg::Vector3D{0., 0.5*wS, -0.5*h}, kBlack, kSolid));
    canvas->add(drawLine(Amg::Vector3D{0., -0.5*wL, 0.5*h}, Amg::Vector3D{0., -0.5*wS, -0.5*h}, kBlack, kSolid));
    canvas->add(drawLine(Amg::Vector3D{0.,  0.5*wL, 0.5*h}, Amg::Vector3D{0., 0.5*wS, -0.5*h}, kBlack, kSolid));
    for (unsigned int tube = 1 ; tube <= reEle.numTubesInLay(); ++ tube) {
        const IdentifierHash measHash = reEle.measurementHash(layer, tube);
        if (!reEle.isValid(measHash)) {
            continue;
        }
        canvas->add(drawBox(reEle.localTubePos(measHash), reEle.tubeLength(measHash), reEle.tubeRadius()));
    } 
    canvas->add(drawLabel(std::format("{:}, layer: {:}", m_idHelperSvc->toStringDetEl(detId), layer), 0.2, 0.05));
}
 
StatusCode GeoModelMdtTest::visualizeTubeStaggering(const EventContext& ctx, const ActsTrk::GeometryContext& gctx) const {
    const MuonMDT_CablingMap* cabling{nullptr};
    ATH_CHECK(SG::get(cabling, m_cablingKey,ctx));
    if (!cabling){
        return StatusCode::SUCCESS;
    }
 
    std::unordered_map<Identifier, std::shared_ptr<MuonValR4::IRootVisualizationService::ICanvasObject>> canvases{};
    const MdtIdHelper& idHelper{m_idHelperSvc->mdtIdHelper()};
    for (const Identifier& test_me : m_testStations) {
        Identifier chamberId = idHelper.elementID(test_me);
        ATH_MSG_DEBUG("Preparing canvases to draw the tube mapping per chamber");
        const MdtReadoutElement* readEle = m_detMgr->getMdtReadoutElement(test_me);
        if (!readEle) {
            ATH_MSG_DEBUG("Detector element does not exist. ");
            continue;
        }
        if(canvases.count(chamberId)) continue;
        const std::string canName = std::format("{:}{:}{:}{:}",
                                                m_idHelperSvc->stationNameString(readEle->identify()),
                                                std::abs(m_idHelperSvc->stationEta(readEle->identify())),
                                                m_idHelperSvc->stationEta(readEle->identify()) > 0 ? 'A' : 'C',
                                                m_idHelperSvc->stationPhi(readEle->identify()));
        canvases.insert(std::make_pair(chamberId,
                                        m_visualSvc->prepareCanvas(ctx, m_clientTokenStaggeringVis, canName)));
        
        canvases[chamberId] -> add(MuonValR4::drawLabelNDC(m_idHelperSvc->toStringChamber(test_me),0.4,0.85,30));
 
    }
 
    for (const Identifier& test_me : m_testStations) {
        Identifier chamberId = idHelper.elementID(test_me);
        const MdtReadoutElement* readEle = m_detMgr->getMdtReadoutElement(test_me);
        if (!readEle) {
            ATH_MSG_DEBUG("Detector element does not exist. ");
            continue;
        }
        auto canvas = canvases[chamberId];
        const Chamber* chamber = m_detMgr->getChamber(readEle->identify());
        const double r = readEle->tubeRadius();
        for(uint tube = 1 ; tube <= readEle->numTubesInLay(); ++tube){
            for(uint layer = 1 ; layer <= readEle->numLayers() ; ++ layer){
                bool is_valid{false};
                const Identifier tube_id = idHelper.channelID(test_me, readEle->multilayer(), 
                                                              layer, tube, is_valid);
                const IdentifierHash tube_hash = readEle->measurementHash(tube_id);
                if (!readEle->isValid(tube_hash)) {
                    continue;
                }
                
                MdtCablingData cabling_data{};
                cabling->convert(tube_id,cabling_data);
                if (!cabling->getOnlineId(cabling_data, msgStream())) {
                    // to be uncommented once BIS1 has the correct number of tubes
                    //ATH_MSG_WARNING("Could no retrieve a valid online channel for "<<m_idHelperSvc->toString(tube_id));
                    continue;
                    return StatusCode::FAILURE;
                }
 
 
                const Amg::Vector3D tubePos = chamber->globalToLocalTrans(gctx) * readEle->readOutPos(gctx, tube_id);
                const Amg::Vector3D locDir{Amg::Vector3D::UnitY()};
 
                if(tube == 1){
                    Identifier tube2 = idHelper.channelID(test_me, readEle->multilayer(), layer, 2);
                    ATH_MSG_VERBOSE("processing tube" << m_idHelperSvc->toString(tube_id)
                                    << " " << Amg::toString(tubePos) << " " << Amg::toString(chamber->globalToLocalTrans(gctx) * readEle->readOutPos(gctx, tube2)) << " " << 
                                    Amg::toString(chamber->localToGlobalTrans(gctx).linear() * locDir));
                }
                if(tube==1 && layer ==1 && idHelper.stationPhi(tube_id)==1){
                    ATH_MSG_VERBOSE(Amg::toString(chamber->localToGlobalTrans(gctx)));
                }
 
                // In the chamber coordinate system z points to the next tube layer, x points along the tube and y points to the next tube in the same layer 
                // If one looks from the readout side and the readout side is at positive local x one looks in negative x direction therefore one needs to invert the sign of y coordinate, e.g. the tube staggering flips.
                const double x{tubePos.y() * -1. * tubePos.x()/std::abs(tubePos.x()) * std::copysign(1.0, idHelper.stationEta(tube_id))   };
                const double y{tubePos.z()};
                canvas->expandPad(x -r , y -r);
                canvas->expandPad(x +r , y +r);
                if (tube == 1 && layer == 1) {
                    canvas->add(MuonValR4::drawDriftCircle(Amg::Vector3D{0., x, y}, r, kRed, MuonValR4::fullFilling));
                } else {
                    canvas->add(MuonValR4::drawDriftCircle(Amg::Vector3D{0., x, y}, r, kBlack, MuonValR4::hollowFilling));
                }
                canvas->add(MuonValR4::drawLabel(std::to_string(cabling_data.channelId),x ,y, 3));
            }
        }
 
    }
    return StatusCode::SUCCESS;
}
     
StatusCode GeoModelMdtTest::finalize() {
    ATH_CHECK(m_tree.write());
    return StatusCode::SUCCESS;
}
StatusCode GeoModelMdtTest::execute() {
    const EventContext& ctx{Gaudi::Hive::currentContext()};
    
    const ActsTrk::GeometryContext* geoContextHandle{nullptr};
    ATH_CHECK(SG::get(geoContextHandle, m_geoCtxKey, ctx));
 
    const ActsTrk::GeometryContext& gctx{*geoContextHandle};
 
    if(m_visualStaggering)ATH_CHECK(visualizeTubeStaggering(ctx, gctx));
 
    const MdtIdHelper& id_helper{m_idHelperSvc->mdtIdHelper()};
    for (const Identifier& test_me : m_testStations) {
      const int ml = id_helper.multilayer(test_me);
      const std::string detStr = m_idHelperSvc->toStringDetEl(test_me);
      ATH_MSG_DEBUG("Test retrieval of Mdt detector element "<<detStr);
      const MdtReadoutElement* reElement = m_detMgr->getMdtReadoutElement(test_me);
      if (!reElement) {
         continue;
      }
      if (reElement->identify() != test_me) {
         ATH_MSG_FATAL("Expected to retrieve "<<detStr<<". But got instead "<<m_idHelperSvc->toStringDetEl(reElement->identify()));
         return StatusCode::FAILURE;
      }
      ATH_CHECK(dumpToTree(ctx,gctx,reElement));
      const Amg::Transform3D globToLocal{reElement->globalToLocalTrans(gctx)};
      const Amg::Transform3D& localToGlob{reElement->localToGlobalTrans(gctx)};
      const Amg::Transform3D transClosure = globToLocal * localToGlob;
      if (!Amg::doesNotDeform(transClosure)) {
            ATH_MSG_FATAL("Closure test failed for "<<detStr<<". Ended up with "<< Amg::toString(transClosure) );
            return StatusCode::FAILURE;         
      }
      for (unsigned int lay = 1 ; lay <= reElement->numLayers() ; ++lay ) {
        visualizeTubeLayer(ctx, *reElement, lay);
         for (unsigned int tube = 1; tube <=reElement->numTubesInLay(); ++tube ){
            const Identifier tube_id = id_helper.channelID(test_me,ml,lay,tube);                 
            const IdentifierHash measHash = reElement->measurementHash(tube_id);
            const Identifier cnv_tube_id  = reElement->measurementId(measHash);
            if (tube_id != cnv_tube_id) {
               ATH_MSG_FATAL("Failed to convert "<<m_idHelperSvc->toString(tube_id)<<" back and forth "<<m_idHelperSvc->toString(cnv_tube_id));
               return StatusCode::FAILURE;
            }
         }
      }      
   }
   dumpReadoutSideXML();
   return StatusCode::SUCCESS;
}
void GeoModelMdtTest::dumpReadoutSideXML() const {
   if (m_swapRead.empty()) return;
   std::ofstream swapReadXML{m_swapRead};
   if (!swapReadXML.good()) {
      ATH_MSG_ERROR("Failed to create "<<m_swapRead);
      return;
   }
   std::set<Identifier> chamberIDs{};
   const MdtIdHelper& idHelper{m_idHelperSvc->mdtIdHelper()};
   swapReadXML<<"<Table name=\"MdtTubeROSides\">"<<std::endl;
   unsigned int counter{1};
   for (MdtIdHelper::const_id_iterator itr = idHelper.detectorElement_begin();
                                       itr != idHelper.detectorElement_end(); 
                                       ++itr){
      const Identifier swap{*itr};
      const MdtReadoutElement* readoutEle = m_detMgr->getMdtReadoutElement(swap);
      if(!readoutEle) continue;
      if (!chamberIDs.insert(idHelper.elementID(swap)).second) continue;
      const int side = readoutEle->getParameters().readoutSide;
      swapReadXML<<"    <Row   ";
      swapReadXML<<"MDTTUBEROSIDES_DATA_ID=\""<<counter<<"\" ";
      swapReadXML<<"stationName=\""<<m_idHelperSvc->stationNameString(swap)<<"\" ";
      swapReadXML<<"stationEta=\""<<m_idHelperSvc->stationEta(swap)<<"\" ";
      swapReadXML<<"stationPhi=\""<<m_idHelperSvc->stationPhi(swap)<<"\" ";
      swapReadXML<<"side=\""<<side<<"\" ";
      swapReadXML<<"/>"<<std::endl;
      ++counter;
   }
   swapReadXML<<"</Table>"<<std::endl;
 
}
StatusCode GeoModelMdtTest::dumpToTree(const EventContext& ctx,
                                       const ActsTrk::GeometryContext& gctx, 
                                       const MdtReadoutElement* readoutEle) {
 
                                  
   m_stIndex = readoutEle->stationName();
   m_stEta = readoutEle->stationEta();
   m_stPhi = readoutEle->stationPhi();
   m_stML = readoutEle->multilayer();
   m_chamberDesign = readoutEle->chamberDesign();
   
   m_numLayers = readoutEle->numLayers();
   m_numTubes = readoutEle->numTubesInLay();
   
   m_tubeRad = readoutEle->innerTubeRadius();
   m_tubePitch = readoutEle->tubePitch();

   const Amg::Transform3D& transform {readoutEle->localToGlobalTrans(gctx)};
   m_readoutTransform = transform;
   m_alignableNode  = readoutEle->alignableTransform()->getDefTransform();
 
   const MuonMDT_CablingMap* cabling{nullptr};
   ATH_CHECK(SG::get(cabling, m_cablingKey, ctx)); 
   for (unsigned int lay = 1; lay <= readoutEle->numLayers(); ++lay) {
        for (unsigned int tube = 1; tube <= readoutEle->numTubesInLay(); ++tube) {
            const IdentifierHash measHash{readoutEle->measurementHash(lay,tube)};
            if (!readoutEle->isValid(measHash)) continue;
            const Amg::Transform3D& tubeTransform{readoutEle->localToGlobalTrans(gctx,measHash)};
            m_tubeLay.push_back(lay);
            m_tubeNum.push_back(tube);         
            m_tubeTransform.push_back(tubeTransform);
            m_tubePosInCh.push_back(readoutEle->msSector()->globalToLocalTrans(gctx) * 
                                    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));
            m_wireLength.push_back(readoutEle->wireLength(measHash));
            if (cabling) {
                MdtCablingData translation{};
                if (!cabling->convert(readoutEle->measurementId(measHash), translation) ||
                    !cabling->getOnlineId(translation, msgStream())){
                    ATH_MSG_FATAL("Cabling translation failed");
                    return StatusCode::FAILURE;
                }
                m_cablingCSM.push_back(translation.csm);
                m_cablingMROD.push_back(translation.mrod);
                m_cablingTdcId.push_back(translation.tdcId);
                m_cablingTdcCh.push_back(translation.channelId);
            }
        }
   }
 
   return m_tree.fill(ctx) ? StatusCode::SUCCESS : StatusCode::FAILURE;
}
 
}