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