MTDGeometryBuilder/test/DD4hep_TestBTLPixelTopology.cc

0001 //#define EDM_ML_DEBUG
0002
0003 #include <iostream>
0004 #include <fstream>
0005 #include <string>
0006 #include <utility>
0007 #include <vector>
0008 #include <algorithm>
0009
0010 #include "FWCore/Framework/interface/Frameworkfwd.h"
0011 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/Framework/interface/ESTransientHandle.h"
0015 #include "FWCore/Framework/interface/MakerMacros.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/Utilities/interface/Exception.h"
0018
0019 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0020 #include "Geometry/Records/interface/DDSpecParRegistryRcd.h"
0021
0022 #include "DetectorDescription/DDCMS/interface/DDDetector.h"
0023 #include "DetectorDescription/DDCMS/interface/DDSolidShapes.h"
0024 #include "DetectorDescription/DDCMS/interface/DDFilteredView.h"
0025 #include "DetectorDescription/DDCMS/interface/DDSpecParRegistry.h"
0026
0027 #include "Geometry/MTDCommonData/interface/MTDTopologyMode.h"
0028 #include "Geometry/MTDCommonData/interface/MTDBaseNumber.h"
0029 #include "Geometry/MTDCommonData/interface/BTLNumberingScheme.h"
0030 #include "Geometry/MTDCommonData/interface/ETLNumberingScheme.h"
0031
0032 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0033 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0034
0035 #include "Geometry/Records/interface/MTDTopologyRcd.h"
0036 #include "Geometry/MTDNumberingBuilder/interface/MTDTopology.h"
0037 #include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h"
0038 #include "Geometry/Records/interface/MTDDigiGeometryRecord.h"
0039 #include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h"
0040 #include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
0041
0042 #include "DataFormats/Math/interface/GeantUnits.h"
0043 #include "DataFormats/Math/interface/Rounding.h"
0044 #include <DD4hep/DD4hepUnits.h>
0045
0046 using namespace cms;
0047 using namespace geant_units::operators;
0048
0049 class DD4hep_TestBTLPixelTopology : public edm::one::EDAnalyzer<> {
0050 public:
0051   explicit DD4hep_TestBTLPixelTopology(const edm::ParameterSet&);
0052   ~DD4hep_TestBTLPixelTopology() = default;
0053
0054   void beginJob() override {}
0055   void analyze(edm::Event const&, edm::EventSetup const&) override;
0056   void endJob() override {}
0057
0058   void theBaseNumber(cms::DDFilteredView& fv);
0059
0060 private:
0061   const edm::ESInputTag tag_;
0062
0063   MTDBaseNumber thisN_;
0064   BTLNumberingScheme btlNS_;
0065
0066   edm::ESGetToken<DDDetector, IdealGeometryRecord> dddetToken_;
0067   edm::ESGetToken<DDSpecParRegistry, DDSpecParRegistryRcd> dspecToken_;
0068   edm::ESGetToken<MTDTopology, MTDTopologyRcd> mtdtopoToken_;
0069   edm::ESGetToken<MTDGeometry, MTDDigiGeometryRecord> mtdgeoToken_;
0070
0071   std::stringstream sunitt;
0072   constexpr static double tolerance{0.5e-3_mm};
0073 };
0074
0075 using DD3Vector = ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double>>;
0076 using cms_rounding::roundIfNear0;
0077
0078 DD4hep_TestBTLPixelTopology::DD4hep_TestBTLPixelTopology(const edm::ParameterSet& iConfig)
0079     : tag_(iConfig.getParameter<edm::ESInputTag>("DDDetector")), thisN_(), btlNS_() {
0080   dddetToken_ = esConsumes<DDDetector, IdealGeometryRecord>(tag_);
0081   dspecToken_ = esConsumes<DDSpecParRegistry, DDSpecParRegistryRcd>(tag_);
0082   mtdtopoToken_ = esConsumes<MTDTopology, MTDTopologyRcd>(tag_);
0083   mtdgeoToken_ = esConsumes<MTDGeometry, MTDDigiGeometryRecord>(tag_);
0084 }
0085
0086 void DD4hep_TestBTLPixelTopology::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0087   auto pDD = iSetup.getTransientHandle(dddetToken_);
0088
0089   auto pSP = iSetup.getTransientHandle(dspecToken_);
0090
0091   if (!pDD.isValid()) {
0092     edm::LogError("DD4hep_TestBTLPixelTopology") << "ESTransientHandle<DDCompactView> pDD is not valid!";
0093     return;
0094   }
0095   if (pDD.description()) {
0096     edm::LogInfo("DD4hep_TestBTLPixelTopology") << pDD.description()->type_ << " label: " << pDD.description()->label_;
0097   } else {
0098     edm::LogWarning("DD4hep_TestBTLPixelTopology") << "NO label found pDD.description() returned false.";
0099   }
0100
0101   if (!pSP.isValid()) {
0102     edm::LogError("DD4hep_TestBTLPixelTopology") << "ESTransientHandle<DDSpecParRegistry> pSP is not valid!";
0103     return;
0104   }
0105
0106   auto pTP = iSetup.getTransientHandle(mtdtopoToken_);
0107   if (!pTP.isValid()) {
0108     edm::LogError("DD4hep_TestBTLPixelTopology") << "ESTransientHandle<MTDTopology> pTP is not valid!";
0109     return;
0110   } else {
0111     edm::LogInfo("DD4hep_TestBTLPixelTopology")
0112         << "MTD topology mode = " << pTP.product()->getMTDTopologyMode() << " BTLDetId:CrysLayout = "
0113         << static_cast<int>(MTDTopologyMode::crysLayoutFromTopoMode(pTP.product()->getMTDTopologyMode()));
0114   }
0115
0116   auto pDG = iSetup.getTransientHandle(mtdgeoToken_);
0117   if (!pDG.isValid()) {
0118     edm::LogError("DD4hep_TestBTLPixelTopology") << "ESTransientHandle<MTDGeometry> pDG is not valid!";
0119     return;
0120   } else {
0121     edm::LogInfo("DD4hep_TestBTLPixelTopology")
0122         << "Geometry node for MTDGeom is  " << &(*pDG) << "\n"
0123         << " # detectors = " << pDG.product()->detUnits().size() << "\n"
0124         << " # types     = " << pDG.product()->detTypes().size() << "\n"
0125         << " # BTL dets  = " << pDG.product()->detsBTL().size() << "\n"
0126         << " # ETL dets  = " << pDG.product()->detsETL().size() << "\n"
0127         << " # layers " << pDG.product()->geomDetSubDetector(1) << "  = " << pDG.product()->numberOfLayers(1) << "\n"
0128         << " # layers " << pDG.product()->geomDetSubDetector(2) << "  = " << pDG.product()->numberOfLayers(2) << "\n";
0129   }
0130
0131   DDFilteredView fv(pDD.product(), pDD.product()->description()->worldVolume());
0132   fv.next(0);
0133
0134   DDSpecParRefs specs;
0135   std::string attribute("ReadOutName"), name("FastTimerHitsBarrel");
0136   pSP.product()->filter(specs, attribute, name);
0137
0138   edm::LogVerbatim("DD4hep_TestBTLPixelTopology").log([&specs](auto& log) {
0139     log << "Filtered DD SpecPar Registry size: " << specs.size() << "\n";
0140     for (const auto& t : specs) {
0141       log << "\nSpecPar " << t.first << ":\nRegExps { ";
0142       for (const auto& ki : t.second->paths)
0143         log << ki << " ";
0144       log << "};\n ";
0145       for (const auto& kl : t.second->spars) {
0146         log << kl.first << " = ";
0147         for (const auto& kil : kl.second) {
0148           log << kil << " ";
0149         }
0150         log << "\n ";
0151       }
0152     }
0153   });
0154
0155   bool isBarrel = true;
0156   bool exitLoop = false;
0157   uint32_t level(0);
0158
0159   do {
0160     if (dd4hep::dd::noNamespace(fv.name()) == "BarrelTimingLayer") {
0161       isBarrel = true;
0162       edm::LogInfo("DD4hep_TestBTLPixelTopology") << "isBarrel = " << isBarrel;
0163     } else if (dd4hep::dd::noNamespace(fv.name()) == "EndcapTimingLayer") {
0164       isBarrel = false;
0165       edm::LogInfo("DD4hep_TestBTLPixelTopology") << "isBarrel = " << isBarrel;
0166     }
0167
0168     theBaseNumber(fv);
0169
0170     auto print_path = [&]() {
0171       std::stringstream ss;
0172       ss << " - OCMS[0]/";
0173       for (int ii = thisN_.getLevels() - 1; ii-- > 0;) {
0174         ss << thisN_.getLevelName(ii);
0175         ss << "[";
0176         ss << thisN_.getCopyNumber(ii);
0177         ss << "]/";
0178       }
0179       return ss.str();
0180     };
0181
0182     if (level > 0 && fv.navPos().size() < level) {
0183       level = 0;
0184       exitLoop = true;
0185     }
0186     if (dd4hep::dd::noNamespace(fv.name()) == "BarrelTimingLayer") {
0187       level = fv.navPos().size();
0188     }
0189
0190     // Test only the desired subdetector
0191
0192     if (exitLoop && isBarrel) {
0193       break;
0194     }
0195
0196     // Actions for MTD volumes: search for sensitive detectors
0197
0198     bool isSens = false;
0199
0200     for (auto const& t : specs) {
0201       for (auto const& it : t.second->paths) {
0202         if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(fv.name()), dd4hep::dd::realTopName(it))) {
0203           isSens = true;
0204           break;
0205         }
0206       }
0207     }
0208
0209     if (isSens && isBarrel) {
0210       std::stringstream spix;
0211       spix << print_path() << "\n\n";
0212
0213       BTLDetId theId(btlNS_.getUnitID(thisN_));
0214
0215       DetId geoId = theId.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(pTP.product()->getMTDTopologyMode()));
0216       const MTDGeomDet* thedet = pDG.product()->idToDet(geoId);
0217
0218       if (thedet == nullptr) {
0219         throw cms::Exception("BTLDeviceSim") << "GeographicalID: " << std::hex << geoId.rawId() << " (" << theId.rawId()
0220                                              << ") is invalid!" << std::dec << std::endl;
0221       }
0222       const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(thedet->topology());
0223       const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0224
0225       int origRow = theId.row(topo.nrows());
0226       int origCol = theId.column(topo.nrows());
0227       spix << "rawId= " << theId.rawId() << " geoId= " << geoId.rawId() << " side/rod= " << theId.mtdSide() << " / "
0228            << theId.mtdRR() << " type/RU= " << theId.modType() << " / " << theId.runit()
0229            << " module/geomodule= " << theId.module() << " / " << static_cast<BTLDetId>(geoId).module()
0230            << " crys= " << theId.crystal() << " BTLDetId row/col= " << origRow << " / " << origCol;
0231       spix << "\n";
0232
0233       //
0234       // Test of positions for sensitive detectors
0235       //
0236
0237       auto fround = [&](double in) {
0238         std::stringstream ss;
0239         ss << std::fixed << std::setw(14) << roundIfNear0(in);
0240         return ss.str();
0241       };
0242
0243       if (!dd4hep::isA<dd4hep::Box>(fv.solid())) {
0244         throw cms::Exception("DD4hep_TestBTLPixelTopology") << "MTD sensitive element not a DDBox";
0245         break;
0246       }
0247       dd4hep::Box mySens(fv.solid());
0248
0249       DD3Vector zeroLocal(0., 0., 0.);
0250       DD3Vector cn1Local(mySens.x(), mySens.y(), mySens.z());
0251       DD3Vector cn2Local(-mySens.x(), -mySens.y(), -mySens.z());
0252       DD3Vector zeroGlobal = (fv.rotation())(zeroLocal) + fv.translation();
0253       DD3Vector cn1Global = (fv.rotation())(cn1Local) + fv.translation();
0254       DD3Vector cn2Global = (fv.rotation())(cn2Local) + fv.translation();
0255
0256       const size_t nTest(3);
0257       std::array<Local3DPoint, nTest> refLocalPoints{{Local3DPoint(zeroLocal.x(), zeroLocal.y(), zeroLocal.z()),
0258                                                       Local3DPoint(cn1Local.x(), cn1Local.y(), cn1Local.z()),
0259                                                       Local3DPoint(cn2Local.x(), cn2Local.y(), cn2Local.z())}};
0260       std::array<DD3Vector, nTest> refGlobalPoints{{zeroGlobal, cn1Global, cn2Global}};
0261
0262       for (size_t iloop = 0; iloop < nTest; iloop++) {
0263         // translate from crystal-local coordinates to module-local coordinates to get the row and column
0264
0265         Local3DPoint cmRefLocal(convertMmToCm(refLocalPoints[iloop].x() / dd4hep::mm),
0266                                 convertMmToCm(refLocalPoints[iloop].y() / dd4hep::mm),
0267                                 convertMmToCm(refLocalPoints[iloop].z() / dd4hep::mm));
0268         Local3DPoint modLocal = topo.pixelToModuleLocalPoint(cmRefLocal, origRow, origCol);
0269         const auto& thepixel = topo.pixel(modLocal);
0270         uint8_t recoRow(thepixel.first), recoCol(thepixel.second);
0271
0272         if (origRow != recoRow || origCol != recoCol) {
0273           std::stringstream warnmsg;
0274           warnmsg << "DIFFERENCE row/col, orig= " << origRow << " " << origCol
0275                   << " reco= " << static_cast<uint32_t>(recoRow) << " " << static_cast<uint32_t>(recoCol) << "\n";
0276           spix << warnmsg.str();
0277           sunitt << warnmsg.str();
0278           recoRow = origRow;
0279           recoCol = origCol;
0280         }
0281
0282         // reconstructed global position from reco geometry and rectangluar MTD topology
0283
0284         const auto& modGlobal = thedet->toGlobal(modLocal);
0285
0286         const double deltax = convertCmToMm(modGlobal.x()) - (refGlobalPoints[iloop].x() / dd4hep::mm);
0287         const double deltay = convertCmToMm(modGlobal.y()) - (refGlobalPoints[iloop].y() / dd4hep::mm);
0288         const double deltaz = convertCmToMm(modGlobal.z()) - (refGlobalPoints[iloop].z() / dd4hep::mm);
0289
0290         spix << "Ref#" << iloop << " local= " << fround(refLocalPoints[iloop].x() / dd4hep::mm)
0291              << fround(refLocalPoints[iloop].y() / dd4hep::mm) << fround(refLocalPoints[iloop].z() / dd4hep::mm)
0292              << " Orig global= " << fround(refGlobalPoints[iloop].x() / dd4hep::mm)
0293              << fround(refGlobalPoints[iloop].y() / dd4hep::mm) << fround(refGlobalPoints[iloop].z() / dd4hep::mm)
0294              << " Reco global= " << fround(convertCmToMm(modGlobal.x())) << fround(convertCmToMm(modGlobal.y()))
0295              << fround(convertCmToMm(modGlobal.z())) << " Delta= " << fround(deltax) << fround(deltay) << fround(deltaz)
0296              << "\n";
0297         if (std::abs(deltax) > tolerance || std::abs(deltay) > tolerance || std::abs(deltaz) > tolerance) {
0298           std::stringstream warnmsg;
0299           warnmsg << "DIFFERENCE detId/ref# " << theId.rawId() << " " << iloop << " dx/dy/dz= " << fround(deltax)
0300                   << fround(deltay) << fround(deltaz) << "\n";
0301           spix << warnmsg.str();
0302           sunitt << warnmsg.str();
0303         }
0304       }
0305
0306       spix << "\n";
0307       edm::LogVerbatim("DD4hep_TestBTLPixelTopology") << spix.str();
0308     }
0309   } while (fv.next(0));
0310
0311   if (!sunitt.str().empty()) {
0312     edm::LogVerbatim("MTDUnitTest") << sunitt.str();
0313   }
0314 }
0315
0316 void DD4hep_TestBTLPixelTopology::theBaseNumber(cms::DDFilteredView& fv) {
0317   thisN_.reset();
0318   thisN_.setSize(fv.navPos().size());
0319
0320   for (uint ii = 0; ii < fv.navPos().size(); ii++) {
0321     std::string_view name((fv.geoHistory()[ii])->GetName());
0322     size_t ipos = name.rfind('_');
0323     thisN_.addLevel(name.substr(0, ipos), fv.copyNos()[ii]);
0324 #ifdef EDM_ML_DEBUG
0325     edm::LogVerbatim("DD4hep_TestBTLPixelTopology") << name.substr(0, ipos) << " " << fv.copyNos()[ii];
0326 #endif
0327   }
0328 }
0329
0330 DEFINE_FWK_MODULE(DD4hep_TestBTLPixelTopology);