SiStripClusterizer/src/ThreeThresholdAlgorithm.cc

0001 #include "RecoLocalTracker/SiStripClusterizer/interface/ThreeThresholdAlgorithm.h"
0002 #include "DataFormats/SiStripDigi/interface/SiStripDigi.h"
0003 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0004 #include <cmath>
0005 #include <numeric>
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008 #include "DataFormats/SiStripCluster/interface/SiStripClusterTools.h"
0009
0010 ThreeThresholdAlgorithm::ThreeThresholdAlgorithm(
0011     const edm::ESGetToken<SiStripClusterizerConditions, SiStripClusterizerConditionsRcd>& conditionsToken,
0012     float chan,
0013     float seed,
0014     float cluster,
0015     unsigned holes,
0016     unsigned bad,
0017     unsigned adj,
0018     unsigned maxClusterSize,
0019     bool removeApvShots,
0020     float minGoodCharge)
0021     : StripClusterizerAlgorithm(conditionsToken),
0022       ChannelThreshold(chan),
0023       SeedThreshold(seed),
0024       ClusterThresholdSquared(cluster * cluster),
0025       MaxSequentialHoles(holes),
0026       MaxSequentialBad(bad),
0027       MaxAdjacentBad(adj),
0028       MaxClusterSize(maxClusterSize),
0029       RemoveApvShots(removeApvShots),
0030       minGoodCharge(minGoodCharge) {}
0031
0032 template <class digiDetSet>
0033 inline void ThreeThresholdAlgorithm::clusterizeDetUnit_(const digiDetSet& digis, output_t::TSFastFiller& output) const {
0034   const auto& cond = conditions();
0035   if (cond.isModuleBad(digis.detId()))
0036     return;
0037
0038   auto const& det = cond.findDetId(digis.detId());
0039   if (!det.valid())
0040     return;
0041
0042 #ifdef EDM_ML_DEBUG
0043   if (!cond.isModuleUsable(digis.detId()))
0044     edm::LogWarning("ThreeThresholdAlgorithm") << " id " << digis.detId() << " not usable???" << std::endl;
0045 #endif
0046
0047   typename digiDetSet::const_iterator scan(digis.begin()), end(digis.end());
0048
0049   SiStripApvShotCleaner ApvCleaner;
0050   if (RemoveApvShots) {
0051     ApvCleaner.clean(digis, scan, end);
0052   }
0053
0054   output.reserve(16);
0055   State state(det);
0056   while (scan != end) {
0057     while (scan != end && !candidateEnded(state, scan->strip()))
0058       addToCandidate(state, *scan++);
0059     endCandidate(state, output);
0060   }
0061 }
0062
0063 inline bool ThreeThresholdAlgorithm::candidateEnded(State const& state, const uint16_t& testStrip) const {
0064   uint16_t holes = testStrip - state.lastStrip - 1;
0065   return (((!state.ADCs.empty()) &       // a candidate exists, and
0066            (holes > MaxSequentialHoles)  // too many holes if not all are bad strips, and
0067            ) &&
0068           (holes > MaxSequentialBad ||                             // (too many bad strips anyway, or
0069            !state.det().allBadBetween(state.lastStrip, testStrip)  // not all holes are bad strips)
0070            ));
0071 }
0072
0073 inline void ThreeThresholdAlgorithm::addToCandidate(State& state, uint16_t strip, uint8_t adc) const {
0074   float Noise = state.det().noise(strip);
0075   if (adc < static_cast<uint8_t>(Noise * ChannelThreshold) || state.det().bad(strip))
0076     return;
0077
0078   if (state.candidateLacksSeed)
0079     state.candidateLacksSeed = adc < static_cast<uint8_t>(Noise * SeedThreshold);
0080   if (state.ADCs.empty())
0081     state.lastStrip = strip - 1;  // begin candidate
0082   while (++state.lastStrip < strip)
0083     state.ADCs.push_back(0);  // pad holes
0084
0085   if (state.ADCs.size() <= MaxClusterSize)
0086     state.ADCs.push_back(adc);
0087   state.noiseSquared += Noise * Noise;
0088 }
0089
0090 template <class T>
0091 inline void ThreeThresholdAlgorithm::endCandidate(State& state, T& out) const {
0092   if (candidateAccepted(state)) {
0093     applyGains(state);
0094     if (MaxAdjacentBad > 0)
0095       appendBadNeighbors(state);
0096     if (minGoodCharge <= 0 ||
0097         siStripClusterTools::chargePerCM(state.det().detId, state.ADCs.begin(), state.ADCs.end()) > minGoodCharge)
0098       out.push_back(std::move(SiStripCluster(firstStrip(state), state.ADCs.begin(), state.ADCs.end())));
0099   }
0100   clearCandidate(state);
0101 }
0102
0103 inline bool ThreeThresholdAlgorithm::candidateAccepted(State const& state) const {
0104   return (!state.candidateLacksSeed && state.ADCs.size() <= MaxClusterSize &&
0105           state.noiseSquared * ClusterThresholdSquared <=
0106               std::pow(float(std::accumulate(state.ADCs.begin(), state.ADCs.end(), int(0))), 2.f));
0107 }
0108
0109 inline void ThreeThresholdAlgorithm::applyGains(State& state) const {
0110   uint16_t strip = firstStrip(state);
0111   for (auto& adc : state.ADCs) {
0112 #ifdef EDM_ML_DEBUG
0113     // if(adc > 255) throw InvalidChargeException( SiStripDigi(strip,adc) );
0114 #endif
0115     // if(adc > 253) continue; //saturated, do not scale
0116     auto charge = int(float(adc) * state.det().weight(strip++) + 0.5f);  //adding 0.5 turns truncation into rounding
0117     if (adc < 254)
0118       adc = (charge > 1022 ? 255 : (charge > 253 ? 254 : charge));
0119   }
0120 }
0121
0122 inline void ThreeThresholdAlgorithm::appendBadNeighbors(State& state) const {
0123   uint8_t max = MaxAdjacentBad;
0124   while (0 < max--) {
0125     if (state.det().bad(firstStrip(state) - 1)) {
0126       state.ADCs.insert(state.ADCs.begin(), 0);
0127     }
0128     if (state.det().bad(state.lastStrip + 1)) {
0129       state.ADCs.push_back(0);
0130       state.lastStrip++;
0131     }
0132   }
0133 }
0134
0135 void ThreeThresholdAlgorithm::clusterizeDetUnit(const edm::DetSet<SiStripDigi>& digis,
0136                                                 output_t::TSFastFiller& output) const {
0137   clusterizeDetUnit_(digis, output);
0138 }
0139 void ThreeThresholdAlgorithm::clusterizeDetUnit(const edmNew::DetSet<SiStripDigi>& digis,
0140                                                 output_t::TSFastFiller& output) const {
0141   clusterizeDetUnit_(digis, output);
0142 }
0143
0144 void ThreeThresholdAlgorithm::stripByStripAdd(State& state,
0145                                               uint16_t strip,
0146                                               uint8_t adc,
0147                                               std::vector<SiStripCluster>& out) const {
0148   if (candidateEnded(state, strip))
0149     endCandidate(state, out);
0150   addToCandidate(state, SiStripDigi(strip, adc));
0151 }
0152
0153 void ThreeThresholdAlgorithm::stripByStripEnd(State& state, std::vector<SiStripCluster>& out) const {
0154   endCandidate(state, out);
0155 }