Modification of existing rawprime version of SiStripApproxCluster

DataFormats/SiStripCluster/interface/SiStripApproximateClusterCollection_v1.h

@@ -0,0 +1,107 @@
+#ifndef DataFormats_SiStripCluster_SiStripApproximateClusterCollection_v1_h
+#define DataFormats_SiStripCluster_SiStripApproximateClusterCollection_v1_h
+
+#include <vector>
+
+#include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster_v1.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.
+**/
+namespace v1 {
+  class SiStripApproximateClusterCollection {
+  public:
+    // Helper classes to make creation and iteration easier
+    class Filler {
+    public:
+      void push_back(SiStripApproximateCluster_v1 const& cluster) { clusters_.push_back(cluster); }
+
+    private:
+      friend SiStripApproximateClusterCollection;
+      Filler(std::vector<SiStripApproximateCluster_v1>& clusters) : clusters_(clusters) {}
+
+      std::vector<SiStripApproximateCluster_v1>& clusters_;
+    };
+
+    class const_iterator;
+    class DetSet {
+    public:
+      using const_iterator = std::vector<SiStripApproximateCluster_v1>::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 SiStripApproximateClusterCollection::const_iterator;
+      DetSet(SiStripApproximateClusterCollection const* coll, unsigned int detIndex)
+          : coll_(coll), detIndex_(detIndex), clusEnd_(coll->clusters_.size()) {}
+
+      SiStripApproximateClusterCollection 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 SiStripApproximateClusterCollection;
+      // default-constructed object acts as the sentinel
+      const_iterator() = default;
+      const_iterator(SiStripApproximateClusterCollection const* coll) : coll_(coll) {}
+
+      SiStripApproximateClusterCollection const* coll_ = nullptr;
+      unsigned int index_ = 0;
+    };
+
+    // Actual public interface
+    SiStripApproximateClusterCollection() = 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_v1> clusters_;
+  };
+}  // namespace v1
+#endif

DataFormats/SiStripCluster/interface/SiStripApproximateCluster_v1.h

@@ -0,0 +1,67 @@
+#ifndef DataFormats_SiStripCluster_SiStripApproximateCluster_v1_h
+#define DataFormats_SiStripCluster_SiStripApproximateCluster_v1_h
+
+#include "FWCore/Utilities/interface/typedefs.h"
+#include "assert.h"
+
+class SiStripCluster;
+class SiStripApproximateCluster_v1 {
+public:
+  SiStripApproximateCluster_v1() {}
+
+  explicit SiStripApproximateCluster_v1(cms_uint16_t compBarycenter, cms_uint8_t width, cms_uint8_t compavgCharge)
+      : compBarycenter_(compBarycenter), width_(width), compavgCharge_(compavgCharge) {}
+
+  explicit SiStripApproximateCluster_v1(const SiStripCluster& cluster,
+                                        unsigned int maxNSat,
+                                        float hitPredPos,
+                                        float& previous_cluster,
+                                        unsigned int& module_length,
+                                        unsigned int& previous_module_length,
+                                        bool peakFilter);
+
+  const cms_uint16_t compBarycenter() const { return compBarycenter_; }
+
+  float barycenter(float previous_barycenter = 0,
+                   unsigned int module_length = 0,
+                   unsigned int previous_module_length = 0) const {
+    float _barycenter;
+    cms_uint16_t compBarycenter = (compBarycenter_ & 0x7FFF);
+    if (previous_barycenter == -999)
+      _barycenter = compBarycenter * maxBarycenter_ / maxRange_;
+    else {
+      _barycenter = ((compBarycenter * maxBarycenter_ / maxRange_) - (module_length - previous_module_length)) +
+                    previous_barycenter;
+    }
+    assert(_barycenter <= maxBarycenter_ && "Returning barycenter > maxBarycenter");
+    return _barycenter;
+  }
+  cms_uint8_t width() const { return width_; }
+  float avgCharge() const {
+    cms_uint8_t compavgCharge = (compavgCharge_ & 0x3F);
+    float avgCharge_ = compavgCharge * maxavgCharge_ / maxavgChargeRange_;
+    assert(avgCharge_ <= maxavgCharge_ && "Returning avgCharge > maxavgCharge");
+    return avgCharge_;
+  }
+  bool filter() const { return (compavgCharge_ & (1 << kfilterMask)); }
+  bool isSaturated() const { return (compavgCharge_ & (1 << kSaturatedMask)); }
+  bool peakFilter() const { return (compBarycenter_ & (1 << kpeakFilterMask)); }
+
+private:
+  cms_uint16_t compBarycenter_ = 0;
+  cms_uint8_t width_ = 0;
+  cms_uint8_t compavgCharge_ = 0;
+  static constexpr double maxRange_ = 32767;
+  static constexpr double maxBarycenter_ = 1536.;
+  static constexpr double maxavgChargeRange_ = 63;
+  static constexpr double maxavgCharge_ = 255.;
+  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.;
+  static constexpr int kfilterMask = 6;
+  static constexpr int kpeakFilterMask = 7;
+  static constexpr int kSaturatedMask = 15;
+};
+#endif  // DataFormats_SiStripCluster_SiStripApproximateCluster_v1_h

