Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-01-18 03:42: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/MTDTopologyMode.h"
0026 #include "Geometry/MTDCommonData/interface/MTDBaseNumber.h"
0027 #include "Geometry/MTDCommonData/interface/BTLNumberingScheme.h"
0028 #include "Geometry/MTDCommonData/interface/ETLNumberingScheme.h"
0029 
0030 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0031 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0032 #include "DataFormats/GeometrySurface/interface/MediumProperties.h"
0033 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0034 
0035 #include "Geometry/Records/interface/MTDTopologyRcd.h"
0036 #include "Geometry/MTDGeometryBuilder/interface/MTDTopology.h"
0037 #include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h"
0038 #include "Geometry/MTDGeometryBuilder/interface/MTDGeomDetUnit.h"
0039 #include "Geometry/Records/interface/MTDDigiGeometryRecord.h"
0040 #include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h"
0041 #include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
0042 
0043 #include "DataFormats/Math/interface/GeantUnits.h"
0044 #include "DataFormats/Math/interface/Rounding.h"
0045 #include "DataFormats/Math/interface/angle_units.h"
0046 #include <DD4hep/DD4hepUnits.h>
0047 
0048 using namespace cms;
0049 using namespace geant_units::operators;
0050 using namespace cms_rounding;
0051 
0052 class DD4hep_TestPixelTopology : public edm::one::EDAnalyzer<> {
0053 public:
0054   explicit DD4hep_TestPixelTopology(const edm::ParameterSet&);
0055   ~DD4hep_TestPixelTopology() = default;
0056 
0057   void beginJob() override {}
0058   void analyze(edm::Event const&, edm::EventSetup const&) override;
0059   void endJob() override {}
0060 
0061   void theBaseNumber(cms::DDFilteredView& fv);
0062 
0063 private:
0064   inline std::string fround(const double in, const size_t prec) const {
0065     std::stringstream ss;
0066     ss << std::setprecision(prec) << std::fixed << std::setw(14) << roundIfNear0(in);
0067     return ss.str();
0068   }
0069 
0070   inline std::string fvecround(const auto& vecin, const size_t prec) const {
0071     std::stringstream ss;
0072     ss << std::setprecision(prec) << std::fixed << std::setw(14) << roundVecIfNear0(vecin);
0073     return ss.str();
0074   }
0075 
0076   void analyseRectangle(const GeomDetUnit& det);
0077   void checkRotation(const GeomDetUnit& det);
0078 
0079   const edm::ESInputTag tag_;
0080   std::string ddTopNodeName_;
0081 
0082   MTDBaseNumber thisN_;
0083   BTLNumberingScheme btlNS_;
0084   ETLNumberingScheme etlNS_;
0085 
0086   edm::ESGetToken<DDDetector, IdealGeometryRecord> dddetToken_;
0087   edm::ESGetToken<DDSpecParRegistry, DDSpecParRegistryRcd> dspecToken_;
0088   edm::ESGetToken<MTDTopology, MTDTopologyRcd> mtdtopoToken_;
0089   edm::ESGetToken<MTDGeometry, MTDDigiGeometryRecord> mtdgeoToken_;
0090 
0091   std::stringstream sunitt_;
0092   constexpr static double tolerance{0.5e-3_mm};
0093 };
0094 
0095 using DD3Vector = ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double>>;
0096 using angle_units::operators::convertRadToDeg;
0097 
0098 DD4hep_TestPixelTopology::DD4hep_TestPixelTopology(const edm::ParameterSet& iConfig)
0099     : tag_(iConfig.getParameter<edm::ESInputTag>("DDDetector")),
0100       ddTopNodeName_(iConfig.getUntrackedParameter<std::string>("ddTopNodeName", "BarrelTimingLayer")),
0101       thisN_(),
0102       btlNS_(),
0103       etlNS_() {
0104   dddetToken_ = esConsumes<DDDetector, IdealGeometryRecord>(tag_);
0105   dspecToken_ = esConsumes<DDSpecParRegistry, DDSpecParRegistryRcd>(tag_);
0106   mtdtopoToken_ = esConsumes<MTDTopology, MTDTopologyRcd>(tag_);
0107   mtdgeoToken_ = esConsumes<MTDGeometry, MTDDigiGeometryRecord>(tag_);
0108 }
0109 
0110 void DD4hep_TestPixelTopology::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0111   auto pDD = iSetup.getTransientHandle(dddetToken_);
0112 
0113   auto pSP = iSetup.getTransientHandle(dspecToken_);
0114 
0115   if (ddTopNodeName_ != "BarrelTimingLayer" && ddTopNodeName_ != "EndcapTimingLayer") {
0116     edm::LogWarning("DD4hep_TestPixelTopology") << ddTopNodeName_ << "Not valid top MTD volume";
0117     return;
0118   }
0119 
0120   if (!pDD.isValid()) {
0121     edm::LogError("DD4hep_TestPixelTopology") << "ESTransientHandle<DDCompactView> pDD is not valid!";
0122     return;
0123   }
0124   if (pDD.description()) {
0125     edm::LogVerbatim("DD4hep_TestPixelTopology") << pDD.description()->type_ << " label: " << pDD.description()->label_;
0126   } else {
0127     edm::LogPrint("DD4hep_TestPixelTopology") << "NO label found pDD.description() returned false.";
0128   }
0129 
0130   if (!pSP.isValid()) {
0131     edm::LogError("DD4hep_TestPixelTopology") << "ESTransientHandle<DDSpecParRegistry> pSP is not valid!";
0132     return;
0133   }
0134 
0135   auto pTP = iSetup.getTransientHandle(mtdtopoToken_);
0136   if (!pTP.isValid()) {
0137     edm::LogError("DD4hep_TestPixelTopology") << "ESTransientHandle<MTDTopology> pTP is not valid!";
0138     return;
0139   } else {
0140     edm::LogVerbatim("DD4hep_TestPixelTopology")
0141         << "MTD topology mode = " << pTP.product()->getMTDTopologyMode() << " BtlLayout = "
0142         << static_cast<int>(MTDTopologyMode::crysLayoutFromTopoMode(pTP.product()->getMTDTopologyMode()))
0143         << " EtlLayout = "
0144         << static_cast<int>(MTDTopologyMode::etlLayoutFromTopoMode(pTP.product()->getMTDTopologyMode()));
0145     sunitt_ << "MTD topology mode = " << pTP.product()->getMTDTopologyMode() << " BtlLayout = "
0146             << static_cast<int>(MTDTopologyMode::crysLayoutFromTopoMode(pTP.product()->getMTDTopologyMode()))
0147             << " EtlLayout = "
0148             << static_cast<int>(MTDTopologyMode::etlLayoutFromTopoMode(pTP.product()->getMTDTopologyMode()));
0149   }
0150 
0151   auto pDG = iSetup.getTransientHandle(mtdgeoToken_);
0152   if (!pDG.isValid()) {
0153     edm::LogError("DD4hep_TestPixelTopology") << "ESTransientHandle<MTDGeometry> pDG is not valid!";
0154     return;
0155   }
0156 
0157   DDSpecParRefs specs;
0158   std::string attribute("MtdDDStructure"), name;
0159   bool isBarrel = false;
0160   if (ddTopNodeName_ == "BarrelTimingLayer") {
0161     edm::LogVerbatim("DD4hep_TestPixelTopology") << "  BTL MTDGeometry:\n";
0162     sunitt_ << "  BTL MTDGeometry:\n";
0163     name = "FastTimerHitsBarrel";
0164     isBarrel = true;
0165   } else if (ddTopNodeName_ == "EndcapTimingLayer") {
0166     edm::LogVerbatim("DD4hep_TestPixelTopology") << "  ETL MTDGeometry:\n";
0167     sunitt_ << "  ETL MTDGeometry:\n";
0168     name = "FastTimerHitsEndcap";
0169   } else {
0170     edm::LogError("DD4hep_TestPixelTopology") << "No correct sensitive detector provided, abort" << ddTopNodeName_;
0171     return;
0172   }
0173   pSP.product()->filter(specs, attribute, ddTopNodeName_);
0174   attribute = "ReadOutName";
0175   pSP.product()->filter(specs, attribute, name);
0176 
0177   edm::LogVerbatim("DD4hep_TestPixelTopology").log([&specs](auto& log) {
0178     log << "Filtered DD SpecPar Registry size: " << specs.size() << "\n";
0179     for (const auto& t : specs) {
0180       log << "\nSpecPar " << t.first << ":\nRegExps { ";
0181       for (const auto& ki : t.second->paths)
0182         log << ki << " ";
0183       log << "};\n ";
0184       for (const auto& kl : t.second->spars) {
0185         log << kl.first << " = ";
0186         for (const auto& kil : kl.second) {
0187           log << kil << " ";
0188         }
0189         log << "\n ";
0190       }
0191     }
0192   });
0193 
0194   std::vector<std::string_view> filterName;
0195   for (auto const& t : specs) {
0196     for (auto const& kl : t.second->spars) {
0197       if (kl.first == attribute) {
0198         for (auto const& it : t.second->paths) {
0199           filterName.emplace_back(it);
0200         }
0201       }
0202     }
0203   }
0204 
0205   DDFilteredView fv(pDD.product(), pDD.product()->description()->worldVolume());
0206   fv.mergedSpecifics(specs);
0207   fv.firstChild();
0208 
0209   bool write = false;
0210   bool exitLoop = false;
0211   uint32_t level(0);
0212   uint32_t count(0);
0213   uint32_t nSensBTL(0);
0214   uint32_t nSensETL(0);
0215   uint32_t oldgeoId(0);
0216 
0217   do {
0218     theBaseNumber(fv);
0219 
0220     auto print_path = [&]() {
0221       std::stringstream ss;
0222       ss << " - OCMS[0]/";
0223       for (int ii = thisN_.getLevels() - 1; ii-- > 0;) {
0224         ss << thisN_.getLevelName(ii);
0225         ss << "[";
0226         ss << thisN_.getCopyNumber(ii);
0227         ss << "]/";
0228       }
0229       return ss.str();
0230     };
0231 
0232     if (level > 0 && fv.navPos().size() < level && count == 2) {
0233       exitLoop = true;
0234     }
0235     if (dd4hep::dd::noNamespace(fv.name()) == ddTopNodeName_) {
0236       write = true;
0237       level = fv.navPos().size();
0238       count++;
0239     }
0240 
0241     // Test only the desired subdetector
0242 
0243     if (exitLoop) {
0244       break;
0245     }
0246 
0247     if (write) {
0248       // Actions for MTD volumes: search for sensitive detectors
0249 
0250       bool isSens = false;
0251 
0252       for (auto const& it : filterName) {
0253         if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(fv.name()), dd4hep::dd::realTopName(it))) {
0254           isSens = true;
0255           break;
0256         }
0257       }
0258 
0259       if (isSens) {
0260         DetId theId, geoId;
0261         BTLDetId theIdBTL, modIdBTL;
0262         ETLDetId theIdETL, modIdETL;
0263         if (isBarrel) {
0264           theIdBTL = btlNS_.getUnitID(thisN_);
0265           theId = theIdBTL;
0266           geoId = theIdBTL.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(pTP.product()->getMTDTopologyMode()));
0267           modIdBTL = geoId;
0268         } else {
0269           theIdETL = etlNS_.getUnitID(thisN_);
0270           theId = theIdETL;
0271           geoId = theIdETL.geographicalId();
0272           modIdETL = geoId;
0273         }
0274 
0275         const MTDGeomDet* thedet = pDG.product()->idToDet(geoId);
0276 
0277         if (dynamic_cast<const MTDGeomDetUnit*>((thedet)) == nullptr) {
0278           throw cms::Exception("DD4hep_TestPixelTopology")
0279               << "GeographicalID: " << std::hex << geoId.rawId() << " (" << theId.rawId()
0280               << ") with invalid MTDGeomDetUnit!" << std::dec << std::endl;
0281         }
0282 
0283         bool isNewId(false);
0284         if (geoId != oldgeoId) {
0285           oldgeoId = geoId;
0286           isNewId = true;
0287           if (isBarrel) {
0288             nSensBTL++;
0289           } else {
0290             nSensETL++;
0291           }
0292           const GeomDetUnit theDetUnit = *(dynamic_cast<const MTDGeomDetUnit*>(thedet));
0293 
0294           if (isBarrel) {
0295             edm::LogVerbatim("DD4hep_TestPixelTopology")
0296                 << "geoId= " << modIdBTL.rawId() << " side= " << modIdBTL.mtdSide()
0297                 << " RU/mod= " << modIdBTL.globalRunit() << " / " << modIdBTL.module();
0298             sunitt_ << "geoId= " << modIdBTL.rawId() << " side= " << modIdBTL.mtdSide()
0299                     << " RU/mod= " << modIdBTL.globalRunit() << " / " << modIdBTL.module();
0300           } else {
0301             edm::LogVerbatim("DD4hep_TestPixelTopology")
0302                 << "geoId= " << modIdETL.rawId() << " side= " << modIdETL.mtdSide()
0303                 << " disc/face/sec= " << modIdETL.nDisc() << " / " << modIdETL.discSide() << " / " << modIdETL.sector()
0304                 << " mod/typ/sens= " << modIdETL.module() << " / " << modIdETL.modType() << " / " << modIdETL.sensor();
0305             sunitt_ << "geoId= " << modIdETL.rawId() << " side= " << modIdETL.mtdSide()
0306                     << " disc/face/sec= " << modIdETL.nDisc() << " / " << modIdETL.discSide() << " / "
0307                     << modIdETL.sector() << " mod/typ/sens= " << modIdETL.module() << " / " << modIdETL.modType()
0308                     << " / " << modIdETL.sensor();
0309           }
0310           analyseRectangle(theDetUnit);
0311         }
0312 
0313         if (thedet == nullptr) {
0314           throw cms::Exception("DD4hep_TestPixelTopology") << "GeographicalID: " << std::hex << geoId.rawId() << " ("
0315                                                            << theId.rawId() << ") is invalid!" << std::dec << std::endl;
0316         }
0317         const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(thedet->topology());
0318         const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0319 
0320         int origRow(-1), origCol(-1), recoRow(-1), recoCol(-1);
0321         if (isBarrel) {
0322           origRow = theIdBTL.row(topo.nrows());
0323           origCol = theIdBTL.column(topo.nrows());
0324         }
0325 
0326         //
0327         // Test of positions for sensitive detectors
0328         //
0329 
0330         if (!dd4hep::isA<dd4hep::Box>(fv.solid())) {
0331           throw cms::Exception("DD4hep_TestPixelTopology") << "MTD sensitive element not a DDBox";
0332           break;
0333         }
0334         dd4hep::Box mySens(fv.solid());
0335 
0336         double xoffset(0.);
0337         double yoffset(0.);
0338         if (!isBarrel) {
0339           xoffset = topo.gapxBorder() + 0.5 * topo.pitch().first;
0340           yoffset = topo.gapyBorder() + 0.5 * topo.pitch().second;
0341         }
0342         DD3Vector zeroLocal(0., 0., 0.);
0343         DD3Vector cn1Local(mySens.x() - xoffset, mySens.y() - yoffset, mySens.z());
0344         DD3Vector cn2Local(-mySens.x() + xoffset, -mySens.y() + yoffset, -mySens.z());
0345         DD3Vector zeroGlobal = (fv.rotation())(zeroLocal) + fv.translation();
0346         DD3Vector cn1Global = (fv.rotation())(cn1Local) + fv.translation();
0347         DD3Vector cn2Global = (fv.rotation())(cn2Local) + fv.translation();
0348 
0349         const size_t nTest(3);
0350         std::array<Local3DPoint, nTest> refLocalPoints{{Local3DPoint(zeroLocal.x(), zeroLocal.y(), zeroLocal.z()),
0351                                                         Local3DPoint(cn1Local.x(), cn1Local.y(), cn1Local.z()),
0352                                                         Local3DPoint(cn2Local.x(), cn2Local.y(), cn2Local.z())}};
0353         std::array<DD3Vector, nTest> refGlobalPoints{{zeroGlobal, cn1Global, cn2Global}};
0354 
0355         for (size_t iloop = 0; iloop < nTest; iloop++) {
0356           Local3DPoint cmRefLocal(convertMmToCm(refLocalPoints[iloop].x() / dd4hep::mm),
0357                                   convertMmToCm(refLocalPoints[iloop].y() / dd4hep::mm),
0358                                   convertMmToCm(refLocalPoints[iloop].z() / dd4hep::mm));
0359 
0360           Local3DPoint modLocal, recoRefLocal;
0361           if (isBarrel) {
0362             // if BTL translate from crystal-local coordinates to module-local coordinates to get the row and column
0363             modLocal = topo.pixelToModuleLocalPoint(cmRefLocal, origRow, origCol);
0364             recoRefLocal = topo.moduleToPixelLocalPoint(modLocal);
0365             const auto& thepixel = topo.pixelIndex(modLocal);
0366             recoRow = thepixel.first;
0367             recoCol = thepixel.second;
0368 
0369             if (origRow != recoRow || origCol != recoCol) {
0370               edm::LogVerbatim("DD4hep_TestPixelTopology") << "DIFFERENCE row/col, orig= " << origRow << " " << origCol
0371                                                            << " reco= " << recoRow << " " << recoCol << "\n";
0372               sunitt_ << "DIFFERENCE row/col, orig= " << origRow << " " << origCol << " reco= " << recoRow << " "
0373                       << recoCol << "\n";
0374               recoRow = origRow;
0375               recoCol = origCol;
0376             }
0377           } else {
0378             // if ETL find the pixel corresponding to the referemce point, compute the pixel coordinate and convert back for check
0379             modLocal = cmRefLocal;
0380             const auto& thepixel = topo.pixelIndex(modLocal);
0381             Local3DPoint pixLocal = topo.moduleToPixelLocalPoint(modLocal);
0382             recoRefLocal = topo.pixelToModuleLocalPoint(pixLocal, thepixel.first, thepixel.second);
0383             recoRow = thepixel.first;
0384             recoCol = thepixel.second;
0385           }
0386 
0387           // reconstructed global position from reco geometry and rectangluar MTD topology
0388 
0389           const auto& modGlobal = thedet->toGlobal(modLocal);
0390 
0391           if (isNewId && iloop == nTest - 1) {
0392             edm::LogVerbatim("DD4hep_TestPixelTopology")
0393                 << "row/col= " << recoRow << " / " << recoCol << " local pos= " << fvecround(modLocal, 4)
0394                 << " global pos= " << fvecround(modGlobal, 4) << "\n";
0395             sunitt_ << "row/col= " << recoRow << " / " << recoCol << " local pos= " << fvecround(modLocal, 2)
0396                     << " global pos= " << fvecround(modGlobal, 2) << "\n";
0397           }
0398 
0399           const double deltax = convertCmToMm(modGlobal.x()) - (refGlobalPoints[iloop].x() / dd4hep::mm);
0400           const double deltay = convertCmToMm(modGlobal.y()) - (refGlobalPoints[iloop].y() / dd4hep::mm);
0401           const double deltaz = convertCmToMm(modGlobal.z()) - (refGlobalPoints[iloop].z() / dd4hep::mm);
0402 
0403           const double local_deltax = recoRefLocal.x() - cmRefLocal.x();
0404           const double local_deltay = recoRefLocal.y() - cmRefLocal.y();
0405           const double local_deltaz = recoRefLocal.z() - cmRefLocal.z();
0406 
0407           if (std::abs(deltax) > tolerance || std::abs(deltay) > tolerance || std::abs(deltaz) > tolerance ||
0408               std::abs(local_deltax) > tolerance || std::abs(local_deltay) > tolerance ||
0409               std::abs(local_deltaz) > tolerance) {
0410             edm::LogVerbatim("DD4hep_TestPixelTopology") << print_path() << "\n";
0411             sunitt_ << print_path() << "\n";
0412             if (isBarrel) {
0413               edm::LogVerbatim("DD4hep_TestPixelTopology")
0414                   << "rawId= " << theIdBTL.rawId() << " geoId= " << geoId.rawId() << " side/rod= " << theIdBTL.mtdSide()
0415                   << " / " << theIdBTL.mtdRR() << " RU= " << theIdBTL.globalRunit()
0416                   << " module/geomodule= " << theIdBTL.module() << " / " << static_cast<BTLDetId>(geoId).module()
0417                   << " crys= " << theIdBTL.crystal() << " BTLDetId row/col= " << origRow << " / " << origCol << "\n";
0418               sunitt_ << "rawId= " << theIdBTL.rawId() << " geoId= " << geoId.rawId()
0419                       << " side/rod= " << theIdBTL.mtdSide() << " / " << theIdBTL.mtdRR()
0420                       << " RU= " << theIdBTL.globalRunit() << " module/geomodule= " << theIdBTL.module() << " / "
0421                       << static_cast<BTLDetId>(geoId).module() << " crys= " << theIdBTL.crystal()
0422                       << " BTLDetId row/col= " << origRow << " / " << origCol << "\n";
0423             } else {
0424               edm::LogVerbatim("DD4hep_TestPixelTopology")
0425                   << "geoId= " << modIdETL.rawId() << " side= " << modIdETL.mtdSide()
0426                   << " disc/face/sec= " << modIdETL.nDisc() << " / " << modIdETL.discSide() << " / "
0427                   << modIdETL.sector() << " mod/typ/sens= " << modIdETL.module() << " / " << modIdETL.modType() << " / "
0428                   << modIdETL.sensor() << "\n";
0429               sunitt_ << "geoId= " << modIdETL.rawId() << " side= " << modIdETL.mtdSide()
0430                       << " disc/face/sec= " << modIdETL.nDisc() << " / " << modIdETL.discSide() << " / "
0431                       << modIdETL.sector() << " mod/typ/sens= " << modIdETL.module() << " / " << modIdETL.modType()
0432                       << " / " << modIdETL.sensor() << "\n";
0433             }
0434 
0435             edm::LogVerbatim("DD4hep_TestPixelTopology")
0436                 << "Ref#" << iloop << " local= " << fround(refLocalPoints[iloop].x() / dd4hep::mm, 4)
0437                 << fround(refLocalPoints[iloop].y() / dd4hep::mm, 4)
0438                 << fround(refLocalPoints[iloop].z() / dd4hep::mm, 4)
0439                 << " Orig global= " << fround(refGlobalPoints[iloop].x() / dd4hep::mm, 4)
0440                 << fround(refGlobalPoints[iloop].y() / dd4hep::mm, 4)
0441                 << fround(refGlobalPoints[iloop].z() / dd4hep::mm, 4)
0442                 << " Reco global= " << fround(convertCmToMm(modGlobal.x()), 4)
0443                 << fround(convertCmToMm(modGlobal.y()), 4) << fround(convertCmToMm(modGlobal.z()), 4)
0444                 << " Delta= " << fround(deltax, 4) << fround(deltay, 4) << fround(deltaz, 4)
0445                 << " Local Delta= " << fround(local_deltax, 4) << fround(local_deltay, 4) << fround(local_deltaz, 4)
0446                 << "\n";
0447             sunitt_ << "Ref#" << iloop << " local= " << fround(refLocalPoints[iloop].x() / dd4hep::mm, 2)
0448                     << fround(refLocalPoints[iloop].y() / dd4hep::mm, 2)
0449                     << fround(refLocalPoints[iloop].z() / dd4hep::mm, 2)
0450                     << " Orig global= " << fround(refGlobalPoints[iloop].x() / dd4hep::mm, 2)
0451                     << fround(refGlobalPoints[iloop].y() / dd4hep::mm, 2)
0452                     << fround(refGlobalPoints[iloop].z() / dd4hep::mm, 2)
0453                     << " Reco global= " << fround(convertCmToMm(modGlobal.x()), 2)
0454                     << fround(convertCmToMm(modGlobal.y()), 2) << fround(convertCmToMm(modGlobal.z()), 2)
0455                     << " Delta= " << fround(deltax, 2) << fround(deltay, 2) << fround(deltaz, 2)
0456                     << " Local Delta= " << fround(local_deltax, 2) << fround(local_deltay, 2) << fround(local_deltaz, 2)
0457                     << "\n";
0458 
0459             if (std::abs(deltax) > tolerance || std::abs(deltay) > tolerance || std::abs(deltaz) > tolerance) {
0460               edm::LogVerbatim("DD4hep_TestPixelTopology")
0461                   << "DIFFERENCE detId/ref# " << theId.rawId() << " " << iloop << " dx/dy/dz= " << fround(deltax, 4)
0462                   << fround(deltay, 4) << fround(deltaz, 4) << "\n";
0463               sunitt_ << "DIFFERENCE detId/ref# " << theId.rawId() << " " << iloop << " dx/dy/dz= " << fround(deltax, 2)
0464                       << fround(deltay, 2) << fround(deltaz, 2) << "\n";
0465             }
0466             if (std::abs(local_deltax) > tolerance || std::abs(local_deltay) > tolerance ||
0467                 std::abs(local_deltaz) > tolerance) {
0468               edm::LogVerbatim("DD4hep_TestPixelTopology")
0469                   << "DIFFERENCE detId/ref# " << theId.rawId() << " " << iloop
0470                   << " local dx/dy/dz= " << fround(local_deltax, 4) << fround(local_deltay, 4)
0471                   << fround(local_deltaz, 4) << "\n";
0472               sunitt_ << "DIFFERENCE detId/ref# " << theId.rawId() << " " << iloop
0473                       << " local dx/dy/dz= " << fround(local_deltax, 2) << fround(local_deltay, 2)
0474                       << fround(local_deltaz, 2) << "\n";
0475             }
0476           }
0477         }
0478       }
0479     }
0480   } while (fv.next(0));
0481 
0482   if (isBarrel && nSensBTL != pDG.product()->detsBTL().size()) {
0483     edm::LogVerbatim("DD4hep_TestPixelTopology")
0484         << "DIFFERENCE #ideal = " << nSensBTL << " #reco = " << pDG.product()->detsBTL().size()
0485         << " BTL module numbers are not matching!";
0486     sunitt_ << "DIFFERENCE #ideal = " << nSensBTL << " #reco = " << pDG.product()->detsBTL().size()
0487             << " BTL module numbers are not matching!";
0488   }
0489 
0490   if (!isBarrel && nSensETL != pDG.product()->detsETL().size()) {
0491     edm::LogVerbatim("DD4hep_TestPixelTopology")
0492         << "DIFFERENCE #ideal = " << nSensETL << " #reco = " << pDG.product()->detsBTL().size()
0493         << " ETL module numbers are not matching!";
0494     sunitt_ << "DIFFERENCE #ideal = " << nSensETL << " #reco = " << pDG.product()->detsBTL().size()
0495             << " ETL module numbers are not matching!";
0496   }
0497 
0498   if (!sunitt_.str().empty()) {
0499     edm::LogVerbatim("MTDUnitTest") << sunitt_.str();
0500   }
0501 }
0502 
0503 void DD4hep_TestPixelTopology::theBaseNumber(cms::DDFilteredView& fv) {
0504   thisN_.reset();
0505   thisN_.setSize(fv.navPos().size());
0506 
0507   for (uint ii = 0; ii < fv.navPos().size(); ii++) {
0508     std::string_view name((fv.geoHistory()[ii])->GetName());
0509     size_t ipos = name.rfind('_');
0510     thisN_.addLevel(name.substr(0, ipos), fv.copyNos()[ii]);
0511 #ifdef EDM_ML_DEBUG
0512     edm::LogVerbatim("DD4hep_TestPixelTopology") << ii << " " << name.substr(0, ipos) << " " << fv.copyNos()[ii];
0513 #endif
0514   }
0515 }
0516 
0517 void DD4hep_TestPixelTopology::analyseRectangle(const GeomDetUnit& det) {
0518   const double safety = 0.9999;
0519 
0520   const BoundPlane& p = det.specificSurface();
0521   const Bounds& bounds = det.surface().bounds();
0522   const RectangularPlaneBounds* tb = dynamic_cast<const RectangularPlaneBounds*>(&bounds);
0523   if (tb == nullptr)
0524     return;  // not trapezoidal
0525 
0526   const GlobalPoint& pos = det.position();
0527   double length = tb->length();
0528   double width = tb->width();
0529   double thickness = tb->thickness();
0530 
0531   GlobalVector yShift = det.surface().toGlobal(LocalVector(0, 0, safety * length / 2.));
0532   GlobalPoint outerMiddle = pos + yShift;
0533   GlobalPoint innerMiddle = pos + (-1. * yShift);
0534   if (outerMiddle.perp() < innerMiddle.perp())
0535     std::swap(outerMiddle, innerMiddle);
0536 
0537   edm::LogVerbatim("DD4hep_TestPixelTopology")
0538       << " " << fvecround(pos, 4) << " R= " << fround(std::sqrt(pos.x() * pos.x() + pos.y() * pos.y()), 4)
0539       << " phi= " << fround(convertRadToDeg(pos.phi()), 4) << " outerMiddle " << fvecround(outerMiddle, 4) << "\n"
0540       << " l/w/t " << fround(length, 4) << " / " << fround(width, 4) << " / " << fround(thickness, 4)
0541       << " RadLeng= " << p.mediumProperties().radLen() << " Xi= " << p.mediumProperties().xi()
0542       << " det center inside bounds? " << tb->inside(det.surface().toLocal(pos)) << "\n";
0543   sunitt_ << " " << fvecround(pos, 2) << " R= " << fround(std::sqrt(pos.x() * pos.x() + pos.y() * pos.y()), 2)
0544           << " phi= " << fround(convertRadToDeg(pos.phi()), 2) << " outerMiddle " << fvecround(outerMiddle, 2) << "\n"
0545           << " l/w/t " << fround(length, 2) << " / " << fround(width, 2) << " / " << fround(thickness, 2)
0546           << " RadLeng= " << p.mediumProperties().radLen() << " Xi= " << p.mediumProperties().xi()
0547           << " det center inside bounds? " << tb->inside(det.surface().toLocal(pos)) << "\n";
0548 
0549   checkRotation(det);
0550 }
0551 
0552 void DD4hep_TestPixelTopology::checkRotation(const GeomDetUnit& det) {
0553   const double eps = 10. * std::numeric_limits<float>::epsilon();
0554   static int first = 0;
0555   if (first == 0) {
0556     edm::LogVerbatim("DD4hep_TestPixelTopology") << "factor x numeric_limits<float>::epsilon() " << eps;
0557     first = 1;
0558   }
0559 
0560   const Surface::RotationType& rot(det.surface().rotation());
0561   GlobalVector a(rot.xx(), rot.xy(), rot.xz());
0562   GlobalVector b(rot.yx(), rot.yy(), rot.yz());
0563   GlobalVector c(rot.zx(), rot.zy(), rot.zz());
0564   GlobalVector cref = a.cross(b);
0565   GlobalVector aref = b.cross(c);
0566   GlobalVector bref = c.cross(a);
0567   if ((a - aref).mag() > eps || (b - bref).mag() > eps || (c - cref).mag() > eps) {
0568     edm::LogVerbatim("DD4hep_TestPixelTopology")
0569         << " DIFFERENCE Rotation not good by cross product: " << (a - aref).mag() << ", " << (b - bref).mag() << ", "
0570         << (c - cref).mag() << " for det at pos " << det.surface().position() << "\n";
0571     sunitt_ << " DIFFERENCE Rotation not good by cross product: " << (a - aref).mag() << ", " << (b - bref).mag()
0572             << ", " << (c - cref).mag() << " for det at pos " << det.surface().position() << "\n";
0573   }
0574   if (fabs(a.mag() - 1.) > eps || fabs(b.mag() - 1.) > eps || fabs(c.mag() - 1.) > eps) {
0575     edm::LogVerbatim("DD4hep_TestPixelTopology")
0576         << " DIFFERENCE Rotation not good by vector mag: " << (a).mag() << ", " << (b).mag() << ", " << (c).mag()
0577         << " for det at pos " << det.surface().position() << "\n";
0578     sunitt_ << " DIFFERENCE Rotation not good by vector mag: " << (a).mag() << ", " << (b).mag() << ", " << (c).mag()
0579             << " for det at pos " << det.surface().position() << "\n";
0580   }
0581 }
0582 
0583 DEFINE_FWK_MODULE(DD4hep_TestPixelTopology);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!