File indexing completed on 2025-06-03 00:12:16
0001 #include "L1Trigger/TrackerTFP/interface/HoughTransform.h"
0002
0003 #include <numeric>
0004 #include <algorithm>
0005 #include <iterator>
0006 #include <deque>
0007 #include <vector>
0008 #include <cmath>
0009
0010 namespace trackerTFP {
0011
0012 HoughTransform::HoughTransform(const tt::Setup* setup,
0013 const DataFormats* dataFormats,
0014 const LayerEncoding* layerEncoding,
0015 std::vector<StubHT>& stubs)
0016 : setup_(setup),
0017 dataFormats_(dataFormats),
0018 layerEncoding_(layerEncoding),
0019 inv2R_(&dataFormats_->format(Variable::inv2R, Process::ht)),
0020 phiT_(&dataFormats_->format(Variable::phiT, Process::ht)),
0021 zT_(&dataFormats_->format(Variable::zT, Process::gp)),
0022 phi_(&dataFormats_->format(Variable::phi, Process::ht)),
0023 z_(&dataFormats_->format(Variable::z, Process::gp)),
0024 stubs_(stubs) {
0025 numChannelIn_ = dataFormats_->numChannel(Process::gp);
0026 numChannelOut_ = dataFormats_->numChannel(Process::ht);
0027 chan_ = setup_->kfNumWorker();
0028 mux_ = numChannelOut_ / chan_;
0029 }
0030
0031
0032 void HoughTransform::produce(const std::vector<std::vector<StubGP*>>& streamsIn,
0033 std::vector<std::deque<StubHT*>>& streamsOut) {
0034
0035 auto multiplicity = [](int sum, StubGP* s) { return sum + (s ? 1 + s->inv2RMax() - s->inv2RMin() : 0); };
0036 int nStubs(0);
0037 for (const std::vector<StubGP*>& input : streamsIn)
0038 nStubs += std::accumulate(input.begin(), input.end(), 0, multiplicity);
0039 stubs_.reserve(nStubs);
0040 for (int channelOut = 0; channelOut < numChannelOut_; channelOut++) {
0041 const int inv2Ru = mux_ * (channelOut % chan_) + channelOut / chan_;
0042 const int inv2R = inv2R_->toSigned(inv2Ru);
0043 std::deque<StubHT*>& output = streamsOut[channelOut];
0044 for (int channelIn = numChannelIn_ - 1; channelIn >= 0; channelIn--) {
0045 const std::vector<StubGP*>& input = streamsIn[channelIn];
0046 std::vector<StubHT*> stubs;
0047 stubs.reserve(2 * input.size());
0048
0049 fillIn(inv2R, channelIn, input, stubs);
0050
0051 if (setup_->enableTruncation() && static_cast<int>(stubs.size()) > setup_->numFramesHigh())
0052 stubs.resize(setup_->numFramesHigh());
0053
0054 stubs.erase(std::remove(stubs.begin(), stubs.end(), nullptr), stubs.end());
0055
0056 readOut(stubs, output);
0057 }
0058
0059 if (setup_->enableTruncation() && static_cast<int>(output.size()) > setup_->numFramesHigh())
0060 output.resize(setup_->numFramesHigh());
0061
0062 for (auto it = output.end(); it != output.begin();)
0063 it = (*--it) ? output.begin() : output.erase(it);
0064 }
0065 }
0066
0067
0068 void HoughTransform::fillIn(int inv2R, int sector, const std::vector<StubGP*>& input, std::vector<StubHT*>& output) {
0069 const DataFormat& gp = dataFormats_->format(Variable::phiT, Process::gp);
0070 auto inv2RrangeCheck = [inv2R](StubGP* stub) {
0071 return (stub && stub->inv2RMin() <= inv2R && stub->inv2RMax() >= inv2R) ? stub : nullptr;
0072 };
0073 const int gpPhiT = gp.toSigned(sector % setup_->gpNumBinsPhiT());
0074 const int zT = sector / setup_->gpNumBinsPhiT() - setup_->gpNumBinsZT() / 2;
0075 const double inv2Rf = inv2R_->floating(inv2R);
0076 auto convert = [this, inv2Rf, gpPhiT, zT, gp](StubGP* stub, int phiTht, double phi, double z) {
0077 const double phiTf = phiT_->floating(phiTht);
0078 const int phiT = phiT_->integer(gp.floating(gpPhiT) + phiTf);
0079 const double htPhi = phi - (inv2Rf * stub->r() + phiTf);
0080 stubs_.emplace_back(*stub, stub->r(), htPhi, z, stub->layer(), phiT, zT);
0081 return &stubs_.back();
0082 };
0083
0084 static constexpr int latency = 1;
0085
0086 std::deque<StubHT*> delay(latency, nullptr);
0087
0088 std::deque<StubHT*> stack;
0089
0090 std::deque<StubGP*> stream;
0091 std::transform(input.begin(), input.end(), std::back_inserter(stream), inv2RrangeCheck);
0092
0093 while (!stream.empty() || !stack.empty() ||
0094 !std::all_of(delay.begin(), delay.end(), [](const StubHT* stub) { return !stub; })) {
0095 StubHT* stubHT = nullptr;
0096 StubGP* stubGP = pop_front(stream);
0097 if (stubGP) {
0098 const double phi = stubGP->phi();
0099 const double z = stubGP->z();
0100 const double phiT = phi - inv2Rf * stubGP->r();
0101 const int major = phiT_->integer(phiT);
0102 if (major >= -setup_->htNumBinsPhiT() / 2 && major < setup_->htNumBinsPhiT() / 2) {
0103
0104 stubHT = convert(stubGP, major, phi, z);
0105 }
0106 const double chi = phi_->digi(phiT - phiT_->floating(major));
0107 if (abs(stubGP->r() * inv2R_->base()) + 2. * abs(chi) >= phiT_->base()) {
0108
0109 const int minor = chi >= 0. ? major + 1 : major - 1;
0110 if (minor >= -setup_->htNumBinsPhiT() / 2 && minor < setup_->htNumBinsPhiT() / 2) {
0111
0112 StubHT* stub = convert(stubGP, minor, phi, z);
0113 delay.push_back(stub);
0114 }
0115 }
0116 }
0117
0118 if (static_cast<int>(delay.size()) == latency)
0119 delay.push_back(nullptr);
0120
0121 StubHT* stub = pop_front(delay);
0122 if (stub) {
0123
0124 if (setup_->enableTruncation() && static_cast<int>(stack.size()) == setup_->htDepthMemory() - 1)
0125 pop_front(stack);
0126
0127 stack.push_back(stub);
0128 }
0129
0130 output.push_back(stubHT ? stubHT : pop_front(stack));
0131 }
0132 }
0133
0134
0135 void HoughTransform::readOut(const std::vector<StubHT*>& input, std::deque<StubHT*>& output) const {
0136 auto toBinPhiT = [this](StubHT* stub) {
0137 const DataFormat& gp = dataFormats_->format(Variable::phiT, Process::gp);
0138 const double phiT = phiT_->floating(stub->phiT());
0139 const double local = phiT - gp.digi(phiT);
0140 return phiT_->integer(local) + setup_->htNumBinsPhiT() / 2;
0141 };
0142 auto toLayerId = [this](StubHT* stub) {
0143 const DataFormat& layer = dataFormats_->format(Variable::layer, Process::ctb);
0144 return stub->layer().val(layer.width());
0145 };
0146
0147 TTBV trkFoundPhiTs(0, setup_->htNumBinsPhiT());
0148
0149 std::vector<TTBV> patternHits(setup_->htNumBinsPhiT(), TTBV(0, setup_->numLayers()));
0150
0151 std::vector<int> phiTs;
0152 phiTs.reserve(setup_->htNumBinsPhiT());
0153 for (StubHT* stub : input) {
0154 const int binPhiT = toBinPhiT(stub);
0155 const int layerId = toLayerId(stub);
0156 TTBV& pattern = patternHits[binPhiT];
0157 pattern.set(layerId);
0158 if (trkFoundPhiTs[binPhiT] || noTrack(pattern, stub->zT()))
0159 continue;
0160
0161 trkFoundPhiTs.set(binPhiT);
0162 phiTs.push_back(binPhiT);
0163 }
0164
0165 for (int phiT : phiTs) {
0166 auto samePhiT = [phiT, toBinPhiT](StubHT* stub) { return toBinPhiT(stub) == phiT; };
0167
0168 std::copy_if(input.rbegin(), input.rend(), std::back_inserter(output), samePhiT);
0169 }
0170 }
0171
0172
0173 bool HoughTransform::noTrack(const TTBV& pattern, int zT) const {
0174
0175 if (pattern.count(0, setup_->kfMaxSeedingLayer()) < 2)
0176 return true;
0177
0178 const int minLayers =
0179 setup_->htMinLayers() - (((zT == -4 || zT == 3) && (!pattern.test(5) && !pattern.test(7))) ? 1 : 0);
0180
0181 const TTBV& maybePattern = layerEncoding_->maybePattern(zT);
0182 int nHits(0);
0183 int nGaps(0);
0184 bool doubleGap = false;
0185 for (int layer = 0; layer < setup_->numLayers(); layer++) {
0186 if (pattern.test(layer)) {
0187 doubleGap = false;
0188 if (++nHits == minLayers)
0189 return false;
0190 } else if (nHits < setup_->kfMinLayers() && !maybePattern.test(layer)) {
0191 if (++nGaps == setup_->kfMaxGaps() || doubleGap)
0192 break;
0193 doubleGap = true;
0194 }
0195 }
0196 return true;
0197 }
0198
0199
0200 template <class T>
0201 T* HoughTransform::pop_front(std::deque<T*>& ts) const {
0202 T* t = nullptr;
0203 if (!ts.empty()) {
0204 t = ts.front();
0205 ts.pop_front();
0206 }
0207 return t;
0208 }
0209
0210 }