TrackFindingTracklet/src/DuplicateRemoval.cc

0001 #include "L1Trigger/TrackFindingTracklet/interface/DuplicateRemoval.h"
0002
0003 #include <vector>
0004 #include <numeric>
0005 #include <algorithm>
0006
0007 namespace trklet {
0008
0009   DuplicateRemoval::DuplicateRemoval(const tt::Setup* setup,
0010                                      const trackerTFP::LayerEncoding* layerEncoding,
0011                                      const DataFormats* dataFormats,
0012                                      const ChannelAssignment* channelAssignment,
0013                                      int region)
0014       : setup_(setup),
0015         layerEncoding_(layerEncoding),
0016         dataFormats_(dataFormats),
0017         channelAssignment_(channelAssignment),
0018         region_(region) {
0019     const DataFormat& r = dataFormats_->format(Variable::r, Process::dr);
0020     const int width = setup_->widthAddrBRAM18() - 1;
0021     const double base = r.base() * pow(2., r.width() - width);
0022     const double range = r.range();
0023     r_ = DataFormat(true, width, base, range);
0024     tmNumLayers_ = channelAssignment_->tmNumLayers();
0025     phi_ = dataFormats_->format(Variable::phi, Process::dr);
0026   }
0027
0028   // read in and organize input tracks and stubs
0029   void DuplicateRemoval::consume(const tt::StreamsTrack& streamsTrack, const tt::StreamsStub& streamsStub) {
0030     auto nonNullTrack = [](int sum, const tt::FrameTrack& frame) { return sum + (frame.first.isNonnull() ? 1 : 0); };
0031     auto nonNullStub = [](int sum, const tt::FrameStub& frame) { return sum + (frame.first.isNonnull() ? 1 : 0); };
0032     // count tracks and stubs and reserve corresponding vectors
0033     int sizeStubs(0);
0034     const int offset = region_ * tmNumLayers_;
0035     const tt::StreamTrack& streamTrack = streamsTrack[region_];
0036     input_.reserve(streamTrack.size());
0037     const int sizeTracks = std::accumulate(streamTrack.begin(), streamTrack.end(), 0, nonNullTrack);
0038     for (int layer = 0; layer < tmNumLayers_; layer++) {
0039       const tt::StreamStub& streamStub = streamsStub[offset + layer];
0040       sizeStubs += std::accumulate(streamStub.begin(), streamStub.end(), 0, nonNullStub);
0041     }
0042     tracks_.reserve(sizeTracks);
0043     stubs_.reserve(sizeStubs);
0044     // transform input data into handy structs
0045     for (int frame = 0; frame < static_cast<int>(streamTrack.size()); frame++) {
0046       const tt::FrameTrack& frameTrack = streamTrack[frame];
0047       if (frameTrack.first.isNull()) {
0048         input_.push_back(nullptr);
0049         continue;
0050       }
0051       // lookup layerEncoding
0052       const TrackTM track(frameTrack, dataFormats_);
0053       const double inv2R = abs(track.inv2R());
0054       const double zT = abs(track.zT());
0055       const double cot = zT / setup_->chosenRofZ();
0056       const std::vector<int>& layerEncoding = layerEncoding_->layerEncoding(zT);
0057       std::vector<Stub*> stubs(tmNumLayers_, nullptr);
0058       TTBV hitPattern(0, setup_->numLayers());
0059       for (int layer = 0; layer < tmNumLayers_; layer++) {
0060         const tt::FrameStub& frameStub = streamsStub[offset + layer][frame];
0061         const TTStubRef& ttStubRef = frameStub.first;
0062         if (ttStubRef.isNull())
0063           continue;
0064         // encode layerId
0065         const int decodedLayerId =
0066             layer + setup_->offsetLayerId() +
0067             (layer < setup_->numBarrelLayer() ? 0 : setup_->offsetLayerDisks() - setup_->numBarrelLayer());
0068         const auto it = std::find(layerEncoding.begin(), layerEncoding.end(), decodedLayerId);
0069         const int encodedLayerId =
0070             std::min(static_cast<int>(std::distance(layerEncoding.begin(), it)), setup_->numLayers() - 1);
0071         // kill stub on already occupied layer
0072         if (hitPattern.test(encodedLayerId))
0073           continue;
0074         hitPattern.set(encodedLayerId);
0075         const StubTM stubTM(frameStub, dataFormats_);
0076         const int stubId = stubTM.stubId() / 2;
0077         const bool psTilt = stubTM.stubId() % 2 == 1;
0078         // calculate stub uncertainties
0079         static constexpr int numBarrelPSLayer = 3;
0080         const bool barrel = layer < setup_->numBarrelLayer();
0081         const bool ps = barrel ? layer < numBarrelPSLayer : psTilt;
0082         const bool tilt = barrel && psTilt;
0083         const double length = .5 * (ps ? setup_->pitchColPS() : setup_->pitchCol2S());
0084         const double pitch = .5 * (ps ? setup_->pitchRowPS() : setup_->pitchRow2S());
0085         const double pitchOverR = phi_.digi(pitch / (r_.digi(stubTM.r()) + setup_->chosenRofPhi()));
0086         double lengthZ = length;
0087         double lengthR = 0.;
0088         if (!barrel) {
0089           lengthZ = length * cot;
0090           lengthR = length;
0091         } else if (tilt) {
0092           lengthZ = length * (setup_->tiltApproxSlope() * cot + setup_->tiltApproxIntercept());
0093           lengthR = .5 * setup_->tiltUncertaintyR();
0094         }
0095         const double dR = lengthR + .5 * setup_->scattering();
0096         const double dZ = lengthZ;
0097         const double dPhi = phi_.digi(dR * inv2R) + pitchOverR;
0098         const StubDR stubDR(stubTM, stubTM.r(), stubTM.phi(), stubTM.z(), dPhi, dZ);
0099         stubs_.emplace_back(stubDR.frame(), stubId, encodedLayerId);
0100         stubs[layer] = &stubs_.back();
0101       }
0102       // kill tracks with not enough layers
0103       if (hitPattern.count() < setup_->kfMinLayers()) {
0104         input_.push_back(nullptr);
0105         continue;
0106       }
0107       tracks_.emplace_back(frameTrack, stubs);
0108       input_.push_back(&tracks_.back());
0109     }
0110     // remove all gaps between end and last track
0111     for (auto it = input_.end(); it != input_.begin();)
0112       it = (*--it) ? input_.begin() : input_.erase(it);
0113   }
0114
0115   // fill output products
0116   void DuplicateRemoval::produce(tt::StreamsTrack& streamsTrack, tt::StreamsStub& streamsStub) {
0117     const int offset = region_ * setup_->numLayers();
0118     // remove duplicated tracks, no merge of stubs, one stub per layer expected
0119     std::vector<Track*> cms(channelAssignment_->numComparisonModules(), nullptr);
0120     for (Track*& track : input_) {
0121       if (!track)
0122         // gaps propagate through chain and appear in output stream
0123         continue;
0124       for (Track*& trackCM : cms) {
0125         if (!trackCM) {
0126           // tracks used in CMs propagate through chain and do appear in output stream unaltered
0127           trackCM = track;
0128           break;
0129         }
0130         if (equalEnough(track, trackCM)) {
0131           // tracks compared in CMs propagate through chain and appear in output stream as gap if identified as duplicate or unaltered elsewise
0132           track = nullptr;
0133           break;
0134         }
0135       }
0136     }
0137     // remove all gaps between end and last track
0138     for (auto it = input_.end(); it != input_.begin();)
0139       it = (*--it) ? input_.begin() : input_.erase(it);
0140     // store output
0141     tt::StreamTrack& streamTrack = streamsTrack[region_];
0142     streamTrack.reserve(input_.size());
0143     for (int layer = 0; layer < setup_->numLayers(); layer++)
0144       streamsStub[offset + layer].reserve(input_.size());
0145     for (Track* track : input_) {
0146       if (!track) {
0147         streamTrack.emplace_back(tt::FrameTrack());
0148         for (int layer = 0; layer < setup_->numLayers(); layer++)
0149           streamsStub[offset + layer].emplace_back(tt::FrameStub());
0150         continue;
0151       }
0152       streamTrack.push_back(track->frame_);
0153       TTBV hitPattern(0, setup_->numLayers());
0154       for (Stub* stub : track->stubs_) {
0155         if (!stub)
0156           continue;
0157         hitPattern.set(stub->layer_);
0158         streamsStub[offset + stub->layer_].emplace_back(stub->frame_);
0159       }
0160       for (int layer : hitPattern.ids(false))
0161         streamsStub[offset + layer].emplace_back(tt::FrameStub());
0162     }
0163   }
0164
0165   // compares two tracks, returns true if those are considered duplicates
0166   bool DuplicateRemoval::equalEnough(Track* t0, Track* t1) const {
0167     int same(0);
0168     for (int layer = 0; layer < channelAssignment_->tmNumLayers(); layer++) {
0169       Stub* s0 = t0->stubs_[layer];
0170       Stub* s1 = t1->stubs_[layer];
0171       if (s0 && s1 && s0->stubId_ == s1->stubId_)
0172         same++;
0173     }
0174     return same >= channelAssignment_->minIdenticalStubs();
0175   }
0176
0177 }  // namespace trklet