Add a new method for HGCalGeometry and add a test for corners

Geometry/HGCalCommonData/python/testGeometryExtended_cff.py

@@ -1,6 +1,6 @@
 import FWCore.ParameterSet.Config as cms
 
-from Geometry.HGCalCommonData.testHGCalXML_cfi import *
+from Geometry.CMSCommonData.cmsExtendedGeometry2023D28XML_cfi import *
 from Geometry.TrackerNumberingBuilder.trackerNumberingGeometry_cfi import *
 from Geometry.HcalCommonData.hcalDDConstants_cff import *
 from Geometry.HGCalCommonData.hgcalParametersInitialization_cfi import *

Geometry/HGCalCommonData/test/dumpExtended2023HGCalGeometry_cfg.py

@@ -2,16 +2,12 @@
 
 process = cms.Process("GEODUMP")
 process.load("Geometry.HGCalCommonData.testGeometryExtended_cff")
+process.load('FWCore.MessageService.MessageLogger_cfi')
 
-process.load("FWCore.MessageService.MessageLogger_cfi")
-
-process.MessageLogger = cms.Service("MessageLogger",
-                                    debugModules = cms.untracked.vstring('*'),
-                                    destinations = cms.untracked.vstring('cout'),
-                                    cout = cms.untracked.PSet( threshold = cms.untracked.string('DEBUG'),
-                                                               noLineBreaks = cms.untracked.bool(True)
-                                                               )
-                                    )
+if 'MessageLogger' in process.__dict__:
+    process.MessageLogger.categories.append('G4cerr')
+    process.MessageLogger.categories.append('G4cout')
+    process.MessageLogger.categories.append('HGCalGeom')
 
 process.source = cms.Source("EmptySource")
 
@@ -22,6 +18,7 @@
                             level = cms.untracked.int32(14)
                             ))
 
-process.dump = cms.EDAnalyzer("DumpSimGeometry")
+process.dump = cms.EDAnalyzer("DumpSimGeometry",
+                              outputFileName = cms.untracked.string('CMS2023D28.root'))
 
 process.p = cms.Path(process.dump)

Geometry/HGCalGeometry/interface/HGCalGeometry.h

