Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:30:44

0001 #include "L1Trigger/TrackFindingTMTT/interface/HTrphi.h"
0002 #include "L1Trigger/TrackFindingTMTT/interface/InputData.h"
0003 #include "L1Trigger/TrackFindingTMTT/interface/TP.h"
0004 #include "L1Trigger/TrackFindingTMTT/interface/L1fittedTrack.h"
0005 #include "L1Trigger/TrackFindingTMTT/interface/Settings.h"
0006 #include "L1Trigger/TrackFindingTMTT/interface/PrintL1trk.h"
0007 #include "FWCore/ServiceRegistry/interface/Service.h"
0008 
0009 #include "DataFormats/Math/interface/deltaPhi.h"
0010 #include "FWCore/Utilities/interface/Exception.h"
0011 
0012 #include <vector>
0013 #include <limits>
0014 #include <atomic>
0015 #include <sstream>
0016 #include <mutex>
0017 
0018 using namespace std;
0019 
0020 namespace tmtt {
0021 
0022   //=== The r-phi Hough Transform array for a single (eta,phi) sector.
0023   //===
0024   //=== Its axes are (q/Pt, phiTrk), where phiTrk is the phi at which the track crosses a
0025   //=== user-configurable radius from the beam-line.
0026 
0027   //=== Initialise
0028 
0029   HTrphi::HTrphi(const Settings* settings,
0030                  unsigned int iPhiSec,
0031                  unsigned int iEtaReg,
0032                  float etaMinSector,
0033                  float etaMaxSector,
0034                  float phiCentreSector,
0035                  HTrphi::ErrorMonitor* errMon)
0036       : HTbase(settings, iPhiSec, iEtaReg, settings->houghNbinsPt(), settings->houghNbinsPhi()),
0037         invPtToDphi_((settings->invPtToDphi())),
0038         shape_(static_cast<HTshape>(settings->shape())),  // shape of HT cells
0039 
0040         //--- Specification of HT q/Pt axis.
0041         maxAbsQoverPtAxis_(1. / settings->houghMinPt()),  // Max. |q/Pt| covered by  HT array.
0042         nBinsQoverPtAxis_(settings->houghNbinsPt()),      // No. of bins in HT array in q/Pt.
0043         binSizeQoverPtAxis_(2 * maxAbsQoverPtAxis_ / nBinsQoverPtAxis_),
0044 
0045         //--- Specification of HT phiTrk axis
0046         // (phiTrk corresponds to phi where track crosses radius = chosenRofPhi_).
0047         chosenRofPhi_(settings->chosenRofPhi()),
0048         phiCentreSector_(phiCentreSector),                           // Centre of phiTrk sector.
0049         maxAbsPhiTrkAxis_(M_PI / float(settings->numPhiSectors())),  // Half-width of phiTrk axis in HT array.
0050         nBinsPhiTrkAxis_(settings->houghNbinsPhi()),                 // No. of bins in HT array phiTrk
0051         binSizePhiTrkAxis_(2 * maxAbsPhiTrkAxis_ / nBinsPhiTrkAxis_),
0052         errMon_(errMon) {
0053     // Deal with unusually shaped HT cells.
0054     if (shape_ != HTshape::square)
0055       binSizeQoverPtAxis_ = 2. * maxAbsQoverPtAxis_ / (nBinsQoverPtAxis_ - 1.);
0056     if (shape_ == HTshape::hexagon)
0057       binSizePhiTrkAxis_ = 2. * maxAbsPhiTrkAxis_ / (nBinsPhiTrkAxis_ - 1. / 6.);
0058     else if (shape_ == HTshape::diamond)
0059       binSizePhiTrkAxis_ = 2. * maxAbsPhiTrkAxis_ / (nBinsPhiTrkAxis_ - 1. / 2.);
0060 
0061     // Optionally merge 2x2 neighbouring cells into a single cell at low Pt, to reduce efficiency loss due to
0062     // scattering. (Do this if either of options EnableMerge2x2 or MiniHTstage are enabled.
0063     // N.B These two options are never both enabled).
0064     enableMerge2x2_ = (settings->enableMerge2x2() || settings->miniHTstage());
0065     if (settings->miniHTstage()) {
0066       // Mini-HT stage cfg: Merge all bins, irrespective of Pt.
0067       minInvPtToMerge2x2_ = 0.;
0068     } else {
0069       // Merged cells cfg: Merge bins below specified Pt threshold.
0070       minInvPtToMerge2x2_ = 1. / (settings->maxPtToMerge2x2());
0071       if (minInvPtToMerge2x2_ > maxAbsQoverPtAxis_)
0072         enableMerge2x2_ = false;
0073     }
0074 
0075     // Merging 2x2 cells into 1 merged cell is only allowed if HT array dimensions are even.
0076     // (This restriction could be removed along q/Pt axis, since there are also unmerged cells there. But this
0077     // would require correcting the code after each called to mergedCell() below, since
0078     //  "if (i%2 == 1) iStore = i - 1" not correct in this case).
0079     if (enableMerge2x2_ && (nBinsQoverPtAxis_ % 2 != 0 || nBinsPhiTrkAxis_ % 2 != 0))
0080       throw cms::Exception("BadConfig") << "HTrphi: You are not allowed to set EnableMerge2x2 or MiniHTstage = True if "
0081                                            "you have an odd number of bins "
0082                                            "in r-phi HT array "
0083                                         << nBinsQoverPtAxis_ << " " << nBinsPhiTrkAxis_;
0084 
0085     //--- Other options used when filling the HT.
0086 
0087     // Don't fill all HT cells nominally crossed by line corresponding to stub.
0088     killSomeHTCellsRphi_ = settings->killSomeHTCellsRphi();
0089 
0090     // Used to kill excess stubs or tracks that can't be transmitted within time-multiplexed period.
0091     nReceivedStubs_ = 0;
0092     busyInputSectorKill_ = settings_->busyInputSectorKill();  // Kill excess stubs going fron GP to HT?
0093     busySectorKill_ = settings_->busySectorKill();            // Kill excess tracks flowing out of HT?
0094     // Max. num. of stubs that can be sent from GP to HT within TM period
0095     busyInputSectorNumStubs_ = settings_->busyInputSectorNumStubs();
0096     // Max. num. of stubs that can be sent out of HT within TM period
0097     busySectorNumStubs_ = settings_->busySectorNumStubs();
0098     // or individual m bin (=q/Pt) ranges to be output to optical links.
0099     busySectorMbinRanges_ = settings_->busySectorMbinRanges();
0100     // Specifies which m bins should be grouped together by BusySectorMbinRanges. If empty, then they are grouped in order 0,1,2,3,4,5 ...
0101     busySectorMbinOrder_ = settings_->busySectorMbinOrder();
0102     // m bin ranges option disabled if vector empty
0103     busySectorUseMbinRanges_ = (not busySectorMbinRanges_.empty());
0104     busySectorUseMbinOrder_ = (not busySectorMbinOrder_.empty());
0105 
0106     bool rescaleMbins = false;
0107     if (busySectorUseMbinRanges_) {
0108       // Check if the total number of bins specified in cfg option BusySectorMbinRanges corresponds
0109       // to the number of m bins (q/Pt) in the HT. If not, determine how much the ranges must be scaled
0110       // to make this true.
0111       unsigned int nTotalBins = 0;
0112       for (unsigned int j = 0; j < busySectorMbinRanges_.size(); j++) {
0113         nTotalBins += busySectorMbinRanges_[j];
0114       }
0115       rescaleMbins = (nTotalBins != nBinsQoverPtAxis_);
0116       // No rescaling allowed with MBinOrder option.
0117       if (rescaleMbins && busySectorUseMbinOrder_)
0118         throw cms::Exception("BadConfig") << "HTrphi: BusySectorUserMbin error";
0119       float rescaleFactor = rescaleMbins ? float(nBinsQoverPtAxis_) / float(nTotalBins) : 1.;
0120       // Find lower and upper inclusive limits of each m bin range to be sent to a separate optical link.
0121       busySectorMbinLow_.resize(busySectorMbinRanges_.size());
0122       busySectorMbinHigh_.resize(busySectorMbinRanges_.size());
0123       float mBinSum = 0.;
0124       for (unsigned int i = 0; i < busySectorMbinRanges_.size(); i++) {
0125         busySectorMbinLow_[i] = std::round(mBinSum);
0126         busySectorMbinHigh_[i] = std::round(mBinSum + rescaleFactor * busySectorMbinRanges_[i]) - 1;
0127         mBinSum += rescaleFactor * busySectorMbinRanges_[i];
0128       }
0129     }
0130     //
0131     for (unsigned int i = 0; i < nBinsQoverPtAxis_; i++) {
0132       for (unsigned int j = 0; j < nBinsPhiTrkAxis_; j++) {
0133         pair<float, float> helix = this->helix2Dconventional(i, j);  // Get track params at centre of cell.
0134         float qOverPt = helix.first;
0135         // Check if this cell is merged with its neighbours (as in low Pt region).
0136         bool mergedCell = false;
0137         if (enableMerge2x2_ && this->mergedCell(i, j))
0138           mergedCell = true;
0139         // Initialize each cell in HT array.
0140         HTbase::htArray_(i, j) =
0141             std::make_unique<HTcell>(settings, iPhiSec, iEtaReg, etaMinSector, etaMaxSector, qOverPt, i, mergedCell);
0142       }
0143     }
0144 
0145     std::stringstream text;
0146     text << "\n";
0147     text << "=== R-PHI HOUGH TRANSFORM AXES RANGES: abs(q/Pt) < " << maxAbsQoverPtAxis_ << " & abs(track-phi) < "
0148          << maxAbsPhiTrkAxis_ << " ===\n";
0149     text << "=== R-PHI HOUGH TRANSFORM ARRAY SIZE: q/Pt bins = " << nBinsQoverPtAxis_
0150          << " & track-phi bins = " << nBinsPhiTrkAxis_ << " ===\n";
0151     text << "=== R-PHI HOUGH TRANSFORM BIN SIZE: BIN(q/Pt) = " << binSizeQoverPtAxis_
0152          << " & BIN(track-phi) = " << binSizePhiTrkAxis_ << " ===\n\n";
0153     if (busySectorKill_ && busySectorUseMbinRanges_ && rescaleMbins) {
0154       text << "=== R-PHI HOUGH TRANSFORM WARNING: Rescaled m bin ranges specified by cfg parameter "
0155               "BusySectorMbinRanges, as they were inconsistent with total number of m bins in HT.\n";
0156       text << "=== Rescaled values for BusySectorMbinRanges =";
0157       for (unsigned int i = 0; i < busySectorMbinRanges_.size(); i++) {
0158         text << " " << (busySectorMbinHigh_[i] - busySectorMbinLow_[i] + 1);
0159       }
0160     }
0161     text << "\n";
0162     static std::once_flag printOnce;
0163     std::call_once(
0164         printOnce, [](string t) { PrintL1trk() << t; }, text.str());
0165 
0166     // Note helix parameters at the centre of each HT cell.
0167     cellCenters_.clear();
0168     for (unsigned int m = 0; m < nBinsQoverPtAxis_; m++) {
0169       std::vector<std::pair<float, float> > binCenters;
0170       for (unsigned int c = 0; c < nBinsPhiTrkAxis_; c++)
0171         binCenters.push_back(this->helix2Dhough(m, c));
0172       cellCenters_.push_back(binCenters);
0173     }
0174   }
0175 
0176   //=== Add stub to HT array.
0177   //=== If eta subsectors are being used within each sector, specify which ones the stub is compatible with.
0178 
0179   void HTrphi::store(Stub* stub, const vector<bool>& inEtaSubSecs) {
0180     // Optionally, only store stubs that can be sent from GP to HT within TM period.
0181     if ((!busyInputSectorKill_) || (nReceivedStubs_ < busyInputSectorNumStubs_)) {
0182       nReceivedStubs_++;
0183 
0184       unsigned int jPhiTrkBinMinLast = 0;  // Used for error checking
0185       unsigned int jPhiTrkBinMaxLast = 99999;
0186 
0187       // Loop over q/Pt related bins in HT array.
0188       for (unsigned int i = 0; i < nBinsQoverPtAxis_; i++) {
0189         if (shape_ == HTshape::square) {
0190           //--- This is a traditional HT with square cells.
0191 
0192           // In this q/Pt bin, find the range of phi bins that this stub is consistent with.
0193           pair<unsigned int, unsigned int> jRange = this->iPhiRange(stub, i);
0194           unsigned int jPhiTrkBinMin = jRange.first;
0195           unsigned int jPhiTrkBinMax = jRange.second;
0196 
0197           // Store stubs in these cells.
0198           for (unsigned int j = jPhiTrkBinMin; j <= jPhiTrkBinMax; j++) {
0199             bool canStoreStub = true;
0200             unsigned int iStore = i;
0201             unsigned int jStore = j;
0202 
0203             // Optionally merge 2x2 neighbouring cells into a single cell at low Pt, to reduce efficiency loss
0204             // due to scattering.
0205             if (enableMerge2x2_) {
0206               // Check if this cell is merged with its neighbours (as in low Pt region).
0207               if (this->mergedCell(i, j)) {
0208                 // Get location of cell that this cell is merged into (iStore, jStore).
0209                 // Calculation assumes HT array has even number of bins in both dimensions.
0210                 if (i % 2 == 1)
0211                   iStore = i - 1;
0212                 if (j % 2 == 1)
0213                   jStore = j - 1;
0214                 // If this stub was already stored in this merged 2x2 cell, then don't store it again.
0215                 if (HTbase::htArray_(iStore, jStore)->stubStoredInCell(stub))
0216                   canStoreStub = false;
0217               }
0218             }
0219 
0220             if (canStoreStub)
0221               HTbase::htArray_(iStore, jStore)->store(stub, inEtaSubSecs);  // Calls HTcell::store()
0222           }
0223 
0224           // Check that limitations of firmware would not prevent stub being stored correctly in this HT column.
0225           if (errMon_ != nullptr) {
0226             this->countFirmwareErrors(i, jPhiTrkBinMin, jPhiTrkBinMax, jPhiTrkBinMinLast, jPhiTrkBinMaxLast);
0227             jPhiTrkBinMinLast = jPhiTrkBinMin;
0228             jPhiTrkBinMaxLast = jPhiTrkBinMax;
0229           }
0230 
0231         } else {
0232           //--- This is are novel HT with unusual shaped cells.
0233 
0234           if (shape_ == HTshape::diamond) {
0235             //--- This HT has diamond shaped cells.
0236 
0237             float qOverPtBin = -maxAbsQoverPtAxis_ + i * binSizeQoverPtAxis_;
0238             float phiTrk = reco::deltaPhi(stub->phi(), phiCentreSector_) +
0239                            invPtToDphi_ * qOverPtBin * (stub->r() - chosenRofPhi_) + maxAbsPhiTrkAxis_;
0240             if (i % 2 == 0)
0241               phiTrk += binSizePhiTrkAxis_ / 2.;
0242             unsigned int binCenter = std::floor(phiTrk / binSizePhiTrkAxis_);
0243             if (binCenter < nBinsPhiTrkAxis_)
0244               HTbase::htArray_(i, binCenter)->store(stub, inEtaSubSecs);
0245 
0246           } else if (shape_ == HTshape::hexagon) {
0247             //--- This HT has hexagonal cells (with two of its sides parallel to the phi axis).
0248 
0249             float qOverPtBin = -maxAbsQoverPtAxis_ + i * binSizeQoverPtAxis_;
0250             float qOverPtBinVar = binSizeQoverPtAxis_;
0251             float phiTrk = reco::deltaPhi(stub->phi(), phiCentreSector_) +
0252                            invPtToDphi_ * qOverPtBin * (stub->r() - chosenRofPhi_) + maxAbsPhiTrkAxis_;
0253             float phiTrkVar = invPtToDphi_ * qOverPtBinVar * std::abs(stub->r() - chosenRofPhi_);
0254             float phiTrkMin = phiTrk - phiTrkVar;
0255             float phiTrkMax = phiTrk + phiTrkVar;
0256             if (i % 2 == 0)
0257               phiTrk += binSizePhiTrkAxis_ / 6.;
0258             else {
0259               phiTrk -= binSizePhiTrkAxis_ / 3.;
0260               phiTrkMin -= binSizePhiTrkAxis_ / 2.;
0261               phiTrkMax -= binSizePhiTrkAxis_ / 2.;
0262             }
0263             unsigned int iCenter = std::floor(phiTrk / binSizePhiTrkAxis_ * 3.);
0264             unsigned int iMin = std::floor(phiTrkMin / binSizePhiTrkAxis_ * 3.);
0265             unsigned int iMax = std::floor(phiTrkMax / binSizePhiTrkAxis_ * 3.);
0266             std::pair<bool, unsigned int> binCenter;
0267             std::pair<bool, unsigned int> binMin;
0268             std::pair<bool, unsigned int> binMax;
0269             binCenter.second = iCenter / 3;
0270             binMin.second = iMin / 3;
0271             binMax.second = iMax / 3;
0272             binCenter.first = !(iCenter % 3 == 2);
0273             binMin.first = (iMin % 3 == 0);
0274             binMax.first = (iMax % 3 == 0);
0275             if (binCenter.first && binCenter.second < nBinsPhiTrkAxis_)
0276               HTbase::htArray_(i, binCenter.second)->store(stub, inEtaSubSecs);
0277             else if (binMin.first && binMin.second < nBinsPhiTrkAxis_)
0278               HTbase::htArray_(i, binMin.second)->store(stub, inEtaSubSecs);
0279             else if (binMax.first && binMax.second < nBinsPhiTrkAxis_)
0280               HTbase::htArray_(i, binMax.second)->store(stub, inEtaSubSecs);
0281 
0282           } else if (shape_ == HTshape::brick) {
0283             //--- This HT has square cells with alternate rows shifted horizontally by 0.5*cell_width.
0284 
0285             float qOverPtBin = -maxAbsQoverPtAxis_ + i * binSizeQoverPtAxis_;
0286             float qOverPtBinVar = binSizeQoverPtAxis_;
0287             float phiTrk = reco::deltaPhi(stub->phi(), phiCentreSector_) +
0288                            invPtToDphi_ * qOverPtBin * (stub->r() - chosenRofPhi_) + maxAbsPhiTrkAxis_;
0289             float phiTrkVar = invPtToDphi_ * qOverPtBinVar * std::abs(stub->r() - chosenRofPhi_);
0290             float phiTrkMin = phiTrk - phiTrkVar;
0291             float phiTrkMax = phiTrk + phiTrkVar;
0292             unsigned int iMin = std::floor(phiTrkMin / binSizePhiTrkAxis_ * 2.);
0293             unsigned int iMax = std::floor(phiTrkMax / binSizePhiTrkAxis_ * 2.);
0294             std::pair<bool, unsigned int> binMin;
0295             std::pair<bool, unsigned int> binMax;
0296             binMin.second = iMin / 2;
0297             binMax.second = iMax / 2;
0298             binMin.first = (iMin % 2 == i % 2);
0299             binMax.first = (iMax % 2 == i % 2);
0300             if (binMin.first && binMin.second < nBinsPhiTrkAxis_)
0301               HTbase::htArray_(i, binMin.second)->store(stub, inEtaSubSecs);
0302             else if (binMax.first && binMax.second < nBinsPhiTrkAxis_)
0303               HTbase::htArray_(i, binMax.second)->store(stub, inEtaSubSecs);
0304           }
0305         }
0306       }
0307       // Note max. |gradient| that the line corresponding to any stub in any of the r-phi HT arrays could have.
0308       // Firmware assumes this should not exceed 1.0;
0309       if (errMon_ != nullptr) {
0310         errMon_->maxLineGradient = max(errMon_->maxLineGradient.load(), this->calcLineGradArray(stub->r()));
0311       }
0312     }
0313   }
0314 
0315   //=== Determine the m-bin (q/pt) range the specified track is in. (Used if outputting each m bin range on a different opto-link).
0316 
0317   unsigned int HTrphi::getMbinRange(const L1track2D& trk) const {
0318     if (busySectorUseMbinRanges_) {
0319       unsigned int mBin = trk.cellLocationHT().first;
0320       unsigned int mBinOrder;
0321       if (busySectorUseMbinOrder_) {
0322         // User wants to group bins in a wierd order.
0323         mBinOrder = 99999;
0324         for (unsigned int k = 0; k < busySectorMbinOrder_.size(); k++) {
0325           if (mBin == busySectorMbinOrder_[k])
0326             mBinOrder = k;
0327         }
0328         if (mBinOrder == 99999)
0329           throw cms::Exception("LogicError") << "HTrphi::getMbinRange() mBinOrder calculation wrong.";
0330       } else {
0331         // User grouping bins in numerical order 0,1,2,3,4,5...
0332         mBinOrder = mBin;
0333       }
0334       for (unsigned int i = 0; i < busySectorMbinRanges_.size(); i++) {
0335         if (mBinOrder >= busySectorMbinLow_[i] && mBinOrder <= busySectorMbinHigh_[i])
0336           return i;
0337       }
0338       throw cms::Exception("LogicError") << "HTrphi::getMbinRange() messed up";
0339     } else {
0340       return 0;
0341     }
0342   }
0343 
0344   //=== For a given Q/Pt bin, find the range of phi bins that a given stub is consistent with.
0345   //=== Return as a pair (min bin, max bin)
0346   //=== If it range lies outside the HT array, then the min bin will be set larger than the max bin.
0347 
0348   pair<unsigned int, unsigned int> HTrphi::iPhiRange(const Stub* stub, unsigned int iQoverPtBin, bool debug) const {
0349     // Note q/Pt value corresponding to centre of this bin.
0350     float qOverPtBin = -maxAbsQoverPtAxis_ + (iQoverPtBin + 0.5) * binSizeQoverPtAxis_;
0351     // Note change in this q/Pt value needed to reach either edge of the bin.
0352     float qOverPtBinVar = 0.5 * binSizeQoverPtAxis_;
0353 
0354     // Reducing effective bin width can reduce fake rate.
0355     //qOverPtVar = 0.4*binSizeQoverPtAxis_;
0356 
0357     // Calculate range of track-phi that would allow a track in this q/Pt range to pass through the stub.
0358     float phiTrk = stub->phi() + invPtToDphi_ * qOverPtBin * (stub->r() - chosenRofPhi_);
0359     // The next line does the phiTrk calculation without the usual approximation, but it doesn't
0360     // improve performance.
0361     //float phiTrk    = stub->phi() + asin(invPtToDphi_ * qOverPtBin * stub->r()) - asin(invPtToDphi_ * qOverPtBin * chosenRofPhi_);
0362     float phiTrkVar = invPtToDphi_ * qOverPtBinVar * std::abs(stub->r() - chosenRofPhi_);
0363     float phiTrkMin = phiTrk - phiTrkVar;
0364     float phiTrkMax = phiTrk + phiTrkVar;
0365 
0366     float deltaPhiMin = reco::deltaPhi(phiTrkMin, phiCentreSector_);  // Offset to centre of sector.
0367     float deltaPhiMax = reco::deltaPhi(phiTrkMax, phiCentreSector_);
0368     pair<float, float> phiTrkRange(deltaPhiMin, deltaPhiMax);
0369 
0370     // Determine which HT array cell range in track-phi this range "phiTrkRange" corresponds to.
0371     pair<unsigned int, unsigned int> iPhiTrkBinRange = this->HTbase::convertCoordRangeToBinRange(
0372         phiTrkRange, nBinsPhiTrkAxis_, (-maxAbsPhiTrkAxis_), binSizePhiTrkAxis_, killSomeHTCellsRphi_);
0373 
0374     return iPhiTrkBinRange;
0375   }
0376 
0377   //=== Check that limitations of firmware would not prevent stub being stored correctly in this HT column.
0378 
0379   void HTrphi::countFirmwareErrors(unsigned int iQoverPtBin,
0380                                    unsigned int jPhiTrkBinMin,
0381                                    unsigned int jPhiTrkBinMax,
0382                                    unsigned int jPhiTrkBinMinLast,
0383                                    unsigned int jPhiTrkBinMaxLast) {
0384     // Only do check if stub is being stored somewhere in this HT column.
0385     if (jPhiTrkBinMax >= jPhiTrkBinMin) {
0386       //--- Remaining code below checks that firmware could successfully store this stub in this column.
0387       //   (a) Does cell lie NE, E or SE of cell filled in previous column?
0388       bool OK_a = (jPhiTrkBinMin + 1 >= jPhiTrkBinMinLast) && (jPhiTrkBinMax <= jPhiTrkBinMaxLast + 1);
0389       //   (b) Are no more than 2 cells filled in this column
0390       bool OK_b = (jPhiTrkBinMax - jPhiTrkBinMin + 1 <= 2);
0391 
0392       if (!OK_a)
0393         errMon_->numErrorsTypeA++;
0394       if (!OK_b)
0395         errMon_->numErrorsTypeB++;
0396       errMon_->numErrorsNorm++;  // No. of times a stub is added to an HT column.
0397     }
0398   }
0399 
0400   //=== Get the values of the track helix params corresponding to middle of a specified HT cell (i,j).
0401   //=== The helix parameters returned will be those corresponding to the two axes of the HT array.
0402   //=== So they might be (q/pt, phi0) or (q/pt, phi65) etc. depending on the configuration.
0403 
0404   pair<float, float> HTrphi::helix2Dhough(unsigned int i, unsigned int j) const {
0405     unsigned int qOverPtBin = i;
0406     unsigned int phiTrkBin = j;
0407 
0408     // If using merged 2x2 cells in low Pt parts of array, must correct for this.
0409     bool merged = false;
0410     if (enableMerge2x2_) {
0411       // Check if this cell is merged with its neighbours (as in low Pt region).
0412       if (this->mergedCell(i, j)) {
0413         merged = true;
0414         // Get location of cell that this cell is merged into (iStore, jStore).
0415         // Calculation assumes HT array has even number of bins in both dimensions.
0416         if (i % 2 == 1)
0417           qOverPtBin = i - 1;
0418         if (j % 2 == 1)
0419           phiTrkBin = j - 1;
0420       }
0421     }
0422 
0423     float qOverPtBinCenter = .5;
0424     float phiTrkBinCenter = .5;
0425 
0426     if (shape_ != HTshape::square) {
0427       qOverPtBinCenter = 0.;
0428 
0429       float evenPhiPos = 0., oddPhiPos = 0.;
0430       if (shape_ == HTshape::hexagon) {
0431         evenPhiPos = 1. / 6.;
0432         oddPhiPos = 2. / 3.;
0433       } else if (shape_ == HTshape::diamond) {
0434         evenPhiPos = 0.;
0435         oddPhiPos = 0.5;
0436       } else if (shape_ == HTshape::brick) {
0437         evenPhiPos = 0.25;
0438         oddPhiPos = 0.75;
0439       }
0440       phiTrkBinCenter = (qOverPtBin % 2 == 0) ? evenPhiPos : oddPhiPos;
0441     }
0442 
0443     float qOverPt = -maxAbsQoverPtAxis_ + (qOverPtBin + qOverPtBinCenter) * binSizeQoverPtAxis_;
0444     float phiTrk = -maxAbsPhiTrkAxis_ + (phiTrkBin + phiTrkBinCenter) * binSizePhiTrkAxis_;
0445 
0446     if (merged) {
0447       qOverPt += 0.5 * binSizeQoverPtAxis_;
0448       phiTrk += 0.5 * binSizePhiTrkAxis_;
0449     }
0450 
0451     // Correct phiTrk to centre of sector, taking care of 2*pi wrapping
0452     phiTrk = reco::deltaPhi(phiTrk + phiCentreSector_, 0.);
0453     return pair<float, float>(qOverPt, phiTrk);
0454   }
0455 
0456   //=== Get the values of the track helix params corresponding to middle of a specified HT cell (i,j).
0457   //=== The helix parameters returned will be always be (q/pt, phi0), irrespective of how the axes
0458   //=== of the HT array are defined.
0459 
0460   pair<float, float> HTrphi::helix2Dconventional(unsigned int i, unsigned int j) const {
0461     // Get the helix parameters corresponding to the axes definitions of the HT.
0462     pair<float, float> helix2Dht = this->helix2Dhough(i, j);
0463     // Convert to the conventionally agreed pair of helix parameters, (q/pt, phi0).
0464     float qOverPt = helix2Dht.first;  // easy
0465     // If HT defined track phi other than at r=0, must correct to get phi0. Allow for 2*pi wrapping of phi.
0466     float phi0 = reco::deltaPhi(helix2Dht.second + invPtToDphi_ * chosenRofPhi_ * qOverPt, 0.);
0467     return pair<float, float>(qOverPt, phi0);
0468   }
0469 
0470   //=== Which cell in HT array should this TP be in, based on its true trajectory?
0471   //=== (If TP is outside HT array, it it put in the closest bin inside it).
0472 
0473   pair<unsigned int, unsigned int> HTrphi::trueCell(const TP* tp) const {
0474     // Get HT axis variables corresponding to this TP.
0475     float qOverPt = tp->qOverPt();
0476     float phiTrk = tp->trkPhiAtR(chosenRofPhi_);
0477     // Measure phi relative to centre of sector.
0478     float deltaPhi = reco::deltaPhi(phiTrk, phiCentreSector_);
0479     // Convert to bin numbers inside HT array.
0480     int iQoverPt = floor((qOverPt - (-maxAbsQoverPtAxis_)) / binSizeQoverPtAxis_);
0481     int iPhiTrk = floor((deltaPhi - (-maxAbsPhiTrkAxis_)) / binSizePhiTrkAxis_);
0482     // Check if this cell was within the HT array.
0483     if (iQoverPt >= 0 && iQoverPt < int(nBinsQoverPtAxis_) && iPhiTrk >= 0 && iPhiTrk < int(nBinsPhiTrkAxis_)) {
0484       // Check if this cell is merged with its neighbours (as in low Pt region), and if so return merged cell location.
0485       // New: because 2nd stage mini HT may recreate tracks from merged cells with finer cell granularity, one can't predict
0486       //      if a merged cell was used to create a track merely by looking at its cell location.
0487       //      So instead ask L1track3D, which knows if it was created from a merged HT cell or not.
0488       ;
0489     } else {
0490       // TP is not in this HT array at all. Flag this by setting "outside" bin index to 0 (Nbins) if outside array below (above).
0491       if (iQoverPt < 0)
0492         iQoverPt = 0;
0493       if (iQoverPt >= int(nBinsQoverPtAxis_))
0494         iQoverPt = nBinsQoverPtAxis_ - 1;
0495       if (iPhiTrk < 0)
0496         iPhiTrk = 0;
0497       if (iPhiTrk >= int(nBinsPhiTrkAxis_))
0498         iPhiTrk = nBinsPhiTrkAxis_ - 1;
0499     }
0500     return pair<unsigned int, unsigned int>(iQoverPt, iPhiTrk);
0501   }
0502 
0503   //=== Which cell in HT array should this fitted track be in, based on its fitted trajectory?
0504   //=== Always uses beam-spot constrained trajectory if available.
0505   //=== (If fitted track is outside HT array, it it put in the closest bin inside it).
0506 
0507   pair<unsigned int, unsigned int> HTrphi::cell(const L1fittedTrack* fitTrk) const {
0508     bool beamConstraint = fitTrk->done_bcon();  // Is beam-spot constraint available? (e.g. 5 param helix fit)
0509     // Get HT axis variables corresponding to this fitted track.
0510     float qOverPt = beamConstraint ? fitTrk->qOverPt_bcon() : fitTrk->qOverPt();
0511     // Convert phi0 to phi at chosen radius used in HT.
0512     float phiTrk = fitTrk->phiAtChosenR(beamConstraint);
0513     // Measure phi relative to centre of sector.
0514     float deltaPhi = reco::deltaPhi(phiTrk, phiCentreSector_);
0515     // Convert to bin numbers inside HT array.
0516     int iQoverPt = 999999;
0517     int iPhiTrk = 999999;
0518 
0519     if (shape_ == HTshape::square) {
0520       //--- This is a traditional HT with square cells.
0521 
0522       iQoverPt = floor((qOverPt - (-maxAbsQoverPtAxis_)) / binSizeQoverPtAxis_);
0523       iPhiTrk = floor((deltaPhi - (-maxAbsPhiTrkAxis_)) / binSizePhiTrkAxis_);
0524 
0525       // Check if this cell was within the HT array.
0526       if (iQoverPt >= 0 && iQoverPt < int(nBinsQoverPtAxis_) && iPhiTrk >= 0 && iPhiTrk < int(nBinsPhiTrkAxis_)) {
0527         // Check if this cell is merged with its neighbours (as in low Pt region), and if so return merged cell location.
0528         // New: because 2nd stage mini HT may recreate tracks from merged cells with finer cell granularity, one can't predict
0529         //      if a merged cell was used to create a track merely by looking at its cell location.
0530         //      So instead ask L1track3D, which knows if it was created from a merged HT cell or not.
0531         ;
0532       } else {
0533         // Fitted track is not in this HT array at all. Flag this by setting "outside" bin index to 0 (Nbins-1) if outside array below (above).
0534         if (iQoverPt < 0)
0535           iQoverPt = 0;
0536         if (iQoverPt >= int(nBinsQoverPtAxis_))
0537           iQoverPt = nBinsQoverPtAxis_ - 1;
0538         if (iPhiTrk < 0)
0539           iPhiTrk = 0;
0540         if (iPhiTrk >= int(nBinsPhiTrkAxis_))
0541           iPhiTrk = nBinsPhiTrkAxis_ - 1;
0542       }
0543 
0544     } else {
0545       //--- This is are novel HT with unusual shaped cells.
0546 
0547       float minD = std::numeric_limits<float>::infinity();
0548       float d(0);
0549       unsigned int m(0);
0550       for (const auto& binCenters : cellCenters_) {
0551         unsigned int c(0);
0552         for (auto cellCenter : binCenters) {
0553           d = std::pow((cellCenter.first - qOverPt) / (float)binSizeQoverPtAxis_, 2) +
0554               std::pow((cellCenter.second - phiTrk) / (float)binSizePhiTrkAxis_, 2);
0555           if (d < minD) {
0556             minD = d;
0557             iQoverPt = m;
0558             iPhiTrk = c;
0559           }
0560           c++;
0561         }
0562         m++;
0563       }
0564       // Fitted track is not in this HT array at all. Flag this by setting "outside" bin index to 0 (Nbins-1) if outside array below (above).
0565       if (iQoverPt < 0)
0566         iQoverPt = 0;
0567       if (iQoverPt >= int(nBinsQoverPtAxis_))
0568         iQoverPt = nBinsQoverPtAxis_ - 1;
0569       if (iPhiTrk < 0)
0570         iPhiTrk = 0;
0571       if (iPhiTrk >= int(nBinsPhiTrkAxis_))
0572         iPhiTrk = nBinsPhiTrkAxis_ - 1;
0573     }
0574 
0575     return pair<unsigned int, unsigned int>(iQoverPt, iPhiTrk);
0576   }
0577 
0578   //=== Check if specified cell is merged with its 2x2 neighbours into a single cell,
0579   //=== as it is in low Pt region.
0580 
0581   bool HTrphi::mergedCell(unsigned int iQoverPtBin, unsigned int jPhiTrkBin) const {
0582     bool merge = false;
0583 
0584     if (enableMerge2x2_) {
0585       unsigned int i = iQoverPtBin;
0586       //unsigned int j = jPhiTrkBin;
0587 
0588       // Calculate number of merged bins on each q/Pt side of array.
0589       float fMergeBins = (maxAbsQoverPtAxis_ - minInvPtToMerge2x2_) / (2. * binSizeQoverPtAxis_);
0590       // Number of unmerged bins this corresponds to, which must be even, since each merged bin comprises two normal q/pt bins.
0591       unsigned int numQoverPtBinsToMerge = 2 * min((unsigned int)(std::round(fMergeBins)), (nBinsQoverPtAxis_ / 4));
0592       const float small = 0.001;
0593       if (minInvPtToMerge2x2_ < small && (unsigned int)(std::round(2. * fMergeBins)) % 2 == 1)
0594         numQoverPtBinsToMerge++;
0595       unsigned int iB = (nBinsQoverPtAxis_ - 1) - i;  // Count backwards across array.
0596       if (min(i, iB) < numQoverPtBinsToMerge)
0597         merge = true;
0598     }
0599 
0600     return merge;
0601   }
0602 
0603   //=== Calculate line |gradient| of stubs in HT array, so can check it doesn't exceed 1.
0604 
0605   float HTrphi::calcLineGradArray(float r) const {
0606     float grad = std::abs(invPtToDphi_ * (r - chosenRofPhi_));
0607     // Convert it to units of bin width.
0608     grad *= binSizeQoverPtAxis_ / binSizePhiTrkAxis_;
0609     if (shape_ == HTshape::hexagon)
0610       grad *= 3.;
0611     else if (shape_ == HTshape::diamond)
0612       grad *= 2.;
0613     else if (shape_ == HTshape::brick)
0614       grad *= 4.;
0615     return grad;
0616   }
0617 
0618   //=== If requested, kill those tracks in this sector that can't be read out during the time-multiplexed period, because
0619   //=== the HT has associated too many stubs to tracks.
0620 
0621   list<L1track2D> HTrphi::killTracksBusySec(const list<L1track2D>& tracks) const {
0622     list<L1track2D> outTracks;
0623 
0624     if (busySectorKill_) {
0625       unsigned int nStubsOut = 0;  // #stubs assigned to tracks in this sector.
0626       // #stubs assigned to each m bin range in this sector.
0627       vector<unsigned int> nStubsOutInRange(busySectorMbinRanges_.size(), 0);
0628 
0629       for (const L1track2D& trk : tracks) {
0630         bool keep = true;
0631         unsigned int nStubs = trk.numStubs();  // #stubs on this track.
0632         if (busySectorUseMbinRanges_) {  // Are tracks from different m bin ranges output seperately to increase bandwidth?
0633           unsigned int mBinRange = this->getMbinRange(trk);  // Which m bin range is this track in?
0634           nStubsOutInRange[mBinRange] += nStubs;
0635           if (nStubsOutInRange[mBinRange] > busySectorNumStubs_)
0636             keep = false;
0637         } else {
0638           nStubsOut += nStubs;
0639           if (nStubsOut > busySectorNumStubs_)
0640             keep = false;
0641         }
0642 
0643         if (keep)
0644           outTracks.push_back(trk);
0645       }
0646 
0647     } else {
0648       outTracks = tracks;
0649     }
0650 
0651     return outTracks;
0652   }
0653 
0654   //=== Define the order in which the hardware processes rows of the HT array when it outputs track candidates.
0655   //=== Currently corresponds to highest Pt tracks first.
0656   //=== If two tracks have the same Pt, the -ve charge one is output before the +ve charge one.
0657 
0658   vector<unsigned int> HTrphi::rowOrder(unsigned int numRows) const {
0659     vector<unsigned int> iOrder;
0660 
0661     // Logic slightly different depending on whether HT array has even or odd number of rows.
0662     const bool oddNumRows = (numRows % 2 == 1);
0663 
0664     // This selects middle rows first before moving to exterior ones.
0665     if (oddNumRows) {
0666       unsigned int middleRow = (numRows - 1) / 2;
0667       iOrder.push_back(middleRow);
0668       for (unsigned int i = 1; i <= (numRows - 1) / 2; i++) {
0669         iOrder.push_back(middleRow - i);  // -ve charge
0670         iOrder.push_back(middleRow + i);  // +ve charge
0671       }
0672     } else {
0673       unsigned int startRowPos = numRows / 2;
0674       unsigned int startRowNeg = startRowPos - 1;
0675       for (unsigned int i = 0; i < numRows / 2; i++) {
0676         iOrder.push_back(startRowNeg - i);  // -ve charge
0677         iOrder.push_back(startRowPos + i);  // +ve charge
0678       }
0679     }
0680 
0681     return iOrder;
0682   }
0683 
0684 }  // namespace tmtt