PixelLowPtUtilities/interface/SlidingPeakFinder.h

0001 #ifndef RecoTracker_PixelLowPtUtilities_SlidingPeakFinder_h
0002 #define RecoTracker_PixelLowPtUtilities_SlidingPeakFinder_h
0003
0004 #include <algorithm>
0005 #include <vector>
0006 #include <cmath>
0007 #include <cstdint>
0008 #include "CondFormats/SiStripObjects/interface/SiStripNoises.h"
0009 #include "CondFormats/DataRecord/interface/SiStripNoisesRcd.h"
0010
0011 class SlidingPeakFinder {
0012 public:
0013   SlidingPeakFinder(unsigned int size) : size_(size), half_((size + 1) / 2) {}
0014
0015   template <typename Test>
0016   bool apply(const uint8_t *x,
0017              const uint8_t *begin,
0018              const uint8_t *end,
0019              const Test &test,
0020              bool verbose = false,
0021              int firststrip = 0) {
0022     const uint8_t *ileft = (x != begin) ? std::min_element(x - 1, x + half_) : begin - 1;
0023     const uint8_t *iright = ((x + size_) < end) ? std::min_element(x + half_, std::min(x + size_ + 1, end)) : end;
0024     uint8_t left = (ileft < begin ? 0 : *ileft);
0025     uint8_t right = (iright >= end ? 0 : *iright);
0026     uint8_t center = *std::max_element(x, std::min(x + size_, end));
0027     uint8_t maxmin = std::max(left, right);
0028     if (maxmin < center) {
0029       bool ret = test(center, maxmin);
0030       if (ret) {
0031         ret = test(ileft, iright, begin, end);
0032       }
0033       return ret;
0034     } else {
0035       return false;
0036     }
0037   }
0038
0039   template <typename V, typename Test>
0040   bool apply(const V &ampls, const Test &test, bool verbose = false, int firststrip = 0) {
0041     const uint8_t *begin = &*ampls.begin();
0042     const uint8_t *end = &*ampls.end();
0043     for (const uint8_t *x = begin; x < end - (half_ - 1); ++x) {
0044       if (apply(x, begin, end, test, verbose, firststrip)) {
0045         return true;
0046       }
0047     }
0048     return false;
0049   }
0050
0051 private:
0052   unsigned int size_, half_;
0053 };
0054
0055 struct PeakFinderTest {
0056   PeakFinderTest(float mip,
0057                  uint32_t detid,
0058                  uint32_t firstStrip,
0059                  const SiStripNoises *theNoise,
0060                  float seedCutMIPs,
0061                  float seedCutSN,
0062                  float subclusterCutMIPs,
0063                  float subclusterCutSN)
0064       : mip_(mip),
0065         detid_(detid),
0066         firstStrip_(firstStrip),
0067         noiseObj_(theNoise),
0068         noises_(theNoise->getRange(detid)),
0069         subclusterCutMIPs_(subclusterCutMIPs),
0070         sumCut_(subclusterCutMIPs_ * mip_),
0071         subclusterCutSN2_(subclusterCutSN * subclusterCutSN) {
0072     cut_ = std::min<float>(seedCutMIPs * mip, seedCutSN * noiseObj_->getNoise(firstStrip + 1, noises_));
0073   }
0074
0075   bool operator()(uint8_t max, uint8_t min) const { return max - min > cut_; }
0076   bool operator()(const uint8_t *left, const uint8_t *right, const uint8_t *begin, const uint8_t *end) const {
0077     int yleft = (left < begin ? 0 : *left);
0078     int yright = (right >= end ? 0 : *right);
0079     float sum = 0.0;
0080     int maxval = 0;
0081     float noise = 0;
0082     for (const uint8_t *x = left + 1; x < right; ++x) {
0083       int baseline = (yleft * int(right - x) + yright * int(x - left)) / int(right - left);
0084       sum += int(*x) - baseline;
0085       noise += std::pow(noiseObj_->getNoise(firstStrip_ + int(x - begin), noises_), 2);
0086       maxval = std::max(maxval, int(*x) - baseline);
0087     }
0088     if (sum > sumCut_ && sum * sum > noise * subclusterCutSN2_)
0089       return true;
0090     return false;
0091   }
0092
0093 private:
0094   float mip_;
0095   unsigned int detid_;
0096   int firstStrip_;
0097   const SiStripNoises *noiseObj_;
0098   SiStripNoises::Range noises_;
0099   uint8_t cut_;
0100   float subclusterCutMIPs_, sumCut_, subclusterCutSN2_;
0101 };
0102 #endif