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