Phase2-hgx321 Make some fixes to the issues in storing hits for the scintillator part of HGCal

Geometry/HGCalCommonData/plugins/DDHGCalMixLayer.cc

@@ -448,8 +448,8 @@ void DDHGCalMixLayer::positionMix(const DDLogicalPart& glog,
       double phi2 = dphi * (fimax - fimin + 1);
 #ifdef EDM_ML_DEBUG
       edm::LogVerbatim("HGCalGeom") << "DDHGCalMixLayer: Layer " << copy << " iR "
-                                    << std::get<1>(HGCalTileIndex::tileUnpack(tileIndex_[ly])) << ":"
-                                    << std::get<2>(HGCalTileIndex::tileUnpack(tileIndex_[ly])) << " R " << r1 << ":"
+                                    << std::get<1>(HGCalTileIndex::tileUnpack(tileIndex_[ti])) << ":"
+                                    << std::get<2>(HGCalTileIndex::tileUnpack(tileIndex_[ti])) << " R " << r1 << ":"
                                     << r2 << " Thick " << (2.0 * hthickl) << " phi " << fimin << ":" << fimax << ":"
                                     << convertRadToDeg(phi1) << ":" << convertRadToDeg(phi2);
       ;

Geometry/HGCalCommonData/plugins/DDHGCalMixRotatedLayer.cc

@@ -470,8 +470,8 @@ void DDHGCalMixRotatedLayer::positionMix(const DDLogicalPart& glog,
       auto cshift = cassette_.getShift(layer + 1, 1, cassette);
 #ifdef EDM_ML_DEBUG
       edm::LogVerbatim("HGCalGeom") << "DDHGCalMixRotatedLayer: Layer " << copy << " iR "
-                                    << std::get<1>(HGCalTileIndex::tileUnpack(tileIndex_[ly])) << ":"
-                                    << std::get<2>(HGCalTileIndex::tileUnpack(tileIndex_[ly])) << " R " << r1 << ":"
+                                    << std::get<1>(HGCalTileIndex::tileUnpack(tileIndex_[ti])) << ":"
+                                    << std::get<2>(HGCalTileIndex::tileUnpack(tileIndex_[ti])) << " R " << r1 << ":"
                                     << r2 << " Thick " << (2.0 * hthickl) << " phi " << fimin << ":" << fimax << ":"
                                     << convertRadToDeg(phi1) << ":" << convertRadToDeg(phi2) << " cassette " << cassette
                                     << " Shift " << cshift.first << ":" << cshift.second;

Geometry/HGCalCommonData/plugins/dd4hep/DDHGCalMixLayer.cc

@@ -377,8 +377,8 @@ struct HGCalMixLayer {
         double phi2 = (forFireworks_ == 1) ? (dphi * (fimax - fimin + 1)) : (dphi * fimax);
 #ifdef EDM_ML_DEBUG
         edm::LogVerbatim("HGCalGeom") << "DDHGCalMixLayer: Layer " << copy << " iR "
-                                      << std::get<1>(HGCalTileIndex::tileUnpack(tileIndex_[ly])) << ":"
-                                      << std::get<2>(HGCalTileIndex::tileUnpack(tileIndex_[ly])) << " R "
+                                      << std::get<1>(HGCalTileIndex::tileUnpack(tileIndex_[ti])) << ":"
+                                      << std::get<2>(HGCalTileIndex::tileUnpack(tileIndex_[ti])) << " R "
                                       << cms::convert2mm(r1) << ":" << cms::convert2mm(r2) << " Thick "
                                       << cms::convert2mm((2.0 * hthickl)) << " phi " << fimin << ":" << fimax << ":"
                                       << convertRadToDeg(phi1) << ":" << convertRadToDeg(phi2);

Geometry/HGCalCommonData/plugins/dd4hep/DDHGCalMixRotatedLayer.cc

@@ -397,8 +397,8 @@ struct HGCalMixRotatedLayer {
         auto cshift = cassette_.getShift(layer + 1, 1, cassette);
 #ifdef EDM_ML_DEBUG
         edm::LogVerbatim("HGCalGeom") << "DDHGCalMixRotatedLayer: Layer " << copy << " iR "
-                                      << std::get<1>(HGCalTileIndex::tileUnpack(tileIndex_[ly])) << ":"
-                                      << std::get<2>(HGCalTileIndex::tileUnpack(tileIndex_[ly])) << " R "
+                                      << std::get<1>(HGCalTileIndex::tileUnpack(tileIndex_[ti])) << ":"
+                                      << std::get<2>(HGCalTileIndex::tileUnpack(tileIndex_[ti])) << " R "
                                       << cms::convert2mm(r1) << ":" << cms::convert2mm(r2) << " Thick "
                                       << cms::convert2mm((2.0 * hthickl)) << " phi " << fimin << ":" << fimax << ":"
                                       << convertRadToDeg(phi1) << ":" << convertRadToDeg(phi2) << " cassette "

Geometry/HGCalCommonData/src/HGCalDDDConstants.cc

@@ -183,35 +183,35 @@ std::array<int, 3> HGCalDDDConstants::assignCellTrap(float x, float y, float z,
   const auto& indx = getIndex(layer, reco);
   if (indx.first < 0)
     return std::array<int, 3>{{irad, iphi, type}};
-  double xx = (z > 0) ? x : -x;
-  double phi = (((y == 0.0) && (x == 0.0)) ? 0. : std::atan2(y, xx));
-  if (phi < 0)
-    phi += (2.0 * M_PI);
-  if (indx.second != 0)
-    iphi = 1 + static_cast<int>(phi / indx.second);
+  double xx = (reco) ? ((z > 0) ? x : -x)
+                     : ((z > 0) ? (HGCalParameters::k_ScaleFromDDD * x) : -(HGCalParameters::k_ScaleFromDDD * x));
+  double yy = (reco) ? y : HGCalParameters::k_ScaleFromDDD * y;
+  int ll = layer - hgpar_->firstLayer_;
+  xx -= hgpar_->xLayerHex_[ll];
+  yy -= hgpar_->yLayerHex_[ll];
   if (mode_ == HGCalGeometryMode::TrapezoidCassette) {
     int cassette = HGCalTileIndex::tileCassette(iphi, hgpar_->phiOffset_, hgpar_->nphiCassette_, hgpar_->cassettes_);
     auto cshift = hgcassette_.getShift(layer, 1, cassette);
 #ifdef EDM_ML_DEBUG
     edm::LogVerbatim("HGCalGeom") << "Cassette " << cassette << " Shift " << cshift.first << ":" << cshift.second;
 #endif
-    if (reco) {
-      x -= cshift.first;
-      y -= cshift.second;
-    } else {
-      x = HGCalParameters::k_ScaleFromDDD * x - cshift.first;
-      y = HGCalParameters::k_ScaleFromDDD * y - cshift.second;
-    }
+    xx -= cshift.first;
+    yy -= cshift.second;
   }
+  double phi = (((yy == 0.0) && (xx == 0.0)) ? 0. : std::atan2(yy, xx));
+  if (phi < 0)
+    phi += (2.0 * M_PI);
+  if (indx.second != 0)
+    iphi = 1 + static_cast<int>(phi / indx.second);
   type = hgpar_->scintType(layer);
-  double r = std::sqrt(x * x + y * y);
+  double r = std::sqrt(xx * xx + yy * yy);
   auto ir = std::lower_bound(hgpar_->radiusLayer_[type].begin(), hgpar_->radiusLayer_[type].end(), r);
   irad = static_cast<int>(ir - hgpar_->radiusLayer_[type].begin());
   irad = std::clamp(irad, hgpar_->iradMinBH_[indx.first], hgpar_->iradMaxBH_[indx.first]);
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("HGCalGeom") << "assignCellTrap Input " << x << ":" << y << ":" << z << ":" << layer << ":" << reco
-                                << " x|r " << xx << ":" << r << " phi " << phi << " o/p " << irad << ":" << iphi << ":"
-                                << type;
+                                << " x|y|r " << xx << ":" << yy << ":" << r << " phi " << phi << ":"
+                                << convertRadToDeg(phi) << " o/p " << irad << ":" << iphi << ":" << type;
 #endif
   return std::array<int, 3>{{irad, iphi, type}};
 }
@@ -803,12 +803,16 @@ std::pair<float, float> HGCalDDDConstants::locateCellTrap(int lay, int irad, int
       edm::LogVerbatim("HGCalGeom") << "locateCellTrap:: Input " << lay << ":" << irad << ":" << iphi << ":" << reco
                                     << " IR " << ir << ":" << hgpar_->iradMinBH_[indx.first] << ":"
                                     << hgpar_->iradMaxBH_[indx.first] << " Type " << type << " Z " << indx.first << ":"
-                                    << z << " phi " << phi << " R " << r << ":" << range.first << ":" << range.second;
+                                    << z << " phi " << phi << ":" << convertRadToDeg(phi) << " R " << r << ":"
+                                    << range.first << ":" << range.second;
     if ((mode_ != HGCalGeometryMode::TrapezoidFile) && (mode_ != HGCalGeometryMode::TrapezoidModule) &&
         (mode_ != HGCalGeometryMode::TrapezoidCassette))
       r = std::max(range.first, std::min(r, range.second));
     x = r * std::cos(phi);
     y = r * std::sin(phi);
+    int ll = lay - hgpar_->firstLayer_;
+    x += hgpar_->xLayerHex_[ll];
+    y += hgpar_->yLayerHex_[ll];
     if (mode_ == HGCalGeometryMode::TrapezoidCassette) {
       int cassette = HGCalTileIndex::tileCassette(iphi, hgpar_->phiOffset_, hgpar_->nphiCassette_, hgpar_->cassettes_);
       auto cshift = hgcassette_.getShift(lay, 1, cassette);

Geometry/HGCalCommonData/test/python/testHGCalParametersDDD_cfg.py

@@ -3,7 +3,7 @@
 
 process = cms.Process("HGCalParametersTest",Phase2C11)
 process.load("SimGeneral.HepPDTESSource.pdt_cfi")
-process.load("Configuration.Geometry.GeometryExtended2026D82Reco_cff")
+process.load("Configuration.Geometry.GeometryExtended2026D88Reco_cff")
 process.load('FWCore.MessageService.MessageLogger_cfi')
 
 if hasattr(process,'MessageLogger'):

Geometry/HGCalGeometry/src/HGCalGeometry.cc

@@ -231,6 +231,9 @@ GlobalPoint HGCalGeometry::getPosition(const DetId& detid, bool debug) const {
                                       << xy.second << " ID " << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":"
                                       << id.iCell1 << ":" << id.iCell2 << " Global " << glob;
     }
+  } else {
+    edm::LogVerbatim("HGCalGeom") << "Cannot recognize " << std::hex << detid.rawId() << " cellIndex " << cellIndex
+                                  << ":" << maxSize << " Type " << m_topology.tileTrapezoid();
   }
   return glob;
 }

Geometry/HGCalGeometry/src/HGCalGeometryLoader.cc

@@ -107,7 +107,7 @@ HGCalGeometry* HGCalGeometryLoader::build(const HGCalTopology& topology) {
 #endif
           if (ok) {
             DetId detId = static_cast<DetId>(id);
-            const auto& w = topology.dddConstants().locateCellTrap(layer, ring, iphi, true);
+            const auto& w = topology.dddConstants().locateCellTrap(layer, ring, iphi, true, false);
             double xx = (zside > 0) ? w.first : -w.first;
             CLHEP::Hep3Vector h3v(xx, w.second, mytr.h3v.z());
             const HepGeom::Transform3D ht3d(mytr.hr, h3v);

SimG4CMS/Calo/interface/HGCScintSD.h

@@ -44,10 +44,12 @@ class HGCScintSD : public CaloSD, public Observer<const BeginOfJob *> {
   std::string nameX_;
   HGCalGeometryMode::GeometryMode geom_mode_;
   double eminHit_, slopeMin_, distanceFromEdge_;
-  int levelT1_, levelT2_;
+  int levelT1_, levelT2_, firstLayer_;
   bool storeAllG4Hits_, fiducialCut_;
   bool useBirk_;
   double birk1_, birk2_, birk3_, weight_;
+  std::string fileName_;
+  std::vector<int> tiles_;
 };
 
 #endif  // HGCScintSD_h

SimG4CMS/Calo/interface/HGCalNumberingScheme.h

@@ -38,7 +38,7 @@ class HGCalNumberingScheme {
   const DetId::Detector det_;
   const std::string name_;
   int firstLayer_;
-  std::vector<int> wafers_;
+  std::vector<int> indices_;
 };
 
 #endif

SimG4CMS/Calo/plugins/HGCalHitPartial.cc → SimG4CMS/Calo/plugins/HGCalTestPartialWaferHits.cc

@@ -28,10 +28,10 @@
 #include <string>
 #include <vector>
 
-class HGCalHitPartial : public edm::one::EDAnalyzer<> {
+class HGCalTestPartialWaferHits : public edm::one::EDAnalyzer<> {
 public:
-  HGCalHitPartial(const edm::ParameterSet& ps);
-  ~HGCalHitPartial() override = default;
+  HGCalTestPartialWaferHits(const edm::ParameterSet& ps);
+  ~HGCalTestPartialWaferHits() override = default;
   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
 
 protected:
@@ -46,7 +46,7 @@ class HGCalHitPartial : public edm::one::EDAnalyzer<> {
   std::vector<int> wafers_;
 };
 
-HGCalHitPartial::HGCalHitPartial(const edm::ParameterSet& ps)
+HGCalTestPartialWaferHits::HGCalTestPartialWaferHits(const edm::ParameterSet& ps)
     : g4Label_(ps.getParameter<std::string>("moduleLabel")),
       caloHitSource_(ps.getParameter<std::string>("caloHitSource")),
       nameSense_(ps.getParameter<std::string>("nameSense")),
@@ -78,7 +78,7 @@ HGCalHitPartial::HGCalHitPartial(const edm::ParameterSet& ps)
   }
 }
 
-void HGCalHitPartial::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+void HGCalTestPartialWaferHits::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
   edm::ParameterSetDescription desc;
   desc.add<std::string>("moduleLabel", "g4SimHits");
   desc.add<std::string>("caloHitSource", "HGCHitsEE");
@@ -87,7 +87,7 @@ void HGCalHitPartial::fillDescriptions(edm::ConfigurationDescriptions& descripti
   descriptions.add("hgcalHitPartialEE", desc);
 }
 
-void HGCalHitPartial::analyze(const edm::Event& e, const edm::EventSetup& iS) {
+void HGCalTestPartialWaferHits::analyze(const edm::Event& e, const edm::EventSetup& iS) {
   // get hcalGeometry
   const HGCalGeometry* geom = &iS.getData(geomToken_);
   const HGCalDDDConstants& hgc = geom->topology().dddConstants();
@@ -97,7 +97,7 @@ void HGCalHitPartial::analyze(const edm::Event& e, const edm::EventSetup& iS) {
   bool getHits = (hitsCalo.isValid());
   uint32_t nhits = (getHits) ? hitsCalo->size() : 0;
   uint32_t good(0), allSi(0), all(0);
-  edm::LogVerbatim("HGCalSim") << "HGCalHitPartial: Input flags Hits " << getHits << " with " << nhits
+  edm::LogVerbatim("HGCalSim") << "HGCalTestPartialWaferHits: Input flags Hits " << getHits << " with " << nhits
                                << " hits first Layer " << firstLayer;
 
   if (getHits) {
@@ -139,4 +139,4 @@ void HGCalHitPartial::analyze(const edm::Event& e, const edm::EventSetup& iS) {
 }
 
 //define this as a plug-in
-DEFINE_FWK_MODULE(HGCalHitPartial);
+DEFINE_FWK_MODULE(HGCalTestPartialWaferHits);

SimG4CMS/Calo/plugins/HGCalTestScintHits.cc

@@ -0,0 +1,139 @@
+#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/FileInPath.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+#include "FWCore/Utilities/interface/Exception.h"
+
+#include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
+
+#include "SimDataFormats/CaloHit/interface/PCaloHit.h"
+#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
+#include "SimG4CMS/Calo/interface/CaloSimUtils.h"
+
+#include "Geometry/CaloTopology/interface/HGCalTopology.h"
+#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
+#include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
+#include "Geometry/Records/interface/IdealGeometryRecord.h"
+
+#include <fstream>
+#include <string>
+#include <vector>
+
+class HGCalTestScintHits : public edm::one::EDAnalyzer<> {
+public:
+  HGCalTestScintHits(const edm::ParameterSet& ps);
+  ~HGCalTestScintHits() override = default;
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+protected:
+  void analyze(edm::Event const&, edm::EventSetup const&) override;
+  void beginJob() override {}
+  void endJob() override {}
+
+private:
+  const std::string g4Label_, caloHitSource_, nameSense_, tileFileName_;
+  const edm::EDGetTokenT<edm::PCaloHitContainer> tok_calo_;
+  const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
+  std::vector<int> tiles_;
+};
+
+HGCalTestScintHits::HGCalTestScintHits(const edm::ParameterSet& ps)
+    : g4Label_(ps.getParameter<std::string>("moduleLabel")),
+      caloHitSource_(ps.getParameter<std::string>("caloHitSource")),
+      nameSense_(ps.getParameter<std::string>("nameSense")),
+      tileFileName_(ps.getParameter<std::string>("tileFileName")),
+      tok_calo_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, caloHitSource_))),
+      geomToken_(esConsumes<HGCalGeometry, IdealGeometryRecord>(edm::ESInputTag{"", nameSense_})) {
+  edm::LogVerbatim("HGCalSim") << "Test Hit ID using SimHits for " << nameSense_ << " with module Label: " << g4Label_
+                               << "   Hits: " << caloHitSource_ << " Tile file " << tileFileName_;
+  if (!tileFileName_.empty()) {
+    edm::FileInPath filetmp("SimG4CMS/Calo/data/" + tileFileName_);
+    std::string fileName = filetmp.fullPath();
+    std::ifstream fInput(fileName.c_str());
+    if (!fInput.good()) {
+      edm::LogVerbatim("HGCalSim") << "Cannot open file " << fileName;
+    } else {
+      char buffer[80];
+      while (fInput.getline(buffer, 80)) {
+        std::vector<std::string> items = CaloSimUtils::splitString(std::string(buffer));
+        if (items.size() > 2) {
+          int layer = std::atoi(items[0].c_str());
+          int ring = std::atoi(items[1].c_str());
+          int phi = std::atoi(items[2].c_str());
+          tiles_.emplace_back(HGCalTileIndex::tileIndex(layer, ring, phi));
+        }
+      }
+      edm::LogVerbatim("HGCalSim") << "Reads in " << tiles_.size() << " tile information from " << fileName;
+      fInput.close();
+    }
+  }
+}
+
+void HGCalTestScintHits::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<std::string>("moduleLabel", "g4SimHits");
+  desc.add<std::string>("caloHitSource", "HGCHitsHEback");
+  desc.add<std::string>("nameSense", "HGCalHEScintillatorSensitive");
+  desc.add<std::string>("tileFileName", "");
+  descriptions.add("hgcalHitScintillator", desc);
+}
+
+void HGCalTestScintHits::analyze(const edm::Event& e, const edm::EventSetup& iS) {
+  // get hcalGeometry
+  const HGCalGeometry* geom = &iS.getData(geomToken_);
+  const HGCalDDDConstants& hgc = geom->topology().dddConstants();
+  int firstLayer = hgc.firstLayer() - 1;
+  // get the hit collection
+  const edm::Handle<edm::PCaloHitContainer>& hitsCalo = e.getHandle(tok_calo_);
+  bool getHits = (hitsCalo.isValid());
+  uint32_t nhits = (getHits) ? hitsCalo->size() : 0;
+  uint32_t good(0), all(0);
+  edm::LogVerbatim("HGCalSim") << "HGCalTestScintHits: Input flags Hits " << getHits << " with " << nhits
+                               << " hits first Layer " << firstLayer;
+
+  if (getHits) {
+    std::vector<PCaloHit> hits;
+    hits.insert(hits.end(), hitsCalo->begin(), hitsCalo->end());
+    if (!hits.empty()) {
+      // Loop over all hits
+      for (auto hit : hits) {
+        ++all;
+        DetId id(hit.id());
+        HGCScintillatorDetId hid(id);
+        std::pair<int, int> typm = hgc.tileType(hid.layer(), hid.ring(), 0);
+        if (typm.first >= 0) {
+          hid.setType(typm.first);
+          hid.setSiPM(typm.second);
+          id = static_cast<DetId>(hid);
+        }
+        bool toCheck(true);
+        if (!tiles_.empty()) {
+          int indx = HGCalTileIndex::tileIndex(firstLayer + hid.layer(), hid.ring(), hid.iphi());
+          if (std::find(tiles_.begin(), tiles_.end(), indx) != tiles_.end())
+            toCheck = true;
+        }
+        if (toCheck) {
+          ++good;
+          GlobalPoint pos = geom->getPosition(id);
+          bool valid1 = geom->topology().valid(id);
+          bool valid2 = hgc.isValidTrap(hid.layer(), hid.ring(), hid.iphi());
+          edm::LogVerbatim("HGCalSim") << "Hit[" << all << ":" << good << "]" << hid << " at (" << pos.x() << ", "
+                                       << pos.y() << ", " << pos.z() << ") Validity " << valid1 << ":" << valid2
+                                       << " Energy " << hit.energy() << " Time " << hit.time();
+        }
+      }
+    }
+  }
+  edm::LogVerbatim("HGCalSim") << "Total hits = " << all << ":" << nhits << " Good DetIds = " << good;
+}
+
+//define this as a plug-in
+DEFINE_FWK_MODULE(HGCalTestScintHits);