Merge pull request #46730 from mmusich/mm_dev_SiStripZeroSuppression_…

…crash_v2

[14.1.X] `SiStripZeroSupperssion`: discard processing unphysical clusters in `suppressHybridData`

DataFormats/SiStripDetId/interface/SiStripDetId.h

@@ -3,6 +3,7 @@
 
 #include "DataFormats/DetId/interface/DetId.h"
 #include "DataFormats/SiStripDetId/interface/SiStripEnums.h"
+#include "FWCore/Utilities/interface/Exception.h"
 #include <ostream>
 
 class SiStripDetId;
@@ -59,6 +60,9 @@ class SiStripDetId : public DetId {
   /** Returns strip length of strip tracker sensor, otherwise null. */
   inline double stripLength() const;
 
+  /** Returns number of APVs connected to each strip tracker module */
+  inline unsigned int numberOfAPVs() const;
+
   // ---------- Constructors that set "reserved" field ----------
 
   /** Construct from a raw value and set "reserved" field. */
@@ -183,4 +187,25 @@ double SiStripDetId::stripLength() const { return 0.; }
 
 uint16_t SiStripDetId::reserved() const { return static_cast<uint16_t>((id_ >> reservedStartBit_) & reservedMask_); }
 
+unsigned int SiStripDetId::numberOfAPVs() const {
+  const auto &moduleGeometry = this->moduleGeometry();
+
+  // Check for unknown geometry and throw a meaningful exception
+  if (moduleGeometry == SiStripModuleGeometry::UNKNOWNGEOMETRY) {
+    throw cms::Exception("InvalidModuleGeometry")
+        << "Error in SiStripDetId::numberOfAPVs(): Module geometry is UNKNOWNGEOMETRY. "
+        << "This indicates an invalid or unsupported module geometry for the current DetId: " << this->rawId();
+  }
+
+  // Determine the number of APVs based on the module geometry
+  if (moduleGeometry == SiStripModuleGeometry::IB1 || moduleGeometry == SiStripModuleGeometry::OB1 ||
+      moduleGeometry == SiStripModuleGeometry::W1A || moduleGeometry == SiStripModuleGeometry::W1B ||
+      moduleGeometry == SiStripModuleGeometry::W2A || moduleGeometry == SiStripModuleGeometry::W2B ||
+      moduleGeometry == SiStripModuleGeometry::W5) {
+    return 6;
+  } else {
+    return 4;
+  }
+}
+
 #endif  // DataFormats_SiStripDetId_SiStripDetId_h

RecoLocalTracker/SiStripClusterizer/plugins/ClustersFromRawProducer.cc

@@ -12,6 +12,7 @@
 #include "DataFormats/FEDRawData/interface/FEDRawDataCollection.h"
 #include "EventFilter/SiStripRawToDigi/interface/SiStripFEDBuffer.h"
 #include "DataFormats/SiStripCommon/interface/SiStripConstants.h"
+#include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
 
 #include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h"
 
@@ -407,7 +408,9 @@ void ClusterFiller::fill(StripClusterizerAlgorithm::output_t::TSFastFiller& reco
               buffer->channel(fedCh), std::back_inserter(unpDigis), ipair * 256, isNonLite, mode, legacy_, lmode, pCode);
           if (fedchannelunpacker::StatusCode::SUCCESS == st_ch) {
             edm::DetSet<SiStripDigi> suppDigis{id};
-            rawAlgos.suppressHybridData(unpDigis, suppDigis, ipair * 2);
+            unsigned int detId = suppDigis.id;
+            uint16_t maxNStrips = SiStripDetId(detId).numberOfAPVs() * 128;
+            rawAlgos.suppressHybridData(maxNStrips, unpDigis, suppDigis, ipair * 2);
             std::copy(std::begin(suppDigis), std::end(suppDigis), perStripAdder);
           }
         }

RecoLocalTracker/SiStripZeroSuppression/interface/SiStripRawProcessingAlgorithms.h

