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 using namespace cms_rounding;
0032
0033 // class declaration
0034
0035 class MTDDigiGeometryAnalyzer : public edm::one::EDAnalyzer<> {
0036 public:
0037   explicit MTDDigiGeometryAnalyzer(const edm::ParameterSet&);
0038   ~MTDDigiGeometryAnalyzer() override = default;
0039
0040   void beginJob() override {}
0041   void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0042   void endJob() override {}
0043
0044 private:
0045   inline std::string fround(const double in, const size_t prec) const {
0046     std::stringstream ss;
0047     ss << std::setprecision(prec) << std::fixed << std::setw(14) << roundIfNear0(in);
0048     return ss.str();
0049   }
0050
0051   inline std::string fvecround(const auto& vecin, const size_t prec) const {
0052     std::stringstream ss;
0053     ss << std::setprecision(prec) << std::fixed << std::setw(14) << roundVecIfNear0(vecin);
0054     return ss.str();
0055   }
0056
0057   void checkRectangularMTDTopology(const RectangularMTDTopology&);
0058   void checkPixelsAcceptance(const GeomDetUnit& det);
0059
0060   std::stringstream sunitt_;
0061
0062   edm::ESGetToken<MTDGeometry, MTDDigiGeometryRecord> mtdgeoToken_;
0063 };
0064
0065 MTDDigiGeometryAnalyzer::MTDDigiGeometryAnalyzer(const edm::ParameterSet& iConfig) {
0066   mtdgeoToken_ = esConsumes<MTDGeometry, MTDDigiGeometryRecord>();
0067 }
0068
0069 // ------------ method called to produce the data  ------------
0070 void MTDDigiGeometryAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0071   //
0072   // get the MTDGeometry
0073   //
0074   auto pDD = iSetup.getTransientHandle(mtdgeoToken_);
0075   edm::LogVerbatim("MTDDigiGeometryAnalyzer")
0076       << "MTDGeometry:\n"
0077       << " # detectors  = " << pDD->detUnits().size() << "\n"
0078       << " # types      = " << pDD->detTypes().size() << "\n"
0079       << " # BTL dets   = " << pDD->detsBTL().size() << "\n"
0080       << " # ETL dets   = " << pDD->detsETL().size() << "\n"
0081       << " # layers " << pDD->geomDetSubDetector(1) << "   = " << pDD->numberOfLayers(1) << "\n"
0082       << " # layers " << pDD->geomDetSubDetector(2) << "   = " << pDD->numberOfLayers(2) << "\n"
0083       << " # dets       = " << pDD->dets().size() << "\n"
0084       << " # detUnitIds = " << pDD->detUnitIds().size() << "\n"
0085       << " # detIds     = " << pDD->detIds().size() << "\n";
0086
0087   sunitt_ << "MTDGeometry:\n"
0088           << " # detectors  = " << pDD->detUnits().size() << "\n"
0089           << " # types      = " << pDD->detTypes().size() << "\n"
0090           << " # BTL dets   = " << pDD->detsBTL().size() << "\n"
0091           << " # ETL dets   = " << pDD->detsETL().size() << "\n"
0092           << " # layers " << pDD->geomDetSubDetector(1) << "   = " << pDD->numberOfLayers(1) << "\n"
0093           << " # layers " << pDD->geomDetSubDetector(2) << "   = " << pDD->numberOfLayers(2) << "\n"
0094           << " # dets       = " << pDD->dets().size() << "\n"
0095           << " # detUnitIds = " << pDD->detUnitIds().size() << "\n"
0096           << " # detIds     = " << pDD->detIds().size() << "\n";
0097
0098   for (auto const& it : pDD->detTypes()) {
0099     if (dynamic_cast<const MTDGeomDetType*>((it)) != nullptr) {
0100       const PixelTopology& p = (dynamic_cast<const MTDGeomDetType*>((it)))->specificTopology();
0101       const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(p);
0102       auto pitchval = topo.pitch();
0103       edm::LogVerbatim("MTDDigiGeometryAnalyzer")
0104           << "\n Subdetector " << it->subDetector() << " MTD Det " << it->name() << "\n"
0105           << " Rows     " << topo.nrows() << " Columns " << topo.ncolumns() << " ROCS X   " << topo.rocsX()
0106           << " ROCS Y  " << topo.rocsY() << " Rows/ROC " << topo.rowsperroc() << " Cols/ROC " << topo.colsperroc()
0107           << " Pitch X " << fround(pitchval.first, 4) << " Pitch Y " << fround(pitchval.second, 4)
0108           << " Sensor Interpad X " << fround(topo.gapxInterpad(), 4) << " Sensor Interpad Y "
0109           << fround(topo.gapyInterpad(), 4) << " Sensor Border X " << fround(topo.gapxBorder(), 4)
0110           << " Sensor Border Y " << fround(topo.gapyBorder(), 4) << "\n";
0111       sunitt_ << "\n Subdetector " << it->subDetector() << " MTD Det " << it->name() << "\n"
0112               << " Rows     " << topo.nrows() << " Columns " << topo.ncolumns() << " ROCS X   " << topo.rocsX()
0113               << " ROCS Y  " << topo.rocsY() << " Rows/ROC " << topo.rowsperroc() << " Cols/ROC " << topo.colsperroc()
0114               << " Pitch X " << fround(pitchval.first, 2) << " Pitch Y " << fround(pitchval.second, 2)
0115               << " Sensor Interpad X " << fround(topo.gapxInterpad(), 2) << " Sensor Interpad Y "
0116               << fround(topo.gapyInterpad(), 2) << " Sensor Border X " << fround(topo.gapxBorder(), 2)
0117               << " Sensor Border Y " << fround(topo.gapyBorder(), 2) << "\n";
0118       checkRectangularMTDTopology(topo);
0119     }
0120   }
0121
0122   edm::LogVerbatim("MTDDigiGeometryAnalyzer") << "\nAcceptance of BTL module:";
0123   sunitt_ << "\nAcceptance of BTL module:";
0124   auto const& btldet = *(dynamic_cast<const MTDGeomDetUnit*>(pDD->detsBTL().front()));
0125   checkPixelsAcceptance(btldet);
0126   edm::LogVerbatim("MTDDigiGeometryAnalyzer") << "\nAcceptance of ETL module:";
0127   sunitt_ << "\nAcceptance of ETL module:";
0128   auto const& etldet = *(dynamic_cast<const MTDGeomDetUnit*>(pDD->detsETL().front()));
0129   checkPixelsAcceptance(etldet);
0130
0131   edm::LogVerbatim("MTDUnitTest") << sunitt_.str();
0132 }
0133
0134 void MTDDigiGeometryAnalyzer::checkRectangularMTDTopology(const RectangularMTDTopology& topo) {
0135   edm::LogVerbatim("MTDDigiGeometryAnalyzer") << "Pixel center location:\n";
0136   sunitt_ << "Pixel center location:\n";
0137   LocalPoint center(0, 0, 0);
0138   for (int r = 0; r < topo.nrows(); r++) {
0139     for (int c = 0; c < topo.ncolumns(); c++) {
0140       edm::LogVerbatim("MTDDigiGeometryAnalyzer") << std::setw(7) << r << std::setw(7) << c << " "
0141                                                   << fvecround(topo.pixelToModuleLocalPoint(center, r, c), 4) << "\n";
0142       sunitt_ << std::setw(7) << r << std::setw(7) << c << " "
0143               << fvecround(topo.pixelToModuleLocalPoint(center, r, c), 2) << "\n";
0144     }
0145   }
0146 }
0147
0148 void MTDDigiGeometryAnalyzer::checkPixelsAcceptance(const GeomDetUnit& det) {
0149   const Bounds& bounds = det.surface().bounds();
0150   const RectangularPlaneBounds* tb = dynamic_cast<const RectangularPlaneBounds*>(&bounds);
0151   if (tb == nullptr)
0152     return;  // not trapezoidal
0153
0154   double length = tb->length();
0155   double width = tb->width();
0156   edm::LogVerbatim("MTDDigiGeometryAnalyzer")
0157       << " X (width) = " << fround(width, 4) << " Y (length) = " << fround(length, 4);
0158   sunitt_ << " X (width) = " << fround(width, 2) << " Y (length) = " << fround(length, 2);
0159
0160   const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(det.topology());
0161   const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0162
0163   const size_t maxindex = 100000;
0164   size_t inpixel(0);
0165   for (size_t index = 0; index < maxindex; index++) {
0166     double ax = CLHEP::RandFlat::shoot(-width * 0.5, width * 0.5);
0167     double ay = CLHEP::RandFlat::shoot(-length * 0.5, length * 0.5);
0168     LocalPoint hit(ax, ay, 0);
0169     auto const indici = topo.pixelIndex(hit);
0170     assert(indici.first < topo.nrows() && indici.second < topo.ncolumns());  // sanity check on the index definition
0171     if (topo.isInPixel(hit)) {
0172       inpixel++;
0173     }
0174   }
0175   double acc = (double)inpixel / (double)maxindex;
0176   double accerr = std::sqrt(acc * (1. - acc) / (double)maxindex);
0177   edm::LogVerbatim("MTDDigiGeometryAnalyzer") << " Acceptance: " << fround(acc, 3) << " +/- " << fround(accerr, 3);
0178   sunitt_ << " Acceptance: " << fround(acc, 3) << " +/- " << fround(accerr, 3);
0179 }
0180
0181 //define this as a plug-in
0182 DEFINE_FWK_MODULE(MTDDigiGeometryAnalyzer);