From db2346d1a92446905ec3a139aa0658777376cb16 Mon Sep 17 00:00:00 2001 From: Fabio Cossutti Date: Tue, 12 Nov 2024 17:57:40 +0100 Subject: [PATCH 1/2] Update MTDGeometry test, use only DD4hep backend, add ETL pixel position test --- .../MTDGeometryBuilder/test/BuildFile.xml | 15 +- .../test/DD4hep_TestBTLPixelTopology.cc | 344 ------------ .../test/DD4hep_TestPixelTopology.cc | 517 ++++++++++++++++++ .../test/MTDDigiGeometryAnalyzer.cc | 158 +----- .../MTDGeometryBuilder/test/dd4hep_mtd_cfg.py | 64 --- Geometry/MTDGeometryBuilder/test/mtd_cfg.py | 34 +- Geometry/MTDGeometryBuilder/test/runTest.sh | 12 +- 7 files changed, 563 insertions(+), 581 deletions(-) delete mode 100644 Geometry/MTDGeometryBuilder/test/DD4hep_TestBTLPixelTopology.cc create mode 100644 Geometry/MTDGeometryBuilder/test/DD4hep_TestPixelTopology.cc delete mode 100644 Geometry/MTDGeometryBuilder/test/dd4hep_mtd_cfg.py diff --git a/Geometry/MTDGeometryBuilder/test/BuildFile.xml b/Geometry/MTDGeometryBuilder/test/BuildFile.xml index a4e689eabb3c7..e9d7d8d29d77f 100644 --- a/Geometry/MTDGeometryBuilder/test/BuildFile.xml +++ b/Geometry/MTDGeometryBuilder/test/BuildFile.xml @@ -1,17 +1,4 @@ - - - - - - - - - - - - - - + diff --git a/Geometry/MTDGeometryBuilder/test/DD4hep_TestBTLPixelTopology.cc b/Geometry/MTDGeometryBuilder/test/DD4hep_TestBTLPixelTopology.cc deleted file mode 100644 index 5ec8d2157cb76..0000000000000 --- a/Geometry/MTDGeometryBuilder/test/DD4hep_TestBTLPixelTopology.cc +++ /dev/null @@ -1,344 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/one/EDAnalyzer.h" -#include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/EventSetup.h" -#include "FWCore/Framework/interface/ESTransientHandle.h" -#include "FWCore/Framework/interface/MakerMacros.h" -#include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "FWCore/Utilities/interface/Exception.h" - -#include "Geometry/Records/interface/IdealGeometryRecord.h" -#include "Geometry/Records/interface/DDSpecParRegistryRcd.h" - -#include "DetectorDescription/DDCMS/interface/DDDetector.h" -#include "DetectorDescription/DDCMS/interface/DDSolidShapes.h" -#include "DetectorDescription/DDCMS/interface/DDFilteredView.h" -#include "DetectorDescription/DDCMS/interface/DDSpecParRegistry.h" - -#include "Geometry/MTDCommonData/interface/MTDTopologyMode.h" -#include "Geometry/MTDCommonData/interface/MTDBaseNumber.h" -#include "Geometry/MTDCommonData/interface/BTLNumberingScheme.h" -#include "Geometry/MTDCommonData/interface/ETLNumberingScheme.h" - -#include "DataFormats/ForwardDetId/interface/BTLDetId.h" -#include "DataFormats/ForwardDetId/interface/ETLDetId.h" - -#include "Geometry/Records/interface/MTDTopologyRcd.h" -#include "Geometry/MTDNumberingBuilder/interface/MTDTopology.h" -#include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h" -#include "Geometry/Records/interface/MTDDigiGeometryRecord.h" -#include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h" -#include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h" - -#include "DataFormats/Math/interface/GeantUnits.h" -#include "DataFormats/Math/interface/Rounding.h" -#include - -using namespace cms; -using namespace geant_units::operators; - -class DD4hep_TestBTLPixelTopology : public edm::one::EDAnalyzer<> { -public: - explicit DD4hep_TestBTLPixelTopology(const edm::ParameterSet&); - ~DD4hep_TestBTLPixelTopology() = default; - - void beginJob() override {} - void analyze(edm::Event const&, edm::EventSetup const&) override; - void endJob() override {} - - void theBaseNumber(cms::DDFilteredView& fv); - -private: - const edm::ESInputTag tag_; - - MTDBaseNumber thisN_; - BTLNumberingScheme btlNS_; - - edm::ESGetToken dddetToken_; - edm::ESGetToken dspecToken_; - edm::ESGetToken mtdtopoToken_; - edm::ESGetToken mtdgeoToken_; - - std::stringstream sunitt; - constexpr static double tolerance{0.5e-3_mm}; -}; - -using DD3Vector = ROOT::Math::DisplacementVector3D>; -using cms_rounding::roundIfNear0; - -DD4hep_TestBTLPixelTopology::DD4hep_TestBTLPixelTopology(const edm::ParameterSet& iConfig) - : tag_(iConfig.getParameter("DDDetector")), thisN_(), btlNS_() { - dddetToken_ = esConsumes(tag_); - dspecToken_ = esConsumes(tag_); - mtdtopoToken_ = esConsumes(tag_); - mtdgeoToken_ = esConsumes(tag_); -} - -void DD4hep_TestBTLPixelTopology::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) { - auto pDD = iSetup.getTransientHandle(dddetToken_); - - auto pSP = iSetup.getTransientHandle(dspecToken_); - - if (!pDD.isValid()) { - edm::LogError("DD4hep_TestBTLPixelTopology") << "ESTransientHandle pDD is not valid!"; - return; - } - if (pDD.description()) { - edm::LogInfo("DD4hep_TestBTLPixelTopology") << pDD.description()->type_ << " label: " << pDD.description()->label_; - } else { - edm::LogWarning("DD4hep_TestBTLPixelTopology") << "NO label found pDD.description() returned false."; - } - - if (!pSP.isValid()) { - edm::LogError("DD4hep_TestBTLPixelTopology") << "ESTransientHandle pSP is not valid!"; - return; - } - - auto pTP = iSetup.getTransientHandle(mtdtopoToken_); - if (!pTP.isValid()) { - edm::LogError("DD4hep_TestBTLPixelTopology") << "ESTransientHandle pTP is not valid!"; - return; - } else { - edm::LogInfo("DD4hep_TestBTLPixelTopology") - << "MTD topology mode = " << pTP.product()->getMTDTopologyMode() << " BTLDetId:CrysLayout = " - << static_cast(MTDTopologyMode::crysLayoutFromTopoMode(pTP.product()->getMTDTopologyMode())); - } - - auto pDG = iSetup.getTransientHandle(mtdgeoToken_); - if (!pDG.isValid()) { - edm::LogError("DD4hep_TestBTLPixelTopology") << "ESTransientHandle pDG is not valid!"; - return; - } else { - edm::LogInfo("DD4hep_TestBTLPixelTopology") - << "Geometry node for MTDGeom is " << &(*pDG) << "\n" - << " # detectors = " << pDG.product()->detUnits().size() << "\n" - << " # types = " << pDG.product()->detTypes().size() << "\n" - << " # BTL dets = " << pDG.product()->detsBTL().size() << "\n" - << " # ETL dets = " << pDG.product()->detsETL().size() << "\n" - << " # layers " << pDG.product()->geomDetSubDetector(1) << " = " << pDG.product()->numberOfLayers(1) << "\n" - << " # layers " << pDG.product()->geomDetSubDetector(2) << " = " << pDG.product()->numberOfLayers(2) << "\n"; - } - - DDFilteredView fv(pDD.product(), pDD.product()->description()->worldVolume()); - fv.next(0); - - DDSpecParRefs specs; - std::string attribute("ReadOutName"), name("FastTimerHitsBarrel"); - pSP.product()->filter(specs, attribute, name); - - edm::LogVerbatim("DD4hep_TestBTLPixelTopology").log([&specs](auto& log) { - log << "Filtered DD SpecPar Registry size: " << specs.size() << "\n"; - for (const auto& t : specs) { - log << "\nSpecPar " << t.first << ":\nRegExps { "; - for (const auto& ki : t.second->paths) - log << ki << " "; - log << "};\n "; - for (const auto& kl : t.second->spars) { - log << kl.first << " = "; - for (const auto& kil : kl.second) { - log << kil << " "; - } - log << "\n "; - } - } - }); - - bool isBarrel = true; - bool exitLoop = false; - uint32_t level(0); - - do { - if (dd4hep::dd::noNamespace(fv.name()) == "BarrelTimingLayer") { - isBarrel = true; - edm::LogInfo("DD4hep_TestBTLPixelTopology") << "isBarrel = " << isBarrel; - } else if (dd4hep::dd::noNamespace(fv.name()) == "EndcapTimingLayer") { - isBarrel = false; - edm::LogInfo("DD4hep_TestBTLPixelTopology") << "isBarrel = " << isBarrel; - } - - theBaseNumber(fv); - - auto print_path = [&]() { - std::stringstream ss; - ss << " - OCMS[0]/"; - for (int ii = thisN_.getLevels() - 1; ii-- > 0;) { - ss << thisN_.getLevelName(ii); - ss << "["; - ss << thisN_.getCopyNumber(ii); - ss << "]/"; - } - return ss.str(); - }; - - if (level > 0 && fv.navPos().size() < level) { - level = 0; - exitLoop = true; - } - if (dd4hep::dd::noNamespace(fv.name()) == "BarrelTimingLayer") { - level = fv.navPos().size(); - } - - // Test only the desired subdetector - - if (exitLoop && isBarrel) { - break; - } - - // Actions for MTD volumes: search for sensitive detectors - - bool isSens = false; - - for (auto const& t : specs) { - for (auto const& it : t.second->paths) { - if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(fv.name()), dd4hep::dd::realTopName(it))) { - isSens = true; - break; - } - } - } - - if (isSens && isBarrel) { - std::stringstream spix; - spix << print_path() << "\n\n"; - - BTLDetId theId(btlNS_.getUnitID(thisN_)); - - DetId geoId = theId.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(pTP.product()->getMTDTopologyMode())); - const MTDGeomDet* thedet = pDG.product()->idToDet(geoId); - - if (thedet == nullptr) { - throw cms::Exception("BTLDeviceSim") << "GeographicalID: " << std::hex << geoId.rawId() << " (" << theId.rawId() - << ") is invalid!" << std::dec << std::endl; - } - const ProxyMTDTopology& topoproxy = static_cast(thedet->topology()); - const RectangularMTDTopology& topo = static_cast(topoproxy.specificTopology()); - - int origRow = theId.row(topo.nrows()); - int origCol = theId.column(topo.nrows()); - spix << "rawId= " << theId.rawId() << " geoId= " << geoId.rawId() << " side/rod= " << theId.mtdSide() << " / " - << theId.mtdRR() << " type/RU= " << theId.modType() << " / " << theId.runit() - << " module/geomodule= " << theId.module() << " / " << static_cast(geoId).module() - << " crys= " << theId.crystal() << " BTLDetId row/col= " << origRow << " / " << origCol; - spix << "\n"; - - // - // Test of positions for sensitive detectors - // - - auto fround = [&](double in) { - std::stringstream ss; - ss << std::fixed << std::setw(14) << roundIfNear0(in); - return ss.str(); - }; - - if (!dd4hep::isA(fv.solid())) { - throw cms::Exception("DD4hep_TestBTLPixelTopology") << "MTD sensitive element not a DDBox"; - break; - } - dd4hep::Box mySens(fv.solid()); - - DD3Vector zeroLocal(0., 0., 0.); - DD3Vector cn1Local(mySens.x(), mySens.y(), mySens.z()); - DD3Vector cn2Local(-mySens.x(), -mySens.y(), -mySens.z()); - DD3Vector zeroGlobal = (fv.rotation())(zeroLocal) + fv.translation(); - DD3Vector cn1Global = (fv.rotation())(cn1Local) + fv.translation(); - DD3Vector cn2Global = (fv.rotation())(cn2Local) + fv.translation(); - - const size_t nTest(3); - std::array refLocalPoints{{Local3DPoint(zeroLocal.x(), zeroLocal.y(), zeroLocal.z()), - Local3DPoint(cn1Local.x(), cn1Local.y(), cn1Local.z()), - Local3DPoint(cn2Local.x(), cn2Local.y(), cn2Local.z())}}; - std::array refGlobalPoints{{zeroGlobal, cn1Global, cn2Global}}; - - for (size_t iloop = 0; iloop < nTest; iloop++) { - // translate from crystal-local coordinates to module-local coordinates to get the row and column - - Local3DPoint cmRefLocal(convertMmToCm(refLocalPoints[iloop].x() / dd4hep::mm), - convertMmToCm(refLocalPoints[iloop].y() / dd4hep::mm), - convertMmToCm(refLocalPoints[iloop].z() / dd4hep::mm)); - Local3DPoint modLocal = topo.pixelToModuleLocalPoint(cmRefLocal, origRow, origCol); - const auto& thepixel = topo.pixelIndex(modLocal); - uint8_t recoRow = static_cast(thepixel.first); - uint8_t recoCol = static_cast(thepixel.second); - - if (origRow != recoRow || origCol != recoCol) { - std::stringstream warnmsg; - warnmsg << "DIFFERENCE row/col, orig= " << origRow << " " << origCol - << " reco= " << static_cast(recoRow) << " " << static_cast(recoCol) << "\n"; - spix << warnmsg.str(); - sunitt << warnmsg.str(); - recoRow = origRow; - recoCol = origCol; - } - - Local3DPoint recoRefLocal = topo.moduleToPixelLocalPoint(modLocal); - - // reconstructed global position from reco geometry and rectangluar MTD topology - - const auto& modGlobal = thedet->toGlobal(modLocal); - - const double deltax = convertCmToMm(modGlobal.x()) - (refGlobalPoints[iloop].x() / dd4hep::mm); - const double deltay = convertCmToMm(modGlobal.y()) - (refGlobalPoints[iloop].y() / dd4hep::mm); - const double deltaz = convertCmToMm(modGlobal.z()) - (refGlobalPoints[iloop].z() / dd4hep::mm); - - const double local_deltax = recoRefLocal.x() - cmRefLocal.x(); - const double local_deltay = recoRefLocal.y() - cmRefLocal.y(); - const double local_deltaz = recoRefLocal.z() - cmRefLocal.z(); - - spix << "Ref#" << iloop << " local= " << fround(refLocalPoints[iloop].x() / dd4hep::mm) - << fround(refLocalPoints[iloop].y() / dd4hep::mm) << fround(refLocalPoints[iloop].z() / dd4hep::mm) - << " Orig global= " << fround(refGlobalPoints[iloop].x() / dd4hep::mm) - << fround(refGlobalPoints[iloop].y() / dd4hep::mm) << fround(refGlobalPoints[iloop].z() / dd4hep::mm) - << " Reco global= " << fround(convertCmToMm(modGlobal.x())) << fround(convertCmToMm(modGlobal.y())) - << fround(convertCmToMm(modGlobal.z())) << " Delta= " << fround(deltax) << fround(deltay) << fround(deltaz) - << " Local Delta= " << fround(local_deltax) << fround(local_deltay) << fround(local_deltaz) << "\n"; - if (std::abs(deltax) > tolerance || std::abs(deltay) > tolerance || std::abs(deltaz) > tolerance) { - std::stringstream warnmsg; - warnmsg << "DIFFERENCE detId/ref# " << theId.rawId() << " " << iloop << " dx/dy/dz= " << fround(deltax) - << fround(deltay) << fround(deltaz) << "\n"; - spix << warnmsg.str(); - sunitt << warnmsg.str(); - } - if (std::abs(local_deltax) > tolerance || std::abs(local_deltay) > tolerance || - std::abs(local_deltaz) > tolerance) { - std::stringstream warnmsg; - warnmsg << "DIFFERENCE detId/ref# " << theId.rawId() << " " << iloop - << " local dx/dy/dz= " << fround(local_deltax) << fround(local_deltay) << fround(local_deltaz) - << "\n"; - spix << warnmsg.str(); - sunitt << warnmsg.str(); - } - } - - spix << "\n"; - edm::LogVerbatim("DD4hep_TestBTLPixelTopology") << spix.str(); - } - } while (fv.next(0)); - - if (!sunitt.str().empty()) { - edm::LogVerbatim("MTDUnitTest") << sunitt.str(); - } -} - -void DD4hep_TestBTLPixelTopology::theBaseNumber(cms::DDFilteredView& fv) { - thisN_.reset(); - thisN_.setSize(fv.navPos().size()); - - for (uint ii = 0; ii < fv.navPos().size(); ii++) { - std::string_view name((fv.geoHistory()[ii])->GetName()); - size_t ipos = name.rfind('_'); - thisN_.addLevel(name.substr(0, ipos), fv.copyNos()[ii]); -#ifdef EDM_ML_DEBUG - edm::LogVerbatim("DD4hep_TestBTLPixelTopology") << ii << " " << name.substr(0, ipos) << " " << fv.copyNos()[ii]; -#endif - } -} - -DEFINE_FWK_MODULE(DD4hep_TestBTLPixelTopology); diff --git a/Geometry/MTDGeometryBuilder/test/DD4hep_TestPixelTopology.cc b/Geometry/MTDGeometryBuilder/test/DD4hep_TestPixelTopology.cc new file mode 100644 index 0000000000000..3cbd0f3adbd36 --- /dev/null +++ b/Geometry/MTDGeometryBuilder/test/DD4hep_TestPixelTopology.cc @@ -0,0 +1,517 @@ +#include +#include +#include +#include +#include +#include + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/one/EDAnalyzer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESTransientHandle.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/Exception.h" + +#include "Geometry/Records/interface/IdealGeometryRecord.h" +#include "Geometry/Records/interface/DDSpecParRegistryRcd.h" + +#include "DetectorDescription/DDCMS/interface/DDDetector.h" +#include "DetectorDescription/DDCMS/interface/DDSolidShapes.h" +#include "DetectorDescription/DDCMS/interface/DDFilteredView.h" +#include "DetectorDescription/DDCMS/interface/DDSpecParRegistry.h" + +#include "Geometry/MTDCommonData/interface/MTDTopologyMode.h" +#include "Geometry/MTDCommonData/interface/MTDBaseNumber.h" +#include "Geometry/MTDCommonData/interface/BTLNumberingScheme.h" +#include "Geometry/MTDCommonData/interface/ETLNumberingScheme.h" + +#include "DataFormats/ForwardDetId/interface/BTLDetId.h" +#include "DataFormats/ForwardDetId/interface/ETLDetId.h" +#include "DataFormats/GeometrySurface/interface/MediumProperties.h" +#include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h" + +#include "Geometry/Records/interface/MTDTopologyRcd.h" +#include "Geometry/MTDNumberingBuilder/interface/MTDTopology.h" +#include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h" +#include "Geometry/MTDGeometryBuilder/interface/MTDGeomDetUnit.h" +#include "Geometry/Records/interface/MTDDigiGeometryRecord.h" +#include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h" +#include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h" + +#include "DataFormats/Math/interface/GeantUnits.h" +#include "DataFormats/Math/interface/Rounding.h" +#include "DataFormats/Math/interface/angle_units.h" +#include + +using namespace cms; +using namespace geant_units::operators; + +class DD4hep_TestPixelTopology : public edm::one::EDAnalyzer<> { +public: + explicit DD4hep_TestPixelTopology(const edm::ParameterSet&); + ~DD4hep_TestPixelTopology() = default; + + void beginJob() override {} + void analyze(edm::Event const&, edm::EventSetup const&) override; + void endJob() override {} + + void theBaseNumber(cms::DDFilteredView& fv); + +private: + void analyseRectangle(const GeomDetUnit& det); + void checkRotation(const GeomDetUnit& det); + + const edm::ESInputTag tag_; + std::string ddTopNodeName_; + + MTDBaseNumber thisN_; + BTLNumberingScheme btlNS_; + ETLNumberingScheme etlNS_; + + edm::ESGetToken dddetToken_; + edm::ESGetToken dspecToken_; + edm::ESGetToken mtdtopoToken_; + edm::ESGetToken mtdgeoToken_; + + std::stringstream sunitt; + constexpr static double tolerance{0.5e-3_mm}; +}; + +using DD3Vector = ROOT::Math::DisplacementVector3D>; +using angle_units::operators::convertRadToDeg; +using cms_rounding::roundIfNear0, cms_rounding::roundVecIfNear0; + +DD4hep_TestPixelTopology::DD4hep_TestPixelTopology(const edm::ParameterSet& iConfig) + : tag_(iConfig.getParameter("DDDetector")), + ddTopNodeName_(iConfig.getUntrackedParameter("ddTopNodeName", "BarrelTimingLayer")), + thisN_(), + btlNS_(), + etlNS_() { + dddetToken_ = esConsumes(tag_); + dspecToken_ = esConsumes(tag_); + mtdtopoToken_ = esConsumes(tag_); + mtdgeoToken_ = esConsumes(tag_); +} + +void DD4hep_TestPixelTopology::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) { + auto pDD = iSetup.getTransientHandle(dddetToken_); + + auto pSP = iSetup.getTransientHandle(dspecToken_); + + if (ddTopNodeName_ != "BarrelTimingLayer" && ddTopNodeName_ != "EndcapTimingLayer") { + edm::LogWarning("DD4hep_TestPixelTopology") << ddTopNodeName_ << "Not valid top MTD volume"; + return; + } + + if (!pDD.isValid()) { + edm::LogError("DD4hep_TestPixelTopology") << "ESTransientHandle pDD is not valid!"; + return; + } + if (pDD.description()) { + edm::LogVerbatim("DD4hep_TestPixelTopology") << pDD.description()->type_ << " label: " << pDD.description()->label_; + } else { + edm::LogPrint("DD4hep_TestPixelTopology") << "NO label found pDD.description() returned false."; + } + + if (!pSP.isValid()) { + edm::LogError("DD4hep_TestPixelTopology") << "ESTransientHandle pSP is not valid!"; + return; + } + + auto pTP = iSetup.getTransientHandle(mtdtopoToken_); + if (!pTP.isValid()) { + edm::LogError("DD4hep_TestPixelTopology") << "ESTransientHandle pTP is not valid!"; + return; + } else { + sunitt << "MTD topology mode = " << pTP.product()->getMTDTopologyMode() << " BtlLayout = " + << static_cast(MTDTopologyMode::crysLayoutFromTopoMode(pTP.product()->getMTDTopologyMode())) + << " EtlLayout = " + << static_cast(MTDTopologyMode::etlLayoutFromTopoMode(pTP.product()->getMTDTopologyMode())); + } + + auto pDG = iSetup.getTransientHandle(mtdgeoToken_); + if (!pDG.isValid()) { + edm::LogError("DD4hep_TestPixelTopology") << "ESTransientHandle pDG is not valid!"; + return; + } + + DDSpecParRefs specs; + std::string attribute("MtdDDStructure"), name; + bool isBarrel = false; + if (ddTopNodeName_ == "BarrelTimingLayer") { + sunitt << " BTL MTDGeometry:\n"; + name = "FastTimerHitsBarrel"; + isBarrel = true; + } else if (ddTopNodeName_ == "EndcapTimingLayer") { + sunitt << " ETL MTDGeometry:\n"; + name = "FastTimerHitsEndcap"; + } else { + edm::LogError("DD4hep_TestPixelTopology") << "No correct sensitive detector provided, abort" << ddTopNodeName_; + return; + } + pSP.product()->filter(specs, attribute, ddTopNodeName_); + attribute = "ReadOutName"; + pSP.product()->filter(specs, attribute, name); + + edm::LogVerbatim("DD4hep_TestPixelTopology").log([&specs](auto& log) { + log << "Filtered DD SpecPar Registry size: " << specs.size() << "\n"; + for (const auto& t : specs) { + log << "\nSpecPar " << t.first << ":\nRegExps { "; + for (const auto& ki : t.second->paths) + log << ki << " "; + log << "};\n "; + for (const auto& kl : t.second->spars) { + log << kl.first << " = "; + for (const auto& kil : kl.second) { + log << kil << " "; + } + log << "\n "; + } + } + }); + + DDFilteredView fv(pDD.product(), pDD.product()->description()->worldVolume()); + fv.mergedSpecifics(specs); + fv.firstChild(); + + bool write = false; + bool exitLoop = false; + uint32_t level(0); + uint32_t count(0); + uint32_t nSensBTL(0); + uint32_t nSensETL(0); + uint32_t oldgeoId(0); + + do { + theBaseNumber(fv); + + auto print_path = [&]() { + std::stringstream ss; + ss << " - OCMS[0]/"; + for (int ii = thisN_.getLevels() - 1; ii-- > 0;) { + ss << thisN_.getLevelName(ii); + ss << "["; + ss << thisN_.getCopyNumber(ii); + ss << "]/"; + } + return ss.str(); + }; + + if (level > 0 && fv.navPos().size() < level && count == 2) { + exitLoop = true; + } + if (dd4hep::dd::noNamespace(fv.name()) == ddTopNodeName_) { + write = true; + level = fv.navPos().size(); + count++; + } + + // Test only the desired subdetector + + if (exitLoop) { + break; + } + + if (write) { + // Actions for MTD volumes: search for sensitive detectors + + bool isSens = false; + + for (auto const& t : specs) { + for (auto const& kl : t.second->spars) { + if (kl.first == attribute) { + for (auto const& it : t.second->paths) { + if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(fv.name()), dd4hep::dd::realTopName(it))) { + isSens = true; + break; + } + } + } + } + } + + if (isSens) { + DetId theId, geoId; + BTLDetId theIdBTL, modIdBTL; + ETLDetId theIdETL, modIdETL; + if (isBarrel) { + theIdBTL = btlNS_.getUnitID(thisN_); + theId = theIdBTL; + geoId = theIdBTL.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(pTP.product()->getMTDTopologyMode())); + modIdBTL = geoId; + } else { + theIdETL = etlNS_.getUnitID(thisN_); + theId = theIdETL; + geoId = theIdETL.geographicalId(); + modIdETL = geoId; + } + + const MTDGeomDet* thedet = pDG.product()->idToDet(geoId); + + if (dynamic_cast((thedet)) == nullptr) { + throw cms::Exception("DD4hep_TestPixelTopology") + << "GeographicalID: " << std::hex << geoId.rawId() << " (" << theId.rawId() + << ") with invalid MTDGeomDetUnit!" << std::dec << std::endl; + } + + bool isNewId(false); + if (geoId != oldgeoId) { + oldgeoId = geoId; + isNewId = true; + if (isBarrel) { + nSensBTL++; + } else { + nSensETL++; + } + const GeomDetUnit theDetUnit = *(dynamic_cast(thedet)); + + if (isBarrel) { + sunitt << "geoId= " << modIdBTL.rawId() << " side= " << modIdBTL.mtdSide() + << " RU/mod= " << modIdBTL.globalRunit() << " / " << modIdBTL.module(); + } else { + sunitt << "geoId= " << modIdETL.rawId() << " side= " << modIdETL.mtdSide() + << " disc/face/sec= " << modIdETL.nDisc() << " / " << modIdETL.discSide() << " / " + << modIdETL.sector() << " mod/typ/sens= " << modIdETL.module() << " / " << modIdETL.modType() + << " / " << modIdETL.sensor(); + } + analyseRectangle(theDetUnit); + } + + if (thedet == nullptr) { + throw cms::Exception("DD4hep_TestPixelTopology") << "GeographicalID: " << std::hex << geoId.rawId() << " (" + << theId.rawId() << ") is invalid!" << std::dec << std::endl; + } + const ProxyMTDTopology& topoproxy = static_cast(thedet->topology()); + const RectangularMTDTopology& topo = static_cast(topoproxy.specificTopology()); + + int origRow(-1), origCol(-1), recoRow(-1), recoCol(-1); + if (isBarrel) { + origRow = theIdBTL.row(topo.nrows()); + origCol = theIdBTL.column(topo.nrows()); + } + + // + // Test of positions for sensitive detectors + // + + auto fround = [&](double in) { + std::stringstream ss; + ss << std::fixed << std::setw(14) << roundIfNear0(in); + return ss.str(); + }; + + if (!dd4hep::isA(fv.solid())) { + throw cms::Exception("DD4hep_TestPixelTopology") << "MTD sensitive element not a DDBox"; + break; + } + dd4hep::Box mySens(fv.solid()); + + double xoffset(0.); + double yoffset(0.); + if (!isBarrel) { + xoffset = topo.gapxBorder() + 0.5 * topo.pitch().first; + yoffset = topo.gapyBorder() + 0.5 * topo.pitch().second; + } + DD3Vector zeroLocal(0., 0., 0.); + DD3Vector cn1Local(mySens.x() - xoffset, mySens.y() - yoffset, mySens.z()); + DD3Vector cn2Local(-mySens.x() + xoffset, -mySens.y() + yoffset, -mySens.z()); + DD3Vector zeroGlobal = (fv.rotation())(zeroLocal) + fv.translation(); + DD3Vector cn1Global = (fv.rotation())(cn1Local) + fv.translation(); + DD3Vector cn2Global = (fv.rotation())(cn2Local) + fv.translation(); + + const size_t nTest(3); + std::array refLocalPoints{{Local3DPoint(zeroLocal.x(), zeroLocal.y(), zeroLocal.z()), + Local3DPoint(cn1Local.x(), cn1Local.y(), cn1Local.z()), + Local3DPoint(cn2Local.x(), cn2Local.y(), cn2Local.z())}}; + std::array refGlobalPoints{{zeroGlobal, cn1Global, cn2Global}}; + + for (size_t iloop = 0; iloop < nTest; iloop++) { + Local3DPoint cmRefLocal(convertMmToCm(refLocalPoints[iloop].x() / dd4hep::mm), + convertMmToCm(refLocalPoints[iloop].y() / dd4hep::mm), + convertMmToCm(refLocalPoints[iloop].z() / dd4hep::mm)); + + Local3DPoint modLocal, recoRefLocal; + if (isBarrel) { + // if BTL translate from crystal-local coordinates to module-local coordinates to get the row and column + modLocal = topo.pixelToModuleLocalPoint(cmRefLocal, origRow, origCol); + recoRefLocal = topo.moduleToPixelLocalPoint(modLocal); + const auto& thepixel = topo.pixelIndex(modLocal); + recoRow = thepixel.first; + recoCol = thepixel.second; + + if (origRow != recoRow || origCol != recoCol) { + sunitt << "DIFFERENCE row/col, orig= " << origRow << " " << origCol << " reco= " << recoRow << " " + << recoCol << "\n"; + recoRow = origRow; + recoCol = origCol; + } + } else { + // if ETL find the pixel corresponding to the referemce point, compute the pixel coordinate and convert back for check + modLocal = cmRefLocal; + const auto& thepixel = topo.pixelIndex(modLocal); + Local3DPoint pixLocal = topo.moduleToPixelLocalPoint(modLocal); + recoRefLocal = topo.pixelToModuleLocalPoint(pixLocal, thepixel.first, thepixel.second); + recoRow = thepixel.first; + recoCol = thepixel.second; + } + + // reconstructed global position from reco geometry and rectangluar MTD topology + + const auto& modGlobal = thedet->toGlobal(modLocal); + + if (isNewId && iloop == nTest - 1) { + sunitt << "row/col= " << recoRow << " / " << recoCol << " local pos= " << std::fixed << std::setw(14) + << roundVecIfNear0(modLocal) << " global pos= " << std::setw(14) << roundVecIfNear0(modGlobal) + << "\n"; + } + + const double deltax = convertCmToMm(modGlobal.x()) - (refGlobalPoints[iloop].x() / dd4hep::mm); + const double deltay = convertCmToMm(modGlobal.y()) - (refGlobalPoints[iloop].y() / dd4hep::mm); + const double deltaz = convertCmToMm(modGlobal.z()) - (refGlobalPoints[iloop].z() / dd4hep::mm); + + const double local_deltax = recoRefLocal.x() - cmRefLocal.x(); + const double local_deltay = recoRefLocal.y() - cmRefLocal.y(); + const double local_deltaz = recoRefLocal.z() - cmRefLocal.z(); + + if (std::abs(deltax) > tolerance || std::abs(deltay) > tolerance || std::abs(deltaz) > tolerance || + std::abs(local_deltax) > tolerance || std::abs(local_deltay) > tolerance || + std::abs(local_deltaz) > tolerance) { + sunitt << print_path() << "\n"; + if (isBarrel) { + sunitt << "rawId= " << theIdBTL.rawId() << " geoId= " << geoId.rawId() + << " side/rod= " << theIdBTL.mtdSide() << " / " << theIdBTL.mtdRR() + << " RU= " << theIdBTL.globalRunit() << " module/geomodule= " << theIdBTL.module() << " / " + << static_cast(geoId).module() << " crys= " << theIdBTL.crystal() + << " BTLDetId row/col= " << origRow << " / " << origCol << "\n"; + } else { + sunitt << "geoId= " << modIdETL.rawId() << " side= " << modIdETL.mtdSide() + << " disc/face/sec= " << modIdETL.nDisc() << " / " << modIdETL.discSide() << " / " + << modIdETL.sector() << " mod/typ/sens= " << modIdETL.module() << " / " << modIdETL.modType() + << " / " << modIdETL.sensor() << "\n"; + } + + sunitt << "Ref#" << iloop << " local= " << fround(refLocalPoints[iloop].x() / dd4hep::mm) + << fround(refLocalPoints[iloop].y() / dd4hep::mm) << fround(refLocalPoints[iloop].z() / dd4hep::mm) + << " Orig global= " << fround(refGlobalPoints[iloop].x() / dd4hep::mm) + << fround(refGlobalPoints[iloop].y() / dd4hep::mm) << fround(refGlobalPoints[iloop].z() / dd4hep::mm) + << " Reco global= " << fround(convertCmToMm(modGlobal.x())) << fround(convertCmToMm(modGlobal.y())) + << fround(convertCmToMm(modGlobal.z())) << " Delta= " << fround(deltax) << fround(deltay) + << fround(deltaz) << " Local Delta= " << fround(local_deltax) << fround(local_deltay) + << fround(local_deltaz) << "\n"; + + if (std::abs(deltax) > tolerance || std::abs(deltay) > tolerance || std::abs(deltaz) > tolerance) { + sunitt << "DIFFERENCE detId/ref# " << theId.rawId() << " " << iloop << " dx/dy/dz= " << fround(deltax) + << fround(deltay) << fround(deltaz) << "\n"; + } + if (std::abs(local_deltax) > tolerance || std::abs(local_deltay) > tolerance || + std::abs(local_deltaz) > tolerance) { + sunitt << "DIFFERENCE detId/ref# " << theId.rawId() << " " << iloop + << " local dx/dy/dz= " << fround(local_deltax) << fround(local_deltay) << fround(local_deltaz) + << "\n"; + } + } + } + } + } + } while (fv.next(0)); + + if (isBarrel && nSensBTL != pDG.product()->detsBTL().size()) { + sunitt << "DIFFERENCE #ideal = " << nSensBTL << " #reco = " << pDG.product()->detsBTL().size() + << " BTL module numbers are not matching!"; + } + + if (!isBarrel && nSensETL != pDG.product()->detsETL().size()) { + sunitt << "DIFFERENCE #ideal = " << nSensETL << " #reco = " << pDG.product()->detsBTL().size() + << " ETL module numbers are not matching!"; + } + + if (!sunitt.str().empty()) { + edm::LogVerbatim("DD4hep_TestPixelTopology") << sunitt.str(); + edm::LogVerbatim("MTDUnitTest") << sunitt.str(); + } +} + +void DD4hep_TestPixelTopology::theBaseNumber(cms::DDFilteredView& fv) { + thisN_.reset(); + thisN_.setSize(fv.navPos().size()); + + for (uint ii = 0; ii < fv.navPos().size(); ii++) { + std::string_view name((fv.geoHistory()[ii])->GetName()); + size_t ipos = name.rfind('_'); + thisN_.addLevel(name.substr(0, ipos), fv.copyNos()[ii]); +#ifdef EDM_ML_DEBUG + edm::LogVerbatim("DD4hep_TestPixelTopology") << ii << " " << name.substr(0, ipos) << " " << fv.copyNos()[ii]; +#endif + } +} + +void DD4hep_TestPixelTopology::analyseRectangle(const GeomDetUnit& det) { + const double safety = 0.9999; + + const BoundPlane& p = det.specificSurface(); + const Bounds& bounds = det.surface().bounds(); + const RectangularPlaneBounds* tb = dynamic_cast(&bounds); + if (tb == nullptr) + return; // not trapezoidal + + const GlobalPoint& pos = det.position(); + double length = tb->length(); + double width = tb->width(); + double thickness = tb->thickness(); + + GlobalVector yShift = det.surface().toGlobal(LocalVector(0, 0, safety * length / 2.)); + GlobalPoint outerMiddle = pos + yShift; + GlobalPoint innerMiddle = pos + (-1. * yShift); + if (outerMiddle.perp() < innerMiddle.perp()) + std::swap(outerMiddle, innerMiddle); + + auto fround = [&](double in) { + std::stringstream ss; + ss << std::fixed << std::setw(14) << roundIfNear0(in); + return ss.str(); + }; + + auto fvecround = [&](GlobalPoint vecin) { + std::stringstream ss; + ss << std::fixed << std::setw(14) << roundVecIfNear0(vecin); + return ss.str(); + }; + + sunitt << " " << fvecround(pos) << " R= " << fround(std::sqrt(pos.x() * pos.x() + pos.y() * pos.y())) + << " phi= " << fround(convertRadToDeg(pos.phi())) << " outerMiddle " << fvecround(outerMiddle) << "\n" + << " l/w/t " << fround(length) << " / " << fround(width) << " / " << fround(thickness) + << " RadLeng= " << p.mediumProperties().radLen() << " Xi= " << p.mediumProperties().xi() + << " det center inside bounds? " << tb->inside(det.surface().toLocal(pos)) << "\n"; + + checkRotation(det); +} + +void DD4hep_TestPixelTopology::checkRotation(const GeomDetUnit& det) { + const double eps = std::numeric_limits::epsilon(); + static int first = 0; + if (first == 0) { + edm::LogVerbatim("DD4hep_TestPixelTopology") + << "numeric_limits::epsilon() " << std::numeric_limits::epsilon(); + first = 1; + } + + const Surface::RotationType& rot(det.surface().rotation()); + GlobalVector a(rot.xx(), rot.xy(), rot.xz()); + GlobalVector b(rot.yx(), rot.yy(), rot.yz()); + GlobalVector c(rot.zx(), rot.zy(), rot.zz()); + GlobalVector cref = a.cross(b); + GlobalVector aref = b.cross(c); + GlobalVector bref = c.cross(a); + if ((a - aref).mag() > eps || (b - bref).mag() > eps || (c - cref).mag() > eps) { + sunitt << " Rotation not good by cross product: " << (a - aref).mag() << ", " << (b - bref).mag() << ", " + << (c - cref).mag() << " for det at pos " << det.surface().position() << "\n"; + } + if (fabs(a.mag() - 1.) > eps || fabs(b.mag() - 1.) > eps || fabs(c.mag() - 1.) > eps) { + sunitt << " Rotation not good by bector mag: " << (a).mag() << ", " << (b).mag() << ", " << (c).mag() + << " for det at pos " << det.surface().position() << "\n"; + } +} + +DEFINE_FWK_MODULE(DD4hep_TestPixelTopology); diff --git a/Geometry/MTDGeometryBuilder/test/MTDDigiGeometryAnalyzer.cc b/Geometry/MTDGeometryBuilder/test/MTDDigiGeometryAnalyzer.cc index dc80d7eeb85d7..cc5852034ac6d 100644 --- a/Geometry/MTDGeometryBuilder/test/MTDDigiGeometryAnalyzer.cc +++ b/Geometry/MTDGeometryBuilder/test/MTDDigiGeometryAnalyzer.cc @@ -40,8 +40,6 @@ class MTDDigiGeometryAnalyzer : public edm::one::EDAnalyzer<> { void endJob() override {} private: - void analyseRectangle(const GeomDetUnit& det); - void checkRotation(const GeomDetUnit& det); void checkRectangularMTDTopology(const RectangularMTDTopology&); void checkPixelsAcceptance(const GeomDetUnit& det); @@ -62,161 +60,51 @@ void MTDDigiGeometryAnalyzer::analyze(const edm::Event& iEvent, const edm::Event // get the MTDGeometry // auto pDD = iSetup.getTransientHandle(mtdgeoToken_); - edm::LogInfo("MTDDigiGeometryAnalyzer") - << "Geometry node for MTDGeom is " << &(*pDD) << "\n" - << " # detectors = " << pDD->detUnits().size() << "\n" - << " # types = " << pDD->detTypes().size() << "\n" - << " # BTL dets = " << pDD->detsBTL().size() << "\n" - << " # ETL dets = " << pDD->detsETL().size() << "\n" - << " # layers " << pDD->geomDetSubDetector(1) << " = " << pDD->numberOfLayers(1) << "\n" - << " # layers " << pDD->geomDetSubDetector(2) << " = " << pDD->numberOfLayers(2) << "\n"; - sunitt << std::fixed << std::setw(7) << pDD->detUnits().size() << std::setw(7) << pDD->detTypes().size() << "\n"; - for (auto const& it : pDD->detUnits()) { - if (dynamic_cast((it)) != nullptr) { - const BoundPlane& p = (dynamic_cast((it)))->specificSurface(); - const MTDDetId mtdId(it->geographicalId()); - std::stringstream moduleLabel; - if (mtdId.mtdSubDetector() == 1) { - const BTLDetId btlId(it->geographicalId()); - moduleLabel << " BTL side " << btlId.mtdSide() << " Rod " << btlId.mtdRR() << " type/RU " << btlId.modType() - << "/" << btlId.runit() << " mod " << btlId.module(); - } else if (mtdId.mtdSubDetector() == 2) { - const ETLDetId etlId(it->geographicalId()); - moduleLabel << " ETL side " << mtdId.mtdSide() << " Disc/Side/Sector " << etlId.nDisc() << " " - << etlId.discSide() << " " << etlId.sector(); - } else { - edm::LogWarning("MTDDigiGeometryanalyzer") << (it->geographicalId()).rawId() << " unknown MTD subdetector!"; - } - edm::LogVerbatim("MTDDigiGeometryAnalyzer") - << "---------------------------------------------------------- \n" - << it->geographicalId().rawId() << moduleLabel.str() << " RadLeng Pixel " << p.mediumProperties().radLen() - << " Xi Pixel " << p.mediumProperties().xi(); - - const GeomDetUnit theDet = *(dynamic_cast(it)); - analyseRectangle(theDet); - } - } + sunitt << "MTDGeometry:\n" + << " # detectors = " << pDD->detUnits().size() << "\n" + << " # types = " << pDD->detTypes().size() << "\n" + << " # BTL dets = " << pDD->detsBTL().size() << "\n" + << " # ETL dets = " << pDD->detsETL().size() << "\n" + << " # layers " << pDD->geomDetSubDetector(1) << " = " << pDD->numberOfLayers(1) << "\n" + << " # layers " << pDD->geomDetSubDetector(2) << " = " << pDD->numberOfLayers(2) << "\n" + << " # dets = " << pDD->dets().size() << "\n" + << " # detUnitIds = " << pDD->detUnitIds().size() << "\n" + << " # detIds = " << pDD->detIds().size() << "\n"; for (auto const& it : pDD->detTypes()) { if (dynamic_cast((it)) != nullptr) { const PixelTopology& p = (dynamic_cast((it)))->specificTopology(); const RectangularMTDTopology& topo = static_cast(p); auto pitchval = topo.pitch(); - edm::LogVerbatim("MTDDigiGeometryAnalyzer") - << "\n Subdetector " << it->subDetector() << " MTD Det " << it->name() << "\n" - << " Rows " << topo.nrows() << " Columns " << topo.ncolumns() << " ROCS X " << topo.rocsX() - << " ROCS Y " << topo.rocsY() << " Rows/ROC " << topo.rowsperroc() << " Cols/ROC " << topo.colsperroc() - << " Pitch X " << pitchval.first << " Pitch Y " << pitchval.second << " Sensor Interpad X " - << topo.gapxInterpad() << " Sensor Interpad Y " << topo.gapyInterpad() << " Sensor Border X " - << topo.gapxBorder() << " Sensor Border Y " << topo.gapyBorder(); - sunitt << std::fixed << std::setw(7) << it->subDetector() << std::setw(4) << topo.nrows() << std::setw(4) - << topo.ncolumns() << std::setw(4) << std::setw(4) << topo.rocsX() << std::setw(4) << topo.rocsY() - << std::setw(4) << topo.rowsperroc() << std::setw(4) << topo.colsperroc() << std::setw(10) - << pitchval.first << std::setw(10) << pitchval.second << std::setw(10) << topo.gapxInterpad() - << std::setw(10) << topo.gapyInterpad() << std::setw(10) << topo.gapxBorder() << std::setw(10) - << topo.gapyBorder() << "\n"; + sunitt << "\n Subdetector " << it->subDetector() << " MTD Det " << it->name() << "\n" + << " Rows " << topo.nrows() << " Columns " << topo.ncolumns() << " ROCS X " << topo.rocsX() + << " ROCS Y " << topo.rocsY() << " Rows/ROC " << topo.rowsperroc() << " Cols/ROC " << topo.colsperroc() + << " Pitch X " << pitchval.first << " Pitch Y " << pitchval.second << " Sensor Interpad X " + << topo.gapxInterpad() << " Sensor Interpad Y " << topo.gapyInterpad() << " Sensor Border X " + << topo.gapxBorder() << " Sensor Border Y " << topo.gapyBorder() << "\n"; checkRectangularMTDTopology(topo); } } - edm::LogInfo("MTDDigiGeometryAnalyzer") << "Acceptance of BTL module:"; + sunitt << "\nAcceptance of BTL module:"; auto const& btldet = *(dynamic_cast(pDD->detsBTL().front())); checkPixelsAcceptance(btldet); - edm::LogInfo("MTDDigiGeometryAnalyzer") << "Acceptance of ETL module:"; + sunitt << "\nAcceptance of ETL module:"; auto const& etldet = *(dynamic_cast(pDD->detsETL().front())); checkPixelsAcceptance(etldet); - edm::LogInfo("MTDDigiGeometryAnalyzer") << "Additional MTD geometry content:" - << "\n" - << " # dets = " << pDD->dets().size() << "\n" - << " # detUnitIds = " << pDD->detUnitIds().size() << "\n" - << " # detIds = " << pDD->detIds().size() << "\n"; - sunitt << std::fixed << std::setw(7) << pDD->dets().size() << std::setw(7) << pDD->detUnitIds().size() << std::setw(7) - << pDD->detIds().size() << "\n"; - + edm::LogVerbatim("MTDDigiGeometryAnalyzer") << sunitt.str(); edm::LogVerbatim("MTDUnitTest") << sunitt.str(); } void MTDDigiGeometryAnalyzer::checkRectangularMTDTopology(const RectangularMTDTopology& topo) { - std::stringstream pixelinfo; - pixelinfo << "Pixel center location:\n"; + sunitt << "Pixel center location:\n"; LocalPoint center(0, 0, 0); for (int r = 0; r < topo.nrows(); r++) { for (int c = 0; c < topo.ncolumns(); c++) { - sunitt << r << " " << c << " " << topo.pixelToModuleLocalPoint(center, r, c) << "\n"; - pixelinfo << r << " " << c << " " << topo.pixelToModuleLocalPoint(center, r, c) << "\n"; + sunitt << std::setw(7) << r << std::setw(7) << c << " " << topo.pixelToModuleLocalPoint(center, r, c) << "\n"; } } - edm::LogVerbatim("MTDDigiGeometryAnalyzer") << pixelinfo.str(); -} - -void MTDDigiGeometryAnalyzer::analyseRectangle(const GeomDetUnit& det) { - const double safety = 0.9999; - - const Bounds& bounds = det.surface().bounds(); - const RectangularPlaneBounds* tb = dynamic_cast(&bounds); - if (tb == nullptr) - return; // not trapezoidal - - const GlobalPoint& pos = det.position(); - double length = tb->length(); - double width = tb->width(); - double thickness = tb->thickness(); - - GlobalVector yShift = det.surface().toGlobal(LocalVector(0, 0, safety * length / 2.)); - GlobalPoint outerMiddle = pos + yShift; - GlobalPoint innerMiddle = pos + (-1. * yShift); - if (outerMiddle.perp() < innerMiddle.perp()) - std::swap(outerMiddle, innerMiddle); - - auto fround = [&](double in) { - std::stringstream ss; - ss << std::fixed << std::setw(14) << roundIfNear0(in); - return ss.str(); - }; - - auto fvecround = [&](GlobalPoint vecin) { - std::stringstream ss; - ss << std::fixed << std::setw(14) << roundVecIfNear0(vecin); - return ss.str(); - }; - - edm::LogVerbatim("MTDDigiGeometryAnalyzer") - << "Det at pos " << fvecround(pos) << " radius " << fround(std::sqrt(pos.x() * pos.x() + pos.y() * pos.y())) - << " has length " << fround(length) << " width " << fround(width) << " thickness " << fround(thickness) << "\n" - << "det center inside bounds? " << tb->inside(det.surface().toLocal(pos)) << "\n" - << "outerMiddle " << fvecround(outerMiddle); - sunitt << det.geographicalId().rawId() << fvecround(pos) << fround(length) << fround(width) << fround(thickness) - << tb->inside(det.surface().toLocal(pos)) << fvecround(outerMiddle) << "\n"; - - checkRotation(det); -} - -void MTDDigiGeometryAnalyzer::checkRotation(const GeomDetUnit& det) { - const double eps = std::numeric_limits::epsilon(); - static int first = 0; - if (first == 0) { - edm::LogVerbatim("MTDDigiGeometryAnalyzer") - << "numeric_limits::epsilon() " << std::numeric_limits::epsilon(); - first = 1; - } - - const Surface::RotationType& rot(det.surface().rotation()); - GlobalVector a(rot.xx(), rot.xy(), rot.xz()); - GlobalVector b(rot.yx(), rot.yy(), rot.yz()); - GlobalVector c(rot.zx(), rot.zy(), rot.zz()); - GlobalVector cref = a.cross(b); - GlobalVector aref = b.cross(c); - GlobalVector bref = c.cross(a); - if ((a - aref).mag() > eps || (b - bref).mag() > eps || (c - cref).mag() > eps) { - edm::LogWarning("MTDDigiGeometryAnalyzer") - << " Rotation not good by cross product: " << (a - aref).mag() << ", " << (b - bref).mag() << ", " - << (c - cref).mag() << " for det at pos " << det.surface().position(); - } - if (fabs(a.mag() - 1.) > eps || fabs(b.mag() - 1.) > eps || fabs(c.mag() - 1.) > eps) { - edm::LogWarning("MTDDigiGeometryAnalyzer") << " Rotation not good by bector mag: " << (a).mag() << ", " << (b).mag() - << ", " << (c).mag() << " for det at pos " << det.surface().position(); - } } void MTDDigiGeometryAnalyzer::checkPixelsAcceptance(const GeomDetUnit& det) { @@ -227,7 +115,7 @@ void MTDDigiGeometryAnalyzer::checkPixelsAcceptance(const GeomDetUnit& det) { double length = tb->length(); double width = tb->width(); - edm::LogVerbatim("MTDDigiGeometryAnalyzer") << "X (width) = " << width << " Y (length) = " << length; + sunitt << " X (width) = " << width << " Y (length) = " << length; const ProxyMTDTopology& topoproxy = static_cast(det.topology()); const RectangularMTDTopology& topo = static_cast(topoproxy.specificTopology()); @@ -246,7 +134,7 @@ void MTDDigiGeometryAnalyzer::checkPixelsAcceptance(const GeomDetUnit& det) { } double acc = (double)inpixel / (double)maxindex; double accerr = std::sqrt(acc * (1. - acc) / (double)maxindex); - edm::LogVerbatim("MTDDigiGeometryAnalyzer") << "Acceptance: " << acc << " +/- " << accerr; + sunitt << " Acceptance: " << acc << " +/- " << accerr; } //define this as a plug-in diff --git a/Geometry/MTDGeometryBuilder/test/dd4hep_mtd_cfg.py b/Geometry/MTDGeometryBuilder/test/dd4hep_mtd_cfg.py deleted file mode 100644 index 42df806bdad16..0000000000000 --- a/Geometry/MTDGeometryBuilder/test/dd4hep_mtd_cfg.py +++ /dev/null @@ -1,64 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -import Configuration.Geometry.defaultPhase2ConditionsEra_cff as _settings -_PH2_GLOBAL_TAG, _PH2_ERA = _settings.get_era_and_conditions(_settings.DEFAULT_VERSION) - -from Configuration.ProcessModifiers.dd4hep_cff import dd4hep - -process = cms.Process("GeometryTest",_PH2_ERA,dd4hep) - -process.source = cms.Source("EmptyIOVSource", - lastValue = cms.uint64(1), - timetype = cms.string('runnumber'), - firstValue = cms.uint64(1), - interval = cms.uint64(1) - ) - -process.maxEvents = cms.untracked.PSet( - input = cms.untracked.int32(1) -) - -process.load("FWCore.MessageLogger.MessageLogger_cfi") -process.MessageLogger.cerr.threshold = cms.untracked.string('INFO') -process.MessageLogger.cerr.INFO = cms.untracked.PSet( - limit = cms.untracked.int32(0) -) -process.MessageLogger.cerr.MTDDigiGeometryAnalyzer = cms.untracked.PSet( - limit = cms.untracked.int32(-1) -) -process.MessageLogger.cerr.DD4hep_TestBTLPixelTopology = cms.untracked.PSet( - limit = cms.untracked.int32(-1) -) -process.MessageLogger.files.mtdGeometryDD4hep = cms.untracked.PSet( - DEBUG = cms.untracked.PSet( - limit = cms.untracked.int32(0) - ), - ERROR = cms.untracked.PSet( - limit = cms.untracked.int32(0) - ), - FWKINFO = cms.untracked.PSet( - limit = cms.untracked.int32(0) - ), - INFO = cms.untracked.PSet( - limit = cms.untracked.int32(0) - ), - MTDUnitTest = cms.untracked.PSet( - limit = cms.untracked.int32(-1) - ), - WARNING = cms.untracked.PSet( - limit = cms.untracked.int32(0) - ), - noLineBreaks = cms.untracked.bool(True), - threshold = cms.untracked.string('INFO') -) - -process.load("Configuration.Geometry.GeometryDD4hepExtendedRun4DefaultReco_cff") - -process.Timing = cms.Service("Timing") - -process.prod = cms.EDAnalyzer("MTDDigiGeometryAnalyzer") -process.prod1 = cms.EDAnalyzer("DD4hep_TestBTLPixelTopology", - DDDetector = cms.ESInputTag('',''), -) - -process.p1 = cms.Path(process.prod+process.prod1) diff --git a/Geometry/MTDGeometryBuilder/test/mtd_cfg.py b/Geometry/MTDGeometryBuilder/test/mtd_cfg.py index 59e7243471ceb..54cb0715bc0d5 100644 --- a/Geometry/MTDGeometryBuilder/test/mtd_cfg.py +++ b/Geometry/MTDGeometryBuilder/test/mtd_cfg.py @@ -3,7 +3,9 @@ import Configuration.Geometry.defaultPhase2ConditionsEra_cff as _settings _PH2_GLOBAL_TAG, _PH2_ERA = _settings.get_era_and_conditions(_settings.DEFAULT_VERSION) -process = cms.Process("GeometryTest", _PH2_ERA) +from Configuration.ProcessModifiers.dd4hep_cff import dd4hep + +process = cms.Process("GeometryTest",_PH2_ERA,dd4hep) process.source = cms.Source("EmptyIOVSource", lastValue = cms.uint64(1), @@ -22,9 +24,15 @@ limit = cms.untracked.int32(0) ) process.MessageLogger.cerr.MTDDigiGeometryAnalyzer = cms.untracked.PSet( - limit = cms.untracked.int32(-1) + limit = cms.untracked.int32(0) +) +process.MessageLogger.cerr.DD4hep_TestPixelTopology = cms.untracked.PSet( + limit = cms.untracked.int32(0) +) +process.MessageLogger.cerr.MTDUnitTest = cms.untracked.PSet( + limit = cms.untracked.int32(0) ) -process.MessageLogger.files.mtdGeometryDDD = cms.untracked.PSet( +process.MessageLogger.files.mtdGeometryDD4hep = cms.untracked.PSet( DEBUG = cms.untracked.PSet( limit = cms.untracked.int32(0) ), @@ -47,18 +55,18 @@ threshold = cms.untracked.string('INFO') ) -process.load("Configuration.Geometry.GeometryExtendedRun4Default_cff") - -process.load("Geometry.MTDNumberingBuilder.mtdNumberingGeometry_cff") - -process.load("Geometry.MTDNumberingBuilder.mtdTopology_cfi") -process.load("Geometry.MTDGeometryBuilder.mtdParameters_cff") - -process.load("Geometry.MTDGeometryBuilder.mtdGeometry_cfi") -process.mtdGeometry.applyAlignment = cms.bool(False) +process.load("Configuration.Geometry.GeometryDD4hepExtendedRun4DefaultReco_cff") process.Timing = cms.Service("Timing") process.prod = cms.EDAnalyzer("MTDDigiGeometryAnalyzer") +process.prod1 = cms.EDAnalyzer("DD4hep_TestPixelTopology", + DDDetector = cms.ESInputTag('',''), + ddTopNodeName = cms.untracked.string('BarrelTimingLayer') +) +process.prod2 = cms.EDAnalyzer("DD4hep_TestPixelTopology", + DDDetector = cms.ESInputTag('',''), + ddTopNodeName = cms.untracked.string('EndcapTimingLayer') +) -process.p1 = cms.Path(process.prod) +process.p1 = cms.Path(cms.wait(process.prod)+cms.wait(process.prod1)+process.prod2) diff --git a/Geometry/MTDGeometryBuilder/test/runTest.sh b/Geometry/MTDGeometryBuilder/test/runTest.sh index bc2127a23aed7..825b1813daa39 100755 --- a/Geometry/MTDGeometryBuilder/test/runTest.sh +++ b/Geometry/MTDGeometryBuilder/test/runTest.sh @@ -14,7 +14,6 @@ function checkDiff { TEST_DIR=$CMSSW_BASE/src/Geometry/MTDGeometryBuilder/test F1=${TEST_DIR}/mtd_cfg.py -F2=${TEST_DIR}/dd4hep_mtd_cfg.py REF_FILE="Geometry/TestReference/data/mtdGeometryRef.log.gz" REF="" @@ -26,8 +25,7 @@ for d in $(echo $CMSSW_SEARCH_PATH | tr ':' '\n') ; do done [ -z $REF ] && exit 1 -FILE1=mtdGeometryDDD.log -FILE2=mtdGeometryDD4hep.log +FILE1=mtdGeometryDD4hep.log LOG=mtdgblog DIF=mtdgbdif @@ -39,11 +37,3 @@ rm -f $LOG $DIF $FILE1 cmsRun $F1 >& $LOG || die "Failure using cmsRun $F1" $? gzip -f $FILE1 || die "$FILE1 compression fail" $? (zdiff $FILE1.gz $REF >& $DIF || [ -s $DIF ] && checkDiff $DIF || echo "OK") || die "Failure in comparison for $FILE1" $? - -rm -f $LOG $DIF $FILE2 -echo "===== Test \"cmsRun dd4hep_mtd_cfg.py\" ====" - -cmsRun $F2 >& $LOG || die "Failure using cmsRun $F2" $? -gzip -f $FILE2 || die "$FILE2 compression fail" $? -(zdiff $FILE2.gz $REF >& $DIF || [ -s $DIF ] && checkDiff $DIF || echo "OK") || die "Failure in comparison for $FILE2" $? - From a07ea53e1386725b09a250d3aea94022041102e3 Mon Sep 17 00:00:00 2001 From: Fabio Cossutti Date: Tue, 12 Nov 2024 23:23:38 +0100 Subject: [PATCH 2/2] Speed up filtering of sensitive volume --- .../test/DD4hep_TestPixelTopology.cc | 25 +++++++++++-------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/Geometry/MTDGeometryBuilder/test/DD4hep_TestPixelTopology.cc b/Geometry/MTDGeometryBuilder/test/DD4hep_TestPixelTopology.cc index 3cbd0f3adbd36..156503f34650d 100644 --- a/Geometry/MTDGeometryBuilder/test/DD4hep_TestPixelTopology.cc +++ b/Geometry/MTDGeometryBuilder/test/DD4hep_TestPixelTopology.cc @@ -172,6 +172,17 @@ void DD4hep_TestPixelTopology::analyze(const edm::Event& iEvent, const edm::Even } }); + std::vector filterName; + for (auto const& t : specs) { + for (auto const& kl : t.second->spars) { + if (kl.first == attribute) { + for (auto const& it : t.second->paths) { + filterName.emplace_back(it); + } + } + } + } + DDFilteredView fv(pDD.product(), pDD.product()->description()->worldVolume()); fv.mergedSpecifics(specs); fv.firstChild(); @@ -219,16 +230,10 @@ void DD4hep_TestPixelTopology::analyze(const edm::Event& iEvent, const edm::Even bool isSens = false; - for (auto const& t : specs) { - for (auto const& kl : t.second->spars) { - if (kl.first == attribute) { - for (auto const& it : t.second->paths) { - if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(fv.name()), dd4hep::dd::realTopName(it))) { - isSens = true; - break; - } - } - } + for (auto const& it : filterName) { + if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(fv.name()), dd4hep::dd::realTopName(it))) { + isSens = true; + break; } }