Phase2-hgx291 Bug fix to address rotated layers in D86

Geometry/HGCalCommonData/interface/HGCalDDDConstants.h

@@ -30,7 +30,7 @@ class HGCalDDDConstants {
   ~HGCalDDDConstants();
 
   std::pair<int, int> assignCell(float x, float y, int lay, int subSec, bool reco) const;
-  std::array<int, 5> assignCellHex(float x, float y, int lay, bool reco) const;
+  std::array<int, 5> assignCellHex(float x, float y, int lay, bool reco, bool extend = false, bool debug = false) const;
   std::array<int, 3> assignCellTrap(float x, float y, float z, int lay, bool reco) const;
   std::pair<double, double> cellEtaPhiTrap(int type, int irad) const;
   bool cellInLayer(int waferU, int waferV, int cellU, int cellV, int lay, bool reco) const;
@@ -136,6 +136,7 @@ class HGCalDDDConstants {
                          int& cellV,
                          int& celltype,
                          double& wt,
+                         bool extend = false,
                          bool debug = false) const;
   bool waferHexagon6() const {
     return ((mode_ == HGCalGeometryMode::Hexagon) || (mode_ == HGCalGeometryMode::HexagonFull));
@@ -212,13 +213,15 @@ class HGCalDDDConstants {
               const double& cellR,
               const std::vector<double>& posX,
               const std::vector<double>& posY) const;
-  void cellHex(double xloc, double yloc, int cellType, int& cellU, int& cellV, bool debug = false) const;
+  void cellHex(
+      double xloc, double yloc, int cellType, int& cellU, int& cellV, bool extend = false, bool debug = false) const;
   std::pair<int, float> getIndex(int lay, bool reco) const;
   int layerFromIndex(int index, bool reco) const;
   bool isValidCell(int layindex, int wafer, int cell) const;
   bool isValidCell8(int lay, int waferU, int waferV, int cellU, int cellV, int type) const;
   int32_t waferIndex(int wafer, int index) const;
   bool waferInLayerTest(int wafer, int lay, bool full) const;
+  std::pair<double, double> waferPositionNoRot(int lay, int waferU, int waferV, bool reco, bool debug = false) const;
   std::pair<double, double> waferPosition(int waferU, int waferV, bool reco) const;
 
   HGCalGeomTools geomTools_;
@@ -229,9 +232,10 @@ class HGCalDDDConstants {
   const HGCalParameters* hgpar_;
   constexpr static double tan30deg_ = 0.5773502693;
   const double sqrt3_;
+  const HGCalGeometryMode::GeometryMode mode_;
+  const bool fullAndPart_;
   double rmax_, hexside_;
-  HGCalGeometryMode::GeometryMode mode_;
-  bool fullAndPart_;
+  double rmaxT_, hexsideT_;
   int32_t tot_wafers_, modHalf_;
   std::array<uint32_t, 2> tot_layers_;
   Simrecovecs max_modules_layer_;

Geometry/HGCalCommonData/src/HGCalParametersFromDD.cc

@@ -502,6 +502,7 @@ void HGCalParametersFromDD::getCellPosition(HGCalParameters& php, int type) {
     php.cellCoarseIndex_ = cellIndex;
   else
     php.cellFineIndex_ = cellIndex;
+
 #ifdef EDM_ML_DEBUG
   if (type == 1) {
     edm::LogVerbatim("HGCalGeom") << "CellPosition for  type " << type << " for " << php.cellCoarseX_.size()

Geometry/HGCalGeometry/interface/HGCalGeometry.h

@@ -93,6 +93,7 @@ class HGCalGeometry final : public CaloSubdetectorGeometry {
 
   // Get closest cell, etc...
   DetId getClosestCell(const GlobalPoint& r) const override;
+  DetId getClosestCellHex(const GlobalPoint& r, bool extend) const;
 
   /** \brief Get a list of all cells within a dR of the given cell
 

Geometry/HGCalGeometry/src/HGCalGeometry.cc

@@ -219,12 +219,17 @@ GlobalPoint HGCalGeometry::getPosition(const DetId& detid) const {
                                     << " Global " << glob;
 #endif
     } else {
-      xy = m_topology.dddConstants().locateCell(id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, false, true);
-      const HepGeom::Point3D<float> lcoord(xy.first, xy.second, 0);
-      glob = m_cellVec[cellIndex].getPosition(lcoord);
 #ifdef EDM_ML_DEBUG
-      edm::LogVerbatim("HGCalGeom") << "getPositionWafer:: index " << cellIndex << " Local " << lcoord.x() << ":"
-                                    << lcoord.y() << " ID " << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":"
+      edm::LogVerbatim("HGCalGeom") << "getPosition for " << HGCSiliconDetId(detid) << " Layer " << id.iLay << " Wafer "
+                                    << id.iSec1 << ":" << id.iSec2 << " Cell " << id.iCell1 << ":" << id.iCell2;
+#endif
+      xy = m_topology.dddConstants().locateCell(id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, true, true);
+      double xx = id.zSide * xy.first;
+      double zz = id.zSide * m_topology.dddConstants().waferZ(id.iLay, true);
+      glob = GlobalPoint(xx, xy.second, zz);
+#ifdef EDM_ML_DEBUG
+      edm::LogVerbatim("HGCalGeom") << "getPositionWafer:: index " << cellIndex << " Local " << xy.first << ":"
+                                    << xy.second << " ID " << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":"
                                     << id.iCell1 << ":" << id.iCell2 << " Global " << glob;
 #endif
     }
@@ -497,7 +502,51 @@ DetId HGCalGeometry::getClosestCell(const GlobalPoint& r) const {
       id.iType = kxy[2];
     } else {
       id.iLay = m_topology.dddConstants().getLayer(r.z(), true);
-      const auto& kxy = m_topology.dddConstants().assignCellHex(local.x(), local.y(), id.iLay, true);
+#ifdef EDM_ML_DEBUG
+      edm::LogVerbatim("HGCalGeom") << "ZZ " << r.z() << " Layer " << id.iLay << " Global " << r << "  Local " << local;
+#endif
+      const auto& kxy = m_topology.dddConstants().assignCellHex(local.x(), local.y(), id.iLay, false, true);
+      id.iSec1 = kxy[0];
+      id.iSec2 = kxy[1];
+      id.iType = kxy[2];
+      id.iCell1 = kxy[3];
+      id.iCell2 = kxy[4];
+    }
+#ifdef EDM_ML_DEBUG
+    edm::LogVerbatim("HGCalGeom") << "getClosestCell: local " << local << " Id " << id.det << ":" << id.zSide << ":"
+                                  << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":" << id.iType << ":"
+                                  << id.iCell1 << ":" << id.iCell2;
+#endif
+
+    //check if returned cell is valid
+    if (id.iCell1 >= 0)
+      return m_topology.encode(id);
+  }
+
+  //if not valid or out of bounds return a null DetId
+  return DetId();
+}
+
+DetId HGCalGeometry::getClosestCellHex(const GlobalPoint& r, bool extend) const {
+  unsigned int cellIndex = getClosestCellIndex(r);
+  if (cellIndex < m_cellVec.size() && m_det != DetId::HGCalHSc) {
+    HGCalTopology::DecodedDetId id = m_topology.decode(m_validGeomIds[cellIndex]);
+    if (id.det == 0)
+      id.det = static_cast<int>(m_topology.detector());
+    HepGeom::Point3D<float> local;
+    if (r.z() > 0) {
+      local = HepGeom::Point3D<float>(r.x(), r.y(), 0);
+      id.zSide = 1;
+    } else {
+      local = HepGeom::Point3D<float>(-r.x(), r.y(), 0);
+      id.zSide = -1;
+    }
+    if (m_topology.waferHexagon8()) {
+      id.iLay = m_topology.dddConstants().getLayer(r.z(), true);
+#ifdef EDM_ML_DEBUG
+      edm::LogVerbatim("HGCalGeom") << "ZZ " << r.z() << " Layer " << id.iLay << " Global " << r << "  Local " << local;
+#endif
+      const auto& kxy = m_topology.dddConstants().assignCellHex(local.x(), local.y(), id.iLay, extend, true);
       id.iSec1 = kxy[0];
       id.iSec2 = kxy[1];
       id.iType = kxy[2];

Geometry/HGCalGeometry/src/HGCalGeometryLoader.cc

@@ -137,7 +137,7 @@ HGCalGeometry* HGCalGeometryLoader::build(const HGCalTopology& topology) {
           int type = topology.dddConstants().getTypeHex(layer, u, v);
           DetId detId = (topology.isHFNose() ? (DetId)(HFNoseDetId(zside, type, layer, u, v, 0, 0))
                                              : (DetId)(HGCSiliconDetId(det, zside, type, layer, u, v, 0, 0)));
-          const auto& w = topology.dddConstants().waferPosition(layer, u, v, true);
+          const auto& w = topology.dddConstants().waferPosition(layer, u, v, true, true);
           double xx = (zside > 0) ? w.first : -w.first;
           CLHEP::Hep3Vector h3v(xx, w.second, mytr.h3v.z());
           const HepGeom::Transform3D ht3d(mytr.hr, h3v);

Geometry/HGCalGeometry/test/HGCalGeometryRotCheck.cc

@@ -0,0 +1,87 @@
+#include <cmath>
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "CommonTools/UtilAlgos/interface/TFileService.h"
+
+#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/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/HGCSiliconDetId.h"
+
+class HGCalGeometryRotCheck : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
+public:
+  explicit HGCalGeometryRotCheck(const edm::ParameterSet&);
+  ~HGCalGeometryRotCheck() override = default;
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+  void beginRun(edm::Run const&, edm::EventSetup const&) override;
+  void analyze(edm::Event const& iEvent, edm::EventSetup const&) override {}
+  void endRun(edm::Run const&, edm::EventSetup const&) override {}
+
+private:
+  const std::string nameDetector_;
+  const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
+  const std::vector<int> layers_;
+};
+
+HGCalGeometryRotCheck::HGCalGeometryRotCheck(const edm::ParameterSet& iC)
+    : nameDetector_(iC.getParameter<std::string>("detectorName")),
+      geomToken_(esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(
+          edm::ESInputTag{"", nameDetector_})),
+      layers_(iC.getParameter<std::vector<int>>("layers")) {}
+
+void HGCalGeometryRotCheck::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  std::vector<int> layer = {27, 28, 29, 30, 31, 32, 33};
+  desc.add<std::string>("detectorName", "HGCalHESiliconSensitive");
+  desc.add<std::vector<int>>("layers", layer);
+  descriptions.add("hgcalGeometryRotCheck", desc);
+}
+
+void HGCalGeometryRotCheck::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
+  const auto& geomR = iSetup.getData(geomToken_);
+  const HGCalGeometry* geom = &geomR;
+  int layerF = *(layers_.begin());
+  int layerL = *(--layers_.end());
+  int layerOff = geom->topology().dddConstants().getLayerOffset();
+  edm::LogVerbatim("HGCalGeom") << nameDetector_ << " with layers in the range " << layerF << ":" << layerL
+                                << " Offset " << layerOff;
+
+  auto rrF = geom->topology().dddConstants().rangeRLayer(layerF - layerOff, true);
+  auto rrE = geom->topology().dddConstants().rangeRLayer(layerL - layerOff, true);
+  edm::LogVerbatim("HGCalGeom") << " RFront " << rrF.first << ":" << rrF.second << " RBack " << rrE.first << ":"
+                                << rrE.second;
+  double r = rrE.first + 5.0;
+  const int nPhi = 10;
+  while (r <= rrF.second) {
+    for (int k = 0; k < nPhi; ++k) {
+      double phi = 2.0 * k * M_PI / nPhi;
+      for (auto lay : layers_) {
+        double zz = geom->topology().dddConstants().waferZ(lay - layerOff, true);
+        GlobalPoint global1(r * cos(phi), r * sin(phi), zz);
+        DetId id = geom->getClosestCellHex(global1, true);
+        HGCSiliconDetId detId = HGCSiliconDetId(id);
+        GlobalPoint global2 = geom->getPosition(id);
+        double dx = global1.x() - global2.x();
+        double dy = global1.y() - global2.y();
+        double dR = std::sqrt(dx * dx + dy * dy);
+        edm::LogVerbatim("HGCalGeom") << "Layer: " << lay << " ID " << detId << " I/P " << global1 << " O/P " << global2
+                                      << " dR " << dR;
+      }
+    }
+    r += 100.0;
+  }
+}
+
+//define this as a plug-in
+#include "FWCore/Framework/interface/MakerMacros.h"
+DEFINE_FWK_MODULE(HGCalGeometryRotCheck);

Geometry/HGCalGeometry/test/HGCalGeometryRotTest.cc

@@ -0,0 +1,90 @@
+#include <cmath>
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "CommonTools/UtilAlgos/interface/TFileService.h"
+
+#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/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/HGCSiliconDetId.h"
+
+class HGCalGeometryRotTest : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
+public:
+  explicit HGCalGeometryRotTest(const edm::ParameterSet&);
+  ~HGCalGeometryRotTest() override = default;
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+  void beginRun(edm::Run const&, edm::EventSetup const&) override;
+  void analyze(edm::Event const& iEvent, edm::EventSetup const&) override {}
+  void endRun(edm::Run const&, edm::EventSetup const&) override {}
+
+private:
+  const std::string nameDetector_;
+  const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
+  const std::vector<int> layers_;
+  const std::vector<int> waferU_, waferV_, cellU_, cellV_, types_;
+};
+
+HGCalGeometryRotTest::HGCalGeometryRotTest(const edm::ParameterSet& iC)
+    : nameDetector_(iC.getParameter<std::string>("detectorName")),
+      geomToken_(esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(
+          edm::ESInputTag{"", nameDetector_})),
+      layers_(iC.getParameter<std::vector<int>>("layers")),
+      waferU_(iC.getParameter<std::vector<int>>("waferUs")),
+      waferV_(iC.getParameter<std::vector<int>>("waferVs")),
+      cellU_(iC.getParameter<std::vector<int>>("cellUs")),
+      cellV_(iC.getParameter<std::vector<int>>("cellVs")),
+      types_(iC.getParameter<std::vector<int>>("types")) {}
+
+void HGCalGeometryRotTest::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  std::vector<int> layer = {27, 28, 29, 30, 31, 32};
+  std::vector<int> waferU = {-2, -3, 1, -8, 2, 8};
+  std::vector<int> waferV = {0, -2, 3, 0, 9, 0};
+  std::vector<int> cellU = {16, 4, 8, 11, 11, 5};
+  std::vector<int> cellV = {20, 10, 17, 13, 9, 2};
+  std::vector<int> type = {0, 0, 0, 2, 2, 2};
+  desc.add<std::string>("detectorName", "HGCalHESiliconSensitive");
+  desc.add<std::vector<int>>("layers", layer);
+  desc.add<std::vector<int>>("waferUs", waferU);
+  desc.add<std::vector<int>>("waferVs", waferV);
+  desc.add<std::vector<int>>("cellUs", cellU);
+  desc.add<std::vector<int>>("cellVs", cellV);
+  desc.add<std::vector<int>>("types", type);
+  descriptions.add("hgcalGeometryRotTest", desc);
+}
+
+void HGCalGeometryRotTest::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
+  const auto& geomR = iSetup.getData(geomToken_);
+  const HGCalGeometry* geom = &geomR;
+  DetId::Detector det = (nameDetector_ == "HGCalEESensitive") ? DetId::HGCalEE : DetId::HGCalHSi;
+  int layerF = *(layers_.begin());
+  int layerL = *(--layers_.end());
+  int layerOff = geom->topology().dddConstants().getLayerOffset();
+  edm::LogVerbatim("HGCalGeom") << nameDetector_ << " with layers in the range " << layerF << ":" << layerL
+                                << " Offset " << layerOff << " and for " << waferU_.size() << " wafers and cells";
+
+  for (unsigned int k = 0; k < waferU_.size(); ++k) {
+    for (auto lay : layers_) {
+      HGCSiliconDetId detId(det, 1, types_[k], lay - layerOff, waferU_[k], waferV_[k], cellU_[k], cellV_[k]);
+      GlobalPoint global = geom->getPosition(DetId(detId));
+      auto xy = geom->topology().dddConstants().waferPosition(lay - layerOff, waferU_[k], waferV_[k], true);
+      edm::LogVerbatim("HGCalGeom") << "Layer: " << lay << " U " << waferU_[k] << " V " << waferV_[k] << " Position ("
+                                    << xy.first << ", " << xy.second << ")";
+      edm::LogVerbatim("HGCalGeom") << detId << " Position " << global;
+    }
+  }
+}
+
+//define this as a plug-in
+#include "FWCore/Framework/interface/MakerMacros.h"
+DEFINE_FWK_MODULE(HGCalGeometryRotTest);

Geometry/HGCalGeometry/test/python/testHGCalGeometryRotCheck_cfg.py

@@ -0,0 +1,39 @@
+import FWCore.ParameterSet.Config as cms
+
+from Configuration.Eras.Era_Phase2C11I13M9_cff import Phase2C11I13M9
+process = cms.Process('PROD',Phase2C11I13M9)
+process.load('Configuration.Geometry.GeometryExtended2026D86_cff')
+process.load('Configuration.Geometry.GeometryExtended2026D86Reco_cff')
+
+process.load("SimGeneral.HepPDTESSource.pdt_cfi")
+process.load('Geometry.HGCalGeometry.hgcalGeometryRotCheck_cfi')
+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.p1 = cms.Path(process.generator*process.hgcalGeometryRotCheck)