New version of RawPrime (SiStripApproximateCluster)

DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h

@@ -1,7 +1,9 @@
 #ifndef DataFormats_SiStripCluster_SiStripApproximateCluster_h
 #define DataFormats_SiStripCluster_SiStripApproximateCluster_h
-
+#include "FWCore/Utilities/interface/Exception.h"
 #include "FWCore/Utilities/interface/typedefs.h"
+#include <climits>
+#include <bits/stdc++.h>
 
 class SiStripCluster;
 class SiStripApproximateCluster {
@@ -13,25 +15,125 @@ class SiStripApproximateCluster {
                                      cms_uint8_t avgCharge,
                                      bool filter,
                                      bool isSaturated,
-                                     bool peakFilter = false)
+                                     bool peakFilter = false,
+                                     cms_uint8_t version_ = 1)
       : barycenter_(barycenter),
         width_(width),
         avgCharge_(avgCharge),
         filter_(filter),
         isSaturated_(isSaturated),
-        peakFilter_(peakFilter) {}
+        peakFilter_(peakFilter),
+        version_(version_) {}
 
   explicit SiStripApproximateCluster(const SiStripCluster& cluster,
                                      unsigned int maxNSat,
                                      float hitPredPos,
-                                     bool peakFilter);
+                                     bool peakFilter,
+                                     cms_uint8_t version_ = 1,
+                                     float previous_barycenter = 0,
+                                     unsigned int offset_module_change = 0);
+
+  // getBarycenter returns the barycenter as a *float* in strips (e.g. 1.0 means center of strip 1)
+  float getBarycenter(float previous_barycenter = 0, unsigned int offset_module_change = 0) const {
+    switch (version_) {
+      case 1:
+        return barycenter_ * 0.1;  // in the old format barycenter_ is in tenths of strips
+      case 2: {
+        // Drop the first bit (encoding the saturation info)
+        // barycenterRangeMax_ is 1 in 15 bits (the 16th bit is used to encode isSaturated_ info)
+        double barycenter_decoded = (barycenter_ & barycenterRangeMax_);
+        // rescale to get the original barycenter value (float)
+        return barycenter_decoded / barycenterScale_ - offset_module_change + previous_barycenter;
+      }
+      default:
+        throw cms::Exception("VersionNotSupported")
+            << "Version " << int(version_) << " of SiStripApproximateCluster not supported";
+    }
+  }
+  //   // kept for compatibility with v1, should be removed in the future
+  //   cms_uint16_t barycenter(float previous_barycenter, unsigned int offset_module_change) const {
+  //   switch (version_){
+  //     case 2: {
+  //       return std::round(getBarycenter(previous_barycenter, offset_module_change)*10.); // return barycenter in tenths of strips for compatibility with v1
+  //     }
+  //     default: throw cms::Exception("VersionNotSupported") << "Version " << int(version_) << " of SiStripApproximateCluster not supported for SiStripApproximateCluster::barycenter(float,unsigned int)";
+  //   }
+  //  }
+
+  // getAvgCharge returns the average charge as a *float* in ADC counts (e.g. 1.0 means 1 ADC count)
+  float getAvgCharge() const {
+    switch (version_) {
+      case 1:
+        return avgCharge_;
+      case 2: {
+        // Drop the first two bits (encoding filter_ and peakFilter_ info)
+        // avgChargeRangeMax_ is 1 in 6 bits (the 2 highest bits are used to encode filter_ and peakFilter_ info)
+        double avgCharge_decoded = (avgCharge_ & avgChargeRangeMax_);
+        // Rescale to get the original average charge (float)
+        return (avgCharge_decoded)*avgChargeScale_ + avgChargeOffset_;
+      }
+      default:
+        throw cms::Exception("VersionNotSupported")
+            << "Version " << int(version_) << " of SiStripApproximateCluster not supported";
+    }
+  }
+
+  //barycenter() returns barycenter position in tenths of strip (i.e. 10 means center of strip 1) (0-1536)
+  //avgCharge() returns the average charge in ADC counts (0-255)
+  //width() returns the cluster width (0-255)
+  //version() returns true if the cluster is in the new format (Fall 2025)
 
   cms_uint16_t barycenter() const { return barycenter_; }
   cms_uint8_t width() const { return width_; }
