TrackerTFP/src/ZHoughTransform.cc

0001 #include "L1Trigger/TrackerTFP/interface/ZHoughTransform.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   ZHoughTransform::ZHoughTransform(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         region_(region),
0026         input_(dataFormats->numChannel(Process::mht)),
0027         stage_(0) {}
0028
0029   // read in and organize input product (fill vector input_)
0030   void ZHoughTransform::consume(const StreamsStub& streams) {
0031     auto valid = [](int sum, const FrameStub& frame) { return sum + (frame.first.isNonnull() ? 1 : 0); };
0032     const int offset = region_ * dataFormats_->numChannel(Process::mht);
0033     int nStubsMHT(0);
0034     for (int channel = 0; channel < dataFormats_->numChannel(Process::mht); channel++) {
0035       const StreamStub& stream = streams[offset + channel];
0036       nStubsMHT += accumulate(stream.begin(), stream.end(), 0, valid);
0037     }
0038     stubsZHT_.reserve(nStubsMHT * (setup_->zhtNumCells() * setup_->zhtNumStages()));
0039     for (int channel = 0; channel < dataFormats_->numChannel(Process::mht); channel++) {
0040       const StreamStub& stream = streams[offset + channel];
0041       vector<StubZHT*>& stubs = input_[channel];
0042       stubs.reserve(stream.size());
0043       // Store input stubs in vector, so rest of ZHT algo can work with pointers to them (saves CPU)
0044       for (const FrameStub& frame : stream) {
0045         StubZHT* stub = nullptr;
0046         if (frame.first.isNonnull()) {
0047           StubMHT stubMHT(frame, dataFormats_);
0048           stubsZHT_.emplace_back(stubMHT);
0049           stub = &stubsZHT_.back();
0050         }
0051         stubs.push_back(stub);
0052       }
0053     }
0054   }
0055
0056   // fill output products
0057   void ZHoughTransform::produce(StreamsStub& accepted, StreamsStub& lost) {
0058     vector<deque<StubZHT*>> streams(dataFormats_->numChannel(Process::mht));
0059     for (int channel = 0; channel < dataFormats_->numChannel(Process::mht); channel++)
0060       streams[channel] = deque<StubZHT*>(input_[channel].begin(), input_[channel].end());
0061     vector<deque<StubZHT*>> stubsCells(dataFormats_->numChannel(Process::mht) * setup_->zhtNumCells());
0062     for (stage_ = 0; stage_ < setup_->zhtNumStages(); stage_++) {
0063       // fill ZHT cells
0064       for (int channel = 0; channel < dataFormats_->numChannel(Process::mht); channel++)
0065         fill(channel, streams[channel], stubsCells);
0066       // perform static load balancing
0067       for (int channel = 0; channel < dataFormats_->numChannel(Process::mht); channel++) {
0068         vector<deque<StubZHT*>> tmp(setup_->zhtNumCells());
0069         // gather streams to mux together: same ZHT cell of 4 adjacent ZHT input streams
0070         for (int k = 0; k < setup_->zhtNumCells(); k++)
0071           //swap(tmp[k], stubsCells[(channel / setup_->zhtNumCells()) * dataFormats_->numChannel(Process::mht) + channel % setup_->zhtNumCells() + k * setup_->zhtNumCells()]);
0072           swap(tmp[k], stubsCells[channel * setup_->zhtNumCells() + k]);
0073         slb(tmp, streams[channel], lost[channel]);
0074       }
0075     }
0076     // fill output product
0077     for (int channel = 0; channel < dataFormats_->numChannel(Process::mht); channel++) {
0078       deque<StubZHT*>& stubs = streams[channel];
0079       StreamStub& stream = accepted[region_ * dataFormats_->numChannel(Process::mht) + channel];
0080       merge(stubs, stream);
0081     }
0082   }
0083
0084   // perform finer pattern recognition per track
0085   void ZHoughTransform::fill(int channel, const deque<StubZHT*>& stubs, vector<deque<StubZHT*>>& streams) {
0086     if (stubs.empty())
0087       return;
0088     const double baseZT =
0089         dataFormats_->format(Variable::zT, Process::zht).base() * pow(2, setup_->zhtNumStages() - stage_);
0090     const double baseCot =
0091         dataFormats_->format(Variable::cot, Process::zht).base() * pow(2, setup_->zhtNumStages() - stage_);
0092     int id;
0093     auto different = [&id](StubZHT* stub) { return !stub || id != stub->trackId(); };
0094     for (auto it = stubs.begin(); it != stubs.end();) {
0095       if (!*it) {
0096         const auto begin = find_if(it, stubs.end(), [](StubZHT* stub) { return stub; });
0097         const int nGaps = distance(it, begin);
0098         for (deque<StubZHT*>& stream : streams)
0099           stream.insert(stream.end(), nGaps, nullptr);
0100         it = begin;
0101         continue;
0102       }
0103       const auto start = it;
0104       const double cotGlobal = (*start)->cotf() + setup_->sectorCot((*start)->sectorEta());
0105       id = (*it)->trackId();
0106       it = find_if(it, stubs.end(), different);
0107       const int size = distance(start, it);
0108       // create finer track candidates stub container
0109       vector<vector<StubZHT*>> mhtCells(setup_->zhtNumCells());
0110       for (vector<StubZHT*>& mhtCell : mhtCells)
0111         mhtCell.reserve(size);
0112       // fill finer track candidates stub container
0113       for (auto stub = start; stub != it; stub++) {
0114         const double r = (*stub)->r() + setup_->chosenRofPhi() - setup_->chosenRofZ();
0115         const double chi = (*stub)->chi();
0116         const double dChi = setup_->dZ((*stub)->ttStubRef(), cotGlobal);
0117         // identify finer track candidates for this stub
0118         // 0 and 1 belong to the ZHT cells with smaller cot; 0 and 2 belong to those with smaller zT
0119         vector<int> cells;
0120         cells.reserve(setup_->zhtNumCells());
0121         const bool compA = 2. * abs(chi) < baseZT + dChi;
0122         const bool compB = 2. * abs(chi) < abs(r) * baseCot + dChi;
0123         const bool compC = 2. * abs(chi) < dChi;
0124         if (chi >= 0. && r >= 0.) {
0125           cells.push_back(1);
0126           if (compA)
0127             cells.push_back(3);
0128           if (compB)
0129             cells.push_back(0);
0130           if (compC)
0131             cells.push_back(2);
0132         }
0133         if (chi >= 0. && r < 0.) {
0134           cells.push_back(3);
0135           if (compA)
0136             cells.push_back(1);
0137           if (compB)
0138             cells.push_back(2);
0139           if (compC)
0140             cells.push_back(0);
0141         }
0142         if (chi < 0. && r >= 0.) {
0143           cells.push_back(2);
0144           if (compA)
0145             cells.push_back(0);
0146           if (compB)
0147             cells.push_back(3);
0148           if (compC)
0149             cells.push_back(1);
0150         }
0151         if (chi < 0. && r < 0.) {
0152           cells.push_back(0);
0153           if (compA)
0154             cells.push_back(2);
0155           if (compB)
0156             cells.push_back(1);
0157           if (compC)
0158             cells.push_back(3);
0159         }
0160         for (int cell : cells) {
0161           const double cot = (cell / setup_->zhtNumBinsZT() - .5) * baseCot / 2.;
0162           const double zT = (cell % setup_->zhtNumBinsZT() - .5) * baseZT / 2.;
0163           stubsZHT_.emplace_back(**stub, zT, cot, cell);
0164           mhtCells[cell].push_back(&stubsZHT_.back());
0165         }
0166       }
0167       // perform pattern recognition
0168       for (int sel = 0; sel < setup_->zhtNumCells(); sel++) {
0169         deque<StubZHT*>& stream = streams[channel * setup_->zhtNumCells() + sel];
0170         vector<StubZHT*>& mhtCell = mhtCells[sel];
0171         set<int> layers;
0172         auto toLayer = [](StubZHT* stub) { return stub->layer(); };
0173         transform(mhtCell.begin(), mhtCell.end(), inserter(layers, layers.begin()), toLayer);
0174         if ((int)layers.size() < setup_->mhtMinLayers())
0175           mhtCell.clear();
0176         for (StubZHT* stub : mhtCell)
0177           stream.push_back(stub);
0178         stream.insert(stream.end(), size - (int)mhtCell.size(), nullptr);
0179       }
0180     }
0181     for (int sel = 0; sel < setup_->zhtNumCells(); sel++) {
0182       deque<StubZHT*>& stream = streams[channel * setup_->zhtNumCells() + sel];
0183       // remove all gaps between end and last stub
0184       for (auto it = stream.end(); it != stream.begin();)
0185         it = (*--it) ? stream.begin() : stream.erase(it);
0186       // read out fine track cannot start before rough track has read in completely, add gaps to take this into account
0187       int pos(0);
0188       for (auto it = stream.begin(); it != stream.end();) {
0189         if (!(*it)) {
0190           it = stream.erase(it);
0191           continue;
0192         }
0193         id = (*it)->trackId();
0194         const int s = distance(it, find_if(it, stream.end(), different));
0195         const int d = distance(stream.begin(), it);
0196         pos += s;
0197         if (d < pos) {
0198           const int diff = pos - d;
0199           it = stream.insert(it, diff, nullptr);
0200           it = next(it, diff);
0201         } else
0202           it = stream.erase(remove(next(stream.begin(), pos), it, nullptr), it);
0203         it = next(it, s);
0204       }
0205       // adjust stream start so that first output stub is in first place in case of quickest track
0206       if (!stream.empty())
0207         stream.erase(stream.begin(), next(stream.begin(), setup_->mhtMinLayers()));
0208     }
0209   }
0210
0211   // Static load balancing of inputs: mux 4 streams to 1 stream
0212   void ZHoughTransform::slb(vector<deque<StubZHT*>>& inputs, deque<StubZHT*>& accepted, StreamStub& lost) const {
0213     accepted.clear();
0214     if (all_of(inputs.begin(), inputs.end(), [](const deque<StubZHT*>& stubs) { return stubs.empty(); }))
0215       return;
0216     // input fifos
0217     vector<deque<StubZHT*>> stacks(setup_->zhtNumCells());
0218     // helper for handshake
0219     TTBV empty(-1, setup_->zhtNumCells(), true);
0220     TTBV enable(0, setup_->zhtNumCells());
0221     // clock accurate firmware emulation, each while trip describes one clock tick, one stub in and one stub out per tick
0222     while (!all_of(inputs.begin(), inputs.end(), [](const deque<StubZHT*>& d) { return d.empty(); }) or
0223            !all_of(stacks.begin(), stacks.end(), [](const deque<StubZHT*>& d) { return d.empty(); })) {
0224       // store stub in fifo
0225       for (int channel = 0; channel < setup_->zhtNumCells(); channel++) {
0226         StubZHT* stub = pop_front(inputs[channel]);
0227         if (stub)
0228           stacks[channel].push_back(stub);
0229       }
0230       // identify empty fifos
0231       for (int channel = 0; channel < setup_->zhtNumCells(); channel++)
0232         empty[channel] = stacks[channel].empty();
0233       // chose new fifo to read from if current fifo got empty
0234       const int iEnableOld = enable.plEncode();
0235       if (enable.none() || empty[iEnableOld]) {
0236         enable.reset();
0237         const int iNotEmpty = empty.plEncode(false);
0238         if (iNotEmpty < setup_->zhtNumCells())
0239           enable.set(iNotEmpty);
0240       }
0241       // read from chosen fifo
0242       const int iEnable = enable.plEncode();
0243       if (enable.any())
0244         accepted.push_back(pop_front(stacks[iEnable]));
0245       else
0246         // gap if no fifo has been chosen
0247         accepted.push_back(nullptr);
0248     }
0249     // perform truncation if desired
0250     if (enableTruncation_ && (int)accepted.size() > setup_->numFrames()) {
0251       const auto limit = next(accepted.begin(), setup_->numFrames());
0252       auto valid = [](int sum, StubZHT* stub) { return sum + (stub ? 1 : 0); };
0253       const int nLost = accumulate(limit, accepted.end(), 0, valid);
0254       lost.reserve(nLost);
0255       for (auto it = limit; it != accepted.end(); it++)
0256         if (*it)
0257           lost.emplace_back((*it)->frame());
0258       accepted.erase(limit, accepted.end());
0259     }
0260     // cosmetics -- remove gaps at the end of stream
0261     for (auto it = accepted.end(); it != accepted.begin();)
0262       it = (*--it) == nullptr ? accepted.erase(it) : accepted.begin();
0263   }
0264
0265   //
0266   void ZHoughTransform::merge(deque<StubZHT*>& stubs, StreamStub& stream) const {
0267     stubs.erase(remove(stubs.begin(), stubs.end(), nullptr), stubs.end());
0268     /*stream.reserve(stubs.size());
0269     transform(stubs.begin(), stubs.end(), back_inserter(stream), [](StubZHT* stub){ return stub->frame(); });
0270     return;*/
0271     map<int, set<pair<int, int>>> candidates;
0272     const int weight = setup_->zhtNumCells() * pow(2, setup_->zhtNumStages());
0273     for (const StubZHT* stub : stubs)
0274       candidates[stub->trackId() / weight].emplace(stub->cot(), stub->zT());
0275     vector<deque<FrameStub>> tracks(candidates.size());
0276     for (auto it = stubs.begin(); it != stubs.end();) {
0277       const auto start = it;
0278       const int id = (*it)->trackId();
0279       const int candId = id / weight;
0280       const auto m = candidates.find(candId);
0281       pair<int, int> cotp(9e9, -9e9);
0282       pair<int, int> zTp(9e9, -9e9);
0283       for (const pair<int, int>& para : m->second) {
0284         cotp = {min(cotp.first, para.first), max(cotp.second, para.first)};
0285         zTp = {min(zTp.first, para.second), max(zTp.second, para.second)};
0286       }
0287       const int cot = (cotp.first + cotp.second) / 2;
0288       const int zT = (cotp.first + cotp.second) / 2;
0289       const int pos = distance(candidates.begin(), m);
0290       deque<FrameStub>& track = tracks[pos];
0291       auto different = [id](const StubZHT* stub) { return id != stub->trackId(); };
0292       it = find_if(it, stubs.end(), different);
0293       for (auto s = start; s != it; s++) {
0294         if (find_if(track.begin(), track.end(), [s](const FrameStub& stub) {
0295               return (*s)->ttStubRef() == stub.first;
0296             }) != track.end())
0297           continue;
0298         const StubZHT stub(**s, cot, zT);
0299         track.push_back(stub.frame());
0300       }
0301     }
0302     const int size = accumulate(tracks.begin(), tracks.end(), 0, [](int sum, const deque<FrameStub>& stubs) {
0303       return sum + (int)stubs.size();
0304     });
0305     stream.reserve(size);
0306     for (deque<FrameStub>& track : tracks)
0307       for (const FrameStub& stub : track)
0308         stream.push_back(stub);
0309   }
0310
0311   // remove and return first element of deque, returns nullptr if empty
0312   template <class T>
0313   T* ZHoughTransform::pop_front(deque<T*>& ts) const {
0314     T* t = nullptr;
0315     if (!ts.empty()) {
0316       t = ts.front();
0317       ts.pop_front();
0318     }
0319     return t;
0320   }
0321
0322 }  // namespace trackerTFP