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 <set>
0009 #include <utility>
0010 #include <cmath>
0011
0012 using namespace std;
0013 using namespace edm;
0014 using namespace tt;
0015
0016 namespace trackerTFP {
0017
0018   HoughTransform::HoughTransform(const ParameterSet& iConfig,
0019                                  const Setup* setup,
0020                                  const DataFormats* dataFormats,
0021                                  int region)
0022       : enableTruncation_(iConfig.getParameter<bool>("EnableTruncation")),
0023         setup_(setup),
0024         dataFormats_(dataFormats),
0025         inv2R_(dataFormats_->format(Variable::inv2R, Process::ht)),
0026         phiT_(dataFormats_->format(Variable::phiT, Process::ht)),
0027         region_(region),
0028         input_(dataFormats_->numChannel(Process::ht), vector<deque<StubGP*>>(dataFormats_->numChannel(Process::gp))) {}
0029
0030   // read in and organize input product (fill vector input_)
0031   void HoughTransform::consume(const StreamsStub& streams) {
0032     const int offset = region_ * dataFormats_->numChannel(Process::gp);
0033     auto validFrame = [](int sum, const FrameStub& frame) { return sum + (frame.first.isNonnull() ? 1 : 0); };
0034     int nStubsGP(0);
0035     for (int sector = 0; sector < dataFormats_->numChannel(Process::gp); sector++) {
0036       const StreamStub& stream = streams[offset + sector];
0037       nStubsGP += accumulate(stream.begin(), stream.end(), 0, validFrame);
0038     }
0039     stubsGP_.reserve(nStubsGP);
0040     for (int sector = 0; sector < dataFormats_->numChannel(Process::gp); sector++) {
0041       const int sectorPhi = sector % setup_->numSectorsPhi();
0042       const int sectorEta = sector / setup_->numSectorsPhi();
0043       for (const FrameStub& frame : streams[offset + sector]) {
0044         // Store input stubs in vector, so rest of HT algo can work with pointers to them (saves CPU)
0045         StubGP* stub = nullptr;
0046         if (frame.first.isNonnull()) {
0047           stubsGP_.emplace_back(frame, dataFormats_, sectorPhi, sectorEta);
0048           stub = &stubsGP_.back();
0049         }
0050         for (int binInv2R = 0; binInv2R < dataFormats_->numChannel(Process::ht); binInv2R++)
0051           input_[binInv2R][sector].push_back(stub && stub->inInv2RBin(binInv2R) ? stub : nullptr);
0052       }
0053     }
0054     // remove all gaps between end and last stub
0055     for (vector<deque<StubGP*>>& input : input_)
0056       for (deque<StubGP*>& stubs : input)
0057         for (auto it = stubs.end(); it != stubs.begin();)
0058           it = (*--it) ? stubs.begin() : stubs.erase(it);
0059     auto validStub = [](int sum, StubGP* stub) { return sum + (stub ? 1 : 0); };
0060     int nStubsHT(0);
0061     for (const vector<deque<StubGP*>>& binInv2R : input_)
0062       for (const deque<StubGP*>& sector : binInv2R)
0063         nStubsHT += accumulate(sector.begin(), sector.end(), 0, validStub);
0064     stubsHT_.reserve(nStubsHT);
0065   }
0066
0067   // fill output products
0068   void HoughTransform::produce(StreamsStub& accepted, StreamsStub& lost) {
0069     for (int binInv2R = 0; binInv2R < dataFormats_->numChannel(Process::ht); binInv2R++) {
0070       const int inv2R = inv2R_.toSigned(binInv2R);
0071       deque<StubHT*> acceptedAll;
0072       deque<StubHT*> lostAll;
0073       for (deque<StubGP*>& inputSector : input_[binInv2R]) {
0074         const int size = inputSector.size();
0075         vector<StubHT*> acceptedSector;
0076         vector<StubHT*> lostSector;
0077         acceptedSector.reserve(size);
0078         lostSector.reserve(size);
0079         // associate stubs with inv2R and phiT bins
0080         fillIn(inv2R, inputSector, acceptedSector, lostSector);
0081         // Process::ht collects all stubs before readout starts -> remove all gaps
0082         acceptedSector.erase(remove(acceptedSector.begin(), acceptedSector.end(), nullptr), acceptedSector.end());
0083         // identify tracks
0084         readOut(acceptedSector, lostSector, acceptedAll, lostAll);
0085       }
0086       // truncate accepted stream
0087       const auto limit = enableTruncation_
0088                              ? next(acceptedAll.begin(), min(setup_->numFrames(), (int)acceptedAll.size()))
0089                              : acceptedAll.end();
0090       copy_if(limit, acceptedAll.end(), back_inserter(lostAll), [](StubHT* stub) { return stub; });
0091       acceptedAll.erase(limit, acceptedAll.end());
0092       // store found tracks
0093       auto put = [](const deque<StubHT*>& stubs, StreamStub& stream) {
0094         stream.reserve(stubs.size());
0095         for (StubHT* stub : stubs)
0096           stream.emplace_back(stub ? stub->frame() : FrameStub());
0097       };
0098       const int offset = region_ * dataFormats_->numChannel(Process::ht);
0099       put(acceptedAll, accepted[offset + binInv2R]);
0100       // store lost tracks
0101       put(lostAll, lost[offset + binInv2R]);
0102     }
0103   }
0104
0105   // associate stubs with phiT bins in this inv2R column
0106   void HoughTransform::fillIn(int inv2R,
0107                               deque<StubGP*>& inputSector,
0108                               vector<StubHT*>& acceptedSector,
0109                               vector<StubHT*>& lostSector) {
0110     // fifo, used to store stubs which belongs to a second possible track
0111     deque<StubHT*> stack;
0112     // clock accurate firmware emulation, each while trip describes one clock tick, one stub in and one stub out per tick
0113     while (!inputSector.empty() || !stack.empty()) {
0114       StubHT* stubHT = nullptr;
0115       StubGP* stubGP = pop_front(inputSector);
0116       if (stubGP) {
0117         const double phiT = stubGP->phi() - inv2R_.floating(inv2R) * stubGP->r();
0118         const int major = phiT_.integer(phiT);
0119         if (phiT_.inRange(major)) {
0120           // major candidate has pt > threshold (3 GeV)
0121           // stubHT records which HT bin this stub is added to
0122           stubsHT_.emplace_back(*stubGP, major, inv2R);
0123           stubHT = &stubsHT_.back();
0124         }
0125         const double chi = phiT - phiT_.floating(major);
0126         if (abs(stubGP->r() * inv2R_.base()) + 2. * abs(chi) >= phiT_.base()) {
0127           // stub belongs to two candidates
0128           const int minor = chi >= 0. ? major + 1 : major - 1;
0129           if (phiT_.inRange(minor)) {
0130             // second (minor) candidate has pt > threshold (3 GeV)
0131             stubsHT_.emplace_back(*stubGP, minor, inv2R);
0132             if (enableTruncation_ && (int)stack.size() == setup_->htDepthMemory() - 1)
0133               // buffer overflow
0134               lostSector.push_back(pop_front(stack));
0135             // store minor stub in fifo
0136             stack.push_back(&stubsHT_.back());
0137           }
0138         }
0139       }
0140       // take a minor stub if no major stub available
0141       acceptedSector.push_back(stubHT ? stubHT : pop_front(stack));
0142     }
0143     // truncate to many input stubs
0144     const auto limit = enableTruncation_
0145                            ? next(acceptedSector.begin(), min(setup_->numFrames(), (int)acceptedSector.size()))
0146                            : acceptedSector.end();
0147     copy_if(limit, acceptedSector.end(), back_inserter(lostSector), [](StubHT* stub) { return stub; });
0148     acceptedSector.erase(limit, acceptedSector.end());
0149   }
0150
0151   // identify tracks
0152   void HoughTransform::readOut(const vector<StubHT*>& acceptedSector,
0153                                const vector<StubHT*>& lostSector,
0154                                deque<StubHT*>& acceptedAll,
0155                                deque<StubHT*>& lostAll) const {
0156     // used to recognise in which order tracks are found
0157     TTBV trkFoundPhiTs(0, setup_->htNumBinsPhiT());
0158     // hitPattern for all possible tracks, used to find tracks
0159     vector<TTBV> patternHits(setup_->htNumBinsPhiT(), TTBV(0, setup_->numLayers()));
0160     // found unsigned phiTs, ordered in time
0161     vector<int> binsPhiT;
0162     // stub container for all possible tracks
0163     vector<vector<StubHT*>> tracks(setup_->htNumBinsPhiT());
0164     for (int binPhiT = 0; binPhiT < setup_->htNumBinsPhiT(); binPhiT++) {
0165       const int phiT = phiT_.toSigned(binPhiT);
0166       auto samePhiT = [phiT](int sum, StubHT* stub) { return sum + (stub->phiT() == phiT); };
0167       const int numAccepted = accumulate(acceptedSector.begin(), acceptedSector.end(), 0, samePhiT);
0168       const int numLost = accumulate(lostSector.begin(), lostSector.end(), 0, samePhiT);
0169       tracks[binPhiT].reserve(numAccepted + numLost);
0170     }
0171     for (StubHT* stub : acceptedSector) {
0172       const int binPhiT = phiT_.toUnsigned(stub->phiT());
0173       TTBV& pattern = patternHits[binPhiT];
0174       pattern.set(stub->layer());
0175       tracks[binPhiT].push_back(stub);
0176       if (pattern.count() >= setup_->htMinLayers() && !trkFoundPhiTs[binPhiT]) {
0177         // first time track found
0178         trkFoundPhiTs.set(binPhiT);
0179         binsPhiT.push_back(binPhiT);
0180       }
0181     }
0182     // read out found tracks ordered as found
0183     for (int binPhiT : binsPhiT) {
0184       const vector<StubHT*>& track = tracks[binPhiT];
0185       acceptedAll.insert(acceptedAll.end(), track.begin(), track.end());
0186     }
0187     // look for lost tracks
0188     for (StubHT* stub : lostSector) {
0189       const int binPhiT = phiT_.toUnsigned(stub->phiT());
0190       if (!trkFoundPhiTs[binPhiT])
0191         tracks[binPhiT].push_back(stub);
0192     }
0193     for (int binPhiT : trkFoundPhiTs.ids(false)) {
0194       const vector<StubHT*>& track = tracks[binPhiT];
0195       set<int> layers;
0196       auto toLayer = [](StubHT* stub) { return stub->layer(); };
0197       transform(track.begin(), track.end(), inserter(layers, layers.begin()), toLayer);
0198       if ((int)layers.size() >= setup_->htMinLayers())
0199         lostAll.insert(lostAll.end(), track.begin(), track.end());
0200     }
0201   }
0202
0203   // remove and return first element of deque, returns nullptr if empty
0204   template <class T>
0205   T* HoughTransform::pop_front(deque<T*>& ts) const {
0206     T* t = nullptr;
0207     if (!ts.empty()) {
0208       t = ts.front();
0209       ts.pop_front();
0210     }
0211     return t;
0212   }
0213
0214 }  // namespace trackerTFP