@@ -21,7 +21,8 @@ class SiStripRawProcessingAlgorithms {
   void initialize(const edm::EventSetup&);
   void initialize(const edm::EventSetup&, const edm::Event&);
 
-  uint16_t suppressHybridData(const edm::DetSet<SiStripDigi>& inDigis,
+  uint16_t suppressHybridData(const uint16_t maxStrip,
+                              const edm::DetSet<SiStripDigi>& inDigis,
                               edm::DetSet<SiStripDigi>& suppressedDigis,
                               uint16_t firstAPV = 0);
 

RecoLocalTracker/SiStripZeroSuppression/plugins/SiStripZeroSuppression.cc

@@ -5,10 +5,12 @@
 #include "DataFormats/Common/interface/Handle.h"
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
 #include "DataFormats/Common/interface/DetSetVector.h"
+#include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
 #include "DataFormats/SiStripDigi/interface/SiStripDigi.h"
 #include "DataFormats/SiStripDigi/interface/SiStripRawDigi.h"
 #include "RecoLocalTracker/SiStripZeroSuppression/interface/SiStripRawProcessingFactory.h"
 #include "FWCore/Utilities/interface/transform.h"
+#include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
 #include <memory>
 
 SiStripZeroSuppression::SiStripZeroSuppression(edm::ParameterSet const& conf)
@@ -151,8 +153,11 @@ inline void SiStripZeroSuppression::processHybrid(const edm::DetSetVector<SiStri
   for (const auto& inDigis : input) {
     edm::DetSet<SiStripDigi> suppressedDigis(inDigis.id);
 
+    unsigned int detId = inDigis.id;
+    uint16_t maxNStrips = SiStripDetId(detId).numberOfAPVs() * 128;
+
     uint16_t nAPVflagged = 0;
-    nAPVflagged = algorithms->suppressHybridData(inDigis, suppressedDigis);
+    nAPVflagged = algorithms->suppressHybridData(maxNStrips, inDigis, suppressedDigis);
 
     storeExtraOutput(inDigis.id, nAPVflagged);
     if (!suppressedDigis.empty())

RecoLocalTracker/SiStripZeroSuppression/src/SiStripRawProcessingAlgorithms.cc

@@ -55,44 +55,76 @@ void SiStripRawProcessingAlgorithms::initialize(const edm::EventSetup& es, const
  *
  * @return number of restored APVs
  */
-uint16_t SiStripRawProcessingAlgorithms::suppressHybridData(const edm::DetSet<SiStripDigi>& hybridDigis,
+uint16_t SiStripRawProcessingAlgorithms::suppressHybridData(const uint16_t maxNStrips,
+                                                            const edm::DetSet<SiStripDigi>& hybridDigis,
                                                             edm::DetSet<SiStripDigi>& suppressedDigis,
                                                             uint16_t firstAPV) {
-  uint16_t nAPVFlagged = 0;
-  auto beginAPV = hybridDigis.begin();
-  const auto indigis_end = hybridDigis.end();
-  auto iAPV = firstAPV;
-  while (beginAPV != indigis_end) {
-    const auto endAPV = std::lower_bound(beginAPV, indigis_end, SiStripDigi((iAPV + 1) * 128, 0));
-    const auto nDigisInAPV = std::distance(beginAPV, endAPV);
+  uint16_t nAPVFlagged = 0;  // Count of flagged APVs
+  auto currentDigi = hybridDigis.begin();
+  const auto endDigi = hybridDigis.end();
+  auto currentAPV = firstAPV;
+
+  // Loop through the APVs in the range
+  while (currentDigi != endDigi) {
+    // Determine the range of digis belonging to the current APV
+    const auto nextAPVBoundary = SiStripDigi((currentAPV + 1) * 128, 0);
+
+    // Reject any APV larger than the max possible
+    if (nextAPVBoundary.strip() > maxNStrips) {
+      edm::LogError("SiStripRawProcessingAlgorithms")
+          << "In DetId " << suppressedDigis.id << " encountered APV boundary with strip number "
+          << nextAPVBoundary.strip() << ", which exceeds the maximum allowed value for this module (" << maxNStrips
+          << "). Exiting loop.";
+      break;
+    }
+
+    const auto nextAPVDigi = std::lower_bound(currentDigi, endDigi, nextAPVBoundary);
+    const auto nDigisInAPV = std::distance(currentDigi, nextAPVDigi);
+
+    // Handle based on the number of digis in the current APV
     if (nDigisInAPV > 64) {
+      // Process hybrid data for noisy APV
       digivector_t workDigis(128, -1024);
-      for (auto it = beginAPV; it != endAPV; ++it) {
-        workDigis[it->strip() - 128 * iAPV] = it->adc() * 2 - 1024;
+
+      // Populate `workDigis` with values from `currentDigi`
+      for (auto it = currentDigi; it != nextAPVDigi; ++it) {
+        workDigis[it->strip() - 128 * currentAPV] = it->adc() * 2 - 1024;
       }
+
+      // Perform pedestal subtraction
       digivector_t workDigisPedSubtracted(workDigis);
-      subtractorCMN->subtract(hybridDigis.id, iAPV, workDigis);
+      subtractorCMN->subtract(hybridDigis.id, currentAPV, workDigis);
+
+      // Inspect and restore digis
       const auto apvFlagged = restorer->inspectAndRestore(
-          hybridDigis.id, iAPV, workDigisPedSubtracted, workDigis, subtractorCMN->getAPVsCM());
+          hybridDigis.id, currentAPV, workDigisPedSubtracted, workDigis, subtractorCMN->getAPVsCM());
       nAPVFlagged += apvFlagged;
-      if (getAPVFlags()[iAPV]) {
-        suppressor->suppress(workDigis, iAPV, suppressedDigis);
-      } else {  // bad APV: more than 64 but not flagged
-        for (uint16_t i = 0; i != 128; ++i) {
+
+      // Process based on the APV flag
+      if (getAPVFlags()[currentAPV]) {
+        // Suppress flagged APVs
+        suppressor->suppress(workDigis, currentAPV, suppressedDigis);
+      } else {
+        // Handle bad APV: not flagged but exceeds threshold
+        for (uint16_t i = 0; i < 128; ++i) {
           const auto digi = workDigisPedSubtracted[i];
           if (digi > 0) {
-            suppressedDigis.push_back(SiStripDigi(iAPV * 128 + i, suppressor->truncate(digi)));
+            suppressedDigis.push_back(SiStripDigi(currentAPV * 128 + i, suppressor->truncate(digi)));
           }
         }
       }
-    } else {  // already zero-suppressed, copy and truncate
-      std::transform(beginAPV, endAPV, std::back_inserter(suppressedDigis), [this](const SiStripDigi inDigi) {
+    } else {
+      // Already zero-suppressed: copy and truncate
+      std::transform(currentDigi, nextAPVDigi, std::back_inserter(suppressedDigis), [this](const SiStripDigi& inDigi) {
         return SiStripDigi(inDigi.strip(), suppressor->truncate(inDigi.adc()));
       });
     }
-    beginAPV = endAPV;
-    ++iAPV;
+
+    // Move to the next APV
+    currentDigi = nextAPVDigi;
+    ++currentAPV;
   }
+
   return nAPVFlagged;
 }