SiStripZeroSuppression/src/IteratedMedianCMNSubtractor.cc

0001 #include "RecoLocalTracker/SiStripZeroSuppression/interface/IteratedMedianCMNSubtractor.h"
0002
0003 #include "CondFormats/SiStripObjects/interface/SiStripNoises.h"
0004 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
0005 #include <cmath>
0006
0007 void IteratedMedianCMNSubtractor::init(const edm::EventSetup& es) {
0008   if (noiseWatcher_.check(es)) {
0009     noiseHandle = &es.getData(noiseToken_);
0010   }
0011   if (qualityWatcher_.check(es)) {
0012     qualityHandle = &es.getData(qualityToken_);
0013   }
0014 }
0015
0016 void IteratedMedianCMNSubtractor::subtract(uint32_t detId, uint16_t firstAPV, std::vector<int16_t>& digis) {
0017   subtract_(detId, firstAPV, digis);
0018 }
0019 void IteratedMedianCMNSubtractor::subtract(uint32_t detId, uint16_t firstAPV, std::vector<float>& digis) {
0020   subtract_(detId, firstAPV, digis);
0021 }
0022
0023 template <typename T>
0024 inline void IteratedMedianCMNSubtractor::subtract_(uint32_t detId, uint16_t firstAPV, std::vector<T>& digis) {
0025   SiStripNoises::Range detNoiseRange = noiseHandle->getRange(detId);
0026   SiStripQuality::Range detQualityRange = qualityHandle->getRange(detId);
0027
0028   typename std::vector<T>::iterator fs, ls;
0029   float offset = 0;
0030   std::vector<std::pair<float, float> > subset;
0031   subset.reserve(128);
0032
0033   _vmedians.clear();
0034
0035   uint16_t APV = firstAPV;
0036   for (; APV < digis.size() / 128 + firstAPV; ++APV) {
0037     subset.clear();
0038     // fill subset vector with all good strips and their noises
0039     for (uint16_t istrip = APV * 128; istrip < (APV + 1) * 128; ++istrip) {
0040       if (!qualityHandle->IsStripBad(detQualityRange, istrip)) {
0041         std::pair<float, float> pin((float)digis[istrip - firstAPV * 128],
0042                                     (float)noiseHandle->getNoiseFast(istrip, detNoiseRange));
0043         subset.push_back(pin);
0044       }
0045     }
0046
0047     // caluate offset for all good strips (first iteration)
0048     if (!subset.empty())
0049       offset = pairMedian(subset);
0050
0051     // for second, third... iterations, remove strips over threshold
0052     // and recalculate offset on remaining strips
0053     for (int ii = 0; ii < iterations_ - 1; ++ii) {
0054       std::vector<std::pair<float, float> >::iterator si = subset.begin();
0055       while (si != subset.end()) {
0056         if (si->first - offset > cut_to_avoid_signal_ * si->second)
0057           si = subset.erase(si);
0058         else
0059           ++si;
0060       }
0061       if (subset.empty())
0062         break;
0063       offset = pairMedian(subset);
0064     }
0065
0066     _vmedians.push_back(std::pair<short, float>(APV, offset));
0067
0068     // remove offset
0069     fs = digis.begin() + (APV - firstAPV) * 128;
0070     ls = digis.begin() + (APV - firstAPV + 1) * 128;
0071     while (fs < ls) {
0072       *fs = static_cast<T>(*fs - offset);
0073       fs++;
0074     }
0075   }
0076 }
0077
0078 inline float IteratedMedianCMNSubtractor::pairMedian(std::vector<std::pair<float, float> >& sample) {
0079   std::vector<std::pair<float, float> >::iterator mid = sample.begin() + sample.size() / 2;
0080   std::nth_element(sample.begin(), mid, sample.end());
0081   if (sample.size() & 1)  //odd size
0082     return (*mid).first;
0083   return ((*std::max_element(sample.begin(), mid)).first + (*mid).first) / 2.;
0084 }