MTDGeometryBuilder/test/MTDDigiGeometryAnalyzer.cc

0001 // user include files
0002 #include "FWCore/Framework/interface/Frameworkfwd.h"
0003 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0004
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/Framework/interface/ESHandle.h"
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h"
0012 #include "Geometry/Records/interface/MTDDigiGeometryRecord.h"
0013 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0014 #include "Geometry/MTDGeometryBuilder/interface/MTDGeomDetType.h"
0015 #include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
0016
0017 #include "DataFormats/ForwardDetId/interface/MTDDetId.h"
0018 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0019 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0020 #include "Geometry/MTDGeometryBuilder/interface/MTDGeomDetUnit.h"
0021 #include "DataFormats/GeometrySurface/interface/MediumProperties.h"
0022 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024
0025 #include "DataFormats/Math/interface/Rounding.h"
0026
0027 #include <fstream>
0028
0029 #include "CLHEP/Random/RandFlat.h"
0030
0031 // class declaration
0032
0033 class MTDDigiGeometryAnalyzer : public edm::one::EDAnalyzer<> {
0034 public:
0035   explicit MTDDigiGeometryAnalyzer(const edm::ParameterSet&);
0036   ~MTDDigiGeometryAnalyzer() override = default;
0037
0038   void beginJob() override {}
0039   void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0040   void endJob() override {}
0041
0042 private:
0043   void analyseRectangle(const GeomDetUnit& det);
0044   void checkRotation(const GeomDetUnit& det);
0045   void checkRectangularMTDTopology(const RectangularMTDTopology&);
0046   void checkPixelsAcceptance(const GeomDetUnit& det);
0047
0048   std::stringstream sunitt;
0049
0050   edm::ESGetToken<MTDGeometry, MTDDigiGeometryRecord> mtdgeoToken_;
0051 };
0052
0053 using cms_rounding::roundIfNear0, cms_rounding::roundVecIfNear0;
0054
0055 MTDDigiGeometryAnalyzer::MTDDigiGeometryAnalyzer(const edm::ParameterSet& iConfig) {
0056   mtdgeoToken_ = esConsumes<MTDGeometry, MTDDigiGeometryRecord>();
0057 }
0058
0059 // ------------ method called to produce the data  ------------
0060 void MTDDigiGeometryAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0061   //
0062   // get the MTDGeometry
0063   //
0064   auto pDD = iSetup.getTransientHandle(mtdgeoToken_);
0065   edm::LogInfo("MTDDigiGeometryAnalyzer")
0066       << "Geometry node for MTDGeom is  " << &(*pDD) << "\n"
0067       << " # detectors = " << pDD->detUnits().size() << "\n"
0068       << " # types     = " << pDD->detTypes().size() << "\n"
0069       << " # BTL dets  = " << pDD->detsBTL().size() << "\n"
0070       << " # ETL dets  = " << pDD->detsETL().size() << "\n"
0071       << " # layers " << pDD->geomDetSubDetector(1) << "  = " << pDD->numberOfLayers(1) << "\n"
0072       << " # layers " << pDD->geomDetSubDetector(2) << "  = " << pDD->numberOfLayers(2) << "\n";
0073   sunitt << std::fixed << std::setw(7) << pDD->detUnits().size() << std::setw(7) << pDD->detTypes().size() << "\n";
0074   for (auto const& it : pDD->detUnits()) {
0075     if (dynamic_cast<const MTDGeomDetUnit*>((it)) != nullptr) {
0076       const BoundPlane& p = (dynamic_cast<const MTDGeomDetUnit*>((it)))->specificSurface();
0077       const MTDDetId mtdId(it->geographicalId());
0078       std::stringstream moduleLabel;
0079       if (mtdId.mtdSubDetector() == 1) {
0080         const BTLDetId btlId(it->geographicalId());
0081         moduleLabel << " BTL side " << btlId.mtdSide() << " Rod " << btlId.mtdRR() << " type/RU " << btlId.modType()
0082                     << "/" << btlId.runit() << " mod " << btlId.module();
0083       } else if (mtdId.mtdSubDetector() == 2) {
0084         const ETLDetId etlId(it->geographicalId());
0085         moduleLabel << " ETL side " << mtdId.mtdSide() << " Disc/Side/Sector " << etlId.nDisc() << " "
0086                     << etlId.discSide() << " " << etlId.sector();
0087       } else {
0088         edm::LogWarning("MTDDigiGeometryanalyzer") << (it->geographicalId()).rawId() << " unknown MTD subdetector!";
0089       }
0090       edm::LogVerbatim("MTDDigiGeometryAnalyzer")
0091           << "---------------------------------------------------------- \n"
0092           << it->geographicalId().rawId() << moduleLabel.str() << " RadLeng Pixel " << p.mediumProperties().radLen()
0093           << " Xi Pixel " << p.mediumProperties().xi();
0094
0095       const GeomDetUnit theDet = *(dynamic_cast<const MTDGeomDetUnit*>(it));
0096       analyseRectangle(theDet);
0097     }
0098   }
0099
0100   for (auto const& it : pDD->detTypes()) {
0101     if (dynamic_cast<const MTDGeomDetType*>((it)) != nullptr) {
0102       const PixelTopology& p = (dynamic_cast<const MTDGeomDetType*>((it)))->specificTopology();
0103       const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(p);
0104       auto pitchval = topo.pitch();
0105       edm::LogVerbatim("MTDDigiGeometryAnalyzer")
0106           << "\n Subdetector " << it->subDetector() << " MTD Det " << it->name() << "\n"
0107           << " Rows     " << topo.nrows() << " Columns " << topo.ncolumns() << " ROCS X   " << topo.rocsX()
0108           << " ROCS Y  " << topo.rocsY() << " Rows/ROC " << topo.rowsperroc() << " Cols/ROC " << topo.colsperroc()
0109           << " Pitch X " << pitchval.first << " Pitch Y " << pitchval.second << " Sensor Interpad X "
0110           << topo.gapxInterpad() << " Sensor Interpad Y " << topo.gapyInterpad() << " Sensor Border X "
0111           << topo.gapxBorder() << " Sensor Border Y " << topo.gapyBorder();
0112       sunitt << std::fixed << std::setw(7) << it->subDetector() << std::setw(4) << topo.nrows() << std::setw(4)
0113              << topo.ncolumns() << std::setw(4) << std::setw(4) << topo.rocsX() << std::setw(4) << topo.rocsY()
0114              << std::setw(4) << topo.rowsperroc() << std::setw(4) << topo.colsperroc() << std::setw(10)
0115              << pitchval.first << std::setw(10) << pitchval.second << std::setw(10) << topo.gapxInterpad()
0116              << std::setw(10) << topo.gapyInterpad() << std::setw(10) << topo.gapxBorder() << std::setw(10)
0117              << topo.gapyBorder() << "\n";
0118       checkRectangularMTDTopology(topo);
0119     }
0120   }
0121
0122   edm::LogInfo("MTDDigiGeometryAnalyzer") << "Acceptance of BTL module:";
0123   auto const& btldet = *(dynamic_cast<const MTDGeomDetUnit*>(pDD->detsBTL().front()));
0124   checkPixelsAcceptance(btldet);
0125   edm::LogInfo("MTDDigiGeometryAnalyzer") << "Acceptance of ETL module:";
0126   auto const& etldet = *(dynamic_cast<const MTDGeomDetUnit*>(pDD->detsETL().front()));
0127   checkPixelsAcceptance(etldet);
0128
0129   edm::LogInfo("MTDDigiGeometryAnalyzer") << "Additional MTD geometry content:"
0130                                           << "\n"
0131                                           << " # dets            = " << pDD->dets().size() << "\n"
0132                                           << " # detUnitIds      = " << pDD->detUnitIds().size() << "\n"
0133                                           << " # detIds          = " << pDD->detIds().size() << "\n";
0134   sunitt << std::fixed << std::setw(7) << pDD->dets().size() << std::setw(7) << pDD->detUnitIds().size() << std::setw(7)
0135          << pDD->detIds().size() << "\n";
0136
0137   edm::LogVerbatim("MTDUnitTest") << sunitt.str();
0138 }
0139
0140 void MTDDigiGeometryAnalyzer::checkRectangularMTDTopology(const RectangularMTDTopology& topo) {
0141   std::stringstream pixelinfo;
0142   pixelinfo << "Pixel center location:\n";
0143   LocalPoint center(0, 0, 0);
0144   for (int r = 0; r < topo.nrows(); r++) {
0145     for (int c = 0; c < topo.ncolumns(); c++) {
0146       sunitt << r << " " << c << " " << topo.pixelToModuleLocalPoint(center, r, c) << "\n";
0147       pixelinfo << r << " " << c << " " << topo.pixelToModuleLocalPoint(center, r, c) << "\n";
0148     }
0149   }
0150   edm::LogVerbatim("MTDDigiGeometryAnalyzer") << pixelinfo.str();
0151 }
0152
0153 void MTDDigiGeometryAnalyzer::analyseRectangle(const GeomDetUnit& det) {
0154   const double safety = 0.9999;
0155
0156   const Bounds& bounds = det.surface().bounds();
0157   const RectangularPlaneBounds* tb = dynamic_cast<const RectangularPlaneBounds*>(&bounds);
0158   if (tb == nullptr)
0159     return;  // not trapezoidal
0160
0161   const GlobalPoint& pos = det.position();
0162   double length = tb->length();
0163   double width = tb->width();
0164   double thickness = tb->thickness();
0165
0166   GlobalVector yShift = det.surface().toGlobal(LocalVector(0, 0, safety * length / 2.));
0167   GlobalPoint outerMiddle = pos + yShift;
0168   GlobalPoint innerMiddle = pos + (-1. * yShift);
0169   if (outerMiddle.perp() < innerMiddle.perp())
0170     std::swap(outerMiddle, innerMiddle);
0171
0172   auto fround = [&](double in) {
0173     std::stringstream ss;
0174     ss << std::fixed << std::setw(14) << roundIfNear0(in);
0175     return ss.str();
0176   };
0177
0178   auto fvecround = [&](GlobalPoint vecin) {
0179     std::stringstream ss;
0180     ss << std::fixed << std::setw(14) << roundVecIfNear0(vecin);
0181     return ss.str();
0182   };
0183
0184   edm::LogVerbatim("MTDDigiGeometryAnalyzer")
0185       << "Det at pos " << fvecround(pos) << " radius " << fround(std::sqrt(pos.x() * pos.x() + pos.y() * pos.y()))
0186       << " has length " << fround(length) << " width " << fround(width) << " thickness " << fround(thickness) << "\n"
0187       << "det center inside bounds? " << tb->inside(det.surface().toLocal(pos)) << "\n"
0188       << "outerMiddle " << fvecround(outerMiddle);
0189   sunitt << det.geographicalId().rawId() << fvecround(pos) << fround(length) << fround(width) << fround(thickness)
0190          << tb->inside(det.surface().toLocal(pos)) << fvecround(outerMiddle) << "\n";
0191
0192   checkRotation(det);
0193 }
0194
0195 void MTDDigiGeometryAnalyzer::checkRotation(const GeomDetUnit& det) {
0196   const double eps = std::numeric_limits<float>::epsilon();
0197   static int first = 0;
0198   if (first == 0) {
0199     edm::LogVerbatim("MTDDigiGeometryAnalyzer")
0200         << "numeric_limits<float>::epsilon() " << std::numeric_limits<float>::epsilon();
0201     first = 1;
0202   }
0203
0204   const Surface::RotationType& rot(det.surface().rotation());
0205   GlobalVector a(rot.xx(), rot.xy(), rot.xz());
0206   GlobalVector b(rot.yx(), rot.yy(), rot.yz());
0207   GlobalVector c(rot.zx(), rot.zy(), rot.zz());
0208   GlobalVector cref = a.cross(b);
0209   GlobalVector aref = b.cross(c);
0210   GlobalVector bref = c.cross(a);
0211   if ((a - aref).mag() > eps || (b - bref).mag() > eps || (c - cref).mag() > eps) {
0212     edm::LogWarning("MTDDigiGeometryAnalyzer")
0213         << " Rotation not good by cross product: " << (a - aref).mag() << ", " << (b - bref).mag() << ", "
0214         << (c - cref).mag() << " for det at pos " << det.surface().position();
0215   }
0216   if (fabs(a.mag() - 1.) > eps || fabs(b.mag() - 1.) > eps || fabs(c.mag() - 1.) > eps) {
0217     edm::LogWarning("MTDDigiGeometryAnalyzer") << " Rotation not good by bector mag: " << (a).mag() << ", " << (b).mag()
0218                                                << ", " << (c).mag() << " for det at pos " << det.surface().position();
0219   }
0220 }
0221
0222 void MTDDigiGeometryAnalyzer::checkPixelsAcceptance(const GeomDetUnit& det) {
0223   const Bounds& bounds = det.surface().bounds();
0224   const RectangularPlaneBounds* tb = dynamic_cast<const RectangularPlaneBounds*>(&bounds);
0225   if (tb == nullptr)
0226     return;  // not trapezoidal
0227
0228   double length = tb->length();
0229   double width = tb->width();
0230   edm::LogVerbatim("MTDDigiGeometryAnalyzer") << "X (width) = " << width << " Y (length) = " << length;
0231
0232   const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(det.topology());
0233   const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0234
0235   const size_t maxindex = 100000;
0236   size_t inpixel(0);
0237   for (size_t index = 0; index < maxindex; index++) {
0238     double ax = CLHEP::RandFlat::shoot(-width * 0.5, width * 0.5);
0239     double ay = CLHEP::RandFlat::shoot(-length * 0.5, length * 0.5);
0240     LocalPoint hit(ax, ay, 0);
0241     auto const indici = topo.pixelIndex(hit);
0242     assert(indici.first < topo.nrows() && indici.second < topo.ncolumns());  // sanity check on the index definition
0243     if (topo.isInPixel(hit)) {
0244       inpixel++;
0245     }
0246   }
0247   double acc = (double)inpixel / (double)maxindex;
0248   double accerr = std::sqrt(acc * (1. - acc) / (double)maxindex);
0249   edm::LogVerbatim("MTDDigiGeometryAnalyzer") << "Acceptance: " << acc << " +/- " << accerr;
0250 }
0251
0252 //define this as a plug-in
0253 DEFINE_FWK_MODULE(MTDDigiGeometryAnalyzer);