Run3-hcx342 Try to fix the issue of HCAL hit relabeller

Geometry/HcalCommonData/interface/HcalDDDRecConstants.h

@@ -61,6 +61,7 @@ class HcalDDDRecConstants {
       return gcons;
     }
   }
+  int findDepth(const int& det, const int& eta, const int& phi, const int& zside, const int& lay) const;
   std::vector<int> getDepth(const int& det, const int& phi, const int& zside, const unsigned int& eta) const;
   std::vector<int> getDepth(const unsigned int& eta, const bool& extra) const;
   int getDepthEta16(const int& det, const int& iphi, const int& zside) const {

Geometry/HcalCommonData/interface/HcalHitRelabeller.h

@@ -14,10 +14,9 @@ class HcalHitRelabeller {
   void setGeometry(const HcalDDDRecConstants*&);
   DetId relabel(const uint32_t testId) const;
   static DetId relabel(const uint32_t testId, const HcalDDDRecConstants* theRecNumber);
-
-private:
   double energyWt(const uint32_t testId) const;
 
+private:
   const HcalDDDRecConstants* theRecNumber;
   bool neutralDensity_;
 };

Geometry/HcalCommonData/src/HcalDDDRecConstants.cc

@@ -7,6 +7,7 @@
 #include "CLHEP/Units/GlobalSystemOfUnits.h"
 #include <algorithm>
 #include <cmath>
+#include <sstream>
 
 //#define EDM_ML_DEBUG
 using namespace geant_units::operators;
@@ -27,6 +28,24 @@ HcalDDDRecConstants::~HcalDDDRecConstants() {
 #endif
 }
 
+int HcalDDDRecConstants::findDepth(
+    const int& det, const int& eta, const int& phi, const int& zside, const int& lay) const {
+  int depth = hcons.findDepth(det, eta, phi, zside, lay);
+  if (depth < 0) {
+    std::vector<int> depths = getDepth(eta, false);
+    if ((lay > 0) && (lay <= static_cast<int>(depths.size())))
+      depth = depths[lay - 1];
+#ifdef EDM_ML_DEBUG
+    std::ostringstream st1;
+    st1 << depths.size() << " depths ";
+    for (const auto& d : depths)
+      st1 << ": " << d;
+    edm::LogVerbatim("HCalGeom") << "HcalDDDRecConstants:: " << st1.str() << " for eta = " << eta << " Depth " << depth;
+#endif
+  }
+  return depth;
+}
+
 std::vector<int> HcalDDDRecConstants::getDepth(const unsigned int& eta, const bool& extra) const {
   if (!extra) {
     std::vector<HcalParameters::LayerItem>::const_iterator last = hpar->layerGroupEtaRec.begin();

Geometry/HcalCommonData/src/HcalHitRelabeller.cc

@@ -89,7 +89,9 @@ double HcalHitRelabeller::energyWt(const uint32_t testId) const {
   int det, z, depth, eta, phi, layer;
   HcalTestNumbering::unpackHcalIndex(testId, det, z, depth, eta, phi, layer);
   int zside = (z == 0) ? (-1) : (1);
-  double wt = (((det == 1) || (det == 2)) && (depth == 1)) ? theRecNumber->getLayer0Wt(det, phi, zside) : 1.0;
+  double wt = ((((det == 1) && (layer <= 1)) || ((det == 2) && (layer <= 2))) && (depth == 1))
+                  ? theRecNumber->getLayer0Wt(det, phi, zside)
+                  : 1.0;
 #ifdef EDM_ML_DEBUG
   edm::LogVerbatim("HcalSim") << "EnergyWT::det: " << det << " z: " << z << ":" << zside << " depth: " << depth
                               << " ieta: " << eta << " iphi: " << phi << " layer: " << layer << " wt " << wt;

Geometry/HcalCommonData/src/HcalLayerDepthMap.cc

@@ -254,7 +254,7 @@ int HcalLayerDepthMap::getMaxDepthLastHE(const int subdet, const int iphi, const
 double HcalLayerDepthMap::getLayer0Wt(const int subdet, const int iphi, const int zside) const {
   double wt = isValid(subdet, iphi, zside) ? wtl0C_ : -1.0;
 #ifdef EDM_ML_DEBUG
-  edm::LogVerbatim("HCalGeom") << "Debug info -- getLayer0Wt::Input " << subdet << ":" << iphi << ":" << zside
+  edm::LogVerbatim("HCalGeom") << "HcalLayerDepthMap -- getLayer0Wt::Input " << subdet << ":" << iphi << ":" << zside
                                << " Output " << wt;
 #endif
   return wt;

Geometry/HcalCommonData/test/HcalHitRelabellerTester.cc

@@ -0,0 +1,93 @@
+// -*- C++ -*-
+//
+// Package:    HcalHitRelabellerTester
+// Class:      HcalHitRelabellerTester
+//
+/**\class HcalHitRelabellerTester HcalHitRelabellerTester.cc test/HcalHitRelabellerTester.cc
+
+ Description: <one line class summary>
+
+ Implementation:
+     <Notes on implementation>
+*/
+//
+// Original Author:  Sunanda Banerjee
+//         Created:  Mon 2023/06/24
+//
+
+// system include files
+#include <memory>
+#include <iostream>
+#include <fstream>
+
+// user include files
+#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 "DataFormats/HcalDetId/interface/HcalTestNumbering.h"
+#include "Geometry/Records/interface/IdealGeometryRecord.h"
+#include "Geometry/Records/interface/HcalRecNumberingRecord.h"
+#include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
+#include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
+
+class HcalHitRelabellerTester : public edm::one::EDAnalyzer<> {
+public:
+  explicit HcalHitRelabellerTester(const edm::ParameterSet&);
+  ~HcalHitRelabellerTester() override = default;
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+  void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
+
+private:
+  bool nd_;
+  edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> token_;
+};
+
+HcalHitRelabellerTester::HcalHitRelabellerTester(const edm::ParameterSet& ps)
+    : nd_(ps.getParameter<bool>("neutralDensity")),
+      token_{esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord>(edm::ESInputTag{})} {
+  edm::LogVerbatim("HCalGeom") << "Construct HcalHitRelabellerTester with Neutral Density: " << nd_;
+}
+
+void HcalHitRelabellerTester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<bool>("neutralDensity", true);
+  descriptions.add("hcalHitRelabellerTester", desc);
+}
+
+// ------------ method called to produce the data  ------------
+void HcalHitRelabellerTester::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  const HcalDDDRecConstants* theRecNumber = &iSetup.getData(token_);
+  std::unique_ptr<HcalHitRelabeller> relabel = std::make_unique<HcalHitRelabeller>(nd_);
+  relabel->setGeometry(theRecNumber);
+  std::vector<int> etas = {-29, -26, -22, -19, -16, -13, -10, -7, -4, -1, 2, 5, 8, 11, 14, 17, 20, 23, 26, 29};
+  std::vector<int> layers = {1, 2};
+  const int iphi = 63;
+  edm::LogVerbatim("HCalGeom") << "HcalHitRelabellerTester:: Testing " << etas.size() << " eta, " << layers.size()
+                               << " layer values and iphi = " << iphi;
+  for (const auto& eta : etas) {
+    int ieta = std::abs(eta);
+    int det = (ieta <= 16) ? 1 : 2;
+    int zside = (eta >= 0) ? 1 : -1;
+    for (const auto& lay : layers) {
+      if (ieta == 16)
+        det = (lay <= 3) ? 1 : 2;
+      int depth = theRecNumber->findDepth(det, ieta, iphi, zside, lay);
+      if (depth > 0) {
+        uint32_t id = HcalTestNumbering::packHcalIndex(det, zside, depth, ieta, iphi, lay);
+        double wt = relabel->energyWt(id);
+        edm::LogVerbatim("HCalGeom") << "Det " << det << " Z " << zside << " Eta" << ieta << "  Phi " << iphi << " Lay "
+                                     << lay << " Depth " << depth << " Layer0Wt " << wt;
+      }
+    }
+  }
+}
+
+//define this as a plug-in
+DEFINE_FWK_MODULE(HcalHitRelabellerTester);

Geometry/HcalCommonData/test/python/testHitRelabeller_cfg.py

@@ -0,0 +1,209 @@
+###############################################################################
+# Way to use this:
+#   cmsRun testHitRelabeller_cfg.py geometry=2021
+#
+#   Options for geometry 2016, 2017, 2018, 2021, 2023
+#
+###############################################################################
+import FWCore.ParameterSet.Config as cms
+import os, sys, imp, re
+import FWCore.ParameterSet.VarParsing as VarParsing
+
+####################################################################
+### SETUP OPTIONS
+options = VarParsing.VarParsing('standard')
+options.register('geometry',
+                 "2021",
+                  VarParsing.VarParsing.multiplicity.singleton,
+                  VarParsing.VarParsing.varType.string,
+                  "geometry of operations: 2016, 2017, 2018, 2021, 2023")
+
+### get and parse the command line arguments
+options.parseArguments()
+
+print(options)
+
+####################################################################
+# Use the options
+
+geomName = "Configuration.Geometry.GeometryExtended" + options.geometry + "Reco_cff"
+
+if (options.geometry == "2016"):
+    from Configuration.Eras.Era_Run2_2016_cff import Run2_2016
+    process = cms.Process('HitRelabeller',Run2_2016)
+    baseName = "auto:run2_data"
+elif (options.geometry == "2017"):
+    from Configuration.Eras.Era_Run2_2017_cff import Run2_2017
+    process = cms.Process('HitRelabeller',Run2_2017)
+    baseName = "auto:run2_data"
+elif (options.geometry == "2018"):
+    from Configuration.Eras.Era_Run2_2018_cff import Run2_2018
+    process = cms.Process('HitRelabeller',Run2_2018)
+    baseName = "auto:run2_data"
+else:
+    from Configuration.Eras.Era_Run3_DDD_cff import Run3_DDD
+    process = cms.Process('HitRelabeller',Run3_DDD)
+    baseName = "auto:run3_data"
+
+print("Base file Name: ", baseName)
+print("Geom file Name: ", geomName)
+
+
+# import of standard configurations
+process.load('Configuration.StandardSequences.Services_cff')
+process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi')
+process.load('FWCore.MessageService.MessageLogger_cfi')
+process.load('Configuration.EventContent.EventContent_cff')
+process.load('SimGeneral.MixingModule.mixNoPU_cfi')
+process.load('Configuration.StandardSequences.MagneticField_cff')
+process.load('Configuration.StandardSequences.Generator_cff')
+process.load('IOMC.EventVertexGenerators.VtxSmearedRealistic25ns13p6TeVEarly2022Collision_cfi')
+process.load('GeneratorInterface.Core.genFilterSummary_cff')
+process.load('Configuration.StandardSequences.SimIdeal_cff')
+process.load('Configuration.StandardSequences.EndOfProcess_cff')
+process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff')
+process.load(geomName)
+
+if hasattr(process,'MessageLogger'):
+    process.MessageLogger.HCalGeom=dict()
+
+process.maxEvents = cms.untracked.PSet(
+    input = cms.untracked.int32(1),
+    output = cms.optional.untracked.allowed(cms.int32,cms.PSet)
+)
+
+# Input source
+process.source = cms.Source("EmptySource")
+
+process.options = cms.untracked.PSet(
+    FailPath = cms.untracked.vstring(),
+    IgnoreCompletely = cms.untracked.vstring(),
+    Rethrow = cms.untracked.vstring(),
+    SkipEvent = cms.untracked.vstring(),
+    allowUnscheduled = cms.obsolete.untracked.bool,
+    canDeleteEarly = cms.untracked.vstring(),
+    emptyRunLumiMode = cms.obsolete.untracked.string,
+    eventSetup = cms.untracked.PSet(
+        forceNumberOfConcurrentIOVs = cms.untracked.PSet(
+
+        ),
+        numberOfConcurrentIOVs = cms.untracked.uint32(1)
+    ),
+    fileMode = cms.untracked.string('FULLMERGE'),
+    forceEventSetupCacheClearOnNewRun = cms.untracked.bool(False),
+    makeTriggerResults = cms.obsolete.untracked.bool,
+    numberOfConcurrentLuminosityBlocks = cms.untracked.uint32(1),
+    numberOfConcurrentRuns = cms.untracked.uint32(1),
+    numberOfStreams = cms.untracked.uint32(0),
+    numberOfThreads = cms.untracked.uint32(1),
+    printDependencies = cms.untracked.bool(False),
+    sizeOfStackForThreadsInKB = cms.optional.untracked.uint32,
+    throwIfIllegalParameter = cms.untracked.bool(True),
+    wantSummary = cms.untracked.bool(False)
+)
+
+# Production Info
+process.configurationMetadata = cms.untracked.PSet(
+    annotation = cms.untracked.string('ZMM_14TeV_TuneCP5_cfi nevts:10'),
+    name = cms.untracked.string('Applications'),
+    version = cms.untracked.string('$Revision: 1.19 $')
+)
+
+# Other statements
+process.genstepfilter.triggerConditions=cms.vstring("generation_step")
+from Configuration.AlCa.GlobalTag import GlobalTag 
+process.GlobalTag = GlobalTag(process.GlobalTag, baseName, '')
+
+process.generator = cms.EDFilter("Pythia8GeneratorFilter",
+                                 pythiaHepMCVerbosity = cms.untracked.bool(False),
+                                 maxEventsToPrint = cms.untracked.int32(0),
+                                 pythiaPylistVerbosity = cms.untracked.int32(0),
+                                 filterEfficiency = cms.untracked.double(1.0),
+                                 comEnergy = cms.double(14000.0),
+                                 PythiaParameters = cms.PSet(
+                                     pythia8CommonSettings = cms.vstring(
+                                         'Tune:preferLHAPDF = 2',
+                                         'Main:timesAllowErrors = 10000',
+                                         'Check:epTolErr = 0.01',
+                                         'Beams:setProductionScalesFromLHEF = off',
+                                         'SLHA:minMassSM = 1000.',
+                                         'ParticleDecays:limitTau0 = on',
+                                         'ParticleDecays:tau0Max = 10',
+                                        'ParticleDecays:allowPhotonRadiation = on',
+                                     ),
+                                     pythia8CP5Settings = cms.vstring(
+                                         'Tune:pp 14',
+                                         'Tune:ee 7',
+                                         'MultipartonInteractions:ecmPow=0.03344',
+                                         'MultipartonInteractions:bProfile=2',
+                                         'MultipartonInteractions:pT0Ref=1.41',
+                                         'MultipartonInteractions:coreRadius=0.7634',
+                                         'MultipartonInteractions:coreFraction=0.63',
+                                         'ColourReconnection:range=5.176',
+                                         'SigmaTotal:zeroAXB=off',
+                                         'SpaceShower:alphaSorder=2',
+                                         'SpaceShower:alphaSvalue=0.118',
+                                         'SigmaProcess:alphaSvalue=0.118',
+                                         'SigmaProcess:alphaSorder=2',
+                                         'MultipartonInteractions:alphaSvalue=0.118',
+                                         'MultipartonInteractions:alphaSorder=2',
+                                         'TimeShower:alphaSorder=2',
+                                         'TimeShower:alphaSvalue=0.118',
+                                         'SigmaTotal:mode = 0',
+                                         'SigmaTotal:sigmaEl = 21.89',
+                                         'SigmaTotal:sigmaTot = 100.309',
+                                         'PDF:pSet=LHAPDF6:NNPDF31_nnlo_as_0118',
+                                     ),
+                                     processParameters = cms.vstring(
+                                         'WeakSingleBoson:ffbar2gmZ = on',
+                                         '23:onMode = off',
+                                         '23:onIfAny = 13',
+                                         'PhaseSpace:mHatMin = 75.',
+                                     ),
+                                     parameterSets = cms.vstring('pythia8CommonSettings',
+                                                                 'pythia8CP5Settings',
+                                                                 'processParameters',
+                                                             )
+                                 )
+)
+
+process.mumugenfilter = cms.EDFilter("MCParticlePairFilter",
+    MaxEta = cms.untracked.vdouble(4.0, 4.0),
+    MinEta = cms.untracked.vdouble(-4.0, -4.0),
+    MinPt = cms.untracked.vdouble(2.5, 2.5),
+    ParticleCharge = cms.untracked.int32(-1),
+    ParticleID1 = cms.untracked.vint32(13),
+    ParticleID2 = cms.untracked.vint32(13),
+    Status = cms.untracked.vint32(1, 1)
+)
+
+
+process.ProductionFilterSequence = cms.Sequence(process.generator+process.mumugenfilter)
+process.load("Geometry.HcalCommonData.hcalHitRelabellerTester_cfi")
+
+# Path and EndPath definitions
+process.generation_step = cms.Path(process.pgen)
+process.simulation_step = cms.Path(process.psim)
+process.genfiltersummary_step = cms.EndPath(process.genFilterSummary)
+process.endjob_step = cms.EndPath(process.endOfProcess)
+process.analysis_step = cms.EndPath(process.hcalHitRelabellerTester)
+
+# Schedule definition
+process.schedule = cms.Schedule(process.generation_step,
+                                process.genfiltersummary_step,
+                                process.simulation_step,
+                                process.endjob_step,
+                                process.analysis_step)
+from PhysicsTools.PatAlgos.tools.helpers import associatePatAlgosToolsTask
+associatePatAlgosToolsTask(process)
+# filter all path with the production filter sequence
+for path in process.paths:
+	getattr(process,path).insert(0, process.ProductionFilterSequence)
+
+
+# Customisation from command line
+
+# Add early deletion of temporary data products to reduce peak memory need
+from Configuration.StandardSequences.earlyDeleteSettings_cff import customiseEarlyDelete
+process = customiseEarlyDelete(process)
+# End adding early deletion