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