Skip to content

Commit

Permalink
Merge pull request #47061 from slava77/CMSSW_14_1_0_pre7/stripsTiming
Browse files Browse the repository at this point in the history 
speedup SiStripClusterizer(FromRaw) using ThreeThresholdAlgorithm
  • Loading branch information
@cmsbuild
cmsbuild authored Jan 9, 2025
2 parents 6eea5fa + 49e294b commit d288c81
Show file tree
Hide file tree
Showing 4 changed files with 184 additions and 139 deletions.
19 changes: 15 additions & 4 deletions CalibFormats/SiStripObjects/interface/SiStripClusterizerConditions.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<bool, kMaxStrips> qualityBits = {};
SiStripNoises::Range noiseRange;
std::array<uint16_t, kMaxStrips> rawNoises = {};
float m_weight[6];
uint32_t detId = 0;
unsigned short ind = invalidI;
Expand Down Expand Up @@ -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];
}
Expand Down
116 changes: 116 additions & 0 deletions RecoLocalTracker/SiStripClusterizer/interface/ThreeThresholdAlgorithm.h
View file
Original file line number Diff line number Diff line change
@@ -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 <cmath>
#include <numeric>

class ThreeThresholdAlgorithm final : public StripClusterizerAlgorithm {
friend class StripClusterizerAlgorithmFactory;

Expand Down Expand Up @@ -64,4 +72,112 @@ class ThreeThresholdAlgorithm final : public StripClusterizerAlgorithm {
float minGoodCharge;
};

template <class digiDetSet>
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<uint8_t>(Noise * ChannelThreshold) || state.det().bad(strip))
return;

if (state.candidateLacksSeed)
state.candidateLacksSeed = adc < static_cast<uint8_t>(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<SiStripDigi>& digis,
output_t::TSFastFiller& output) const {
clusterizeDetUnit_(digis, output);
}

template <class T>
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
74 changes: 53 additions & 21 deletions RecoLocalTracker/SiStripClusterizer/plugins/ClustersFromRawProducer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -89,10 +90,11 @@ namespace {
return buffer;
}

template <class AlgoT>
class ClusterFiller final : public StripClusterizerAlgorithm::output_t::Getter {
public:
ClusterFiller(const FEDRawDataCollection& irawColl,
StripClusterizerAlgorithm& iclusterizer,
const AlgoT& iclusterizer,
SiStripRawProcessingAlgorithms& irawAlgos,
bool idoAPVEmulatorCheck,
bool legacy,
Expand All @@ -119,7 +121,7 @@ namespace {

const FEDRawDataCollection& rawColl;

StripClusterizerAlgorithm& clusterizer;
const AlgoT& clusterizer;
const SiStripClusterizerConditions& conditions;
SiStripRawProcessingAlgorithms& rawAlgos;

Expand Down Expand Up @@ -181,7 +183,7 @@ class SiStripClusterizerFromRaw final : public edm::stream::EDProducer<> {
legacy_(conf.getParameter<bool>("LegacyUnpacker")),
hybridZeroSuppressed_(conf.getParameter<bool>("HybridZeroSuppressed")) {
productToken_ = consumes<FEDRawDataCollection>(conf.getParameter<edm::InputTag>("ProductLabel"));
produces<edmNew::DetSetVector<SiStripCluster> >();
produces<edmNew::DetSetVector<SiStripCluster>>();
assert(clusterizer_.get());
assert(rawAlgos_.get());
}
Expand All @@ -193,12 +195,23 @@ class SiStripClusterizerFromRaw final : public edm::stream::EDProducer<> {
edm::Handle<FEDRawDataCollection> rawData;
ev.getByToken(productToken_, rawData);

std::unique_ptr<edmNew::DetSetVector<SiStripCluster> > output(
onDemand ? new edmNew::DetSetVector<SiStripCluster>(
std::shared_ptr<edmNew::DetSetVector<SiStripCluster>::Getter>(std::make_shared<ClusterFiller>(
*rawData, *clusterizer_, *rawAlgos_, doAPVEmulatorCheck_, legacy_, hybridZeroSuppressed_)),
clusterizer_->conditions().allDetIds())
: new edmNew::DetSetVector<SiStripCluster>());
const ThreeThresholdAlgorithm* clusterizer3 = dynamic_cast<const ThreeThresholdAlgorithm*>(clusterizer_.get());
std::unique_ptr<edmNew::DetSetVector<SiStripCluster>> output;
if (onDemand) {
if (clusterizer3 == nullptr)
output = std::make_unique<edmNew::DetSetVector<SiStripCluster>>(edmNew::DetSetVector<SiStripCluster>(
std::shared_ptr<edmNew::DetSetVector<SiStripCluster>::Getter>(
std::make_shared<ClusterFiller<StripClusterizerAlgorithm>>(
*rawData, *clusterizer_, *rawAlgos_, doAPVEmulatorCheck_, legacy_, hybridZeroSuppressed_)),
clusterizer_->conditions().allDetIds()));
else
output = std::make_unique<edmNew::DetSetVector<SiStripCluster>>(edmNew::DetSetVector<SiStripCluster>(
std::shared_ptr<edmNew::DetSetVector<SiStripCluster>::Getter>(
std::make_shared<ClusterFiller<ThreeThresholdAlgorithm>>(
*rawData, *clusterizer3, *rawAlgos_, doAPVEmulatorCheck_, legacy_, hybridZeroSuppressed_)),
clusterizer_->conditions().allDetIds()));
} else
output = std::make_unique<edmNew::DetSetVector<SiStripCluster>>(edmNew::DetSetVector<SiStripCluster>());

if (onDemand)
assert(output->onDemand());
Expand Down Expand Up @@ -266,20 +279,38 @@ void SiStripClusterizerFromRaw::initialize(const edm::EventSetup& es) {
}

void SiStripClusterizerFromRaw::run(const FEDRawDataCollection& rawColl, edmNew::DetSetVector<SiStripCluster>& output) {
ClusterFiller filler(rawColl, *clusterizer_, *rawAlgos_, doAPVEmulatorCheck_, legacy_, hybridZeroSuppressed_);
const ThreeThresholdAlgorithm* clusterizer3 = dynamic_cast<const ThreeThresholdAlgorithm*>(clusterizer_.get());
if (clusterizer3 == nullptr) {
ClusterFiller<StripClusterizerAlgorithm> 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<ThreeThresholdAlgorithm> 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 AlgoT>
class StripByStripAdder {
public:
typedef std::output_iterator_tag iterator_category;
Expand All @@ -288,7 +319,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) {}
Expand All @@ -303,7 +334,7 @@ namespace {
StripByStripAdder& operator++(int) { return *this; }

private:
StripClusterizerAlgorithm& clusterizer_;
const AlgoT& clusterizer_;
StripClusterizerAlgorithm::State& state_;
StripClusterizerAlgorithm::output_t::TSFastFiller& record_;
};
Expand All @@ -326,7 +357,8 @@ namespace {
};
} // namespace

void ClusterFiller::fill(StripClusterizerAlgorithm::output_t::TSFastFiller& record) const {
template <class AlgoT>
void ClusterFiller<AlgoT>::fill(StripClusterizerAlgorithm::output_t::TSFastFiller& record) const {
try { // edmNew::CapacityExaustedException
incReady();

Expand Down Expand Up @@ -393,7 +425,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<AlgoT>(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;
Expand Down
Loading

0 comments on commit d288c81

Please sign in to comment.