Phase2-hgx319A Test RecHitTool for partial wafers of HGCal and also add a few utility methods for wafer

Geometry/HGCalCommonData/interface/HGCalWaferType.h

@@ -26,10 +26,14 @@ class HGCalWaferType {
                  int choice,
                  unsigned int cutValue,
                  double cutFracArea);
-  ~HGCalWaferType();
-  int getType(double xpos, double ypos, double zpos);
-  static int getType(int index, const std::vector<int>& indices, const std::vector<int>& types);
+  ~HGCalWaferType() = default;
+
+  static int getCassette(int index, const HGCalParameters::waferInfo_map& wafers);
+  static int getOrient(int index, const HGCalParameters::waferInfo_map& wafers);
+  static int getPartial(int index, const HGCalParameters::waferInfo_map& wafers);
   static int getType(int index, const HGCalParameters::waferInfo_map& wafers);
+  static int getType(int index, const std::vector<int>& indices, const std::vector<int>& types);
+  int getType(double xpos, double ypos, double zpos);
   std::pair<double, double> rLimits(double zpos);
 
 private:

Geometry/HGCalCommonData/src/HGCalDDDConstants.cc

@@ -1336,6 +1336,13 @@ void HGCalDDDConstants::waferFromPosition(const double x,
         if (waferHexagon8File()) {
           int index = HGCalWaferIndex::waferIndex(layer, waferU, waferV);
           celltype = HGCalWaferType::getType(index, hgpar_->waferInfoMap_);
+#ifdef EDM_ML_DEBUG
+          if (debug)
+            edm::LogVerbatim("HGCalGeom") << "Position (" << x << ", " << y << ") Wafer type:partial:orient:cassette "
+                                          << celltype << ":" << HGCalWaferType::getPartial(index, hgpar_->waferInfoMap_)
+                                          << ":" << HGCalWaferType::getOrient(index, hgpar_->waferInfoMap_) << ":"
+                                          << HGCalWaferType::getCassette(index, hgpar_->waferInfoMap_);
+#endif
         } else {
           auto itr = hgpar_->typesInLayers_.find(HGCalWaferIndex::waferIndex(layer, waferU, waferV));
           celltype = ((itr == hgpar_->typesInLayers_.end()) ? HGCSiliconDetId::HGCalCoarseThick

Geometry/HGCalCommonData/src/HGCalWaferType.cc

@@ -31,7 +31,33 @@ HGCalWaferType::HGCalWaferType(const std::vector<double>& rad100,
 #endif
 }
 
-HGCalWaferType::~HGCalWaferType() {}
+int HGCalWaferType::getCassette(int index, const HGCalParameters::waferInfo_map& wafers) {
+  auto itr = wafers.find(index);
+  return ((itr == wafers.end()) ? -1 : ((itr->second).cassette));
+}
+
+int HGCalWaferType::getOrient(int index, const HGCalParameters::waferInfo_map& wafers) {
+  auto itr = wafers.find(index);
+  return ((itr == wafers.end()) ? -1 : ((itr->second).orient));
+}
+
+int HGCalWaferType::getPartial(int index, const HGCalParameters::waferInfo_map& wafers) {
+  auto itr = wafers.find(index);
+  return ((itr == wafers.end()) ? -1 : ((itr->second).part));
+}
+
+int HGCalWaferType::getType(int index, const HGCalParameters::waferInfo_map& wafers) {
+  auto itr = wafers.find(index);
+  return ((itr == wafers.end()) ? -1 : ((itr->second).type));
+}
+
+int HGCalWaferType::getType(int index, const std::vector<int>& indices, const std::vector<int>& properties) {
+  auto itr = std::find(std::begin(indices), std::end(indices), index);
+  int type = (itr == std::end(indices))
+                 ? -1
+                 : HGCalProperty::waferThick(properties[static_cast<unsigned int>(itr - std::begin(indices))]);
+  return type;
+}
 
 int HGCalWaferType::getType(double xpos, double ypos, double zpos) {
   std::vector<double> xc(HGCalParameters::k_CornerSize, 0);
@@ -104,19 +130,6 @@ int HGCalWaferType::getType(double xpos, double ypos, double zpos) {
   return type;
 }
 
-int HGCalWaferType::getType(int index, const std::vector<int>& indices, const std::vector<int>& properties) {
-  auto itr = std::find(std::begin(indices), std::end(indices), index);
-  int type = (itr == std::end(indices))
-                 ? -1
-                 : HGCalProperty::waferThick(properties[static_cast<unsigned int>(itr - std::begin(indices))]);
-  return type;
-}
-
-int HGCalWaferType::getType(int index, const HGCalParameters::waferInfo_map& wafers) {
-  auto itr = wafers.find(index);
-  return ((itr == wafers.end()) ? -1 : ((itr->second).type));
-}
-
 std::pair<double, double> HGCalWaferType::rLimits(double zpos) {
   double zz = std::abs(zpos);
   if (zz < zMin_)

RecoLocalCalo/HGCalRecAlgos/python/hgcalToolTesterPartialWafer_cff.py

@@ -0,0 +1,8 @@
+import FWCore.ParameterSet.Config as cms
+
+from RecoLocalCalo.HGCalRecAlgos.hgcalToolTesterPartialWaferEE_cfi import *
+
+hgcalToolTesterPartialWaferHE = hgcalToolTesterPartialWaferEE.clone(
+    nameSense  = "HGCalHESiliconSensitive",
+    caloHitSource = "HGCHitsHEfront"
+)

RecoLocalCalo/HGCalRecAlgos/test/HGCalToolTesterPartialWafer.cc

@@ -0,0 +1,127 @@
+#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/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/Utilities/interface/Exception.h"
+
+#include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
+
+#include "SimDataFormats/CaloHit/interface/PCaloHit.h"
+#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
+
+#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
+#include "Geometry/CaloTopology/interface/HGCalTopology.h"
+#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
+#include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
+#include "Geometry/HGCalCommonData/interface/HGCalParameters.h"
+#include "Geometry/Records/interface/CaloGeometryRecord.h"
+#include "Geometry/Records/interface/IdealGeometryRecord.h"
+
+#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
+
+#include <string>
+#include <vector>
+
+class HGCalToolTesterPartialWafer : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
+public:
+  HGCalToolTesterPartialWafer(const edm::ParameterSet& ps);
+  ~HGCalToolTesterPartialWafer() override = default;
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+protected:
+  void beginRun(edm::Run const&, edm::EventSetup const&) override;
+  void analyze(edm::Event const&, edm::EventSetup const&) override;
+  void endRun(edm::Run const&, edm::EventSetup const&) override {}
+  void beginJob() override {}
+  void endJob() override {}
+
+private:
+  const std::string g4Label_, caloHitSource_, nameSense_;
+  const edm::EDGetTokenT<edm::PCaloHitContainer> tok_calo_;
+  const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
+  const HGCalGeometry* geom_;
+  CaloGeometry geo_;
+  hgcal::RecHitTools tool_;
+};
+
+HGCalToolTesterPartialWafer::HGCalToolTesterPartialWafer(const edm::ParameterSet& ps)
+    : g4Label_(ps.getParameter<std::string>("moduleLabel")),
+      caloHitSource_(ps.getParameter<std::string>("caloHitSource")),
+      nameSense_(ps.getParameter<std::string>("nameSense")),
+      tok_calo_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, caloHitSource_))),
+      geomToken_(esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()),
+      geom_(nullptr) {
+  edm::LogVerbatim("HGCalSim") << "Test Hit ID using SimHits for " << nameSense_ << " with module Label: " << g4Label_
+                               << "   Hits: " << caloHitSource_;
+}
+
+void HGCalToolTesterPartialWafer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<std::string>("moduleLabel", "g4SimHits");
+  desc.add<std::string>("caloHitSource", "HGCHitsEE");
+  desc.add<std::string>("nameSense", "HGCalEESensitive");
+  descriptions.add("hgcalToolTesterPartialWaferEE", desc);
+}
+
+void HGCalToolTesterPartialWafer::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
+  //Setup the tool
+  geo_ = iSetup.getData(geomToken_);
+  geom_ = (nameSense_ == "HGCalEESensitive") ? static_cast<const HGCalGeometry*>(geo_.getSubdetectorGeometry(
+                                                   DetId::HGCalEE, ForwardSubdetector::ForwardEmpty))
+                                             : static_cast<const HGCalGeometry*>(geo_.getSubdetectorGeometry(
+                                                   DetId::HGCalHSi, ForwardSubdetector::ForwardEmpty));
+  edm::LogVerbatim("HGCalSim") << "HGCalToolTesterPartialWafer: beginRun Called for " << nameSense_;
+}
+
+void HGCalToolTesterPartialWafer::analyze(const edm::Event& e, const edm::EventSetup& iS) {
+  // get HGCal geometry constant
+  tool_.setGeometry(geo_);
+  const HGCalDDDConstants& hgc = geom_->topology().dddConstants();
+
+  // 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), allSi(0), all(0);
+  edm::LogVerbatim("HGCalSim") << "HGCalToolTesterPartialWafer: Input flags Hits " << getHits << " with " << nhits
+                               << " hits";
+
+  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());
+        if ((id.det() == DetId::HGCalEE) || (id.det() == DetId::HGCalHSi)) {
+          ++allSi;
+          HGCSiliconDetId hid(id);
+          const auto& info = hgc.waferInfo(hid.layer(), hid.waferU(), hid.waferV());
+          // Only partial wafers
+          if (info.part != HGCalTypes::WaferFull) {
+            ++good;
+            GlobalPoint pos1 = geom_->getPosition(id);
+            GlobalPoint pos2 = tool_.getPosition(id);
+            edm::LogVerbatim("HGCalSim") << "Hit[" << all << ":" << allSi << ":" << good << "]" << HGCSiliconDetId(id)
+                                         << " Wafer Type:Part:Orient:Cassette " << info.type << ":" << info.part << ":"
+                                         << info.orient << ":" << info.cassette << " at (" << pos1.x() << ", "
+                                         << pos1.y() << ", " << pos1.z() << ") or (" << pos2.x() << ", " << pos2.y()
+                                         << ", " << pos2.z() << ")";
+          }
+        }
+      }
+    }
+  }
+  edm::LogVerbatim("HGCalSim") << "Total hits = " << all << ":" << nhits << " Good DetIds = " << allSi << ":" << good;
+}
+
+//define this as a plug-in
+DEFINE_FWK_MODULE(HGCalToolTesterPartialWafer);

