TrackFindingTMTT/src/Sector.cc

0001 #include "L1Trigger/TrackFindingTMTT/interface/Sector.h"
0002 #include "L1Trigger/TrackFindingTMTT/interface/Stub.h"
0003 #include "L1Trigger/TrackFindingTMTT/interface/Settings.h"
0004
0005 #include "DataFormats/Math/interface/deltaPhi.h"
0006
0007 using namespace std;
0008
0009 namespace tmtt {
0010
0011   //=== Initialise
0012
0013   Sector::Sector(const Settings* settings, unsigned int iPhiSec, unsigned int iEtaReg)
0014       : settings_(settings),
0015         // Sector number
0016         iPhiSec_(iPhiSec),
0017         iEtaReg_(iEtaReg),
0018
0019         beamWindowZ_(settings->beamWindowZ()),  // Assumed half-length of beam-spot
0020
0021         //===  Characteristics of this eta region.
0022         // Using lines of specified rapidity drawn from centre of CMS, determine the z coords at which
0023         // they cross the radius chosenRofZ_.
0024         etaMin_(settings->etaRegions()[iEtaReg]),
0025         etaMax_(settings->etaRegions()[iEtaReg + 1]),
0026         chosenRofZ_(settings->chosenRofZ()),
0027         // Get range in z of tracks covered by this sector at chosen radius from beam-line
0028         zOuterMin_(chosenRofZ_ / tan(2. * atan(exp(-etaMin_)))),
0029         zOuterMax_(chosenRofZ_ / tan(2. * atan(exp(-etaMax_)))),
0030
0031         //=== Characteristics of this phi region.
0032         chosenRofPhi_(settings->chosenRofPhi()),
0033         minPt_(settings->houghMinPt()),  // Min Pt covered by  HT array.
0034         assumedPhiTrkRes_(settings->assumedPhiTrkRes()),
0035         useStubPhi_(settings->useStubPhi()),
0036         useStubPhiTrk_(settings->useStubPhiTrk()),
0037         calcPhiTrkRes_(settings->calcPhiTrkRes()),
0038         //=== Check if subsectors in eta are being used within each sector.
0039         numSubSecsEta_(settings->numSubSecsEta()) {
0040     // Centre of phi (tracking) nonant zero must be along x-axis to be consistent with tracker cabling map.
0041     // Define phi sector zero  to start at lower end of phi range in nonant 0.
0042     float phiCentreSec0 = -M_PI / float(settings->numPhiNonants()) + M_PI / float(settings->numPhiSectors());
0043     // Centre of sector in phi
0044     phiCentre_ = 2. * M_PI * float(iPhiSec) / float(settings->numPhiSectors()) + phiCentreSec0;
0045     sectorHalfWidth_ = M_PI / float(settings->numPhiSectors());  // Sector half width excluding overlaps.
0046
0047     // If eta subsectors have equal width in rapidity, do this.
0048     float subSecWidth = (etaMax_ - etaMin_) / float(numSubSecsEta_);
0049     for (unsigned int i = 0; i < numSubSecsEta_; i++) {
0050       float subSecEtaMin = etaMin_ + i * subSecWidth;
0051       float subSecEtaMax = subSecEtaMin + subSecWidth;
0052       float subSecZmin = chosenRofZ_ / tan(2. * atan(exp(-subSecEtaMin)));
0053       float subSecZmax = chosenRofZ_ / tan(2. * atan(exp(-subSecEtaMax)));
0054       zOuterMinSub_.push_back(subSecZmin);
0055       zOuterMaxSub_.push_back(subSecZmax);
0056     }
0057   }
0058
0059   //=== Check if stub is inside this eta region.
0060
0061   bool Sector::insideEta(const Stub* stub) const {
0062     // Lower edge of this eta region defined by line from (r,z) = (0,-beamWindowZ) to (chosenRofZ_, zOuterMin_).
0063     // Upper edge of this eta region defined by line from (r,z) = (0, beamWindowZ) to (chosenRofZ_, zOuterMax_).
0064
0065     bool inside = this->insideEtaRange(stub, zOuterMin_, zOuterMax_);
0066     return inside;
0067   }
0068
0069   //=== Check if stub is within subsectors in eta that sector may be divided into.
0070
0071   vector<bool> Sector::insideEtaSubSecs(const Stub* stub) const {
0072     if (settings_->enableDigitize() && numSubSecsEta_ == 2) {
0073       // Use (complicated) digitized firmware emulation
0074       return subEtaFwCalc(stub->digitalStub()->iDigi_Rt(), stub->digitalStub()->iDigi_Z());
0075
0076     } else {
0077       // Use (simpler) floating point calculation.
0078
0079       vector<bool> insideVec;
0080
0081       // Loop over subsectors.
0082       for (unsigned int i = 0; i < numSubSecsEta_; i++) {
0083         bool inside = this->insideEtaRange(stub, zOuterMinSub_[i], zOuterMaxSub_[i]);
0084         insideVec.push_back(inside);
0085       }
0086
0087       return insideVec;
0088     }
0089   }
0090
0091   //=== Check if stub is within eta sector or subsector that is delimated by specified zTrk range.
0092
0093   bool Sector::insideEtaRange(const Stub* stub, float zRangeMin, float zRangeMax) const {
0094     // Lower edge of this eta region defined by line from (r,z) = (0,-beamWindowZ) to (chosenRofZ_, zRangeMin).
0095     // Upper edge of this eta region defined by line from (r,z) = (0, beamWindowZ) to (chosenRofZ_, zRangeMax).
0096
0097     float zMin, zMax;
0098     bool inside;
0099
0100     // Calculate z coordinate of lower edge of this eta region, evaluated at radius of stub.
0101     zMin = (zRangeMin * stub->r() - beamWindowZ_ * std::abs(stub->r() - chosenRofZ_)) / chosenRofZ_;
0102     // Calculate z coordinate of upper edge of this eta region, evaluated at radius of stub.
0103     zMax = (zRangeMax * stub->r() + beamWindowZ_ * std::abs(stub->r() - chosenRofZ_)) / chosenRofZ_;
0104
0105     inside = (stub->z() > zMin && stub->z() < zMax);
0106     return inside;
0107   }
0108
0109   //=== Check if stub is inside this phi region.
0110
0111   bool Sector::insidePhi(const Stub* stub) const {
0112     // N.B. The logic here for preventing a stub being assigned to > 2 sectors seems overly agressive.
0113     // But attempts at improving it have failed ...
0114
0115     bool okPhi = true;
0116     bool okPhiTrk = true;
0117
0118     if (useStubPhi_) {
0119       float delPhi =
0120           reco::deltaPhi(stub->phi(), phiCentre_);  // Phi difference between stub & sector in range -PI to +PI.
0121       float tolerancePhi = stub->phiDiff(
0122           chosenRofPhi_, minPt_);  // How much stub phi might differ from track phi because of track curvature.
0123       float outsidePhi = std::abs(delPhi) - sectorHalfWidth_ -
0124                          tolerancePhi;  // If > 0, then stub is not compatible with being inside this sector.
0125       if (outsidePhi > 0)
0126         okPhi = false;
0127     }
0128
0129     if (useStubPhiTrk_) {
0130       // Estimate either phi0 of track from stub info, or phi of the track at radius chosenRofPhi_.
0131       float phiTrk = stub->trkPhiAtR(chosenRofPhi_);
0132       // Phi difference between stub & sector in range -PI to +PI.
0133       float delPhiTrk = reco::deltaPhi(phiTrk, phiCentre_);
0134       // Set tolerance equal to nominal resolution assumed in phiTrk
0135       float tolerancePhiTrk = assumedPhiTrkRes_ * (2 * sectorHalfWidth_);
0136       if (calcPhiTrkRes_) {
0137         // Calculate uncertainty in phiTrk due to poor resolution in stub bend
0138         float phiTrkRes = stub->trkPhiAtRcut(chosenRofPhi_);
0139         // Reduce tolerance if this is smaller than the nominal assumed resolution.
0140         tolerancePhiTrk = min(tolerancePhiTrk, phiTrkRes);
0141       }
0142       // If following > 0, then stub is not compatible with being inside this sector.
0143       float outsidePhiTrk = std::abs(delPhiTrk) - sectorHalfWidth_ - tolerancePhiTrk;
0144
0145       if (outsidePhiTrk > 0)
0146         okPhiTrk = false;
0147     }
0148
0149     return (okPhi && okPhiTrk);
0150   }
0151
0152   //=== For performance studies, note which stubs on given tracking particle are inside the sector.
0153   //=== Returns two booleans for each stub, indicating if they are in phi & eta sectors respectively.
0154   //=== AND them together to get (eta,phi) sector decision.
0155
0156   unordered_map<const Stub*, pair<bool, bool> > Sector::stubsInside(const TP& tp) const {
0157     unordered_map<const Stub*, pair<bool, bool> > inside;
0158     // Loop over stubs produced by tracking particle
0159     const vector<const Stub*>& assStubs = tp.assocStubs();
0160     for (const Stub* stub : assStubs) {
0161       // Check if this stub is inside sector
0162       inside[stub] = pair<bool, bool>(this->insidePhi(stub), this->insideEta(stub));
0163     }
0164     return inside;
0165   }
0166
0167   //=== Count number of stubs in given tracking particle which are inside this (phi,eta) sector;
0168   //=== or inside it if only the eta cuts are applied; or inside it if only the phi cuts are applied.
0169   //=== The results are returned as the 3 last arguments of the function.
0170
0171   void Sector::numStubsInside(const TP& tp,
0172                               unsigned int& nStubsInsideEtaPhi,
0173                               unsigned int& nStubsInsideEta,
0174                               unsigned int& nStubsInsidePhi) const {
0175     nStubsInsideEtaPhi = 0;
0176     nStubsInsideEta = 0;
0177     nStubsInsidePhi = 0;
0178     for (const auto& iter : this->stubsInside(tp)) {
0179       bool insidePhi = iter.second.first;
0180       bool insideEta = iter.second.second;
0181       if (insidePhi && insideEta)
0182         nStubsInsideEtaPhi++;
0183       if (insideEta)
0184         nStubsInsideEta++;
0185       if (insidePhi)
0186         nStubsInsidePhi++;
0187     }
0188   }
0189
0190   // Digitize a floating point number to 2s complement integer, dropping anything after the decimal point. (Kristian Harder)
0191
0192   int64_t Sector::forceBitWidth(const float value, const UInt_t nBits) const {
0193     // slightly hand-waving treatment of 2s complement
0194     int64_t sign = 1;
0195     if (value < 0)
0196       sign = -1;
0197     int64_t iValue = int64_t(std::abs(value));
0198     int64_t mask = (int64_t(1) << nBits) - int64_t(1);
0199     int64_t result = sign * (iValue & mask);
0200     if (std::abs(result - value) > 1)
0201       throw cms::Exception("LogicError")
0202           << "Sector::forceBitWidth is messing up by using too few bits to digitize number"
0203           << " nBits=" << nBits << " Input float=" << value << " Output digi = " << result;
0204     return result;
0205     // Check that result is compatible with value. Throw error if not.
0206   }
0207
0208   //=== Check if stub is within subsectors in eta that sector may be divided into. Uses digitized calculation corresponding to GP firmware. (Kristian Harder)
0209   //=== Modified to configurable number of rT and z digisation bits by Ian, with advice from Luis.
0210
0211   vector<bool> Sector::subEtaFwCalc(const int rT, const int z) const {
0212     // Note number of reference bits used to digitize rT and z, used when GP authors determined some constants below.
0213     unsigned int rtBitsRef = 10;
0214     unsigned int zBitsRef = 12;
0215
0216     // This replaces Kristian's hard-wired constants with configurable ones.
0217     unsigned int rtBits = settings_->rtBits();
0218     unsigned int zBits = settings_->zBits();
0219     float rtRange = settings_->rtRange();
0220     float zRange = settings_->zRange();
0221     constexpr float cm_to_mm = 10.;  // firwmare is in mm and CMSSW in cm.
0222     float zBase = cm_to_mm / (pow(2, zBits) / zRange);
0223     float rTBase = cm_to_mm / (pow(2, rtBits) / rtRange);
0224
0225     // Number of bits used by DSP in UltraScale-Plus FPGA (where DSP does D = A*B + C)
0226     constexpr unsigned int nDSPa = 27;
0227     //constexpr unsigned int nDSPb = 18;
0228     constexpr unsigned int nDSPc = 48;
0229     constexpr unsigned int nDSPd = 48;
0230
0231     // unit transformations: firmware uses mm, software uses cm
0232     float BeamWindow = cm_to_mm * beamWindowZ_;
0233     float T_rphi = cm_to_mm * chosenRofPhi_;
0234     float T_rz = cm_to_mm * chosenRofZ_;
0235
0236     // actual algorithm as used in firmware, mostly using same variable names
0237     float Beam_over_T = BeamWindow / T_rz;
0238     // Value chosen so that number digitized below when calculating "bot" uses most of the nDSPa bits, without overflowing them. This is done assuming reference number of bits for rT and z mentioned above.
0239     unsigned int nShiftA = 24;
0240     // Guess from to keep "bot" in correct range (nDSPa) if number of digitsation bits are changed.
0241     nShiftA += (rtBits - rtBitsRef) - (zBits - zBitsRef);
0242     float Beam_over_T_base = 1. / (1 << nShiftA);
0243     int64_t bot = forceBitWidth(Beam_over_T * rTBase / zBase / Beam_over_T_base, nDSPa);
0244     int64_t bw = forceBitWidth(BeamWindow / zBase / Beam_over_T_base, nDSPc);
0245     float etaSecMid = (settings_->etaRegions()[iEtaReg_] + settings_->etaRegions()[iEtaReg_ + 1]) / 2.0;
0246     float tanlSecMid = 1.0 / tan(2.0 * atan(exp(-etaSecMid)));
0247     // Value chosen so that number digitized below when calculating "tanlSec_Mid" uses most of the nDSPa bits, without overflowing them. This is done assuming reference number of bits for rT and z mentioned above.
0248     unsigned int nShiftB = 16;
0249     // Guess to keep "tanlSec_Mid" in correct range (nDSPa) if number of digitsation bits are changed.
0250     nShiftB += (rtBits - rtBitsRef) - (zBits - zBitsRef);
0251     float tanlSecBase = 1. / (1 << nShiftB);
0252     int64_t tanlSec_Mid = forceBitWidth(int(tanlSecMid * rTBase / zBase / tanlSecBase), nDSPa);
0253     // Number of extra bits used to digitise r instead of rT within GP code, if both encoded as signed int.
0254     constexpr unsigned int nExtraBitsR = 2;
0255     unsigned int rBits = rtBits + nExtraBitsR;
0256     int64_t r = forceBitWidth(rT + T_rphi / rTBase, rBits);
0257     int64_t g = forceBitWidth(bot * r - bw, nDSPd);
0258     int64_t absg = abs(g);
0259     // Number of useful bits left of the nDSPd assigned to "absg" after right-shifting by nShiftA bits.
0260     const unsigned nBitsRemainingA = nDSPd - nShiftA;
0261     int64_t shift_g = forceBitWidth((absg >> nShiftA), nBitsRemainingA);
0262     // Number of bits is sum of those in two numbers being multiplied.
0263     int64_t tlsr = forceBitWidth(tanlSec_Mid * r, nDSPa + rBits);
0264     // Number of useful bits left of (nDSPa + rBits) assigned to "tlsr" after right-shifting by nShiftB bits.
0265     const unsigned nBitsRemainingB = (nDSPa + rBits) - nShiftB;
0266     int64_t shift_tlsr = forceBitWidth((tlsr >> nShiftB), nBitsRemainingB);
0267
0268     vector<bool> insideVec;
0269     insideVec.push_back(z <= (shift_tlsr + shift_g));
0270     insideVec.push_back(z >= (shift_tlsr - shift_g));
0271     return insideVec;
0272   }
0273
0274 }  // namespace tmtt