Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-04-14 02:22:42

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 // class declaration
0030 
0031 class MTDDigiGeometryAnalyzer : public edm::one::EDAnalyzer<> {
0032 public:
0033   explicit MTDDigiGeometryAnalyzer(const edm::ParameterSet&);
0034   ~MTDDigiGeometryAnalyzer() override = default;
0035 
0036   void beginJob() override {}
0037   void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0038   void endJob() override {}
0039 
0040 private:
0041   void analyseRectangle(const GeomDetUnit& det);
0042   void checkRotation(const GeomDetUnit& det);
0043   void checkRectangularMTDTopology(const RectangularMTDTopology&);
0044 
0045   std::stringstream sunitt;
0046 
0047   edm::ESGetToken<MTDGeometry, MTDDigiGeometryRecord> mtdgeoToken_;
0048 };
0049 
0050 using cms_rounding::roundIfNear0, cms_rounding::roundVecIfNear0;
0051 
0052 MTDDigiGeometryAnalyzer::MTDDigiGeometryAnalyzer(const edm::ParameterSet& iConfig) {
0053   mtdgeoToken_ = esConsumes<MTDGeometry, MTDDigiGeometryRecord>();
0054 }
0055 
0056 // ------------ method called to produce the data  ------------
0057 void MTDDigiGeometryAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0058   //
0059   // get the MTDGeometry
0060   //
0061   auto pDD = iSetup.getTransientHandle(mtdgeoToken_);
0062   edm::LogInfo("MTDDigiGeometryAnalyzer")
0063       << "Geometry node for MTDGeom is  " << &(*pDD) << "\n"
0064       << " # detectors = " << pDD->detUnits().size() << "\n"
0065       << " # types     = " << pDD->detTypes().size() << "\n"
0066       << " # BTL dets  = " << pDD->detsBTL().size() << "\n"
0067       << " # ETL dets  = " << pDD->detsETL().size() << "\n"
0068       << " # layers " << pDD->geomDetSubDetector(1) << "  = " << pDD->numberOfLayers(1) << "\n"
0069       << " # layers " << pDD->geomDetSubDetector(2) << "  = " << pDD->numberOfLayers(2) << "\n";
0070   sunitt << std::fixed << std::setw(7) << pDD->detUnits().size() << std::setw(7) << pDD->detTypes().size() << "\n";
0071   for (auto const& it : pDD->detUnits()) {
0072     if (dynamic_cast<const MTDGeomDetUnit*>((it)) != nullptr) {
0073       const BoundPlane& p = (dynamic_cast<const MTDGeomDetUnit*>((it)))->specificSurface();
0074       const MTDDetId mtdId(it->geographicalId());
0075       std::stringstream moduleLabel;
0076       if (mtdId.mtdSubDetector() == 1) {
0077         moduleLabel << " BTL side " << mtdId.mtdSide() << " Rod " << mtdId.mtdRR() << " mod "
0078                     << (static_cast<const BTLDetId>(mtdId)).module();
0079       } else if (mtdId.mtdSubDetector() == 2) {
0080         const ETLDetId etlId(it->geographicalId());
0081         moduleLabel << " ETL side " << mtdId.mtdSide() << " Disc/Side/Sector " << etlId.nDisc() << " "
0082                     << etlId.discSide() << " " << etlId.sector();
0083       } else {
0084         edm::LogWarning("MTDDigiGeometryanalyzer") << (it->geographicalId()).rawId() << " unknown MTD subdetector!";
0085       }
0086       edm::LogVerbatim("MTDDigiGeometryAnalyzer")
0087           << "---------------------------------------------------------- \n"
0088           << it->geographicalId().rawId() << moduleLabel.str() << " RadLeng Pixel " << p.mediumProperties().radLen()
0089           << " Xi Pixel " << p.mediumProperties().xi();
0090 
0091       const GeomDetUnit theDet = *(dynamic_cast<const MTDGeomDetUnit*>(it));
0092       analyseRectangle(theDet);
0093     }
0094   }
0095 
0096   for (auto const& it : pDD->detTypes()) {
0097     if (dynamic_cast<const MTDGeomDetType*>((it)) != nullptr) {
0098       const PixelTopology& p = (dynamic_cast<const MTDGeomDetType*>((it)))->specificTopology();
0099       const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(p);
0100       auto pitchval = topo.pitch();
0101       edm::LogVerbatim("MTDDigiGeometryAnalyzer")
0102           << "\n Subdetector " << it->subDetector() << " MTD Det " << it->name() << "\n"
0103           << " Rows     " << topo.nrows() << " Columns " << topo.ncolumns() << " ROCS X   " << topo.rocsX()
0104           << " ROCS Y  " << topo.rocsY() << " Rows/ROC " << topo.rowsperroc() << " Cols/ROC " << topo.colsperroc()
0105           << " Pitch X " << pitchval.first << " Pitch Y " << pitchval.second << " Sensor Interpad X "
0106           << topo.gapxInterpad() << " Sensor Interpad Y " << topo.gapyInterpad() << " Sensor Border X "
0107           << topo.gapxBorder() << " Sensor Border Y " << topo.gapyBorder();
0108       sunitt << std::fixed << std::setw(7) << it->subDetector() << std::setw(4) << topo.nrows() << std::setw(4)
0109              << topo.ncolumns() << std::setw(4) << std::setw(4) << topo.rocsX() << std::setw(4) << topo.rocsY()
0110              << std::setw(4) << topo.rowsperroc() << std::setw(4) << topo.colsperroc() << std::setw(10)
0111              << pitchval.first << std::setw(10) << pitchval.second << std::setw(10) << topo.gapxInterpad()
0112              << std::setw(10) << topo.gapyInterpad() << std::setw(10) << topo.gapxBorder() << std::setw(10)
0113              << topo.gapyBorder() << "\n";
0114       checkRectangularMTDTopology(topo);
0115     }
0116   }
0117 
0118   edm::LogInfo("MTDDigiGeometryAnalyzer") << "Additional MTD geometry content:"
0119                                           << "\n"
0120                                           << " # dets            = " << pDD->dets().size() << "\n"
0121                                           << " # detUnitIds      = " << pDD->detUnitIds().size() << "\n"
0122                                           << " # detIds          = " << pDD->detIds().size() << "\n";
0123   sunitt << std::fixed << std::setw(7) << pDD->dets().size() << std::setw(7) << pDD->detUnitIds().size() << std::setw(7)
0124          << pDD->detIds().size() << "\n";
0125 
0126   edm::LogVerbatim("MTDUnitTest") << sunitt.str();
0127 }
0128 
0129 void MTDDigiGeometryAnalyzer::checkRectangularMTDTopology(const RectangularMTDTopology& topo) {
0130   std::stringstream pixelinfo;
0131   pixelinfo << "Pixel center location:\n";
0132   LocalPoint center(0, 0, 0);
0133   for (int r = 0; r < topo.nrows(); r++) {
0134     for (int c = 0; c < topo.ncolumns(); c++) {
0135       sunitt << r << " " << c << " " << topo.pixelToModuleLocalPoint(center, r, c) << "\n";
0136       pixelinfo << r << " " << c << " " << topo.pixelToModuleLocalPoint(center, r, c) << "\n";
0137     }
0138   }
0139   edm::LogVerbatim("MTDDigiGeometryAnalyzer") << pixelinfo.str();
0140 }
0141 
0142 void MTDDigiGeometryAnalyzer::analyseRectangle(const GeomDetUnit& det) {
0143   const double safety = 0.9999;
0144 
0145   const Bounds& bounds = det.surface().bounds();
0146   const RectangularPlaneBounds* tb = dynamic_cast<const RectangularPlaneBounds*>(&bounds);
0147   if (tb == nullptr)
0148     return;  // not trapezoidal
0149 
0150   const GlobalPoint& pos = det.position();
0151   double length = tb->length();
0152   double width = tb->width();
0153   double thickness = tb->thickness();
0154 
0155   GlobalVector yShift = det.surface().toGlobal(LocalVector(0, 0, safety * length / 2.));
0156   GlobalPoint outerMiddle = pos + yShift;
0157   GlobalPoint innerMiddle = pos + (-1. * yShift);
0158   if (outerMiddle.perp() < innerMiddle.perp())
0159     std::swap(outerMiddle, innerMiddle);
0160 
0161   auto fround = [&](double in) {
0162     std::stringstream ss;
0163     ss << std::fixed << std::setw(14) << roundIfNear0(in);
0164     return ss.str();
0165   };
0166 
0167   auto fvecround = [&](GlobalPoint vecin) {
0168     std::stringstream ss;
0169     ss << std::fixed << std::setw(14) << roundVecIfNear0(vecin);
0170     return ss.str();
0171   };
0172 
0173   edm::LogVerbatim("MTDDigiGeometryAnalyzer")
0174       << "Det at pos " << fvecround(pos) << " radius " << fround(std::sqrt(pos.x() * pos.x() + pos.y() * pos.y()))
0175       << " has length " << fround(length) << " width " << fround(width) << " thickness " << fround(thickness) << "\n"
0176       << "det center inside bounds? " << tb->inside(det.surface().toLocal(pos)) << "\n"
0177       << "outerMiddle " << fvecround(outerMiddle);
0178   sunitt << det.geographicalId().rawId() << fvecround(pos) << fround(length) << fround(width) << fround(thickness)
0179          << tb->inside(det.surface().toLocal(pos)) << fvecround(outerMiddle) << "\n";
0180 
0181   checkRotation(det);
0182 }
0183 
0184 void MTDDigiGeometryAnalyzer::checkRotation(const GeomDetUnit& det) {
0185   const double eps = std::numeric_limits<float>::epsilon();
0186   static int first = 0;
0187   if (first == 0) {
0188     edm::LogVerbatim("MTDDigiGeometryAnalyzer")
0189         << "numeric_limits<float>::epsilon() " << std::numeric_limits<float>::epsilon();
0190     first = 1;
0191   }
0192 
0193   const Surface::RotationType& rot(det.surface().rotation());
0194   GlobalVector a(rot.xx(), rot.xy(), rot.xz());
0195   GlobalVector b(rot.yx(), rot.yy(), rot.yz());
0196   GlobalVector c(rot.zx(), rot.zy(), rot.zz());
0197   GlobalVector cref = a.cross(b);
0198   GlobalVector aref = b.cross(c);
0199   GlobalVector bref = c.cross(a);
0200   if ((a - aref).mag() > eps || (b - bref).mag() > eps || (c - cref).mag() > eps) {
0201     edm::LogWarning("MTDDigiGeometryAnalyzer")
0202         << " Rotation not good by cross product: " << (a - aref).mag() << ", " << (b - bref).mag() << ", "
0203         << (c - cref).mag() << " for det at pos " << det.surface().position();
0204   }
0205   if (fabs(a.mag() - 1.) > eps || fabs(b.mag() - 1.) > eps || fabs(c.mag() - 1.) > eps) {
0206     edm::LogWarning("MTDDigiGeometryAnalyzer") << " Rotation not good by bector mag: " << (a).mag() << ", " << (b).mag()
0207                                                << ", " << (c).mag() << " for det at pos " << det.surface().position();
0208   }
0209 }
0210 
0211 //define this as a plug-in
0212 DEFINE_FWK_MODULE(MTDDigiGeometryAnalyzer);