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