-  cms_uint8_t avgCharge() const { return avgCharge_; }
-  bool filter() const { return filter_; }
-  bool isSaturated() const { return isSaturated_; }
-  bool peakFilter() const { return peakFilter_; }
+  cms_uint8_t avgCharge() const {  // should be a float in the future instead of an int
+    switch (version_) {
+      case 1:
+        return avgCharge_;
+      default:
+        throw cms::Exception("VersionNotSupported")
+            << "Version " << int(version_)
+            << " of SiStripApproximateCluster not supported for SiStripApproximateCluster::avgCharge()";
+    }
+  }
+
+  bool filter() const {
+    switch (version_) {
+      case 1:
+        return filter_;
+      // In v2, filter_ info are encoded in avgCharge_
+      case 2:
+        return (avgCharge_ & (1 << kFilterBit));
+      default:
+        throw cms::Exception("VersionNotSupported")
+            << "Version " << int(version_) << " of SiStripApproximateCluster not supported";
+    }
+  }
+  bool peakFilter() const {
+    switch (version_) {
+      case 1:
+        return peakFilter_;
+      // In v2, filter_ and peakFilter_ info are encoded in avgCharge_
+      case 2:
+        return (avgCharge_ & (1 << kPeakFilterBit));
+      default:
+        throw cms::Exception("VersionNotSupported")
+            << "Version " << int(version_) << " of SiStripApproximateCluster not supported";
+    }
+  }
+
+  bool isSaturated() const {
+    switch (version_) {
+      case 1:
+        return isSaturated_;
+      // In v2, isSaturated_ info is encoded in barycenter_
+      case 2:
+        return (barycenter_ & (1 << kSaturatedBit));
+      default:
+        throw cms::Exception("VersionNotSupported")
+            << "Version " << int(version_) << " of SiStripApproximateCluster not supported";
+    }
+  }
+  char version() const { return version_; }
 
 private:
   cms_uint16_t barycenter_ = 0;
@@ -40,10 +142,43 @@ class SiStripApproximateCluster {
   bool filter_ = false;
   bool isSaturated_ = false;
   bool peakFilter_ = false;
+  // v2 --> new version
+  cms_uint8_t version_ = 1;
   static constexpr double trimMaxADC_ = 30.;
   static constexpr double trimMaxFracTotal_ = .15;
   static constexpr double trimMaxFracNeigh_ = .25;
   static constexpr double maxTrimmedSizeDiffNeg_ = .7;
   static constexpr double maxTrimmedSizeDiffPos_ = 1.;
+
+  ////// Encoding constants for v2 ///////////
+  // maximum value of barycenter_ is 768 strips (128 strips/APV * 6 APVs)
+  // multiplied by a factor 2 as we save the distance from the previous cluster, which can be in another module
+  static constexpr double barycenterMax_ = 768. * 2;
+  // get the total number of bits in barycenter_ (16 bits for cms_uint16_t)
+  static constexpr int nbits_barycenter_ = sizeof(barycenter_) * CHAR_BIT;
+  // position of the bit used to encode isSaturated_ in barycenter_
+  static constexpr int kSaturatedBit = nbits_barycenter_ - 1;
+  // get the largest number storable in barycenter_ with the remaining bits (2^15 -1 = 32767)
+  static constexpr int barycenterRangeMax_ = (1 << (nbits_barycenter_ - 1)) - 1;
+
+  // maximum value of avgCharge_ is 255 ADC counts
+  static constexpr double avgChargeMax_ = 255.;
+  // get the number of bits in avgCharge_ (8 bits for cms_uint8_t)
+  static constexpr int nbits_avgCharge_ = sizeof(avgCharge_) * CHAR_BIT;
+  // positions of the bit used to encode filter_ and peakFilter_ in avgCharge_
+  static constexpr int kPeakFilterBit = nbits_avgCharge_ - 1;
+  static constexpr int kFilterBit = nbits_avgCharge_ - 2;
+  // get the largest number storable in avgCharge_ with the remaining bits (2^6 -1 = 63)
+  static constexpr int avgChargeRangeMax_ = (1 << (nbits_avgCharge_ - 2)) - 1;
+
+public:
+  // constants used in compression and decompression,
+  // floor approximate to integer to get no approximation error for integers
+  static constexpr int barycenterScale_ = barycenterRangeMax_ / barycenterMax_;
+  static constexpr int avgChargeScale_ = avgChargeMax_ / avgChargeRangeMax_;
+
+  static constexpr float barycenterOffset_ = +0.5;  // to get no approximation error for integers
+  static constexpr float avgChargeOffset_ = +2.0;   // to minimize bias
+  ////////////////////////////////////////
 };
 #endif  // DataFormats_SiStripCluster_SiStripApproximateCluster_h

