File indexing completed on 2021-11-23 04:29:05
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 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
0177
0178 if (exitLoop && isBarrel) {
0179 break;
0180 }
0181
0182
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
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
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);