Merge pull request #42739 from fabiocos/fc-pixel

MTD geometry: add test for pixel acceptance in module
cms-sw · Sep 11, 2023 · d6ec86c · d6ec86c
2 parents 99f2750 + a390194
commit d6ec86c
Showing 1 changed file with 38 additions and 0 deletions.
diff --git a/Geometry/MTDGeometryBuilder/test/MTDDigiGeometryAnalyzer.cc b/Geometry/MTDGeometryBuilder/test/MTDDigiGeometryAnalyzer.cc
@@ -26,6 +26,8 @@
 
 #include <fstream>
 
+#include "CLHEP/Random/RandFlat.h"
+
 // class declaration
 
 class MTDDigiGeometryAnalyzer : public edm::one::EDAnalyzer<> {
@@ -41,6 +43,7 @@ class MTDDigiGeometryAnalyzer : public edm::one::EDAnalyzer<> {
   void analyseRectangle(const GeomDetUnit& det);
   void checkRotation(const GeomDetUnit& det);
   void checkRectangularMTDTopology(const RectangularMTDTopology&);
+  void checkPixelsAcceptance(const GeomDetUnit& det);
 
   std::stringstream sunitt;
 
@@ -116,6 +119,13 @@ void MTDDigiGeometryAnalyzer::analyze(const edm::Event& iEvent, const edm::Event
     }
   }
 
+  edm::LogInfo("MTDDigiGeometryAnalyzer") << "Acceptance of BTL module:";
+  auto const& btldet = *(dynamic_cast<const MTDGeomDetUnit*>(pDD->detsBTL().front()));
+  checkPixelsAcceptance(btldet);
+  edm::LogInfo("MTDDigiGeometryAnalyzer") << "Acceptance of ETL module:";
+  auto const& etldet = *(dynamic_cast<const MTDGeomDetUnit*>(pDD->detsETL().front()));
+  checkPixelsAcceptance(etldet);
+
   edm::LogInfo("MTDDigiGeometryAnalyzer") << "Additional MTD geometry content:"
                                           << "\n"
                                           << " # dets            = " << pDD->dets().size() << "\n"
@@ -209,5 +219,33 @@ void MTDDigiGeometryAnalyzer::checkRotation(const GeomDetUnit& det) {
   }
 }
 
+void MTDDigiGeometryAnalyzer::checkPixelsAcceptance(const GeomDetUnit& det) {
+  const Bounds& bounds = det.surface().bounds();
+  const RectangularPlaneBounds* tb = dynamic_cast<const RectangularPlaneBounds*>(&bounds);
+  if (tb == nullptr)
+    return;  // not trapezoidal
+
+  double length = tb->length();
+  double width = tb->width();
+  edm::LogVerbatim("MTDDigiGeometryAnalyzer") << "X (width) = " << width << " Y (length) = " << length;
+
+  const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(det.topology());
+  const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
+
+  const size_t maxindex = 100000;
+  size_t inpixel(0);
+  for (size_t index = 0; index < maxindex; index++) {
+    double ax = CLHEP::RandFlat::shoot(-width * 0.5, width * 0.5);
+    double ay = CLHEP::RandFlat::shoot(-length * 0.5, length * 0.5);
+    LocalPoint hit(ax, ay, 0);
+    if (topo.isInPixel(hit)) {
+      inpixel++;
+    }
+  }
+  double acc = (double)inpixel / (double)maxindex;
+  double accerr = std::sqrt(acc * (1. - acc) / (double)maxindex);
+  edm::LogVerbatim("MTDDigiGeometryAnalyzer") << "Acceptance: " << acc << " +/- " << accerr;
+}
+
 //define this as a plug-in
 DEFINE_FWK_MODULE(MTDDigiGeometryAnalyzer);