MTDGeometryBuilder/test/DD4hep_TestBTLPixelTopology.cc

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