MuonPhaseII/MuonDetDescr/MuonGeoModelTestR4/src/GeoModelTgcTest.cxx

Go to the documentation of this file.

 
/*
  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
*/
#include "GeoModelTgcTest.h"
#include <ActsGeometryInterfaces/ActsGeometryContext.h>
#include <MuonReadoutGeometryR4/TgcReadoutElement.h>
#include <EventPrimitives/EventPrimitivesToStringConverter.h>
#include <fstream>
 
using namespace ActsTrk;
 
namespace MuonGMR4{
 
GeoModelTgcTest::GeoModelTgcTest(const std::string& name, ISvcLocator* pSvcLocator):
    AthHistogramAlgorithm(name,pSvcLocator) {}
 
StatusCode GeoModelTgcTest::initialize() {
    ATH_CHECK(m_idHelperSvc.retrieve());
    ATH_CHECK(m_geoCtxKey.initialize());
    ATH_CHECK(m_tree.init(this));
 
    const TgcIdHelper& id_helper{m_idHelperSvc->tgcIdHelper()};
    for (const std::string& testCham : m_selectStat) {
        if (std::find(m_excludeStat.begin(), m_excludeStat.end(), testCham) != m_excludeStat.end()) {
            continue;
        }
        if (testCham.size() != 7) {
            ATH_MSG_FATAL("Wrong format given " << testCham);
            return StatusCode::FAILURE;
        }
        std::string statName = testCham.substr(0, 3);
        unsigned int statEta = std::atoi(testCham.substr(3, 1).c_str()) *
                               (testCham[4] == 'A' ? 1 : -1);
        unsigned int statPhi = std::atoi(testCham.substr(5, 2).c_str());
        bool is_valid{false};
        const Identifier eleId{id_helper.elementID(statName, statEta, statPhi, is_valid)};
        if (!is_valid) {
            ATH_MSG_FATAL("Failed to deduce a station name for " << testCham);
            return StatusCode::FAILURE;
        }
        m_testStations.insert(eleId);
    }
    if (m_testStations.empty()){
        std::set<Identifier> excludedStations{};
        for (const std::string& testCham : m_excludeStat) {
            if (testCham.size() != 7) {
                ATH_MSG_FATAL("Wrong format given " << testCham);
                return StatusCode::FAILURE;
            }
            std::string statName = testCham.substr(0, 3);
            unsigned int statEta = std::atoi(testCham.substr(3, 1).c_str()) *
                                   (testCham[4] == 'A' ? 1 : -1);
            unsigned int statPhi = std::atoi(testCham.substr(5, 2).c_str());
            bool is_valid{false};
            const Identifier eleId{id_helper.elementID(statName, statEta, statPhi, is_valid)};
            if (!is_valid) {
                ATH_MSG_FATAL("Failed to deduce a station name for " << testCham);
                return StatusCode::FAILURE;
            }
            excludedStations.insert(eleId);
        }
        std::copy_if(id_helper.detectorElement_begin(), 
                     id_helper.detectorElement_end(), 
                     std::inserter(m_testStations, m_testStations.end()),
                     [&](const Identifier& id) {
                        return excludedStations.count(id) == 0;
                     });
        if (!excludedStations.empty()) {
            std::stringstream excluded_report{};
            for (const Identifier& id : excludedStations){
                excluded_report << " *** " << m_idHelperSvc->toString(id) << std::endl;
            }
            ATH_MSG_INFO("Test all station except the following excluded ones " << std::endl << excluded_report.str());
        }
    } else {
        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());
    }
    ATH_CHECK(detStore()->retrieve(m_detMgr));
    return StatusCode::SUCCESS;
}
StatusCode GeoModelTgcTest::finalize() {
    ATH_CHECK(m_tree.write());
    return StatusCode::SUCCESS;
}
StatusCode GeoModelTgcTest::execute() {
    const EventContext& ctx{Gaudi::Hive::currentContext()};
 
    SG::ReadHandle<ActsGeometryContext> geoContextHandle{m_geoCtxKey, ctx};
    ATH_CHECK(geoContextHandle.isPresent());
    const ActsGeometryContext& gctx{*geoContextHandle};
 
    for (const Identifier& test_me : m_testStations) {
      ATH_MSG_DEBUG("Test retrieval of Tgc detector element "<<m_idHelperSvc->toStringDetEl(test_me));
      const TgcReadoutElement* reElement = m_detMgr->getTgcReadoutElement(test_me);
      if (!reElement) {
         continue;
      }
      if (reElement->identify() != test_me) {
         ATH_MSG_FATAL("Expected to retrieve "<<m_idHelperSvc->toStringDetEl(test_me)
                      <<". But got instead "<<m_idHelperSvc->toStringDetEl(reElement->identify()));
         return StatusCode::FAILURE;
      }
      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 "<<m_idHelperSvc->toStringDetEl(test_me)
                        <<". Ended up with "<< Amg::toString(transClosure) );
            return StatusCode::FAILURE;                  
      }
      const TgcIdHelper& id_helper{m_idHelperSvc->tgcIdHelper()};
      for (unsigned int gasGap = 1; gasGap <= reElement->nGasGaps(); ++gasGap) {
        for (bool isStrip : {false, true}) {
            const unsigned int nChan = isStrip ? reElement->numStrips(gasGap)
                                               : reElement->numWireGangs(gasGap);
            for (unsigned int chan = 1; chan <= nChan ; ++chan) {
                bool isValid{false};
                const Identifier channelId = id_helper.channelID(reElement->identify(),
                                                                 gasGap, isStrip, chan, isValid);
                if (!isValid) {
                    ATH_MSG_DEBUG("No valid Identifier constructed from the fields "
                                <<m_idHelperSvc->toStringDetEl(reElement->identify())
                                <<"isStrip: "<<(isStrip ? "yay" : "nay")<<" gasGap: "<<gasGap<<
                                " channel: "<<chan);
                    continue;
                }
                const IdentifierHash measHash{reElement->measurementHash(channelId)};
                const Identifier backCnv = reElement->measurementId(measHash);
                if (backCnv != channelId) {
                    ATH_MSG_FATAL("Forward-backward conversion of the Identifier "<<m_idHelperSvc->toString(channelId)
                                <<"failed. Got instead "<<m_idHelperSvc->toString(backCnv));
                    return StatusCode::FAILURE;
                }
                if (reElement->layerHash(channelId) != reElement->layerHash(measHash)) {
                    ATH_MSG_FATAL("The cosntruction of the layer hash from the Identifier "<<m_idHelperSvc->toString(channelId)
                    <<" gave something else than doing it from the measurement hash "<<measHash<<". "<<
                    reElement->layerHash(channelId)<<" vs. "<<reElement->layerHash(measHash));
                }
            }
        }        
    } 
    ATH_CHECK(dumpToTree(ctx, gctx, reElement));  
   }
   return StatusCode::SUCCESS;
}
StatusCode GeoModelTgcTest::dumpToTree(const EventContext& ctx,
                                       const ActsGeometryContext& gctx, 
                                       const TgcReadoutElement* reElement) {
   
   m_stIndex    = reElement->stationName();
   m_stEta      = reElement->stationEta();
   m_stPhi      = reElement->stationPhi();
   m_stLayout   = reElement->chamberDesign();
   m_nGasGaps   = reElement->nGasGaps();
   m_readoutTransform = reElement->localToGlobalTrans(gctx);
 
  m_alignableNode  = reElement->alignableTransform()->getDefTransform();
 
   m_shortWidth = reElement->moduleWidthS();
   m_longWidth = reElement->moduleWidthL();
   m_height = reElement->moduleHeight();
   m_thickness = reElement->moduleThickness();
 
   const TgcIdHelper& idHelper{m_idHelperSvc->tgcIdHelper()};   
   for (unsigned int gap = 1; gap <= reElement->nGasGaps(); ++gap) {
        for (unsigned int strip = 1 ; strip <= reElement->numStrips(gap); ++strip) {
            const Identifier measId = idHelper.channelID(reElement->identify(),gap, true, strip);
            const RadialStripDesign& layout{reElement->stripLayout(gap)};
 
            const Amg::Transform3D& localToGlobal{reElement->localToGlobalTrans(gctx,measId)};
            if (strip == 1) {
                m_layTans.push_back(localToGlobal);
                m_layMeasPhi.push_back(true);
                m_layNumber.push_back(gap);
                m_layShortWidth.push_back(2.*layout.shortHalfHeight());
                m_layLongWidth.push_back(2.*layout.longHalfHeight());
                m_layHeight.push_back(2.*layout.halfWidth());
                m_layNumWires.push_back(0);
            }
            m_stripGasGap.push_back(gap);
            m_stripNum.push_back(strip);
            m_stripCenter.push_back(reElement->channelPosition(gctx, measId));
            const RadialStripDesign& stripDesign{reElement->stripLayout(gap)};
            const Amg::Vector2D locTop2D{stripDesign.leftEdge(strip).value_or(Amg::Vector2D::Zero())};
            const Amg::Vector2D locBot2D{stripDesign.rightEdge(strip).value_or(Amg::Vector2D::Zero())};
            const Amg::Vector3D globTop{localToGlobal * Amg::Vector3D{locTop2D.x(), locTop2D.y(), 0}};
            const Amg::Vector3D globBot{localToGlobal * Amg::Vector3D{locBot2D.x(), locBot2D.y(), 0}};
            
            m_stripBottom.push_back(globBot);
            m_stripTop.push_back(globTop);
            m_locStripTop.push_back(locTop2D);
            m_locStripCenter.push_back(stripDesign.center(strip).value_or(Amg::Vector2D::Zero()));
            m_locStripBottom.push_back(locBot2D);
 
        }
        for (unsigned int gang = 1; gang <= reElement->numWireGangs(gap); ++gang) {
            const Identifier measId = idHelper.channelID(reElement->identify(),gap, false, gang);
            const WireGroupDesign& layout{reElement->wireGangLayout(gap)};
            if (gang == 1) {
                m_layTans.push_back(reElement->localToGlobalTrans(gctx,measId));
                m_layMeasPhi.push_back(false);
                m_layNumber.push_back(gap);
                m_layShortWidth.push_back(2.*layout.shortHalfHeight());
                m_layLongWidth.push_back(2.*layout.longHalfHeight());
                m_layHeight.push_back(2.*layout.halfWidth());
                m_layNumWires.push_back(layout.nAllWires());
            }
            m_gangNum.push_back(gang);
            m_gangGasGap.push_back(gap);
            m_gangCenter.push_back(reElement->channelPosition(gctx, measId));
            m_gangNumWires.push_back(layout.numWiresInGroup(gang));
            m_locGangPos.push_back(layout.center(gang).value_or(Amg::Vector2D::Zero()));
            m_gangLength.push_back(layout.stripLength(gang));
        }
   }
   return m_tree.fill(ctx) ? StatusCode::SUCCESS : StatusCode::FAILURE;
}
 
}