diff --git a/CalibFormats/SiStripObjects/interface/SiStripClusterizerConditions.h b/CalibFormats/SiStripObjects/interface/SiStripClusterizerConditions.h index 35d07be8c8587..9323e81474f1a 100644 --- a/CalibFormats/SiStripObjects/interface/SiStripClusterizerConditions.h +++ b/CalibFormats/SiStripObjects/interface/SiStripClusterizerConditions.h @@ -14,18 +14,21 @@ class SiStripClusterizerConditions { struct Det { bool valid() const { return ind != invalidI; } - uint16_t rawNoise(const uint16_t strip) const { return SiStripNoises::getRawNoise(strip, noiseRange); } - float noise(const uint16_t strip) const { return SiStripNoises::getNoise(strip, noiseRange); } + uint16_t rawNoise(const uint16_t strip) const { return rawNoises[strip]; } + float noise(const uint16_t strip) const { return 0.1f * float(rawNoises[strip]); } float weight(const uint16_t strip) const { return m_weight[strip / 128]; } - bool bad(const uint16_t strip) const { return quality->IsStripBad(qualityRange, strip); } + bool bad(const uint16_t strip) const { return qualityBits[strip]; } bool allBadBetween(uint16_t L, const uint16_t& R) const { while (++L < R && bad(L)) { }; return L == R; } + static constexpr uint16_t kMaxStrips = 768; SiStripQuality const* quality; - SiStripNoises::Range noiseRange; SiStripQuality::Range qualityRange; + std::array qualityBits = {}; + SiStripNoises::Range noiseRange; + std::array rawNoises = {}; float m_weight[6]; uint32_t detId = 0; unsigned short ind = invalidI; @@ -72,7 +75,15 @@ class SiStripClusterizerConditions { auto& det = m_dets.emplace_back(); det.quality = m_quality; det.qualityRange = qualityRange; + for (uint16_t s = 0U; s < det.kMaxStrips; ++s) + det.qualityBits[s] = m_quality->IsStripBad(qualityRange, s); det.noiseRange = noiseRange; + auto maxRange8 = (noiseRange.second - noiseRange.first) * 8; + for (uint16_t s = 0U; s < det.kMaxStrips; ++s) { + if (9 * s >= maxRange8) + break; + det.rawNoises[s] = SiStripNoises::getRawNoise(s, noiseRange); + } for (uint32_t i = 0; i != invGains.size(); ++i) { det.m_weight[i] = invGains[i]; } diff --git a/RecoLocalTracker/SiStripClusterizer/interface/ThreeThresholdAlgorithm.h b/RecoLocalTracker/SiStripClusterizer/interface/ThreeThresholdAlgorithm.h index 2c4551ef40550..268aac8451306 100644 --- a/RecoLocalTracker/SiStripClusterizer/interface/ThreeThresholdAlgorithm.h +++ b/RecoLocalTracker/SiStripClusterizer/interface/ThreeThresholdAlgorithm.h @@ -1,8 +1,16 @@ #ifndef RecoLocalTracker_SiStripClusterizer_ThreeThresholdAlgorithm_h #define RecoLocalTracker_SiStripClusterizer_ThreeThresholdAlgorithm_h + +#include "DataFormats/SiStripCluster/interface/SiStripCluster.h" +#include "DataFormats/SiStripCluster/interface/SiStripClusterTools.h" +#include "DataFormats/SiStripDigi/interface/SiStripDigi.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" #include "RecoLocalTracker/SiStripClusterizer/interface/StripClusterizerAlgorithm.h" #include "RecoLocalTracker/SiStripClusterizer/interface/SiStripApvShotCleaner.h" +#include +#include + class ThreeThresholdAlgorithm final : public StripClusterizerAlgorithm { friend class StripClusterizerAlgorithmFactory; @@ -64,4 +72,112 @@ class ThreeThresholdAlgorithm final : public StripClusterizerAlgorithm { float minGoodCharge; }; +template +inline void ThreeThresholdAlgorithm::clusterizeDetUnit_(const digiDetSet& digis, output_t::TSFastFiller& output) const { + const auto& cond = conditions(); + if (cond.isModuleBad(digis.detId())) + return; + + auto const& det = cond.findDetId(digis.detId()); + if (!det.valid()) + return; + +#ifdef EDM_ML_DEBUG + if (!cond.isModuleUsable(digis.detId())) + edm::LogWarning("ThreeThresholdAlgorithm") << " id " << digis.detId() << " not usable???" << std::endl; +#endif + + typename digiDetSet::const_iterator scan(digis.begin()), end(digis.end()); + + SiStripApvShotCleaner ApvCleaner; + if (RemoveApvShots) { + ApvCleaner.clean(digis, scan, end); + } + + output.reserve(16); + State state(det); + while (scan != end) { + while (scan != end && !candidateEnded(state, scan->strip())) + addToCandidate(state, *scan++); + endCandidate(state, output); + } +} + +inline bool ThreeThresholdAlgorithm::candidateEnded(State const& state, const uint16_t& testStrip) const { + uint16_t holes = testStrip - state.lastStrip - 1; + return (((!state.ADCs.empty()) & // a candidate exists, and + (holes > MaxSequentialHoles) // too many holes if not all are bad strips, and + ) && + (holes > MaxSequentialBad || // (too many bad strips anyway, or + !state.det().allBadBetween(state.lastStrip, testStrip) // not all holes are bad strips) + )); +} + +inline void ThreeThresholdAlgorithm::addToCandidate(State& state, uint16_t strip, uint8_t adc) const { + float Noise = state.det().noise(strip); + if (adc < static_cast(Noise * ChannelThreshold) || state.det().bad(strip)) + return; + + if (state.candidateLacksSeed) + state.candidateLacksSeed = adc < static_cast(Noise * SeedThreshold); + if (state.ADCs.empty()) + state.lastStrip = strip - 1; // begin candidate + while (++state.lastStrip < strip) + state.ADCs.push_back(0); // pad holes + + if (state.ADCs.size() <= MaxClusterSize) + state.ADCs.push_back(adc); + state.noiseSquared += Noise * Noise; +} + +inline void ThreeThresholdAlgorithm::clusterizeDetUnit(const edmNew::DetSet& digis, + output_t::TSFastFiller& output) const { + clusterizeDetUnit_(digis, output); +} + +template +inline void ThreeThresholdAlgorithm::endCandidate(State& state, T& out) const { + if (candidateAccepted(state)) { + applyGains(state); + if (MaxAdjacentBad > 0) + appendBadNeighbors(state); + if (minGoodCharge <= 0 || + siStripClusterTools::chargePerCM(state.det().detId, state.ADCs.begin(), state.ADCs.end()) > minGoodCharge) + out.push_back(std::move(SiStripCluster(firstStrip(state), state.ADCs.begin(), state.ADCs.end()))); + } + clearCandidate(state); +} + +inline bool ThreeThresholdAlgorithm::candidateAccepted(State const& state) const { + return (!state.candidateLacksSeed && state.ADCs.size() <= MaxClusterSize && + state.noiseSquared * ClusterThresholdSquared <= + std::pow(float(std::accumulate(state.ADCs.begin(), state.ADCs.end(), int(0))), 2.f)); +} + +inline void ThreeThresholdAlgorithm::applyGains(State& state) const { + uint16_t strip = firstStrip(state); + for (auto& adc : state.ADCs) { +#ifdef EDM_ML_DEBUG + // if(adc > 255) throw InvalidChargeException( SiStripDigi(strip,adc) ); +#endif + // if(adc > 253) continue; //saturated, do not scale + auto charge = int(float(adc) * state.det().weight(strip++) + 0.5f); //adding 0.5 turns truncation into rounding + if (adc < 254) + adc = (charge > 1022 ? 255 : (charge > 253 ? 254 : charge)); + } +} + +inline void ThreeThresholdAlgorithm::appendBadNeighbors(State& state) const { + uint8_t max = MaxAdjacentBad; + while (0 < max--) { + if (state.det().bad(firstStrip(state) - 1)) { + state.ADCs.insert(state.ADCs.begin(), 0); + } + if (state.det().bad(state.lastStrip + 1)) { + state.ADCs.push_back(0); + state.lastStrip++; + } + } +} + #endif diff --git a/RecoLocalTracker/SiStripClusterizer/plugins/ClustersFromRawProducer.cc b/RecoLocalTracker/SiStripClusterizer/plugins/ClustersFromRawProducer.cc index 8622f0c32a929..7a02ba9887df7 100644 --- a/RecoLocalTracker/SiStripClusterizer/plugins/ClustersFromRawProducer.cc +++ b/RecoLocalTracker/SiStripClusterizer/plugins/ClustersFromRawProducer.cc @@ -4,6 +4,7 @@ #include "RecoLocalTracker/SiStripZeroSuppression/interface/SiStripRawProcessingFactory.h" #include "RecoLocalTracker/SiStripClusterizer/interface/StripClusterizerAlgorithm.h" +#include "RecoLocalTracker/SiStripClusterizer/interface/ThreeThresholdAlgorithm.h" #include "RecoLocalTracker/SiStripZeroSuppression/interface/SiStripRawProcessingAlgorithms.h" #include "DataFormats/SiStripCluster/interface/SiStripCluster.h" @@ -88,10 +89,11 @@ namespace { return buffer; } + template class ClusterFiller final : public StripClusterizerAlgorithm::output_t::Getter { public: ClusterFiller(const FEDRawDataCollection& irawColl, - StripClusterizerAlgorithm& iclusterizer, + const AlgoT& iclusterizer, SiStripRawProcessingAlgorithms& irawAlgos, bool idoAPVEmulatorCheck, bool legacy, @@ -118,7 +120,7 @@ namespace { const FEDRawDataCollection& rawColl; - StripClusterizerAlgorithm& clusterizer; + const AlgoT& clusterizer; const SiStripClusterizerConditions& conditions; SiStripRawProcessingAlgorithms& rawAlgos; @@ -180,7 +182,7 @@ class SiStripClusterizerFromRaw final : public edm::stream::EDProducer<> { legacy_(conf.getParameter("LegacyUnpacker")), hybridZeroSuppressed_(conf.getParameter("HybridZeroSuppressed")) { productToken_ = consumes(conf.getParameter("ProductLabel")); - produces >(); + produces>(); assert(clusterizer_.get()); assert(rawAlgos_.get()); } @@ -192,12 +194,23 @@ class SiStripClusterizerFromRaw final : public edm::stream::EDProducer<> { edm::Handle rawData; ev.getByToken(productToken_, rawData); - std::unique_ptr > output( - onDemand ? new edmNew::DetSetVector( - std::shared_ptr::Getter>(std::make_shared( - *rawData, *clusterizer_, *rawAlgos_, doAPVEmulatorCheck_, legacy_, hybridZeroSuppressed_)), - clusterizer_->conditions().allDetIds()) - : new edmNew::DetSetVector()); + const ThreeThresholdAlgorithm* clusterizer3 = dynamic_cast(clusterizer_.get()); + std::unique_ptr> output; + if (onDemand) { + if (clusterizer3 == nullptr) + output = std::make_unique>(edmNew::DetSetVector( + std::shared_ptr::Getter>( + std::make_shared>( + *rawData, *clusterizer_, *rawAlgos_, doAPVEmulatorCheck_, legacy_, hybridZeroSuppressed_)), + clusterizer_->conditions().allDetIds())); + else + output = std::make_unique>(edmNew::DetSetVector( + std::shared_ptr::Getter>( + std::make_shared>( + *rawData, *clusterizer3, *rawAlgos_, doAPVEmulatorCheck_, legacy_, hybridZeroSuppressed_)), + clusterizer_->conditions().allDetIds())); + } else + output = std::make_unique>(edmNew::DetSetVector()); if (onDemand) assert(output->onDemand()); @@ -265,20 +278,38 @@ void SiStripClusterizerFromRaw::initialize(const edm::EventSetup& es) { } void SiStripClusterizerFromRaw::run(const FEDRawDataCollection& rawColl, edmNew::DetSetVector& output) { - ClusterFiller filler(rawColl, *clusterizer_, *rawAlgos_, doAPVEmulatorCheck_, legacy_, hybridZeroSuppressed_); + const ThreeThresholdAlgorithm* clusterizer3 = dynamic_cast(clusterizer_.get()); + if (clusterizer3 == nullptr) { + ClusterFiller filler( + rawColl, *clusterizer_, *rawAlgos_, doAPVEmulatorCheck_, legacy_, hybridZeroSuppressed_); - // loop over good det in cabling - for (auto idet : clusterizer_->conditions().allDetIds()) { - StripClusterizerAlgorithm::output_t::TSFastFiller record(output, idet); + // loop over good det in cabling + for (auto idet : clusterizer_->conditions().allDetIds()) { + StripClusterizerAlgorithm::output_t::TSFastFiller record(output, idet); - filler.fill(record); + filler.fill(record); - if (record.empty()) - record.abort(); - } // end loop over dets + if (record.empty()) + record.abort(); + } // end loop over dets + } else { + ClusterFiller filler( + rawColl, *clusterizer3, *rawAlgos_, doAPVEmulatorCheck_, legacy_, hybridZeroSuppressed_); + + // loop over good det in cabling + for (auto idet : clusterizer_->conditions().allDetIds()) { + StripClusterizerAlgorithm::output_t::TSFastFiller record(output, idet); + + filler.fill(record); + + if (record.empty()) + record.abort(); + } // end loop over dets + } } namespace { + template class StripByStripAdder { public: typedef std::output_iterator_tag iterator_category; @@ -287,7 +318,7 @@ namespace { typedef void pointer; typedef void reference; - StripByStripAdder(StripClusterizerAlgorithm& clusterizer, + StripByStripAdder(const AlgoT& clusterizer, StripClusterizerAlgorithm::State& state, StripClusterizerAlgorithm::output_t::TSFastFiller& record) : clusterizer_(clusterizer), state_(state), record_(record) {} @@ -302,7 +333,7 @@ namespace { StripByStripAdder& operator++(int) { return *this; } private: - StripClusterizerAlgorithm& clusterizer_; + const AlgoT& clusterizer_; StripClusterizerAlgorithm::State& state_; StripClusterizerAlgorithm::output_t::TSFastFiller& record_; }; @@ -325,7 +356,8 @@ namespace { }; } // namespace -void ClusterFiller::fill(StripClusterizerAlgorithm::output_t::TSFastFiller& record) const { +template +void ClusterFiller::fill(StripClusterizerAlgorithm::output_t::TSFastFiller& record) const { try { // edmNew::CapacityExaustedException incReady(); @@ -392,7 +424,7 @@ void ClusterFiller::fill(StripClusterizerAlgorithm::output_t::TSFastFiller& reco using namespace sistrip; if LIKELY (fedchannelunpacker::isZeroSuppressed(mode, legacy_, lmode)) { - auto perStripAdder = StripByStripAdder(clusterizer, state, record); + auto perStripAdder = StripByStripAdder(clusterizer, state, record); const auto isNonLite = fedchannelunpacker::isNonLiteZS(mode, legacy_, lmode); const uint8_t pCode = (isNonLite ? buffer->packetCode(legacy_, fedCh) : 0); auto st_ch = fedchannelunpacker::StatusCode::SUCCESS; diff --git a/RecoLocalTracker/SiStripClusterizer/src/ThreeThresholdAlgorithm.cc b/RecoLocalTracker/SiStripClusterizer/src/ThreeThresholdAlgorithm.cc index c2cf5b8962896..329e6310993b5 100644 --- a/RecoLocalTracker/SiStripClusterizer/src/ThreeThresholdAlgorithm.cc +++ b/RecoLocalTracker/SiStripClusterizer/src/ThreeThresholdAlgorithm.cc @@ -1,11 +1,4 @@ #include "RecoLocalTracker/SiStripClusterizer/interface/ThreeThresholdAlgorithm.h" -#include "DataFormats/SiStripDigi/interface/SiStripDigi.h" -#include "DataFormats/SiStripCluster/interface/SiStripCluster.h" -#include -#include -#include "FWCore/MessageLogger/interface/MessageLogger.h" - -#include "DataFormats/SiStripCluster/interface/SiStripClusterTools.h" ThreeThresholdAlgorithm::ThreeThresholdAlgorithm( const edm::ESGetToken& conditionsToken, @@ -29,117 +22,10 @@ ThreeThresholdAlgorithm::ThreeThresholdAlgorithm( RemoveApvShots(removeApvShots), minGoodCharge(minGoodCharge) {} -template -inline void ThreeThresholdAlgorithm::clusterizeDetUnit_(const digiDetSet& digis, output_t::TSFastFiller& output) const { - const auto& cond = conditions(); - if (cond.isModuleBad(digis.detId())) - return; - - auto const& det = cond.findDetId(digis.detId()); - if (!det.valid()) - return; - -#ifdef EDM_ML_DEBUG - if (!cond.isModuleUsable(digis.detId())) - edm::LogWarning("ThreeThresholdAlgorithm") << " id " << digis.detId() << " not usable???" << std::endl; -#endif - - typename digiDetSet::const_iterator scan(digis.begin()), end(digis.end()); - - SiStripApvShotCleaner ApvCleaner; - if (RemoveApvShots) { - ApvCleaner.clean(digis, scan, end); - } - - output.reserve(16); - State state(det); - while (scan != end) { - while (scan != end && !candidateEnded(state, scan->strip())) - addToCandidate(state, *scan++); - endCandidate(state, output); - } -} - -inline bool ThreeThresholdAlgorithm::candidateEnded(State const& state, const uint16_t& testStrip) const { - uint16_t holes = testStrip - state.lastStrip - 1; - return (((!state.ADCs.empty()) & // a candidate exists, and - (holes > MaxSequentialHoles) // too many holes if not all are bad strips, and - ) && - (holes > MaxSequentialBad || // (too many bad strips anyway, or - !state.det().allBadBetween(state.lastStrip, testStrip) // not all holes are bad strips) - )); -} - -inline void ThreeThresholdAlgorithm::addToCandidate(State& state, uint16_t strip, uint8_t adc) const { - float Noise = state.det().noise(strip); - if (adc < static_cast(Noise * ChannelThreshold) || state.det().bad(strip)) - return; - - if (state.candidateLacksSeed) - state.candidateLacksSeed = adc < static_cast(Noise * SeedThreshold); - if (state.ADCs.empty()) - state.lastStrip = strip - 1; // begin candidate - while (++state.lastStrip < strip) - state.ADCs.push_back(0); // pad holes - - if (state.ADCs.size() <= MaxClusterSize) - state.ADCs.push_back(adc); - state.noiseSquared += Noise * Noise; -} - -template -inline void ThreeThresholdAlgorithm::endCandidate(State& state, T& out) const { - if (candidateAccepted(state)) { - applyGains(state); - if (MaxAdjacentBad > 0) - appendBadNeighbors(state); - if (minGoodCharge <= 0 || - siStripClusterTools::chargePerCM(state.det().detId, state.ADCs.begin(), state.ADCs.end()) > minGoodCharge) - out.push_back(std::move(SiStripCluster(firstStrip(state), state.ADCs.begin(), state.ADCs.end()))); - } - clearCandidate(state); -} - -inline bool ThreeThresholdAlgorithm::candidateAccepted(State const& state) const { - return (!state.candidateLacksSeed && state.ADCs.size() <= MaxClusterSize && - state.noiseSquared * ClusterThresholdSquared <= - std::pow(float(std::accumulate(state.ADCs.begin(), state.ADCs.end(), int(0))), 2.f)); -} - -inline void ThreeThresholdAlgorithm::applyGains(State& state) const { - uint16_t strip = firstStrip(state); - for (auto& adc : state.ADCs) { -#ifdef EDM_ML_DEBUG - // if(adc > 255) throw InvalidChargeException( SiStripDigi(strip,adc) ); -#endif - // if(adc > 253) continue; //saturated, do not scale - auto charge = int(float(adc) * state.det().weight(strip++) + 0.5f); //adding 0.5 turns truncation into rounding - if (adc < 254) - adc = (charge > 1022 ? 255 : (charge > 253 ? 254 : charge)); - } -} - -inline void ThreeThresholdAlgorithm::appendBadNeighbors(State& state) const { - uint8_t max = MaxAdjacentBad; - while (0 < max--) { - if (state.det().bad(firstStrip(state) - 1)) { - state.ADCs.insert(state.ADCs.begin(), 0); - } - if (state.det().bad(state.lastStrip + 1)) { - state.ADCs.push_back(0); - state.lastStrip++; - } - } -} - void ThreeThresholdAlgorithm::clusterizeDetUnit(const edm::DetSet& digis, output_t::TSFastFiller& output) const { clusterizeDetUnit_(digis, output); } -void ThreeThresholdAlgorithm::clusterizeDetUnit(const edmNew::DetSet& digis, - output_t::TSFastFiller& output) const { - clusterizeDetUnit_(digis, output); -} void ThreeThresholdAlgorithm::stripByStripAdd(State& state, uint16_t strip,