DataFormats/SiStripCluster/interface/SiStripApproximateClusterCollectionV2.h

@@ -0,0 +1,111 @@
+#ifndef DataFormats_SiStripCluster_SiStripApproximateClusterCollectionV2_v1_h
+#define DataFormats_SiStripCluster_SiStripApproximateClusterCollectionV2_v1_h
+
+#include <vector>
+
+#include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h"
+#include <iostream>
+#include <numeric>
+/**
+ * This class provides a minimal interface that resembles
+ * edmNew::DetSetVector, but is crafted such that we are comfortable
+ * to provide an infinite backwards compatibility guarantee for it
+ * (like all RAW data). Any modifications need to be made with care.
+ * Please consult core software group if in doubt.
+**/
+class SiStripApproximateClusterCollectionV2 {
+public:
+  // Helper classes to make creation and iteration easier
+  class Filler {
+  public:
+    void push_back(SiStripApproximateCluster const& cluster) { clusters_.push_back(cluster); }
+
+  private:
+    friend SiStripApproximateClusterCollectionV2;
+    Filler(std::vector<SiStripApproximateCluster>& clusters) : clusters_(clusters) {}
+
+    std::vector<SiStripApproximateCluster>& clusters_;
+  };
+
+  class const_iterator;
+  class DetSet {
+  public:
+    using const_iterator = std::vector<SiStripApproximateCluster>::const_iterator;
+
+    unsigned int id() const {
+      return std::accumulate(coll_->detIds_.cbegin(), coll_->detIds_.cbegin() + detIndex_ + 1, 0);
+    }
+
+    void move(unsigned int clusBegin) const { clusBegin_ = clusBegin; }
+    const_iterator begin() const { return coll_->clusters_.begin() + clusBegin_; }
+    const_iterator cbegin() const { return begin(); }
+    const_iterator end() const { return coll_->clusters_.begin() + clusEnd_; }
+    const_iterator cend() const { return end(); }
+
+  private:
+    friend SiStripApproximateClusterCollectionV2::const_iterator;
+    DetSet(SiStripApproximateClusterCollectionV2 const* coll, unsigned int detIndex)
+        : coll_(coll),
+          detIndex_(detIndex),
+          clusEnd_(coll->clusters_.size())
+    //clusBegin_(coll_->beginIndices_[detIndex]),
+    //clusEnd_(detIndex == coll_->beginIndices_.size() - 1 ? coll->clusters_.size()
+    //: coll_->beginIndices_[detIndex + 1]) {}
+    {}
+
+    SiStripApproximateClusterCollectionV2 const* const coll_;
+    unsigned int const detIndex_;
+    mutable unsigned int clusBegin_ = 0;
+    unsigned int const clusEnd_;
+  };
+
+  class const_iterator {
+  public:
+    DetSet operator*() const { return DetSet(coll_, index_); }
+
+    const_iterator& operator++() {
+      ++index_;
+      if (index_ == coll_->detIds_.size()) {
+        *this = const_iterator();
+      }
+      return *this;
+    }
+
+    const_iterator operator++(int) {
+      const_iterator clone = *this;
+      ++(*this);
+      return clone;
+    }
+
+    bool operator==(const_iterator const& other) const { return coll_ == other.coll_ and index_ == other.index_; }
+    bool operator!=(const_iterator const& other) const { return not operator==(other); }
+
+  private:
+    friend SiStripApproximateClusterCollectionV2;
+    // default-constructed object acts as the sentinel
+    const_iterator() = default;
+    const_iterator(SiStripApproximateClusterCollectionV2 const* coll) : coll_(coll) {}
+
+    SiStripApproximateClusterCollectionV2 const* coll_ = nullptr;
+    unsigned int index_ = 0;
+  };
+
+  // Actual public interface
+  SiStripApproximateClusterCollectionV2() = default;
+
+  void reserve(std::size_t dets, std::size_t clusters);
+  Filler beginDet(unsigned int detId);
+
+  const_iterator begin() const { return clusters_.empty() ? end() : const_iterator(this); }
+  const_iterator cbegin() const { return begin(); }
+  const_iterator end() const { return const_iterator(); }
+  const_iterator cend() const { return end(); }
+
+private:
+  // The detIds_ and beginIndices_ have one element for each Det. An
+  // element of beginIndices_ points to the first cluster of the Det
+  // in clusters_.
+  std::vector<unsigned int> detIds_;  // DetId for the Det
+  std::vector<SiStripApproximateCluster> clusters_;
+};
+#endif

