TgcReadoutGeomTool.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
#include "TgcReadoutGeomTool.h"
 
#include <GaudiKernel/SystemOfUnits.h>
#include <RDBAccessSvc/IRDBAccessSvc.h>
#include <RDBAccessSvc/IRDBRecordset.h>
 
#include <GeoModelKernel/GeoFullPhysVol.h>
#include <GeoModelKernel/GeoPhysVol.h>
#include <GeoModelKernel/GeoTrd.h>
#include <GeoModelKernel/GeoBox.h>
 
#include <GeoModelRead/ReadGeoModel.h>
#include <MuonReadoutGeometryR4/MuonDetectorManager.h>
#include <MuonReadoutGeometryR4/WireGroupDesign.h>
#include <MuonReadoutGeometryR4/RadialStripDesign.h>
 
#include <format>
#include <RDBAccessSvc/IRDBRecord.h>
 
#include <MuonDetDescrUtils/MuonSectorMapping.h>
 
#ifndef SIMULATIONBASE
#   include "Acts/Utilities/BoundFactory.hpp"
#endif
 
using namespace CxxUtils;
using namespace ActsTrk;
 
namespace MuonGMR4 {
 
using physVolWithTrans = IMuonGeoUtilityTool::physVolWithTrans;
using defineArgs = TgcReadoutElement::defineArgs;
 
 
 
std::unique_ptr<WireGroupDesign> 
    TgcReadoutGeomTool::constructWireDesign(const wTgcTable& table,
                                            const GeoTrd* gapTrd) const {
    if (table.wireGangs.empty()) {
        return nullptr;
    }
    const double halfMinX = std::min(gapTrd->getYHalfLength1(), gapTrd->getYHalfLength2());
    const double halfMaxX = std::max(gapTrd->getYHalfLength1(), gapTrd->getYHalfLength2());
    const double halfY = gapTrd->getZHalfLength();
    auto design = std::make_unique<WireGroupDesign>();
    design->defineTrapezoid(halfMinX, halfMaxX, halfY);
    for (unsigned gang : table.wireGangs) {
        design->declareGroup(gang); 
    }
    const double wireOffSet = -0.5*table.wirePitch * design->nAllWires();
    design->defineStripLayout(Amg::Vector2D{wireOffSet, 0.},
                              table.wirePitch, 0., table.wireGangs.size());
    return design;
}
std::unique_ptr<RadialStripDesign> 
    TgcReadoutGeomTool::constructRadialDesign(const wTgcTable& table,
                                              const GeoTrd* gapTrd) const {
    if (table.bottomStripPos.empty()) {
        return nullptr;
    }
    const double halfMinX = std::min(gapTrd->getYHalfLength1(), gapTrd->getYHalfLength2());
    const double halfMaxX = std::max(gapTrd->getYHalfLength1(), gapTrd->getYHalfLength2());
    const double halfY = gapTrd->getZHalfLength();
 
    auto design = std::make_unique<RadialStripDesign>();
 
    design->defineTrapezoid(halfMinX, halfMaxX, halfY);
    design->defineStripLayout(Amg::Vector2D{-halfY,0.},
                              0.,0.,table.bottomStripPos.size());
    design->flipTrapezoid();
    for (size_t s = 0; s < table.bottomStripPos.size(); ++s) {
        design->addStrip(table.bottomStripPos.at(s),
                         table.topStripPos.at(s));
    }
    return design;
}
 
StatusCode TgcReadoutGeomTool::loadDimensions(TgcReadoutElement::defineArgs& define,
                                              FactoryCache& factoryCache) {
    ATH_MSG_VERBOSE("Load dimensions of "<<m_idHelperSvc->toString(define.detElId)
                     <<std::endl<<std::endl<<m_geoUtilTool->dumpVolume(define.physVol));
    const GeoShape* shape = m_geoUtilTool->extractShape(define.physVol);
    if (!shape) {
        ATH_MSG_FATAL("Failed to deduce a valid shape for "<<m_idHelperSvc->toString(define.detElId));
        return StatusCode::FAILURE;
    }
 
    if (shape->typeID() != GeoTrd::getClassTypeID()) {
        ATH_MSG_FATAL("The shape of "<<m_idHelperSvc->toStringDetEl(define.detElId)
                    <<" is expected to be a trapezoid "<<m_geoUtilTool->dumpShape(shape));
        return StatusCode::FAILURE;
    }
    const GeoTrd* chambTrd = static_cast<const GeoTrd*>(shape);
    define.halfWidthShort = std::min(chambTrd->getYHalfLength1(), chambTrd->getYHalfLength2());
    define.halfWidthLong = std::max(chambTrd->getYHalfLength1(), chambTrd->getYHalfLength2());
    define.halfHeight = chambTrd->getZHalfLength();
    define.halfThickness = chambTrd->getXHalfLength1();
    std::vector<physVolWithTrans> allGasGaps = m_geoUtilTool->findAllLeafNodesByName(define.physVol, "TgcGas");
    if (allGasGaps.empty()) {
        ATH_MSG_FATAL("The volume "<<m_idHelperSvc->toStringDetEl(define.detElId)<<" does not have any childern TgcGas"
            <<std::endl<<m_geoUtilTool->dumpVolume(define.physVol));
        return StatusCode::FAILURE;
    }
    unsigned gasGap{0};
    auto assignReadoutLayer = [&define, & gasGap, this] (const StripLayerPtr& layerReadout) {
        if (!layerReadout) {
            return StatusCode::SUCCESS;
        }
        unsigned layerIdx = static_cast<unsigned>(layerReadout->hash());
        if (layerIdx >= define.sensorLayouts.size()) {
            ATH_MSG_FATAL("The strip index "<<layerIdx<<" is out of range for gasGap "<<gasGap);
            return StatusCode::FAILURE;
        }           
        define.sensorLayouts[layerIdx] = layerReadout;
        return StatusCode::SUCCESS;
    };
    for (const physVolWithTrans& pVolTrans : allGasGaps) {
        const std::string key = std::format("{:}_{:}{:}",define.chambDesign, gasGap+1,
                                            (m_idHelperSvc->stationEta(define.detElId) > 0 ? "A" : "C"));
        StripLayerPtr& wireReadout{factoryCache.wireLayers[key]};
        StripLayerPtr& stripReadout{factoryCache.stripLayers[key]};
        if (!wireReadout || !stripReadout) {
            const wTgcTable& table{factoryCache.parameterBook[key]};
            if (!table.gasGap) {
                ATH_MSG_FATAL("No wTGC table could be found for "<<m_idHelperSvc->toStringDetEl(define.detElId)
                            <<" "<<define.chambDesign<<", gasGap "<<(gasGap+1));
                return StatusCode::FAILURE;
            }
            const GeoShape* gapShape = m_geoUtilTool->extractShape(pVolTrans.volume);
            if (gapShape->typeID() != GeoTrd::getClassTypeID()) {
                ATH_MSG_FATAL("Expected shape "<<m_geoUtilTool->dumpShape(gapShape)<<" to be a trapezoid");
                return StatusCode::FAILURE;
            }
            const GeoTrd* gapTrd = static_cast<const GeoTrd*>(gapShape);
            WireDesignPtr wireDesign = constructWireDesign(table, gapTrd);
            RadialStripDesignPtr radDesign = constructRadialDesign(table, gapTrd);
            if (wireDesign) {
                wireDesign = (*factoryCache.wireLayouts.insert(std::move(wireDesign)).first);
            }
            if (radDesign) {
                radDesign = (*factoryCache.stripReadouts.insert(std::move(radDesign)).first);
            }
            if (!wireReadout && wireDesign) {
                const IdentifierHash layHash = TgcReadoutElement::constructHash(0, gasGap+1, false);
                ATH_MSG_VERBOSE("Wire hash "<<layHash);
                const Amg::Transform3D trans{pVolTrans.transform 
                                             * Amg::getRotateY3D(-90.*Gaudi::Units::deg)
                                             * Amg::getRotateX3D(180.* Gaudi::Units::deg)};
                if (false && radDesign) {
                    wireReadout = std::make_unique<StripLayer>(factoryCache.trfNodeMaker.makeTransform(trans), 
                                                                wireDesign, radDesign, layHash);
                    wireReadout->flipPhiRotation();
                } else {
                    wireReadout = std::make_unique<StripLayer>(factoryCache.trfNodeMaker.makeTransform(trans), 
                                                                wireDesign, layHash);
                }
            } 
            if (!stripReadout && radDesign) {
                const IdentifierHash layHash = TgcReadoutElement::constructHash(0, gasGap+1, true);
                ATH_MSG_VERBOSE("Radial hash "<<layHash);
                const Amg::Transform3D trans{pVolTrans.transform 
                                             * Amg::getRotateZ3D(90.* Gaudi::Units::deg)
                                             * Amg::getRotateX3D(90.*Gaudi::Units::deg)};
                stripReadout = std::make_unique<StripLayer>(factoryCache.trfNodeMaker.makeTransform(trans), 
                                                            radDesign, layHash);
            }
        }
        ATH_CHECK(assignReadoutLayer(stripReadout));
        ATH_CHECK(assignReadoutLayer(wireReadout));
        ++gasGap;
    }
    define.nGasGaps = gasGap;        
    return StatusCode::SUCCESS;
}
StatusCode TgcReadoutGeomTool::buildReadOutElements(MuonDetectorManager& mgr) {
    ATH_CHECK(m_geoDbTagSvc.retrieve());
    ATH_CHECK(m_idHelperSvc.retrieve());
    ATH_CHECK(m_geoUtilTool.retrieve());
    GeoModelIO::ReadGeoModel* sqliteReader = m_geoDbTagSvc->getSqliteReader();
    if (!sqliteReader) {
        ATH_MSG_FATAL("Error, the tool works exclusively from sqlite geometry inputs");
        return StatusCode::FAILURE;
    }
    
    FactoryCache facCache{};
    ATH_CHECK(readParameterBook(facCache));
 
#ifndef SIMULATIONBASE
    auto layerBounds = std::make_shared<Acts::SurfaceBoundFactory>();
#endif    
    const TgcIdHelper& idHelper{m_idHelperSvc->tgcIdHelper()};
    // Get the list of full phys volumes from SQLite, and create detector elements
    physNodeMap mapFPV = sqliteReader->getPublishedNodes<std::string, GeoFullPhysVol*>("Muon");
    for (auto& [key, pv] : mapFPV) {
        std::vector<std::string> key_tokens = tokenize(key, "_");
        if (key_tokens.size() < 4 ||
            !key_tokens[0].starts_with( "TGC"))
            continue;
        
        bool isValid{false};
        const Identifier elementID = idHelper.elementID(key_tokens[1].substr(0,3), atoi(key_tokens[2]), atoi(key_tokens[3]), isValid);
        if (!isValid){
            ATH_MSG_FATAL("Failed to construct the station Identifier from "<<key);
            return StatusCode::FAILURE;
        }        
        defineArgs define{};
        define.physVol = pv;
        define.detElId = elementID;
        define.chambDesign = key_tokens[0];
        define.alignTransform = m_geoUtilTool->findAlignableTransform(define.physVol);
#ifndef SIMULATIONBASE
        define.layerBounds = layerBounds;
#endif
        ATH_MSG_DEBUG("Key  "<<key<<" lead to Identifier "<<m_idHelperSvc->toStringDetEl(define.detElId));
        ATH_CHECK(loadDimensions(define, facCache));
        std::unique_ptr<TgcReadoutElement> readoutEle = std::make_unique<TgcReadoutElement>(std::move(define));
        ATH_CHECK(mgr.addTgcReadoutElement(std::move(readoutEle)));
    }
    ATH_CHECK(writeSectorMapping(mgr));
    return StatusCode::SUCCESS;
}
StatusCode TgcReadoutGeomTool::writeSectorMapping(const MuonDetectorManager& mgr) const {
    std::vector<const TgcReadoutElement*> tgcReadOutEles = mgr.getAllTgcReadoutElements();
    std::unique_ptr<std::vector<int>> tgcSectorMapping = std::make_unique<std::vector<int>>();
    tgcSectorMapping->resize(m_idHelperSvc->tgcIdHelper().module_hash_max());
    Muon::MuonSectorMapping sectorMapping{};
    const ActsGeometryContext gctx{};
    for (const TgcReadoutElement* readoutEle : tgcReadOutEles) {
        int& sectNumb = (*tgcSectorMapping)[m_idHelperSvc->moduleHash(readoutEle->identify())];
        sectNumb = sectorMapping.getSector(readoutEle->center(gctx).phi());
    }
    ATH_CHECK(detStore()->record(std::move(tgcSectorMapping), "TGC_SectorMapping"));
    return StatusCode::SUCCESS;
}
StatusCode TgcReadoutGeomTool::readParameterBook(FactoryCache& cache) {
    ServiceHandle<IRDBAccessSvc> accessSvc(m_geoDbTagSvc->getParamSvcName(), name());
    ATH_CHECK(accessSvc.retrieve());
    IRDBRecordset_ptr paramTable = accessSvc->getRecordsetPtr("TgcSensorLayout", "");
    if (paramTable->size() == 0) {
        ATH_MSG_FATAL("Empty parameter book table found");
        return StatusCode::FAILURE;
    }
    for (const IRDBRecord_ptr& record : *paramTable) {
        const std::string chambType = record->getString("technology");
        const int gasGap = record->getInt("gasGap");
        const std::vector<int> wireGangs{tokenizeInt(record->getString("wireGangs"),",")};
        const std::vector<double> botStrips = tokenizeDouble(record->getString("bottomStrips"), ",");
        const std::vector<double> topStrips = tokenizeDouble(record->getString("topStrips"), ",");
        const std::vector<std::string> sides = tokenize(record->getString("side"), ";");
        const double wirePitch = record->getDouble("wirePitch");
        for (const std::string& side : sides) {
            const std::string key = std::format("{:}_{:}{:}", chambType, gasGap, side);
            wTgcTable& parBook{cache.parameterBook[key]};
            parBook.wireGangs.insert(parBook.wireGangs.end(), wireGangs.begin(), wireGangs.end());
            parBook.bottomStripPos = botStrips;
            parBook.topStripPos = topStrips;
            parBook.wirePitch = wirePitch;
            parBook.gasGap = gasGap;
        }
       
    }
    ATH_MSG_DEBUG("Read in total "<<cache.parameterBook.size()<<" chamber layouts");
    return StatusCode::SUCCESS;
}
}  // namespace MuonGMR4