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