@@ -78,6 +78,7 @@ class HGCalGeometry final: public CaloSubdetectorGeometry {
 
   /// Returns the corner points of this cell's volume.
   CornersVec getCorners( const DetId& id ) const; 
+  CornersVec get8Corners( const DetId& id ) const; 
 
   // avoid sorting set in base class  
   const std::vector<DetId>& getValidDetIds( DetId::Detector det = DetId::Detector(0), int subdet = 0) const override { return m_validIds; }

Geometry/HGCalGeometry/src/HGCalGeometry.cc

@@ -320,7 +320,7 @@ HGCalGeometry::CornersVec HGCalGeometry::getCorners(const DetId& id) const {
       static const int signx[] = {-1,-1,1,1,-1,-1,1,1};
       static const int signy[] = {-1,1,1,-1,-1,1,1,-1};
       static const int signz[] = {-1,-1,-1,-1,1,1,1,1};
-      for (unsigned int i = 0; i != ncorner; ++i) {
+      for (unsigned int i = 0; i < ncorner; ++i) {
 	const HepGeom::Point3D<float> lcoord(xy.first+signx[i]*dx,xy.second+signy[i]*dx,signz[i]*dz);
 	co[i] = m_cellVec2[cellIndex].getPosition(lcoord);
       }
@@ -331,7 +331,7 @@ HGCalGeometry::CornersVec HGCalGeometry::getCorners(const DetId& id) const {
       static const int signx[] = {0,-1,-1,0,1,1,0,-1,-1,0,1,1};
       static const int signy[] = {-2,-1,1,2,1,-1,-2,-1,1,2,1,-1};
       static const int signz[] = {-1,-1,-1,-1,-1,-1,1,1,1,1,1,1};
-      for (unsigned int i = 0; i != ncorner; ++i) {
+      for (unsigned int i = 0; i < ncorner; ++i) {
 	const HepGeom::Point3D<float> lcoord(xy.first+signx[i]*dx,xy.second+signy[i]*dy,signz[i]*dz);
 	co[i] = m_cellVec[cellIndex].getPosition(lcoord);
       }
@@ -340,6 +340,43 @@ HGCalGeometry::CornersVec HGCalGeometry::getCorners(const DetId& id) const {
   return co;
 }
 
+HGCalGeometry::CornersVec HGCalGeometry::get8Corners(const DetId& id) const {
+
+  unsigned int ncorner = FlatTrd::ncorner_;
+  HGCalGeometry::CornersVec co (ncorner, GlobalPoint(0,0,0));
+  unsigned int cellIndex =  indexFor(id);
+  if ((cellIndex <  m_cellVec.size() && m_det != DetId::HGCalHSc) ||
+      (cellIndex <  m_cellVec2.size() && m_det == DetId::HGCalHSc)) {
+    HGCalTopology::DecodedDetId id_ = m_topology.decode(id);
+    std::pair<float,float> xy;
+    if ((mode_ == HGCalGeometryMode::Hexagon) ||
+	(mode_ == HGCalGeometryMode::HexagonFull)) {
+      xy = m_topology.dddConstants().locateCellHex(id_.iCell1,id_.iSec1,true);
+    } else if (mode_ == HGCalGeometryMode::Trapezoid) {
+      xy = m_topology.dddConstants().locateCellTrap(id_.iLay,id_.iSec1,
+						    id_.iCell1,true);
+    } else {
+      xy = m_topology.dddConstants().locateCell(id_.iLay,id_.iSec1,id_.iSec2,
+						id_.iCell1,id_.iCell2,true,false);
+    }
+    float dx = ((m_det == DetId::HGCalHSc) ? 0.5*m_cellVec2[cellIndex].param()[11] :
+		m_cellVec[cellIndex].param()[1]);
+    float dz = ((m_det == DetId::HGCalHSc) ? m_cellVec2[cellIndex].param()[0] :
+		m_cellVec[cellIndex].param()[0]);
+    static const int signx[] = {-1,-1,1,1,-1,-1,1,1};
+    static const int signy[] = {-1,1,1,-1,-1,1,1,-1};
+    static const int signz[] = {-1,-1,-1,-1,1,1,1,1};
+    for (unsigned int i = 0; i < ncorner; ++i) {
+      const HepGeom::Point3D<float> lcoord(xy.first+signx[i]*dx,
+					   xy.second+signy[i]*dx,signz[i]*dz);
+      co[i] = ((m_det == DetId::HGCalHSc) ? 
+	       (m_cellVec2[cellIndex].getPosition(lcoord)) :
+	       (m_cellVec[cellIndex].getPosition(lcoord)));
+    }
+  }
+  return co;
+}
+
 DetId HGCalGeometry::getClosestCell(const GlobalPoint& r) const {
   unsigned int cellIndex = getClosestCellIndex(r);
   if ((cellIndex <  m_cellVec.size() && m_det != DetId::HGCalHSc) ||

Geometry/HGCalGeometry/test/HGCalGeometryCornerTester.cc

@@ -0,0 +1,76 @@
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "FWCore/ParameterSet/interface/ParameterSet.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/Framework/interface/ESTransientHandle.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+#include "Geometry/Records/interface/IdealGeometryRecord.h"
+#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
+#include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
+#include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
+#include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
+
+class HGCalGeometryCornerTester : public edm::one::EDAnalyzer<> {
+public:
+  explicit HGCalGeometryCornerTester(const edm::ParameterSet& );
+  ~HGCalGeometryCornerTester() override;
+
+  void beginJob() override {}
+  void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
+  void endJob() override {}
+
+private:
+  void doTest(const HGCalGeometry* geom);
+  const std::string    name;
+};
+
+HGCalGeometryCornerTester::HGCalGeometryCornerTester(const edm::ParameterSet& iC) : name(iC.getParameter<std::string>("Detector")) { }
+
+
+HGCalGeometryCornerTester::~HGCalGeometryCornerTester() {}
+
+void HGCalGeometryCornerTester::analyze(const edm::Event& , 
+					const edm::EventSetup& iSetup ) {
+
+  edm::ESHandle<HGCalGeometry> geomH;
+  iSetup.get<IdealGeometryRecord>().get(name,geomH);
+  const HGCalGeometry* geom = (geomH.product());
+  if (!geomH.isValid()) {
+    std::cout << "Cannot get valid HGCalGeometry Object for " << name 
+	      << std::endl;
+  } else {
+    doTest(geom);
+  }
+}
+
+void HGCalGeometryCornerTester::doTest(const HGCalGeometry* geom) {
+
+  const std::vector<DetId>& ids = geom->getValidDetIds();
+  std::cout << "doTest: " << ids.size() << " valid ids for " 
+	    << geom->cellElement() << std::endl;
+  unsigned int kount(0);
+  for (auto const& id : ids) {
+    std::cout << "[" << kount << "] ";
+    if      (id.det() == DetId::Forward)  std::cout <<HGCalDetId(id);
+    else if (id.det() == DetId::HGCalHSc) std::cout <<HGCScintillatorDetId(id);
+    else                                  std::cout <<HGCSiliconDetId(id);
+    std::cout << std::endl;
+    const auto cor = geom->getCorners(id);
+    std::cout << cor.size() << " Corners:";
+    for (unsigned k=0; k<cor.size(); ++k)
+      std::cout << " [" << k << "] " << cor[k];
+    std::cout << std::endl;
+    ++kount;
+  }
+}
+
+//define this as a plug-in
+DEFINE_FWK_MODULE(HGCalGeometryCornerTester);

Geometry/HGCalGeometry/test/python/testHGCalCorner_cfg.py

@@ -0,0 +1,57 @@
+import FWCore.ParameterSet.Config as cms
+
+process = cms.Process("PROD")
+process.load("SimGeneral.HepPDTESSource.pdt_cfi")
+
+#process.load("Geometry.CMSCommonData.cmsExtendedGeometry2023D17XML_cfi")
+#process.load("Geometry.HGCalCommonData.hgcalV6NumberingInitialization_cfi")
+#process.load("Geometry.HGCalCommonData.hgcalV6ParametersInitialization_cfi")
+#process.load("Geometry.CaloEventSetup.HGCalV6Topology_cfi")
+process.load("Geometry.HGCalCommonData.testHGCXML_cfi")
+process.load("Geometry.HGCalCommonData.hgcalParametersInitialization_cfi")
+process.load("Geometry.HGCalCommonData.hgcalNumberingInitialization_cfi")
+process.load("Geometry.CaloEventSetup.HGCalV9Topology_cfi")
+process.load("Geometry.HGCalGeometry.HGCalGeometryESProducer_cfi")
+process.load('FWCore.MessageService.MessageLogger_cfi')
+
+if hasattr(process,'MessageLogger'):
+    process.MessageLogger.categories.append('HGCalGeom')
+
+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.prodEE = cms.EDAnalyzer("HGCalGeometryCornerTester",
+                                Detector   = cms.string("HGCalEESensitive"),
+                                )
+
+process.prodHEF = process.prodEE.clone(
+    Detector   = "HGCalHESiliconSensitive",
+)
+
+process.prodHEB = process.prodEE.clone(
+    Detector   = "HGCalHEScintillatorSensitive",
+)
+
+#process.p1 = cms.Path(process.generator*process.prodEE*process.prodHEF)
+process.p1 = cms.Path(process.generator*process.prodEE*process.prodHEF*process.prodHEB)