MTDCommonData/test/TestMTDIdealGeometry.cc

0001 #include <iostream>
0002 #include <fstream>
0003 #include <string>
0004 #include <utility>
0005 #include <vector>
0006
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/Framework/interface/ESTransientHandle.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/Utilities/interface/Exception.h"
0015
0016 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0017
0018 #include "DetectorDescription/Core/interface/DDFilteredView.h"
0019 #include "DetectorDescription/Core/interface/DDSolid.h"
0020 #include "DetectorDescription/Core/interface/DDSolidShapes.h"
0021
0022 #include "Geometry/MTDCommonData/interface/MTDBaseNumber.h"
0023 #include "Geometry/MTDCommonData/interface/BTLNumberingScheme.h"
0024 #include "Geometry/MTDCommonData/interface/ETLNumberingScheme.h"
0025
0026 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0027 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0028
0029 #include "DataFormats/Math/interface/angle_units.h"
0030 #include "DataFormats/Math/interface/Rounding.h"
0031
0032 //#define EDM_ML_DEBUG
0033
0034 class TestMTDIdealGeometry : public edm::one::EDAnalyzer<> {
0035 public:
0036   explicit TestMTDIdealGeometry(const edm::ParameterSet&);
0037   ~TestMTDIdealGeometry() override;
0038
0039   void beginJob() override {}
0040   void analyze(edm::Event const&, edm::EventSetup const&) override;
0041   void endJob() override {}
0042
0043   void theBaseNumber(const DDGeoHistory& gh);
0044
0045 private:
0046   int nNodes_;
0047   std::string ddTopNodeName_;
0048
0049   MTDBaseNumber thisN_;
0050   BTLNumberingScheme btlNS_;
0051   ETLNumberingScheme etlNS_;
0052
0053   edm::ESGetToken<DDCompactView, IdealGeometryRecord> cpvToken_;
0054 };
0055
0056 TestMTDIdealGeometry::TestMTDIdealGeometry(const edm::ParameterSet& iConfig)
0057     : ddTopNodeName_(iConfig.getUntrackedParameter<std::string>("ddTopNodeName", "BarrelTimingLayer")),
0058       thisN_(),
0059       btlNS_(),
0060       etlNS_() {
0061   cpvToken_ = esConsumes<DDCompactView, IdealGeometryRecord>();
0062 }
0063
0064 TestMTDIdealGeometry::~TestMTDIdealGeometry() {}
0065
0066 using angle_units::operators::convertRadToDeg;
0067 using cms_rounding::roundIfNear0;
0068
0069 void TestMTDIdealGeometry::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0070   using namespace edm;
0071
0072   if (ddTopNodeName_ != "BarrelTimingLayer" && ddTopNodeName_ != "EndcapTimingLayer") {
0073     edm::LogWarning("TestMTDIdealGeometry") << ddTopNodeName_ << "Not valid top MTD volume";
0074     return;
0075   }
0076
0077   auto pDD = iSetup.getTransientHandle(cpvToken_);
0078
0079   if (!pDD.isValid()) {
0080     edm::LogError("TestMTDIdealGeometry") << "ESTransientHandle<DDCompactView> pDD is not valid!";
0081     return;
0082   }
0083   if (pDD.description()) {
0084     edm::LogInfo("TestMTDIdealGeometry") << pDD.description()->type_ << " label: " << pDD.description()->label_;
0085   } else {
0086     edm::LogWarning("TestMTDIdealGeometry") << "NO label found pDD.description() returned false.";
0087   }
0088
0089   DDPassAllFilter filter;
0090   DDFilteredView fv(*pDD, filter);
0091
0092   edm::LogInfo("TestMTDIdealGeometry") << "Top Most LogicalPart = " << fv.logicalPart();
0093
0094   using nav_type = DDFilteredView::nav_type;
0095   using id_type = std::map<nav_type, int>;
0096   id_type idMap;
0097   int id = 0;
0098
0099   bool write = false;
0100   bool isBarrel = true;
0101   size_t limit = 0;
0102
0103   do {
0104     nav_type pos = fv.navPos();
0105     idMap[pos] = id;
0106
0107     size_t num = fv.geoHistory().size();
0108
0109     if (num <= limit) {
0110       write = false;
0111     }
0112     if (fv.geoHistory()[num - 1].logicalPart().name() == "btl:BarrelTimingLayer") {
0113       isBarrel = true;
0114       limit = num;
0115       write = true;
0116 #ifdef EDM_ML_DEBUG
0117       edm::LogInfo("TestMTDIdealGeometry") << "isBarrel = " << isBarrel;
0118 #endif
0119     } else if (fv.geoHistory()[num - 1].logicalPart().name() == "etl:EndcapTimingLayer") {
0120       isBarrel = false;
0121       limit = num;
0122       write = true;
0123 #ifdef EDM_ML_DEBUG
0124       edm::LogInfo("TestMTDIdealGeometry") << "isBarrel = " << isBarrel;
0125 #endif
0126     }
0127
0128     // Actions for MTD volumes: searchg for sensitive detectors
0129
0130     std::stringstream ss;
0131     auto print_path = [&]() {
0132       ss << " - OCMS[0]/";
0133       for (uint i = 1; i < fv.geoHistory().size(); i++) {
0134         ss << fv.geoHistory()[i].logicalPart().name().fullname();
0135         ss << "[";
0136         ss << std::to_string(fv.geoHistory()[i].copyno());
0137         ss << "]/";
0138       }
0139     };
0140
0141     if (write && fv.geoHistory()[limit - 1].logicalPart().name().name() == ddTopNodeName_) {
0142       print_path();
0143       edm::LogInfo("TestMTDPath") << ss.str();
0144
0145       bool isSens = false;
0146
0147       if (!fv.geoHistory()[num - 1].logicalPart().specifics().empty()) {
0148         for (auto vec : fv.geoHistory()[num - 1].logicalPart().specifics()) {
0149           for (const auto& elem : *vec) {
0150             if (elem.second.name() == "SensitiveDetector") {
0151               isSens = true;
0152               break;
0153             }
0154           }
0155         }
0156       }
0157
0158       if (isSens) {
0159         theBaseNumber(fv.geoHistory());
0160
0161         //
0162         // Test of numbering scheme for sensitive detectors
0163         //
0164
0165         std::stringstream sunitt;
0166         std::stringstream snum;
0167
0168         if (isBarrel) {
0169           BTLDetId theId(btlNS_.getUnitID(thisN_));
0170           sunitt << theId.rawId();
0171           snum << theId;
0172           snum << "\n";
0173         } else {
0174           ETLDetId theId(etlNS_.getUnitID(thisN_));
0175           sunitt << theId.rawId();
0176           snum << theId;
0177 #ifdef EDM_ML_DEBUG
0178           edm::LogInfo("TestMTDNumbering")
0179               << " ETLDetId = " << theId << "\n geographicalId = " << theId.geographicalId();
0180 #endif
0181         }
0182         edm::LogInfo("TestMTDNumbering") << snum.str();
0183
0184         //
0185         // Test of positions for sensitive detectors
0186         //
0187
0188         std::stringstream spos;
0189
0190         auto fround = [&](double in) {
0191           std::stringstream ss;
0192           ss << std::fixed << std::setw(14) << roundIfNear0(in);
0193           return ss.str();
0194         };
0195
0196         DDBox mySens = fv.geoHistory()[num - 1].logicalPart().solid();
0197         spos << "Solid shape name: " << DDSolidShapesName::name(mySens.shape()) << "\n";
0198         if (static_cast<int>(mySens.shape()) != 1) {
0199           throw cms::Exception("TestMTDPosition") << "MTD sensitive element not a DDBox";
0200           break;
0201         }
0202         spos << "Box dimensions: " << fround(mySens.halfX()) << " " << fround(mySens.halfY()) << " "
0203              << fround(mySens.halfZ()) << "\n";
0204
0205         DD3Vector x, y, z;
0206         fv.rotation().GetComponents(x, y, z);
0207         spos << "Translation vector components: " << fround(fv.translation().x()) << " " << fround(fv.translation().y())
0208              << " " << fround(fv.translation().z()) << " "
0209              << "\n";
0210         spos << "Rotation matrix components: " << fround(x.X()) << " " << fround(x.Y()) << " " << fround(x.Z()) << " "
0211              << fround(y.X()) << " " << fround(y.Y()) << " " << fround(y.Z()) << " " << fround(z.X()) << " "
0212              << fround(z.Y()) << " " << fround(z.Z()) << "\n";
0213
0214         DD3Vector zeroLocal(0., 0., 0.);
0215         DD3Vector cn1Local(mySens.halfX(), mySens.halfY(), mySens.halfZ());
0216         double distLocal = cn1Local.R();
0217         DD3Vector zeroGlobal = (fv.rotation())(zeroLocal) + fv.translation();
0218         DD3Vector cn1Global = (fv.rotation())(cn1Local) + fv.translation();
0219         ;
0220         double distGlobal =
0221             std::sqrt(std::pow(zeroGlobal.X() - cn1Global.X(), 2) + std::pow(zeroGlobal.Y() - cn1Global.Y(), 2) +
0222                       std::pow(zeroGlobal.Z() - cn1Global.Z(), 2));
0223
0224         spos << "Center global   = " << fround(zeroGlobal.X()) << fround(zeroGlobal.Y()) << fround(zeroGlobal.Z())
0225              << " r = " << fround(zeroGlobal.Rho()) << " phi = " << fround(convertRadToDeg(zeroGlobal.Phi())) << "\n";
0226
0227         spos << "Corner 1 local  = " << fround(cn1Local.X()) << fround(cn1Local.Y()) << fround(cn1Local.Z())
0228              << " DeltaR = " << fround(distLocal) << "\n";
0229
0230         spos << "Corner 1 global = " << fround(cn1Global.X()) << fround(cn1Global.Y()) << fround(cn1Global.Z())
0231              << " DeltaR = " << fround(distGlobal) << "\n";
0232
0233         spos << "\n";
0234         if (std::fabs(distGlobal - distLocal) > 1.e-6) {
0235           spos << "DIFFERENCE IN DISTANCE \n";
0236         }
0237         sunitt << fround(zeroGlobal.X()) << fround(zeroGlobal.Y()) << fround(zeroGlobal.Z()) << fround(cn1Global.X())
0238                << fround(cn1Global.Y()) << fround(cn1Global.Z());
0239         edm::LogInfo("TestMTDPosition") << spos.str();
0240
0241         edm::LogVerbatim("MTDUnitTest") << sunitt.str();
0242       }
0243     }
0244     ++id;
0245   } while (fv.next());
0246 }
0247
0248 void TestMTDIdealGeometry::theBaseNumber(const DDGeoHistory& gh) {
0249   thisN_.reset();
0250   thisN_.setSize(gh.size());
0251
0252   for (uint i = gh.size(); i-- > 0;) {
0253     thisN_.addLevel(gh[i].logicalPart().name().name(), gh[i].copyno());
0254 #ifdef EDM_ML_DEBUG
0255     edm::LogInfo("TestMTDIdealGeometry") << gh[i].logicalPart().name().name() << " " << gh[i].copyno();
0256 #endif
0257   }
0258 }
0259
0260 DEFINE_FWK_MODULE(TestMTDIdealGeometry);