DataFormats/SiStripCluster/interface/SiStripCluster.h

@@ -3,6 +3,7 @@
 
 #include "DataFormats/SiStripDigi/interface/SiStripDigi.h"
 #include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h"
+#include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster_v1.h"
 #include <vector>
 #include <numeric>
 #include <iostream>
@@ -45,6 +46,11 @@ class SiStripCluster {
   }
 
   SiStripCluster(const SiStripApproximateCluster cluster, const uint16_t maxStrips);
+  SiStripCluster(const SiStripApproximateCluster_v1 cluster,
+                 const uint16_t maxStrips,
+                 float pc = -999,
+                 unsigned int module_length = 0,
+                 unsigned int previous_module_length = 0);
 
   // extend the cluster
   template <typename Iter>

DataFormats/SiStripCluster/src/SiStripApproximateClusterCollection_v1.cc

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

DataFormats/SiStripCluster/src/SiStripApproximateCluster_v1.cc

@@ -0,0 +1,76 @@
+#include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster_v1.h"
+#include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+
+SiStripApproximateCluster_v1::SiStripApproximateCluster_v1(const SiStripCluster& cluster,
+                                                           unsigned int maxNSat,
+                                                           float hitPredPos,
+                                                           float& previous_cluster,
+                                                           unsigned int& module_length,
+                                                           unsigned int& previous_module_length,
+                                                           bool peakFilter) {
+  bool filter_, isSaturated_, peakFilter_;
+  if (previous_cluster == -999.)
+    compBarycenter_ = std::round(cluster.barycenter() * maxRange_ / maxBarycenter_);
+  else
+    compBarycenter_ =
+        std::round(((cluster.barycenter() - previous_cluster) + (module_length - previous_module_length)) * maxRange_ /
+                   maxBarycenter_);
+  previous_cluster = barycenter(previous_cluster, module_length, previous_module_length);
+  assert(cluster.barycenter() <= maxBarycenter_ && "Got a barycenter > maxBarycenter");
+  assert(compBarycenter_ <= maxRange_ && "Filling compBarycenter > maxRange");
+  width_ = std::min(255, (int)cluster.size());  //cluster.size();
+  float avgCharge_ = cluster.charge() * 1. / width_;
+  assert(avgCharge_ <= maxavgCharge_ && "Got a avgCharge > maxavgCharge");
+  compavgCharge_ = std::round(avgCharge_ * maxavgChargeRange_ / maxavgCharge_);
+  assert(compavgCharge_ <= maxavgChargeRange_ && "Filling compavgCharge > maxavgChargeRange");
+  filter_ = false;
+  isSaturated_ = false;
+  peakFilter_ = peakFilter;
+
+  //mimicing the algorithm used in StripSubClusterShapeTrajectoryFilter...
+  //Looks for 3 adjacent saturated strips (ADC>=254)
+  const auto& ampls = cluster.amplitudes();
+  unsigned int thisSat = (ampls[0] >= 254), maxSat = thisSat;
+  for (unsigned int i = 1, n = ampls.size(); i < n; ++i) {
+    if (ampls[i] >= 254) {
+      thisSat++;
+    } else if (thisSat > 0) {
+      maxSat = std::max<int>(maxSat, thisSat);
+      thisSat = 0;
+    }
+  }
+  if (thisSat > 0) {
+    maxSat = std::max<int>(maxSat, thisSat);
+  }
+  if (maxSat >= maxNSat) {
+    filter_ = true;
+    isSaturated_ = true;
+  }
+
+  unsigned int hitStripsTrim = ampls.size();
+  int sum = std::accumulate(ampls.begin(), ampls.end(), 0);
+  uint8_t trimCut = std::min<uint8_t>(trimMaxADC_, std::floor(trimMaxFracTotal_ * sum));
+  auto begin = ampls.begin();
+  auto last = ampls.end() - 1;
+  while (hitStripsTrim > 1 && (*begin < std::max<uint8_t>(trimCut, trimMaxFracNeigh_ * (*(begin + 1))))) {
+    hitStripsTrim--;
+    ++begin;
+  }
+  while (hitStripsTrim > 1 && (*last < std::max<uint8_t>(trimCut, trimMaxFracNeigh_ * (*(last - 1))))) {
+    hitStripsTrim--;
+    --last;
+  }
+  if (hitStripsTrim < std::floor(std::abs(hitPredPos) - maxTrimmedSizeDiffNeg_)) {
+    filter_ = false;
+  } else if (hitStripsTrim <= std::ceil(std::abs(hitPredPos) + maxTrimmedSizeDiffPos_)) {
+    filter_ = true;
+  } else {
+    filter_ = peakFilter_;
+  }
+  compavgCharge_ = (compavgCharge_ | (filter_ << kfilterMask));
+  compavgCharge_ = (compavgCharge_ | (peakFilter_ << kpeakFilterMask));
+  compBarycenter_ = (compBarycenter_ | (isSaturated_ << kSaturatedMask));
+}

