/*
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;
}