MTDCommonData/test/DD4hep_TestMTDIdealGeometry.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 #include "FWCore/Framework/interface/ConsumesCollector.h"
0017
0018 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0019 #include "Geometry/Records/interface/DDSpecParRegistryRcd.h"
0020
0021 #include "DetectorDescription/DDCMS/interface/DDDetector.h"
0022 #include "DetectorDescription/DDCMS/interface/DDSolidShapes.h"
0023 #include "DetectorDescription/DDCMS/interface/DDFilteredView.h"
0024 #include "DetectorDescription/DDCMS/interface/DDSpecParRegistry.h"
0025
0026 #include "Geometry/MTDCommonData/interface/MTDBaseNumber.h"
0027 #include "Geometry/MTDCommonData/interface/BTLNumberingScheme.h"
0028 #include "Geometry/MTDCommonData/interface/ETLNumberingScheme.h"
0029 #include "Geometry/MTDCommonData/interface/BTLElectronicsMapping.h"
0030 #include "Geometry/MTDCommonData/interface/MTDTopologyMode.h"
0031
0032 #include "Geometry/MTDGeometryBuilder/interface/MTDTopology.h"
0033 #include "Geometry/Records/interface/MTDTopologyRcd.h"
0034
0035 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0036 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0037
0038 #include "DataFormats/Math/interface/angle_units.h"
0039 #include "DataFormats/Math/interface/Rounding.h"
0040 #include <DD4hep/DD4hepUnits.h>
0041
0042 using namespace cms;
0043
0044 class DD4hep_TestMTDIdealGeometry : public edm::one::EDAnalyzer<> {
0045 public:
0046   explicit DD4hep_TestMTDIdealGeometry(const edm::ParameterSet&);
0047   ~DD4hep_TestMTDIdealGeometry() override = default;
0048
0049   void beginJob() override {}
0050   void analyze(edm::Event const&, edm::EventSetup const&) override;
0051   void endJob() override {}
0052
0053   void theBaseNumber(cms::DDFilteredView& fv);
0054
0055 private:
0056   const edm::ESInputTag tag_;
0057   std::string ddTopNodeName_;
0058
0059   MTDBaseNumber thisN_;
0060   BTLNumberingScheme btlNS_;
0061   ETLNumberingScheme etlNS_;
0062
0063   edm::ESGetToken<MTDTopology, MTDTopologyRcd> mtdtopoToken_;
0064   edm::ESGetToken<DDDetector, IdealGeometryRecord> dddetToken_;
0065   edm::ESGetToken<DDSpecParRegistry, DDSpecParRegistryRcd> dspecToken_;
0066 };
0067
0068 using DD3Vector = ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double>>;
0069 using angle_units::operators::convertRadToDeg;
0070 using cms_rounding::roundIfNear0;
0071
0072 DD4hep_TestMTDIdealGeometry::DD4hep_TestMTDIdealGeometry(const edm::ParameterSet& iConfig)
0073     : tag_(iConfig.getParameter<edm::ESInputTag>("DDDetector")),
0074       ddTopNodeName_(iConfig.getUntrackedParameter<std::string>("ddTopNodeName", "BarrelTimingLayer")),
0075       thisN_(),
0076       btlNS_(),
0077       etlNS_() {
0078   mtdtopoToken_ = esConsumes<MTDTopology, MTDTopologyRcd>();
0079   dddetToken_ = esConsumes<DDDetector, IdealGeometryRecord>(tag_);
0080   dspecToken_ = esConsumes<DDSpecParRegistry, DDSpecParRegistryRcd>(tag_);
0081 }
0082
0083 void DD4hep_TestMTDIdealGeometry::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0084   auto pDD = iSetup.getTransientHandle(dddetToken_);
0085
0086   auto pSP = iSetup.getTransientHandle(dspecToken_);
0087
0088   auto topologyHandle = iSetup.getTransientHandle(mtdtopoToken_);
0089   const MTDTopology* topology = topologyHandle.product();
0090   auto btlCrysLayout = MTDTopologyMode::crysLayoutFromTopoMode(topology->getMTDTopologyMode());
0091
0092   if (ddTopNodeName_ != "BarrelTimingLayer" && ddTopNodeName_ != "EndcapTimingLayer") {
0093     edm::LogWarning("DD4hep_TestMTDIdealGeometry") << ddTopNodeName_ << "Not valid top MTD volume";
0094     return;
0095   }
0096
0097   if (!pDD.isValid()) {
0098     edm::LogError("DD4hep_TestMTDIdealGeometry") << "ESTransientHandle<DDCompactView> pDD is not valid!";
0099     return;
0100   }
0101   if (pDD.description()) {
0102     edm::LogInfo("DD4hep_TestMTDIdealGeometry") << pDD.description()->type_ << " label: " << pDD.description()->label_;
0103   } else {
0104     edm::LogWarning("DD4hep_TestMTDIdealGeometry") << "NO label found pDD.description() returned false.";
0105   }
0106
0107   if (!pSP.isValid()) {
0108     edm::LogError("DD4hep_TestMTDIdealGeometry") << "ESTransientHandle<DDSpecParRegistry> pSP is not valid!";
0109     return;
0110   }
0111
0112   DDFilteredView fv(pDD.product(), pDD.product()->description()->worldVolume());
0113   fv.next(0);
0114   edm::LogInfo("DD4hep_TestMTDIdealGeometry") << fv.name();
0115
0116   DDSpecParRefs specs;
0117   std::string attribute("ReadOutName"), name;
0118   if (ddTopNodeName_ == "BarrelTimingLayer") {
0119     name = "FastTimerHitsBarrel";
0120   } else if (ddTopNodeName_ == "EndcapTimingLayer") {
0121     name = "FastTimerHitsEndcap";
0122   }
0123   if (name.empty()) {
0124     edm::LogError("DD4hep_TestMTDIdealGeometry") << "No sensitive detector provided, abort";
0125     return;
0126   }
0127   pSP.product()->filter(specs, attribute, name);
0128
0129   edm::LogVerbatim("Geometry").log([&specs](auto& log) {
0130     log << "Filtered DD SpecPar Registry size: " << specs.size() << "\n";
0131     for (const auto& t : specs) {
0132       log << "\nSpecPar " << t.first << ":\nRegExps { ";
0133       for (const auto& ki : t.second->paths)
0134         log << ki << " ";
0135       log << "};\n ";
0136       for (const auto& kl : t.second->spars) {
0137         log << kl.first << " = ";
0138         for (const auto& kil : kl.second) {
0139           log << kil << " ";
0140         }
0141         log << "\n ";
0142       }
0143     }
0144   });
0145
0146   bool write = false;
0147   bool isBarrel = true;
0148   bool exitLoop = false;
0149   uint32_t level(0);
0150   uint32_t count(0);
0151
0152   do {
0153     if (dd4hep::dd::noNamespace(fv.name()) == "BarrelTimingLayer") {
0154       isBarrel = true;
0155       edm::LogInfo("DD4hep_TestMTDIdealGeometry") << "isBarrel = " << isBarrel;
0156       if (static_cast<int>(btlCrysLayout) < static_cast<int>(BTLDetId::CrysLayout::v4)) {
0157         edm::LogInfo("DD4hep_TestMTDIdealGeometry")
0158             << "BTL electronics mapping not available for BTL crystal layout " << static_cast<int>(btlCrysLayout)
0159             << ", use layout 7 (v4) or later!" << std::endl;
0160       }
0161     } else if (dd4hep::dd::noNamespace(fv.name()) == "EndcapTimingLayer") {
0162       isBarrel = false;
0163       edm::LogInfo("DD4hep_TestMTDIdealGeometry") << "isBarrel = " << isBarrel;
0164     }
0165
0166     if (level > 0 && fv.navPos().size() < level) {
0167       level = 0;
0168       write = false;
0169       if (isBarrel) {
0170         exitLoop = true;
0171       } else if (!isBarrel && count == 2) {
0172         exitLoop = true;
0173       }
0174     }
0175     if (dd4hep::dd::noNamespace(fv.name()) == ddTopNodeName_) {
0176       write = true;
0177       level = fv.navPos().size();
0178       count += 1;
0179     }
0180
0181 #ifdef EDM_ML_DEBUG
0182     edm::LogVerbatim("DD4hep_TestMTDIdealGeometry")
0183         << "level= " << level << " isBarrel= " << isBarrel << " exitLoop= " << exitLoop << " count= " << count << " "
0184         << fv.path();
0185 #endif
0186
0187     // Test only the desired subdetector
0188
0189     if (exitLoop && isBarrel) {
0190       break;
0191     }
0192
0193     // Actions for MTD volumes: searchg for sensitive detectors
0194
0195     if (write) {
0196       std::stringstream ss;
0197
0198       theBaseNumber(fv);
0199
0200       auto print_path = [&]() {
0201         ss << " - OCMS[0]/";
0202         for (int ii = thisN_.getLevels() - 1; ii-- > 0;) {
0203           ss << thisN_.getLevelName(ii);
0204           ss << "[";
0205           ss << thisN_.getCopyNumber(ii);
0206           ss << "]/";
0207         }
0208       };
0209
0210       print_path();
0211
0212       edm::LogInfo("DD4hep_TestMTDPath") << ss.str();
0213
0214       bool isSens = false;
0215
0216       for (auto const& t : specs) {
0217         for (auto const& it : t.second->paths) {
0218           if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(fv.name()), dd4hep::dd::realTopName(it))) {
0219             isSens = true;
0220             break;
0221           }
0222         }
0223       }
0224
0225       if (isSens) {
0226         //
0227         // Test of numbering scheme for sensitive detectors
0228         //
0229
0230         std::stringstream sunitt;
0231         std::stringstream snum;
0232
0233         if (isBarrel) {
0234           BTLDetId theId(btlNS_.getUnitID(thisN_));
0235           sunitt << theId.rawId();
0236           snum << theId;
0237
0238           if (static_cast<int>(btlCrysLayout) >= static_cast<int>(BTLDetId::CrysLayout::v4)) {
0239             BTLElectronicsMapping btlEM = BTLElectronicsMapping(btlCrysLayout);
0240             snum << "\n";
0241             snum << "----------------------------------------------------------------------------" << std::endl;
0242             snum << " CCBoard: " << btlEM.CCBoard(theId) << " FEBoard: " << btlEM.FEBoard(theId)
0243                  << " TOFHIRASIC: " << btlEM.TOFHIRASIC(theId) << "\n SiPMCh   minus: " << btlEM.SiPMCh(theId, 0)
0244                  << " plus: " << btlEM.SiPMCh(theId, 1) << "\n TOFHIRCh minus: " << btlEM.TOFHIRCh(theId, 0)
0245                  << " plus: " << btlEM.TOFHIRCh(theId, 1) << "\n";
0246             snum << "----------------------------------------------------------------------------" << std::endl;
0247           }
0248         } else {
0249           ETLDetId theId(etlNS_.getUnitID(thisN_));
0250           sunitt << theId.rawId();
0251           snum << theId;
0252         }
0253         edm::LogInfo("DD4hep_TestMTDNumbering") << snum.str();
0254
0255         //
0256         // Test of positions for sensitive detectors
0257         //
0258
0259         std::stringstream spos;
0260
0261         auto fround = [&](double in) {
0262           std::stringstream ss;
0263           ss << std::fixed << std::setw(14) << roundIfNear0(in);
0264           return ss.str();
0265         };
0266
0267         if (!dd4hep::isA<dd4hep::Box>(fv.solid())) {
0268           throw cms::Exception("TestMTDIdealGeometry") << "MTD sensitive element not a DDBox";
0269           break;
0270         }
0271         dd4hep::Box mySens(fv.solid());
0272         spos << "Solid shape name: " << DDSolidShapesName::name(fv.legacyShape(fv.shape())) << "\n";
0273         spos << "Box dimensions: " << fround(mySens.x() / dd4hep::mm) << " " << fround(mySens.y() / dd4hep::mm) << " "
0274              << fround(mySens.z() / dd4hep::mm) << "\n";
0275
0276         DD3Vector x, y, z;
0277         fv.rotation().GetComponents(x, y, z);
0278         spos << "Translation vector components: " << fround(fv.translation().x() / dd4hep::mm) << " "
0279              << fround(fv.translation().y() / dd4hep::mm) << " " << fround(fv.translation().z() / dd4hep::mm) << " "
0280              << "\n";
0281         spos << "Rotation matrix components: " << fround(x.X()) << " " << fround(x.Y()) << " " << fround(x.Z()) << " "
0282              << fround(y.X()) << " " << fround(y.Y()) << " " << fround(y.Z()) << " " << fround(z.X()) << " "
0283              << fround(z.Y()) << " " << fround(z.Z()) << "\n";
0284
0285         DD3Vector zeroLocal(0., 0., 0.);
0286         DD3Vector cn1Local(mySens.x(), mySens.y(), mySens.z());
0287         double distLocal = cn1Local.R();
0288         DD3Vector zeroGlobal = (fv.rotation())(zeroLocal) + fv.translation();
0289         DD3Vector cn1Global = (fv.rotation())(cn1Local) + fv.translation();
0290         double distGlobal =
0291             std::sqrt(std::pow(zeroGlobal.X() - cn1Global.X(), 2) + std::pow(zeroGlobal.Y() - cn1Global.Y(), 2) +
0292                       std::pow(zeroGlobal.Z() - cn1Global.Z(), 2));
0293
0294         spos << "Center global   = " << fround(zeroGlobal.X() / dd4hep::mm) << fround(zeroGlobal.Y() / dd4hep::mm)
0295              << fround(zeroGlobal.Z() / dd4hep::mm) << " r = " << fround(zeroGlobal.Rho() / dd4hep::mm)
0296              << " phi = " << fround(convertRadToDeg(zeroGlobal.Phi())) << "\n";
0297
0298         spos << "Corner 1 local  = " << fround(cn1Local.X() / dd4hep::mm) << fround(cn1Local.Y() / dd4hep::mm)
0299              << fround(cn1Local.Z() / dd4hep::mm) << " DeltaR = " << fround(distLocal / dd4hep::mm) << "\n";
0300
0301         spos << "Corner 1 global = " << fround(cn1Global.X() / dd4hep::mm) << fround(cn1Global.Y() / dd4hep::mm)
0302              << fround(cn1Global.Z() / dd4hep::mm) << " DeltaR = " << fround(distGlobal / dd4hep::mm) << "\n";
0303
0304         spos << "\n";
0305         if (std::fabs(distGlobal - distLocal) / dd4hep::mm > 1.e-6) {
0306           spos << "DIFFERENCE IN DISTANCE \n";
0307         }
0308         sunitt << fround(zeroGlobal.X() / dd4hep::mm) << fround(zeroGlobal.Y() / dd4hep::mm)
0309                << fround(zeroGlobal.Z() / dd4hep::mm) << fround(cn1Global.X() / dd4hep::mm)
0310                << fround(cn1Global.Y() / dd4hep::mm) << fround(cn1Global.Z() / dd4hep::mm);
0311         edm::LogInfo("DD4hep_TestMTDPosition") << spos.str();
0312
0313         edm::LogVerbatim("MTDUnitTest") << sunitt.str();
0314       }
0315     }
0316   } while (fv.next(0) && !(exitLoop == 1 && count == 2));
0317 }
0318
0319 void DD4hep_TestMTDIdealGeometry::theBaseNumber(cms::DDFilteredView& fv) {
0320   thisN_.reset();
0321   thisN_.setSize(fv.navPos().size());
0322
0323   for (uint ii = 0; ii < fv.navPos().size(); ii++) {
0324     std::string_view name((fv.geoHistory()[ii])->GetName());
0325     size_t ipos = name.rfind('_');
0326     thisN_.addLevel(name.substr(0, ipos), fv.copyNos()[ii]);
0327 #ifdef EDM_ML_DEBUG
0328     edm::LogVerbatim("DD4hep_TestMTDIdealGeometry") << ii << " " << name.substr(0, ipos) << " " << fv.copyNos()[ii];
0329 #endif
0330   }
0331 }
0332
0333 DEFINE_FWK_MODULE(DD4hep_TestMTDIdealGeometry);