DataFormats/SiStripCluster/src/SiStripCluster.cc

@@ -39,3 +39,25 @@ SiStripCluster::SiStripCluster(const SiStripApproximateCluster cluster, const ui
   }
   firstStrip_ |= approximateMask;
 }
+
+SiStripCluster::SiStripCluster(const SiStripApproximateCluster_v1 cluster,
+                               const uint16_t maxStrips,
+                               float p_bc,
+                               unsigned int module_length,
+                               unsigned int previous_module_length)
+    : error_x(-99999.9) {
+  barycenter_ = cluster.barycenter(p_bc, module_length, previous_module_length);
+  charge_ = cluster.width() * cluster.avgCharge();
+  amplitudes_.resize(cluster.width(), cluster.avgCharge());
+  filter_ = cluster.filter();
+
+  float halfwidth_ = 0.5f * float(cluster.width());
+
+  //initialize firstStrip_
+  firstStrip_ = std::max(barycenter_ - halfwidth_, 0.f);
+
+  if UNLIKELY (firstStrip_ + cluster.width() > maxStrips) {
+    firstStrip_ = maxStrips - cluster.width();
+  }
+  firstStrip_ |= approximateMask;
+}

DataFormats/SiStripCluster/src/classes.h

@@ -7,7 +7,9 @@
 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
 #include "DataFormats/SiStripCluster/interface/SiStripClustersSOA.h"
 #include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h"
+#include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster_v1.h"
 #include "DataFormats/SiStripCluster/interface/SiStripApproximateClusterCollection.h"
+#include "DataFormats/SiStripCluster/interface/SiStripApproximateClusterCollection_v1.h"
 #include "DataFormats/Common/interface/ContainerMask.h"
 
 #endif  // SISTRIPCLUSTER_CLASSES_H

DataFormats/SiStripCluster/src/classes_def.xml

@@ -59,6 +59,36 @@
  <class name="edmNew::DetSetVector<edm::Ref<edmNew::DetSetVector<SiStripApproximateCluster>,SiStripApproximateCluster,edmNew::DetSetVector<SiStripApproximateCluster>::FindForDetSetVector> >" />
  <class name="edm::Wrapper<edmNew::DetSetVector<edm::Ref<edmNew::DetSetVector<SiStripApproximateCluster>,SiStripApproximateCluster,edmNew::DetSetVector<SiStripApproximateCluster>::FindForDetSetVector> > >" />
 
+
+
+					 <class name="SiStripApproximateCluster_v1" ClassVersion="7">
+						 <version ClassVersion="7" checksum="1154754493"/>
+  <version ClassVersion="6" checksum="132211472"/>
+  <version ClassVersion="5" checksum="3495825183"/>
+  <version ClassVersion="4" checksum="2854791577"/>
+  <version ClassVersion="3" checksum="2041370183"/>
+ </class>
+ <class name="v1::SiStripApproximateClusterCollection" ClassVersion="4">
+	<version ClassVersion="4" checksum="2896589077"/>
+  <version ClassVersion="3" checksum="3101417750"/>
+ </class>
+ <class name="edm::Wrapper<v1::SiStripApproximateClusterCollection>"/>
+
+ <class name="edmNew::DetSetVector<SiStripApproximateCluster_v1>"/>
+ <class name="edm::Wrapper<edmNew::DetSetVector<SiStripApproximateCluster_v1>>"/>
+
+ <class name="std::vector<SiStripApproximateCluster_v1>"/>
+
+ <class name="edm::Ref<edmNew::DetSetVector<SiStripApproximateCluster_v1>, SiStripApproximateCluster_v1, edmNew::DetSetVector<SiStripApproximateCluster_v1>::FindForDetSetVector>"/>
+
+
+ <class name="edm::ContainerMask<edmNew::DetSetVector<SiStripApproximateCluster_v1> >"/>
+ <class name="edm::Wrapper<edm::ContainerMask<edmNew::DetSetVector<SiStripApproximateCluster_v1> > >"/>
+
+
+ <class name="std::vector<edm::Ref<edmNew::DetSetVector<SiStripApproximateCluster_v1>,SiStripApproximateCluster_v1,edmNew::DetSetVector<SiStripApproximateCluster_v1>::FindForDetSetVector> >" />
+ <class name="edmNew::DetSetVector<edm::Ref<edmNew::DetSetVector<SiStripApproximateCluster_v1>,SiStripApproximateCluster_v1,edmNew::DetSetVector<SiStripApproximateCluster_v1>::FindForDetSetVector> >" />
+ <class name="edm::Wrapper<edmNew::DetSetVector<edm::Ref<edmNew::DetSetVector<SiStripApproximateCluster_v1>,SiStripApproximateCluster_v1,edmNew::DetSetVector<SiStripApproximateCluster_v1>::FindForDetSetVector> > >" />
  <class name="SiStripClustersSOA" ClassVersion="3">
   <version ClassVersion="3" checksum="2739562998"/>
  </class>