Phase2-hgx294C Make sure that all rings of scinitllators are included in the definitions of HGCalGeometry

Geometry/HGCalCommonData/interface/HGCalDDDConstants.h

@@ -113,12 +113,21 @@ class HGCalDDDConstants {
   std::pair<int, int> rowColumnWafer(const int wafer) const;
   int sectors() const { return hgpar_->nSectors_; }
   std::pair<int, int> simToReco(int cell, int layer, int mod, bool half) const;
+  int tileCount(int layer, int ring) const;
   bool tileExist(int zside, int layer, int ring, int phi) const {
     int indx = HGCalTileIndex::tileIndex(layer, ring, 0);
     auto itr = hgpar_->tileInfoMap_.find(indx);
     bool ok = (itr == hgpar_->tileInfoMap_.end()) ? false : HGCalTileIndex::tileExist(itr->second.hex, zside, phi);
     return ok;
   }
+  std::pair<int, int> tileRings(int layer) const {
+    if (mode_ == HGCalGeometryMode::TrapezoidFile) {
+      int ll = layer - hgpar_->firstLayer_;
+      if (ll >= 0 && ll < static_cast<int>(hgpar_->tileRingRange_.size()))
+        return hgpar_->tileRingRange_[ll];
+    }
+    return std::make_pair(0, 0);
+  }
   int tileSiPM(int sipm) const { return ((sipm > 0) ? HGCalTypes::SiPMSmall : HGCalTypes::SiPMLarge); }
   bool tileTrapezoid() const {
     return ((mode_ == HGCalGeometryMode::Trapezoid) || (mode_ == HGCalGeometryMode::TrapezoidFile) ||
@@ -212,16 +221,17 @@ class HGCalDDDConstants {
   int wafers() const;
   int wafers(int layer, int type) const;
   int waferToCopy(int wafer) const {
-    return ((wafer >= 0) && (wafer < (int)(hgpar_->waferCopy_.size()))) ? hgpar_->waferCopy_[wafer]
-                                                                        : (int)(hgpar_->waferCopy_.size());
+    return ((wafer >= 0) && (wafer < static_cast<int>(hgpar_->waferCopy_.size())))
+               ? hgpar_->waferCopy_[wafer]
+               : static_cast<int>(hgpar_->waferCopy_.size());
   }
   // wafer transverse thickness classification (2 = coarse, 1 = fine)
   int waferTypeT(int wafer) const {
-    return ((wafer >= 0) && (wafer < (int)(hgpar_->waferTypeT_.size()))) ? hgpar_->waferTypeT_[wafer] : 0;
+    return ((wafer >= 0) && (wafer < static_cast<int>(hgpar_->waferTypeT_.size()))) ? hgpar_->waferTypeT_[wafer] : 0;
   }
   // wafer longitudinal thickness classification (1 = 100um, 2 = 200um, 3=300um)
   int waferTypeL(int wafer) const {
-    return ((wafer >= 0) && (wafer < (int)(hgpar_->waferTypeL_.size()))) ? hgpar_->waferTypeL_[wafer] : 0;
+    return ((wafer >= 0) && (wafer < static_cast<int>(hgpar_->waferTypeL_.size()))) ? hgpar_->waferTypeL_[wafer] : 0;
   }
   int waferType(DetId const& id, bool fromFile = false) const;
   int waferType(int layer, int waferU, int waferV, bool fromFile = false) const;

Geometry/HGCalCommonData/interface/HGCalGeomParameters.h

@@ -99,6 +99,10 @@ class HGCalGeomParameters {
         : half(h), wafer(w), xyz(std::move(p)) {}
   };
 
+  constexpr static int siliconFileEE = 2;
+  constexpr static int siliconFileHE = 3;
+  constexpr static int scintillatorFile = 4;
+
 private:
   void loadGeometryHexagon(const std::map<int, HGCalGeomParameters::layerParameters>& layers,
                            std::vector<HGCalParameters::hgtrform>& trforms,
@@ -143,9 +147,6 @@ class HGCalGeomParameters {
   void resetZero(std::vector<double>&);
 
   constexpr static double tan30deg_ = 0.5773502693;
-  constexpr static int siliconFileEE = 2;
-  constexpr static int siliconFileHE = 3;
-  constexpr static int scintillatorFile = 4;
   HGCalGeomTools geomTools_;
   const double sqrt3_;
   double waferSize_;

Geometry/HGCalCommonData/src/HGCalDDDConstants.cc

@@ -7,6 +7,7 @@
 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 #include "FWCore/Utilities/interface/Exception.h"
+#include "Geometry/HGCalCommonData/interface/HGCalGeomParameters.h"
 #include "Geometry/HGCalCommonData/interface/HGCalGeomTools.h"
 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
 #include "Geometry/HGCalCommonData/interface/HGCalTypes.h"
@@ -15,6 +16,7 @@
 #include "Geometry/HGCalCommonData/interface/HGCalWaferType.h"
 
 #include <algorithm>
+#include <bitset>
 #include <iterator>
 #include <functional>
 #include <numeric>
@@ -588,7 +590,7 @@ bool HGCalDDDConstants::isValidTrap(int layer, int irad, int iphi) const {
   const auto& indx = getIndex(layer, true);
   if (indx.first < 0)
     return false;
-  return ((irad >= hgpar_->iradMinBH_[indx.first]) && (irad <= hgpar_->iradMaxBH_[indx.first]) && (iphi > 0) &&
+  return ((irad >= hgpar_->iradMinBH_[indx.first]) && (irad <= (hgpar_->iradMaxBH_[indx.first] + 1)) && (iphi > 0) &&
           (iphi <= hgpar_->scintCells(layer)));
 }
 
@@ -949,11 +951,14 @@ double HGCalDDDConstants::mouseBite(bool reco) const {
 }
 
 int HGCalDDDConstants::numberCells(bool reco) const {
-  int cells(0);
-  unsigned int nlayer = (reco) ? hgpar_->depth_.size() : hgpar_->layer_.size();
-  for (unsigned k = 0; k < nlayer; ++k) {
-    std::vector<int> ncells = numberCells(((reco) ? hgpar_->depth_[k] : hgpar_->layer_[k]), reco);
-    cells = std::accumulate(ncells.begin(), ncells.end(), cells);
+  int cells =
+      (tileTrapezoid() && (hgpar_->waferMaskMode_ == HGCalGeomParameters::scintillatorFile)) ? tileCount(0, -1) : 0;
+  if (cells == 0) {
+    unsigned int nlayer = (reco) ? hgpar_->depth_.size() : hgpar_->layer_.size();
+    for (unsigned k = 0; k < nlayer; ++k) {
+      std::vector<int> ncells = numberCells(((reco) ? hgpar_->depth_[k] : hgpar_->layer_[k]), reco);
+      cells = std::accumulate(ncells.begin(), ncells.end(), cells);
+    }
   }
   return cells;
 }
@@ -1111,6 +1116,32 @@ std::pair<int, int> HGCalDDDConstants::simToReco(int cell, int lay, int mod, boo
   }
 }
 
+int HGCalDDDConstants::tileCount(int layer, int ring) const {
+  int laymin(layer), laymax(layer), ringmin(ring), ringmax(ring), kount(0);
+  if (layer == 0) {
+    laymin = hgpar_->firstLayer_;
+    laymax = lastLayer(true);
+  }
+  for (int lay = laymin; lay <= laymax; ++lay) {
+    if (ring < 0) {
+      int ll = lay - hgpar_->firstLayer_;
+      ringmin = hgpar_->tileRingRange_[ll].first;
+      ringmax = hgpar_->tileRingRange_[ll].second;
+    }
+    for (int rin = ringmin; rin <= ringmax; ++rin) {
+      int indx = HGCalTileIndex::tileIndex(lay, rin + 1, 0);
+      auto itr = hgpar_->tileInfoMap_.find(indx);
+      if (itr != hgpar_->tileInfoMap_.end()) {
+        for (int k = 0; k < 4; ++k) {
+          std::bitset<24> b(itr->second.hex[k]);
+          kount += b.count();
+        }
+      }
+    }
+  }
+  return (3 * kount);
+}
+
 int HGCalDDDConstants::waferFromCopy(int copy) const {
   const int ncopies = hgpar_->waferCopy_.size();
   int wafer(ncopies);

Geometry/HGCalCommonData/src/HGCalGeomParameters.cc

@@ -2109,7 +2109,7 @@ void HGCalGeomParameters::loadCellTrapezoid(HGCalParameters& php) {
     // Minimum and maximum radius index for each layer
     for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
       php.iradMinBH_.emplace_back(1 + php.tileRingRange_[k].first);
-      php.iradMaxBH_.emplace_back(php.tileRingRange_[k].second);
+      php.iradMaxBH_.emplace_back(1 + php.tileRingRange_[k].second);
 #ifdef EDM_ML_DEBUG
       int kk = php.scintType(php.firstLayer_ + (int)(k));
       edm::LogVerbatim("HGCalGeom") << "New Layer " << k << " Type " << kk << " Low edge " << php.iradMinBH_.back()

Geometry/HGCalGeometry/src/HGCalGeometryLoader.cc

@@ -182,7 +182,7 @@ HGCalGeometry* HGCalGeometryLoader::build(const HGCalTopology& topology) {
   geom->sortDetIds();
 
   if (counter != numberExpected) {
-    if (test) {
+    if (topology.tileTrapezoid()) {
       edm::LogVerbatim("HGCalGeom") << "Inconsistent # of cells: expected " << numberExpected << ":" << numberOfCells
                                     << " , inited " << counter;
     } else {

Geometry/HGCalGeometry/test/HGCalValidScintTest.cc

@@ -0,0 +1,108 @@
+#include <iostream>
+#include <sstream>
+#include <map>
+#include <string>
+#include <vector>
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/one/EDAnalyzer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+
+#include "Geometry/Records/interface/IdealGeometryRecord.h"
+#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
+#include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
+#include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
+#include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
+#include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
+
+class HGCalValidScintTest : public edm::one::EDAnalyzer<> {
+public:
+  explicit HGCalValidScintTest(const edm::ParameterSet&);
+  ~HGCalValidScintTest() override = default;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+  void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
+
+private:
+  const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
+  struct layerInfo {
+    int ringMin, ringMax;
+    double rMin, rMax;
+    layerInfo(int minR = 100, double rMn = 0, int maxR = 0, double rMx = 0)
+        : ringMin(minR), ringMax(maxR), rMin(rMn), rMax(rMx){};
+  };
+};
+
+HGCalValidScintTest::HGCalValidScintTest(const edm::ParameterSet& iC)
+    : geomToken_(esConsumes<HGCalGeometry, IdealGeometryRecord>(edm::ESInputTag{"", "HGCalHEScintillatorSensitive"})) {}
+
+void HGCalValidScintTest::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  descriptions.add("hgcalValidScintTest", desc);
+}
+
+void HGCalValidScintTest::analyze(const edm::Event&, const edm::EventSetup& iSetup) {
+  const auto& geom = &iSetup.getData(geomToken_);
+  DetId::Detector det = DetId::HGCalHSc;
+  edm::LogVerbatim("HGCalGeom") << "Perform test for HGCalHEScintillatorSensitive Detector " << det << " Mode "
+                                << geom->topology().dddConstants().geomMode();
+
+  int firstLayer = geom->topology().dddConstants().firstLayer();
+  int lastLayer = geom->topology().dddConstants().lastLayer(true);
+  const std::vector<DetId>& ids = geom->getValidDetIds();
+  edm::LogVerbatim("HGCalGeom") << "doTest: " << ids.size() << " valid ids for " << geom->cellElement();
+
+  int zside(1);
+  std::map<int, layerInfo> layerMap;
+  for (int layer = firstLayer; layer <= lastLayer; ++layer) {
+    std::vector<std::pair<int, int> > done;
+    for (auto const& id : ids) {
+      HGCScintillatorDetId hid(id);
+      if ((hid.zside() != zside) || (hid.layer() != layer))
+        continue;
+      std::pair<int, int> ring = std::make_pair(hid.ring(), 0);
+      if (std::find(done.begin(), done.end(), ring) != done.end())
+        continue;
+      done.emplace_back(ring);
+      edm::LogVerbatim("HGCalGeom") << "Corners for " << hid;
+
+      const auto cor = geom->getNewCorners(id);
+      std::ostringstream st1;
+      for (unsigned int k = 0; k < cor.size(); ++k)
+        st1 << " [" << k << "] " << std::setprecision(4) << cor[k];
+      edm::LogVerbatim("HGCalGeom") << st1.str();
+
+      double r = cor[0].perp();
+      auto itr = layerMap.find(layer);
+      if (itr == layerMap.end()) {
+        layerInfo info;
+        info.ringMin = info.ringMax = ring.first;
+        info.rMin = info.rMax = r;
+        layerMap[layer] = info;
+      } else {
+        layerInfo info = itr->second;
+        if (info.ringMin > ring.first) {
+          info.ringMin = ring.first;
+          info.rMin = r;
+        } else if (info.ringMax < ring.first) {
+          info.ringMax = ring.first;
+          info.rMax = r;
+        }
+        layerMap[layer] = info;
+      }
+    }
+  }
+  edm::LogVerbatim("HGCalGeom") << "\n\nSummary of " << layerMap.size() << " Scintillator layers";
+  for (auto itr = layerMap.begin(); itr != layerMap.end(); ++itr)
+    edm::LogVerbatim("HGCalGeom") << "Layer " << itr->first << " lowest Ring " << (itr->second).ringMin << ":"
+                                  << (itr->second).rMin << " largest Ring " << (itr->second).ringMax << ":"
+                                  << (itr->second).rMax;
+}
+
+//define this as a plug-in
+DEFINE_FWK_MODULE(HGCalValidScintTest);

Geometry/HGCalGeometry/test/python/testHGCalValidScintTest_cfg.py

@@ -0,0 +1,39 @@
+import FWCore.ParameterSet.Config as cms
+
+from Configuration.Eras.Era_Phase2C11_cff import Phase2C11
+process = cms.Process('PROD',Phase2C11)
+
+process.load("SimGeneral.HepPDTESSource.pdt_cfi")
+process.load("Configuration.Geometry.GeometryExtended2026D86Reco_cff")
+process.load('FWCore.MessageService.MessageLogger_cfi')
+
+if hasattr(process,'MessageLogger'):
+    process.MessageLogger.HGCalGeom=dict()
+
+process.load("IOMC.RandomEngine.IOMC_cff")
+process.RandomNumberGeneratorService.generator.initialSeed = 456789
+
+process.source = cms.Source("EmptySource")
+
+process.generator = cms.EDProducer("FlatRandomEGunProducer",
+    PGunParameters = cms.PSet(
+        PartID = cms.vint32(14),
+        MinEta = cms.double(-3.5),
+        MaxEta = cms.double(3.5),
+        MinPhi = cms.double(-3.14159265359),
+        MaxPhi = cms.double(3.14159265359),
+        MinE   = cms.double(9.99),
+        MaxE   = cms.double(10.01)
+    ),
+    AddAntiParticle = cms.bool(False),
+    Verbosity       = cms.untracked.int32(0),
+    firstRun        = cms.untracked.uint32(1)
+)
+
+process.maxEvents = cms.untracked.PSet(
+    input = cms.untracked.int32(1)
+)
+
+process.load("Geometry.HGCalGeometry.hgcalValidScintTest_cfi")
+
+process.p1 = cms.Path(process.generator*process.hgcalValidScintTest)