File indexing completed on 2022-10-25 03:22:15
0001
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
0191
0192 if (exitLoop && isBarrel) {
0193 break;
0194 }
0195
0196
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
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
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
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);