TrackerTFP/src/HoughTransform.cc

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   // fill output products
0032   void HoughTransform::produce(const std::vector<std::vector<StubGP*>>& streamsIn,
0033                                std::vector<std::deque<StubHT*>>& streamsOut) {
0034     // count and reserve ht stubs
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         // associate stubs with inv2R and phiT bins
0049         fillIn(inv2R, channelIn, input, stubs);
0050         // apply truncation
0051         if (setup_->enableTruncation() && static_cast<int>(stubs.size()) > setup_->numFramesHigh())
0052           stubs.resize(setup_->numFramesHigh());
0053         // ht collects all stubs before readout starts -> remove all gaps
0054         stubs.erase(std::remove(stubs.begin(), stubs.end(), nullptr), stubs.end());
0055         // identify tracks
0056         readOut(stubs, output);
0057       }
0058       // apply truncation
0059       if (setup_->enableTruncation() && static_cast<int>(output.size()) > setup_->numFramesHigh())
0060         output.resize(setup_->numFramesHigh());
0061       // remove trailing gaps
0062       for (auto it = output.end(); it != output.begin();)
0063         it = (*--it) ? output.begin() : output.erase(it);
0064     }
0065   }
0066
0067   // associate stubs with phiT bins in this inv2R column
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     // Latency of ht fifo firmware
0084     static constexpr int latency = 1;
0085     // static delay container
0086     std::deque<StubHT*> delay(latency, nullptr);
0087     // fifo, used to store stubs which belongs to a second possible track
0088     std::deque<StubHT*> stack;
0089     // stream of incroming stubs
0090     std::deque<StubGP*> stream;
0091     std::transform(input.begin(), input.end(), std::back_inserter(stream), inv2RrangeCheck);
0092     // clock accurate firmware emulation, each while trip describes one clock tick, one stub in and one stub out per tick
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           // major candidate has pt > threshold (3 GeV)
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           // stub belongs to two candidates
0109           const int minor = chi >= 0. ? major + 1 : major - 1;
0110           if (minor >= -setup_->htNumBinsPhiT() / 2 && minor < setup_->htNumBinsPhiT() / 2) {
0111             // second (minor) candidate has pt > threshold (3 GeV)
0112             StubHT* stub = convert(stubGP, minor, phi, z);
0113             delay.push_back(stub);
0114           }
0115         }
0116       }
0117       // add nullptr to delay pipe if stub didn't fill any cell
0118       if (static_cast<int>(delay.size()) == latency)
0119         delay.push_back(nullptr);
0120       // take fifo latency into account (read before write)
0121       StubHT* stub = pop_front(delay);
0122       if (stub) {
0123         // buffer overflow
0124         if (setup_->enableTruncation() && static_cast<int>(stack.size()) == setup_->htDepthMemory() - 1)
0125           pop_front(stack);
0126         // store minor stub in fifo
0127         stack.push_back(stub);
0128       }
0129       // take a minor stub if no major stub available
0130       output.push_back(stubHT ? stubHT : pop_front(stack));
0131     }
0132   }
0133
0134   // identify tracks
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     // used to recognise in which order tracks are found
0147     TTBV trkFoundPhiTs(0, setup_->htNumBinsPhiT());
0148     // hitPattern for all possible tracks, used to find tracks
0149     std::vector<TTBV> patternHits(setup_->htNumBinsPhiT(), TTBV(0, setup_->numLayers()));
0150     // found phiTs, ordered in time
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       // first time track found
0161       trkFoundPhiTs.set(binPhiT);
0162       phiTs.push_back(binPhiT);
0163     }
0164     // read out found tracks ordered as found
0165     for (int phiT : phiTs) {
0166       auto samePhiT = [phiT, toBinPhiT](StubHT* stub) { return toBinPhiT(stub) == phiT; };
0167       // read out stubs in reverse order to emulate f/w (backtracking linked list)
0168       std::copy_if(input.rbegin(), input.rend(), std::back_inserter(output), samePhiT);
0169     }
0170   }
0171
0172   // track identification
0173   bool HoughTransform::noTrack(const TTBV& pattern, int zT) const {
0174     // not enough seeding layer
0175     if (pattern.count(0, setup_->kfMaxSeedingLayer()) < 2)
0176       return true;
0177     // check min layers req
0178     const int minLayers =
0179         setup_->htMinLayers() - (((zT == -4 || zT == 3) && (!pattern.test(5) && !pattern.test(7))) ? 1 : 0);
0180     // prepare pattern analysis
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   // remove and return first element of deque, returns nullptr if empty
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 }  // namespace trackerTFP