DataFormats/SiStripCluster/interface/SiStripCluster.h

@@ -44,8 +44,10 @@ class SiStripCluster {
     initQB();
   }
 
-  SiStripCluster(const SiStripApproximateCluster cluster, const uint16_t maxStrips);
-
+  SiStripCluster(const SiStripApproximateCluster cluster,
+                 const uint16_t maxStrips,
+                 float previous_barycenter = 0,
+                 unsigned int offset_module_change = 0);
   // extend the cluster
   template <typename Iter>
   void extend(Iter begin, Iter end) {

DataFormats/SiStripCluster/src/SiStripApproximateCluster.cc

@@ -1,18 +1,23 @@
 #include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h"
 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
 #include <algorithm>
+#include <cassert>
 #include <cmath>
 
 SiStripApproximateCluster::SiStripApproximateCluster(const SiStripCluster& cluster,
                                                      unsigned int maxNSat,
                                                      float hitPredPos,
-                                                     bool peakFilter) {
+                                                     bool peakFilter,
+                                                     unsigned char version,
+                                                     float previous_barycenter,
+                                                     unsigned int offset_module_change) {
   barycenter_ = std::round(cluster.barycenter() * 10);
   width_ = cluster.size();
   avgCharge_ = cluster.charge() / cluster.size();
   filter_ = false;
   isSaturated_ = false;
   peakFilter_ = peakFilter;
+  version_ = version;
 
   //mimicing the algorithm used in StripSubClusterShapeTrajectoryFilter...
   //Looks for 3 adjacent saturated strips (ADC>=254)
@@ -54,4 +59,28 @@ SiStripApproximateCluster::SiStripApproximateCluster(const SiStripCluster& clust
   } else {
     filter_ = peakFilter_;
   }
+
+  if (version_ == 2) {
+    // Compression of avgCharge_ to integer
+    avgCharge_ = floor((float(cluster.charge()) / cluster.size()) / avgChargeScale_);
+    // In v2, we encode the filter_ and peakFilter_ info in avgCharge_ as the two highest bits
+    assert(avgCharge_ <= avgChargeRangeMax_ && "Setting avgCharge > avgChargeRangeMax_");
+    // filter_ and peakFilter_ are encoded in the two highest bits of avgCharge_
+    avgCharge_ = (avgCharge_ | (filter_ << kFilterBit));
+    assert(avgCharge_ <= ((1 << (nbits_avgCharge_ - 1)) - 1) && "Setting avgCharge with filter > max single bit range");
+    avgCharge_ = (avgCharge_ | (peakFilter_ << kPeakFilterBit));
+    assert(avgCharge_ <= ((1 << nbits_avgCharge_) - 1) && "Setting avgCharge with peakFilter > max full range");
+
+    // Compression of barycenter_ to integer [note: it contains the distance from the previous cluster]
+    barycenter_ = round(float(cluster.barycenter() - previous_barycenter + (offset_module_change)) * barycenterScale_);
+    assert(barycenter_ <= barycenterRangeMax_ && "Setting barycenter > barycenterRangeMax_");
+    // isSaturated_ is encoded in the highest bit of barycenter_
+    barycenter_ = (barycenter_ | (isSaturated_ << kSaturatedBit));
+    assert(barycenter_ <= ((1 << nbits_barycenter_) - 1) && "Setting barycenter with isSaturated > max full range");
+
+    // Flags set to false to reduce event size (they should be removed when moving to v2 only)
+    filter_ = false;
+    isSaturated_ = false;
+    peakFilter_ = false;
+  }
 }

DataFormats/SiStripCluster/src/SiStripApproximateClusterCollectionV2.cc

@@ -0,0 +1,10 @@
+#include "DataFormats/SiStripCluster/interface/SiStripApproximateClusterCollectionV2.h"
+void SiStripApproximateClusterCollectionV2::reserve(std::size_t dets, std::size_t clusters) {
+  detIds_.reserve(dets);
+  clusters_.reserve(clusters);
+}
+
+SiStripApproximateClusterCollectionV2::Filler SiStripApproximateClusterCollectionV2::beginDet(unsigned int detId) {
+  detIds_.push_back((detIds_.empty()) ? detId : detId - (std::accumulate(detIds_.cbegin(), detIds_.cend(), 0)));
+  return Filler(clusters_);
+}