RecoLocalCalo/HGCalRecAlgos/test/testHGCalPartialWafer_cfg.py

@@ -0,0 +1,74 @@
+import FWCore.ParameterSet.Config as cms
+from Configuration.Eras.Era_Phase2C11_cff import Phase2C11
+
+process = cms.Process("PROD",Phase2C11)
+process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi")
+process.load("IOMC.EventVertexGenerators.VtxSmearedGauss_cfi")
+process.load("Configuration.Geometry.GeometryExtended2026D92Reco_cff")
+process.load("Configuration.StandardSequences.MagneticField_cff")
+process.load("Configuration.EventContent.EventContent_cff")
+process.load('Configuration.StandardSequences.Generator_cff')
+process.load('Configuration.StandardSequences.SimIdeal_cff')
+process.load('FWCore.MessageService.MessageLogger_cfi')
+process.load('RecoLocalCalo.HGCalRecAlgos.hgcalToolTesterPartialWafer_cff')
+process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff")
+from Configuration.AlCa.GlobalTag import GlobalTag
+process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase2_realistic_T21', '')
+
+if hasattr(process,'MessageLogger'):
+    process.MessageLogger.HGCalGeom=dict()
+    process.MessageLogger.HGCalSim=dict()
+    process.MessageLogger.CaloSim=dict()
+
+process.load("IOMC.RandomEngine.IOMC_cff")
+process.RandomNumberGeneratorService.generator.initialSeed = 456789
+process.RandomNumberGeneratorService.g4SimHits.initialSeed = 9876
+process.RandomNumberGeneratorService.VtxSmeared.initialSeed = 123456789
+
+process.Timing = cms.Service("Timing")
+
+process.maxEvents = cms.untracked.PSet(
+    input = cms.untracked.int32(25)
+)
+
+process.source = cms.Source("EmptySource",
+    firstRun        = cms.untracked.uint32(1),
+    firstEvent      = cms.untracked.uint32(1)
+)
+
+process.generator = cms.EDProducer("FlatRandomEGunProducer",
+    PGunParameters = cms.PSet(
+        PartID = cms.vint32(211),
+        MinEta = cms.double(1.50),
+        MaxEta = cms.double(3.00),
+        MinPhi = cms.double(-3.1415926),
+        MaxPhi = cms.double(-1.5707963),
+        MinE   = cms.double(100.00),
+        MaxE   = cms.double(100.00)
+    ),
+    Verbosity       = cms.untracked.int32(0),
+    AddAntiParticle = cms.bool(True)
+)
+
+process.output = cms.OutputModule("PoolOutputModule",
+    process.FEVTSIMEventContent,
+    fileName = cms.untracked.string('hgcV17.root')
+)
+
+process.generation_step = cms.Path(process.pgen)
+process.simulation_step = cms.Path(process.psim)
+process.analysis_step = cms.Path(process.hgcalToolTesterPartialWaferEE+process.hgcalToolTesterPartialWaferHE)
+process.out_step = cms.EndPath(process.output)
+
+process.g4SimHits.Physics.type = 'SimG4Core/Physics/FTFP_BERT_EMM'
+
+# Schedule definition
+process.schedule = cms.Schedule(process.generation_step,
+                                process.simulation_step,
+                                process.analysis_step,
+                                process.out_step
+                                )
+
+# filter all path with the production filter sequence
+for path in process.paths:
+        getattr(process,path)._seq = process.generator * getattr(process,path)._seq

