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