SimG4CMS/Calo/plugins/HGCalHitPartial.cc

@@ -25,10 +25,10 @@
 #include <string>
 #include <vector>
 
-class HGcalHitPartial : public edm::one::EDAnalyzer<> {
+class HGCalHitPartial : public edm::one::EDAnalyzer<> {
 public:
-  HGcalHitPartial(const edm::ParameterSet& ps);
-  ~HGcalHitPartial() override = default;
+  HGCalHitPartial(const edm::ParameterSet& ps);
+  ~HGCalHitPartial() override = default;
   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
 
 protected:
@@ -42,7 +42,7 @@ class HGcalHitPartial : public edm::one::EDAnalyzer<> {
   const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
 };
 
-HGcalHitPartial::HGcalHitPartial(const edm::ParameterSet& ps)
+HGCalHitPartial::HGCalHitPartial(const edm::ParameterSet& ps)
     : g4Label_(ps.getParameter<std::string>("moduleLabel")),
       caloHitSource_(ps.getParameter<std::string>("caloHitSource")),
       nameSense_(ps.getParameter<std::string>("nameSense")),
@@ -52,15 +52,15 @@ HGcalHitPartial::HGcalHitPartial(const edm::ParameterSet& ps)
                                << "   Hits: " << caloHitSource_;
 }
 
-void HGcalHitPartial::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+void HGCalHitPartial::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
   edm::ParameterSetDescription desc;
   desc.add<std::string>("moduleLabel", "g4SimHits");
   desc.add<std::string>("caloHitSource", "HGCHitsEE");
   desc.add<std::string>("nameSense", "HGCalEESensitive");
   descriptions.add("hgcalHitPartialEE", desc);
 }
 
-void HGcalHitPartial::analyze(const edm::Event& e, const edm::EventSetup& iS) {
+void HGCalHitPartial::analyze(const edm::Event& e, const edm::EventSetup& iS) {
   // get hcalGeometry
   const HGCalGeometry* geom = &iS.getData(geomToken_);
   const HGCalDDDConstants& hgc = geom->topology().dddConstants();
@@ -69,7 +69,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 << " hits";
+  edm::LogVerbatim("HGCalSim") << "HGCalHitPartial: Input flags Hits " << getHits << " with " << nhits << " hits";
 
   if (getHits) {
     std::vector<PCaloHit> hits;
@@ -100,4 +100,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(HGCalHitPartial);