Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-12-12 02:57:19

0001 //////////////////////////////////////////////////////////////////
0002 // MatchProcessor
0003 //
0004 // This module is the combined version of the PR+ME+MC
0005 // See more in execute()
0006 //
0007 // Variables such as `best_ideltaphi_barrel` store the "global"
0008 // best value for delta phi, r, z, and r*phi, for instances
0009 // where the same tracklet has multiple stub pairs. This allows
0010 // us to find the truly best match
0011 //////////////////////////////////////////////////////////////////
0012 
0013 #include "L1Trigger/TrackFindingTracklet/interface/MatchProcessor.h"
0014 #include "L1Trigger/TrackFindingTracklet/interface/Globals.h"
0015 #include "L1Trigger/TrackFindingTracklet/interface/Util.h"
0016 #include "L1Trigger/TrackFindingTracklet/interface/HistBase.h"
0017 
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 #include "FWCore/Utilities/interface/Exception.h"
0020 #include "DataFormats/Math/interface/deltaPhi.h"
0021 #include "L1Trigger/TrackFindingTracklet/interface/IMATH_TrackletCalculator.h"
0022 
0023 #include <filesystem>
0024 
0025 using namespace std;
0026 using namespace trklet;
0027 
0028 MatchProcessor::MatchProcessor(string name, Settings const& settings, Globals* global)
0029     : ProcessBase(name, settings, global),
0030       phimatchcuttable_(settings),
0031       zmatchcuttable_(settings),
0032       rphicutPStable_(settings),
0033       rphicut2Stable_(settings),
0034       rcutPStable_(settings),
0035       rcut2Stable_(settings),
0036       alphainner_(settings),
0037       alphaouter_(settings),
0038       rSSinner_(settings),
0039       rSSouter_(settings),
0040       diskRadius_(settings),
0041       fullmatches_(12),
0042       rinvbendlut_(settings),
0043       luttable_(settings),
0044       inputProjBuffer_(3) {
0045   phiregion_ = name[8] - 'A';
0046 
0047   layerdisk_ = initLayerDisk(3);
0048 
0049   barrel_ = layerdisk_ < N_LAYER;
0050 
0051   phishift_ = settings_.nphibitsstub(N_LAYER - 1) - settings_.nphibitsstub(layerdisk_);
0052   dzshift_ = settings_.nzbitsstub(0) - settings_.nzbitsstub(layerdisk_);
0053 
0054   if (barrel_) {
0055     icorrshift_ = ilog2(settings_.kphi(layerdisk_) / (settings_.krbarrel() * settings_.kphider()));
0056     icorzshift_ = ilog2(settings_.kz(layerdisk_) / (settings_.krbarrel() * settings_.kzder()));
0057   } else {
0058     icorrshift_ = ilog2(settings_.kphi(layerdisk_) / (settings_.kz() * settings_.kphiderdisk()));
0059     icorzshift_ = ilog2(settings_.krprojshiftdisk() / (settings_.kz() * settings_.krder()));
0060   }
0061 
0062   luttable_.initBendMatch(layerdisk_);
0063 
0064   nrbits_ = 5;
0065   nphiderbits_ = 6;
0066 
0067   nrprojbits_ = 8;
0068 
0069   if (!barrel_) {
0070     rinvbendlut_.initProjectionBend(
0071         global->ITC_L1L2()->der_phiD_final.K(), layerdisk_ - N_LAYER, nrbits_, nphiderbits_);
0072   }
0073 
0074   nrinv_ = NRINVBITS;
0075 
0076   unsigned int region = getName()[8] - 'A';
0077   assert(region < settings_.nallstubs(layerdisk_));
0078 
0079   if (barrel_) {
0080     phimatchcuttable_.initmatchcut(layerdisk_, TrackletLUT::MatchType::barrelphi, region);
0081     zmatchcuttable_.initmatchcut(layerdisk_, TrackletLUT::MatchType::barrelz, region);
0082   } else {
0083     rphicutPStable_.initmatchcut(layerdisk_, TrackletLUT::MatchType::diskPSphi, region);
0084     rphicut2Stable_.initmatchcut(layerdisk_, TrackletLUT::MatchType::disk2Sphi, region);
0085     rcutPStable_.initmatchcut(layerdisk_, TrackletLUT::MatchType::diskPSr, region);
0086     rcut2Stable_.initmatchcut(layerdisk_, TrackletLUT::MatchType::disk2Sr, region);
0087     alphainner_.initmatchcut(layerdisk_, TrackletLUT::MatchType::alphainner, region);
0088     alphaouter_.initmatchcut(layerdisk_, TrackletLUT::MatchType::alphaouter, region);
0089     rSSinner_.initmatchcut(layerdisk_, TrackletLUT::MatchType::rSSinner, region);
0090     rSSouter_.initmatchcut(layerdisk_, TrackletLUT::MatchType::rSSouter, region);
0091     diskRadius_.initProjectionDiskRadius(nrprojbits_);
0092   }
0093 
0094   for (unsigned int i = 0; i < N_DSS_MOD * 2; i++) {
0095     ialphafactinner_[i] = (1 << settings_.alphashift()) * settings_.krprojshiftdisk() * settings_.half2SmoduleWidth() /
0096                           (1 << (settings_.nbitsalpha() - 1)) / (settings_.rDSSinner(i) * settings_.rDSSinner(i)) /
0097                           settings_.kphi();
0098     ialphafactouter_[i] = (1 << settings_.alphashift()) * settings_.krprojshiftdisk() * settings_.half2SmoduleWidth() /
0099                           (1 << (settings_.nbitsalpha() - 1)) / (settings_.rDSSouter(i) * settings_.rDSSouter(i)) /
0100                           settings_.kphi();
0101   }
0102 
0103   nvm_ = settings_.nvmme(layerdisk_) * settings_.nallstubs(layerdisk_);
0104   nvmbins_ = settings_.nvmme(layerdisk_);
0105   nvmbits_ = settings_.nbitsvmme(layerdisk_) + settings_.nbitsallstubs(layerdisk_);
0106 
0107   nMatchEngines_ = 4;
0108   for (unsigned int iME = 0; iME < nMatchEngines_; iME++) {
0109     MatchEngineUnit tmpME(settings_, barrel_, layerdisk_, luttable_);
0110     tmpME.setimeu(iME);
0111     matchengines_.push_back(tmpME);
0112   }
0113 
0114   // Pick some initial large values
0115   best_ideltaphi_barrel = 0xFFFF;
0116   best_ideltaz_barrel = 0xFFFF;
0117   best_ideltaphi_disk = 0xFFFF;
0118   best_ideltar_disk = 0xFFFF;
0119   curr_tracklet = nullptr;
0120   next_tracklet = nullptr;
0121 }
0122 
0123 void MatchProcessor::addOutput(MemoryBase* memory, string output) {
0124   if (settings_.writetrace()) {
0125     edm::LogVerbatim("Tracklet") << "In " << name_ << " adding output to " << memory->getName() << " to output "
0126                                  << output;
0127   }
0128   if (output.find("matchout") != std::string::npos) {
0129     auto* tmp = dynamic_cast<FullMatchMemory*>(memory);
0130     assert(tmp != nullptr);
0131     unsigned int iSeed = getISeed(tmp->getName());
0132     assert(iSeed < fullmatches_.size());
0133     assert(fullmatches_[iSeed] == nullptr);
0134     fullmatches_[iSeed] = tmp;
0135     return;
0136   }
0137   throw cms::Exception("BadConfig") << __FILE__ << " " << __LINE__ << " could not find output: " << output;
0138 }
0139 
0140 void MatchProcessor::addInput(MemoryBase* memory, string input) {
0141   if (settings_.writetrace()) {
0142     edm::LogVerbatim("Tracklet") << "In " << name_ << " adding input from " << memory->getName() << " to input "
0143                                  << input;
0144   }
0145   if (input == "allstubin") {
0146     auto* tmp = dynamic_cast<AllStubsMemory*>(memory);
0147     assert(tmp != nullptr);
0148     allstubs_ = tmp;
0149     return;
0150   }
0151   if (input == "vmstubin") {
0152     auto* tmp = dynamic_cast<VMStubsMEMemory*>(memory);
0153     assert(tmp != nullptr);
0154     vmstubs_.push_back(tmp);  //to allow more than one stub in?  vmstubs_=tmp;
0155     return;
0156   }
0157   if (input == "projin") {
0158     auto* tmp = dynamic_cast<TrackletProjectionsMemory*>(memory);
0159     assert(tmp != nullptr);
0160     inputprojs_.push_back(tmp);
0161     return;
0162   }
0163   throw cms::Exception("BadConfig") << __FILE__ << " " << __LINE__ << " could not find input: " << input;
0164 }
0165 
0166 void MatchProcessor::execute(unsigned int iSector, double phimin) {
0167   assert(vmstubs_.size() == 1);
0168 
0169   /*
0170     The code is organized in three 'steps' corresponding to the PR, ME, and MC functions. The output from
0171     the PR step is buffered in a 'circular' buffer, and similarly the ME output is put in a circular buffer.     
0172     The implementation is done in steps, emulating what can be done in firmware. On each step we do:
0173     
0174     1) A projection is read and if there is space it is insert into the inputProjBuffer_
0175     
0176     2) Process next match in the ME - if there is an idle ME the next projection is inserted
0177     
0178     3) Readout match from ME and send to match calculator
0179 
0180     However, for the pipelining to work in HLS these blocks are executed in reverse order
0181     
0182   */
0183 
0184   bool print = getName() == "MP_L3PHIC" && iSector == 3;
0185   print = false;
0186 
0187   phimin_ = phimin;
0188 
0189   Tracklet* oldTracklet = nullptr;
0190 
0191   unsigned int countme = 0;
0192   unsigned int countall = 0;
0193   unsigned int countsel = 0;
0194   unsigned int countinputproj = 0;
0195 
0196   unsigned int iprojmem = 0;
0197   while (iprojmem < inputprojs_.size() && inputprojs_[iprojmem]->nTracklets() == 0) {
0198     iprojmem++;
0199   }
0200 
0201   unsigned int iproj = 0;
0202 
0203   inputProjBuffer_.reset();
0204 
0205   for (const auto& inputproj : inputprojs_) {
0206     countinputproj += inputproj->nTracklets();
0207   }
0208 
0209   for (auto& matchengine : matchengines_) {
0210     matchengine.reset();
0211   }
0212 
0213   ProjectionTemp tmpProj_, tmpProj__;
0214   bool good_ = false;
0215   bool good__ = false;
0216 
0217   for (unsigned int istep = 0; istep < settings_.maxStep("MP"); istep++) {
0218     // This print statement is useful for detailed comparison with the HLS code
0219     // It prints out detailed status information for each clock step
0220     /* 
0221     if (print) {
0222       cout << "istep = "<<istep<<" projBuff: "<<inputProjBuffer_.rptr()<<" "<<inputProjBuffer_.wptr()<<" "<<projBuffNearFull;
0223       unsigned int iMEU = 0;
0224       for (auto& matchengine : matchengines_) {
0225     cout <<" MEU"<<iMEU<<": "<<matchengine.rptr()<<" "<<matchengine.wptr()
0226          <<" "<<matchengine.idle()<<" "<<matchengine.empty()
0227          <<" "<<matchengine.TCID();
0228     iMEU++;
0229       }
0230       cout << std::endl;
0231     }
0232     */
0233 
0234     //First do some caching of state at the start of the clock
0235 
0236     bool projdone = false;
0237 
0238     bool projBuffNearFull = inputProjBuffer_.nearfull();
0239 
0240     for (unsigned int iME = 0; iME < nMatchEngines_; iME++) {
0241       matchengines_[iME].setAlmostFull();
0242     }
0243 
0244     //Step 3
0245     //Check if we have candidate match to process
0246 
0247     unsigned int iMEbest = 0;
0248     int bestTCID = matchengines_[0].TCID();
0249     bool meactive = matchengines_[0].active();
0250     for (unsigned int iME = 1; iME < nMatchEngines_; iME++) {
0251       meactive = meactive || matchengines_[iME].active();
0252       int tcid = matchengines_[iME].TCID();
0253       if (tcid < bestTCID) {
0254         bestTCID = tcid;
0255         iMEbest = iME;
0256       }
0257     }
0258 
0259     // check if the matche engine processing the smallest tcid has match
0260 
0261     if (!matchengines_[iMEbest].empty()) {
0262       std::pair<Tracklet*, const Stub*> candmatch = matchengines_[iMEbest].read();
0263 
0264       const Stub* fpgastub = candmatch.second;
0265       Tracklet* tracklet = candmatch.first;
0266 
0267       //Consistency check
0268       if (oldTracklet != nullptr) {
0269         //allow equal here since we can have more than one cadidate match per tracklet projection
0270         //cout << "old new : "<<oldTracklet->TCID()<<" "<<tracklet->TCID()<<" "<<iMEbest<<endl;
0271         assert(oldTracklet->TCID() <= tracklet->TCID());
0272       }
0273       oldTracklet = tracklet;
0274 
0275       bool match = matchCalculator(tracklet, fpgastub, print, istep);
0276 
0277       if (settings_.debugTracklet() && match) {
0278         edm::LogVerbatim("Tracklet") << getName() << " have match";
0279       }
0280 
0281       countall++;
0282       if (match)
0283         countsel++;
0284     }
0285 
0286     //Step 2
0287     //Check if we have ME that can process projection
0288 
0289     bool addedProjection = false;
0290     for (unsigned int iME = 0; iME < nMatchEngines_; iME++) {
0291       if (!matchengines_[iME].idle())
0292         countme++;
0293       //if match engine empty and we have queued projections add to match engine
0294       if ((!addedProjection) && matchengines_[iME].idle() && (!inputProjBuffer_.empty())) {
0295         ProjectionTemp tmpProj = inputProjBuffer_.read();
0296         VMStubsMEMemory* stubmem = vmstubs_[0];
0297 
0298         if (settings_.debugTracklet()) {
0299           edm::LogVerbatim("Tracklet") << getName() << " adding projection to match engine";
0300         }
0301 
0302         int nbins = (1 << N_RZBITS);
0303         if (layerdisk_ >= N_LAYER) {
0304           nbins *= 2;  //twice as many bins in disks (since there are two disks)
0305         }
0306 
0307         matchengines_[iME].init(stubmem,
0308                                 nbins,
0309                                 tmpProj.slot(),
0310                                 tmpProj.iphi(),
0311                                 tmpProj.shift(),
0312                                 tmpProj.projrinv(),
0313                                 tmpProj.projfinerz(),
0314                                 tmpProj.projfinephi(),
0315                                 tmpProj.use(0, 0),
0316                                 tmpProj.use(0, 1),
0317                                 tmpProj.use(1, 0),
0318                                 tmpProj.use(1, 1),
0319                                 tmpProj.isPSseed(),
0320                                 tmpProj.proj());
0321         meactive = true;
0322         addedProjection = true;
0323       } else {
0324         matchengines_[iME].step();
0325       }
0326       matchengines_[iME].processPipeline();
0327     }
0328 
0329     //Step 1
0330     //First step here checks if we have more input projections to put into
0331     //the input puffer for projections
0332 
0333     if (good__) {
0334       inputProjBuffer_.store(tmpProj__);
0335     }
0336 
0337     good__ = good_;
0338     tmpProj__ = tmpProj_;
0339 
0340     good_ = false;
0341 
0342     if (iprojmem < inputprojs_.size()) {
0343       TrackletProjectionsMemory* projMem = inputprojs_[iprojmem];
0344       if (!projBuffNearFull) {
0345         if (settings_.debugTracklet()) {
0346           edm::LogVerbatim("Tracklet") << getName() << " have projection in memory : " << projMem->getName();
0347         }
0348 
0349         Tracklet* proj = projMem->getTracklet(iproj);
0350 
0351         FPGAWord fpgaphi = proj->proj(layerdisk_).fpgaphiproj();
0352 
0353         unsigned int iphi = (fpgaphi.value() >> (fpgaphi.nbits() - nvmbits_)) & (nvmbins_ - 1);
0354 
0355         constexpr int nextrabits = 2;
0356         int overlapbits = nvmbits_ + nextrabits;
0357 
0358         unsigned int extrabits = fpgaphi.bits(fpgaphi.nbits() - overlapbits - nextrabits, nextrabits);
0359 
0360         unsigned int ivmPlus = iphi;
0361 
0362         int shift = 0;
0363 
0364         if (extrabits == ((1U << nextrabits) - 1) && iphi != ((1U << settings_.nbitsvmme(layerdisk_)) - 1)) {
0365           shift = 1;
0366           ivmPlus++;
0367         }
0368         unsigned int ivmMinus = iphi;
0369         if (extrabits == 0 && iphi != 0) {
0370           shift = -1;
0371           ivmMinus--;
0372         }
0373 
0374         int projrinv = -1;
0375         if (barrel_) {
0376           FPGAWord phider = proj->proj(layerdisk_).fpgaphiprojder();
0377           projrinv = (1 << (nrinv_ - 1)) - 1 - (phider.value() >> (phider.nbits() - nrinv_));
0378         } else {
0379           //The next lines looks up the predicted bend based on:
0380           // 1 - r projections
0381           // 2 - phi derivative
0382           // 3 - the sign - i.e. if track is forward or backward
0383 
0384           int rindex =
0385               (proj->proj(layerdisk_).fpgarzproj().value() >> (proj->proj(layerdisk_).fpgarzproj().nbits() - nrbits_)) &
0386               ((1 << nrbits_) - 1);
0387 
0388           int phiprojder = proj->proj(layerdisk_).fpgaphiprojder().value();
0389 
0390           int phiderindex = (phiprojder >> (proj->proj(layerdisk_).fpgaphiprojder().nbits() - nphiderbits_)) &
0391                             ((1 << nphiderbits_) - 1);
0392 
0393           int signindex = proj->proj(layerdisk_).fpgarzprojder().value() < 0;
0394 
0395           int bendindex = (signindex << (nphiderbits_ + nrbits_)) + (rindex << (nphiderbits_)) + phiderindex;
0396 
0397           projrinv = rinvbendlut_.lookup(bendindex);
0398 
0399           proj->proj(layerdisk_).setBendIndex(projrinv);
0400         }
0401         assert(projrinv >= 0);
0402 
0403         unsigned int projfinephi =
0404             (fpgaphi.value() >> (fpgaphi.nbits() - (nvmbits_ + NFINEPHIBITS))) & ((1 << NFINEPHIBITS) - 1);
0405 
0406         unsigned int slot;
0407         bool second;
0408         int projfinerz;
0409 
0410         if (barrel_) {
0411           slot = proj->proj(layerdisk_).fpgarzbin1projvm().value();
0412           second = proj->proj(layerdisk_).fpgarzbin2projvm().value();
0413           projfinerz = proj->proj(layerdisk_).fpgafinerzvm().value();
0414         } else {
0415           //The -1 here is due to not using the full range of bits. Should be fixed.
0416           unsigned int ir = proj->proj(layerdisk_).fpgarzproj().value() >>
0417                             (proj->proj(layerdisk_).fpgarzproj().nbits() - nrprojbits_ - 1);
0418           unsigned int word = diskRadius_.lookup(ir);
0419 
0420           slot = (word >> 1) & ((1 << N_RZBITS) - 1);
0421           if (proj->proj(layerdisk_).fpgarzprojder().value() < 0) {
0422             slot += (1 << N_RZBITS);
0423           }
0424           second = word & 1;
0425           projfinerz = word >> 4;
0426         }
0427 
0428         bool isPSseed = proj->PSseed();
0429 
0430         int nbins = (1 << N_RZBITS);
0431         if (layerdisk_ >= N_LAYER) {
0432           nbins *= 2;  //twice as many bins in disks (since there are two disks)
0433         }
0434 
0435         VMStubsMEMemory* stubmem = vmstubs_[0];
0436         bool usefirstMinus = stubmem->nStubsBin(ivmMinus * nbins + slot) != 0;
0437         bool usesecondMinus = (second && (stubmem->nStubsBin(ivmMinus * nbins + slot + 1) != 0));
0438         bool usefirstPlus = ivmPlus != ivmMinus && (stubmem->nStubsBin(ivmPlus * nbins + slot) != 0);
0439         bool usesecondPlus = ivmPlus != ivmMinus && (second && (stubmem->nStubsBin(ivmPlus * nbins + slot + 1) != 0));
0440 
0441         good_ = usefirstPlus || usesecondPlus || usefirstMinus || usesecondMinus;
0442 
0443         if (good_) {
0444           ProjectionTemp tmpProj(proj,
0445                                  slot,
0446                                  projrinv,
0447                                  projfinerz,
0448                                  projfinephi,
0449                                  ivmMinus,
0450                                  shift,
0451                                  usefirstMinus,
0452                                  usefirstPlus,
0453                                  usesecondMinus,
0454                                  usesecondPlus,
0455                                  isPSseed);
0456           tmpProj_ = tmpProj;
0457         }
0458 
0459         iproj++;
0460         if (iproj == projMem->nTracklets()) {
0461           iproj = 0;
0462           do {
0463             iprojmem++;
0464           } while (iprojmem < inputprojs_.size() && inputprojs_[iprojmem]->nTracklets() == 0);
0465         }
0466 
0467       } else {
0468         projdone = true && !good_ && !good__;
0469       }
0470     }
0471 
0472     //
0473     //  Check if done
0474     //
0475     //
0476     //
0477 
0478     if ((projdone && !meactive && inputProjBuffer_.rptr() == inputProjBuffer_.wptr()) ||
0479         (istep == settings_.maxStep("MP") - 1)) {
0480       if (settings_.writeMonitorData("MP")) {
0481         globals_->ofstream("matchprocessor.txt") << getName() << " " << istep << " " << countall << " " << countsel
0482                                                  << " " << countme << " " << countinputproj << endl;
0483       }
0484       break;
0485     }
0486   }
0487 
0488   if (settings_.writeMonitorData("MC")) {
0489     globals_->ofstream("matchcalculator.txt") << getName() << " " << countall << " " << countsel << endl;
0490   }
0491 }
0492 
0493 bool MatchProcessor::matchCalculator(Tracklet* tracklet, const Stub* fpgastub, bool, unsigned int istep) {
0494   const L1TStub* stub = fpgastub->l1tstub();
0495 
0496   if (layerdisk_ < N_LAYER) {
0497     const Projection& proj = tracklet->proj(layerdisk_);
0498     int ir = fpgastub->r().value();
0499     int iphi = proj.fpgaphiproj().value();
0500     int icorr = (ir * proj.fpgaphiprojder().value()) >> icorrshift_;
0501     iphi += icorr;
0502 
0503     int iz = proj.fpgarzproj().value();
0504     int izcor = (ir * proj.fpgarzprojder().value() + (1 << (icorzshift_ - 1))) >> icorzshift_;
0505     iz += izcor;
0506 
0507     int ideltaz = fpgastub->z().value() - iz;
0508     int ideltaphi = (fpgastub->phi().value() - iphi) << phishift_;
0509 
0510     //Floating point calculations
0511 
0512     double phi = stub->phi();
0513     double r = stub->r();
0514     double z = stub->z();
0515 
0516     if (settings_.useapprox()) {
0517       double dphi = reco::reduceRange(phi - fpgastub->phiapprox(phimin_, 0.0));
0518       assert(std::abs(dphi) < 0.001);
0519       phi = fpgastub->phiapprox(phimin_, 0.0);
0520       z = fpgastub->zapprox();
0521       r = fpgastub->rapprox();
0522     }
0523 
0524     if (phi < 0)
0525       phi += 2 * M_PI;
0526     phi -= phimin_;
0527 
0528     double dr = r - settings_.rmean(layerdisk_);
0529     assert(std::abs(dr) < settings_.drmax());
0530 
0531     double dphi = reco::reduceRange(phi - (proj.phiproj() + dr * proj.phiprojder()));
0532 
0533     double dz = z - (proj.rzproj() + dr * proj.rzprojder());
0534 
0535     double dphiapprox = reco::reduceRange(phi - (proj.phiprojapprox() + dr * proj.phiprojderapprox()));
0536 
0537     double dzapprox = z - (proj.rzprojapprox() + dr * proj.rzprojderapprox());
0538 
0539     int seedindex = tracklet->getISeed();
0540     curr_tracklet = next_tracklet;
0541     next_tracklet = tracklet;
0542 
0543     // Do we have a new tracklet?
0544     bool newtracklet = (istep == 0 || tracklet != curr_tracklet);
0545     if (istep == 0)
0546       best_ideltar_disk = (1 << (fpgastub->r().nbits() - 1));  // Set to the maximum possible
0547     // If so, replace the "best" values with the cut tables
0548     if (newtracklet) {
0549       best_ideltaphi_barrel = (int)phimatchcuttable_.lookup(seedindex);
0550       best_ideltaz_barrel = (int)zmatchcuttable_.lookup(seedindex);
0551     }
0552 
0553     assert(phimatchcuttable_.lookup(seedindex) > 0);
0554     assert(zmatchcuttable_.lookup(seedindex) > 0);
0555 
0556     if (settings_.bookHistos()) {
0557       bool truthmatch = tracklet->stubtruthmatch(stub);
0558 
0559       HistBase* hists = globals_->histograms();
0560       hists->FillLayerResidual(layerdisk_ + 1,
0561                                seedindex,
0562                                dphiapprox * settings_.rmean(layerdisk_),
0563                                ideltaphi * settings_.kphi1() * settings_.rmean(layerdisk_),
0564                                (ideltaz << dzshift_) * settings_.kz(),
0565                                dz,
0566                                truthmatch);
0567     }
0568 
0569     if (settings_.writeMonitorData("Residuals")) {
0570       double pt = 0.01 * settings_.c() * settings_.bfield() / std::abs(tracklet->rinv());
0571 
0572       globals_->ofstream("layerresiduals.txt")
0573           << layerdisk_ + 1 << " " << seedindex << " " << pt << " "
0574           << ideltaphi * settings_.kphi1() * settings_.rmean(layerdisk_) << " "
0575           << dphiapprox * settings_.rmean(layerdisk_) << " "
0576           << phimatchcuttable_.lookup(seedindex) * settings_.kphi1() * settings_.rmean(layerdisk_) << "   "
0577           << (ideltaz << dzshift_) * settings_.kz() << " " << dz << " "
0578           << zmatchcuttable_.lookup(seedindex) * settings_.kz() << endl;
0579     }
0580 
0581     bool imatch = (std::abs(ideltaphi) <= best_ideltaphi_barrel && (ideltaz << dzshift_ < best_ideltaz_barrel) &&
0582                    (ideltaz << dzshift_ >= -best_ideltaz_barrel));
0583     // Update the "best" values
0584     if (imatch) {
0585       best_ideltaphi_barrel = std::abs(ideltaphi);
0586       best_ideltaz_barrel = std::abs(ideltaz << dzshift_);
0587     }
0588 
0589     if (settings_.debugTracklet()) {
0590       edm::LogVerbatim("Tracklet") << getName() << " imatch = " << imatch << " ideltaphi cut " << ideltaphi << " "
0591                                    << phimatchcuttable_.lookup(seedindex) << " ideltaz<<dzshift cut "
0592                                    << (ideltaz << dzshift_) << " " << zmatchcuttable_.lookup(seedindex);
0593     }
0594 
0595     //This would catch significant consistency problems in the configuration - helps to debug if there are problems.
0596     if (std::abs(dphi) > 0.5 * settings_.dphisectorHG() || std::abs(dphiapprox) > 0.5 * settings_.dphisectorHG()) {
0597       throw cms::Exception("LogicError") << "WARNING dphi and/or dphiapprox too large : " << dphi << " " << dphiapprox
0598                                          << endl;
0599     }
0600 
0601     bool keep = true;
0602     if (!settings_.doKF() || !settings_.doMultipleMatches()) {
0603       // Case of allowing only one stub per track per layer (or no KF which implies the same).
0604       if (imatch && tracklet->match(layerdisk_)) {
0605         // Veto match if is not the best one for this tracklet (in given layer)
0606         auto res = tracklet->resid(layerdisk_);
0607         keep = abs(ideltaphi) < abs(res.fpgaphiresid().value());
0608         imatch = keep;
0609       }
0610     }
0611 
0612     if (imatch) {
0613       tracklet->addMatch(layerdisk_,
0614                          ideltaphi,
0615                          ideltaz,
0616                          dphi,
0617                          dz,
0618                          dphiapprox,
0619                          dzapprox,
0620                          (phiregion_ << N_BITSMEMADDRESS) + fpgastub->stubindex().value(),
0621                          fpgastub);
0622 
0623       if (settings_.debugTracklet()) {
0624         edm::LogVerbatim("Tracklet") << "Accepted full match in layer " << getName() << " " << tracklet;
0625       }
0626 
0627       int iSeed = tracklet->getISeed();
0628       assert(fullmatches_[iSeed] != nullptr);
0629       fullmatches_[iSeed]->addMatch(tracklet, fpgastub);
0630 
0631       return true;
0632     } else {
0633       return false;
0634     }
0635   } else {  //disk matches
0636 
0637     //check that stubs and projections in same half of detector
0638     assert(stub->z() * tracklet->t() > 0.0);
0639 
0640     int sign = (tracklet->t() > 0.0) ? 1 : -1;
0641     int disk = sign * (layerdisk_ - N_LAYER + 1);
0642     assert(disk != 0);
0643 
0644     //Perform integer calculations here
0645 
0646     int iz = fpgastub->z().value();
0647 
0648     const Projection& proj = tracklet->proj(layerdisk_);
0649 
0650     int iphi = proj.fpgaphiproj().value();
0651     int iphicorr = (iz * proj.fpgaphiprojder().value()) >> icorrshift_;
0652 
0653     iphi += iphicorr;
0654 
0655     int ir = proj.fpgarzproj().value();
0656     int ircorr = (iz * proj.fpgarzprojder().value()) >> icorzshift_;
0657     ir += ircorr;
0658 
0659     int ideltaphi = fpgastub->phi().value() - iphi;
0660 
0661     int irstub = fpgastub->r().value();
0662     int ialphafact = 0;
0663     if (!stub->isPSmodule()) {
0664       assert(irstub < (int)N_DSS_MOD * 2);
0665       if (layerdisk_ - N_LAYER <= 1) {
0666         ialphafact = ialphafactinner_[irstub];
0667         irstub = settings_.rDSSinner(irstub) / settings_.kr();
0668       } else {
0669         ialphafact = ialphafactouter_[irstub];
0670         irstub = settings_.rDSSouter(irstub) / settings_.kr();
0671       }
0672     }
0673 
0674     constexpr int diff_bits = 1;
0675     int ideltar = (irstub >> diff_bits) - ir;
0676 
0677     if (!stub->isPSmodule()) {
0678       int ialpha = fpgastub->alpha().value();
0679       int iphialphacor = ((ideltar * ialpha * ialphafact) >> settings_.alphashift());
0680       ideltaphi += iphialphacor;
0681     }
0682 
0683     //Perform floating point calculations here
0684 
0685     double phi = stub->phi();
0686     double z = stub->z();
0687     double r = stub->r();
0688 
0689     if (settings_.useapprox()) {
0690       double dphi = reco::reduceRange(phi - fpgastub->phiapprox(phimin_, 0.0));
0691       assert(std::abs(dphi) < 0.001);
0692       phi = fpgastub->phiapprox(phimin_, 0.0);
0693       z = fpgastub->zapprox();
0694       r = fpgastub->rapprox();
0695     }
0696 
0697     if (phi < 0)
0698       phi += 2 * M_PI;
0699     phi -= phimin_;
0700 
0701     double dz = z - sign * settings_.zmean(layerdisk_ - N_LAYER);
0702 
0703     if (std::abs(dz) > settings_.dzmax()) {
0704       edm::LogProblem("Tracklet") << __FILE__ << ":" << __LINE__ << " " << name_ << " " << tracklet->getISeed();
0705       edm::LogProblem("Tracklet") << "stub " << stub->z() << " disk " << disk << " " << dz;
0706       assert(std::abs(dz) < settings_.dzmax());
0707     }
0708 
0709     double phiproj = proj.phiproj() + dz * proj.phiprojder();
0710     double rproj = proj.rzproj() + dz * proj.rzprojder();
0711     double deltar = r - rproj;
0712 
0713     double dr = stub->r() - rproj;
0714     double drapprox = stub->r() - (proj.rzprojapprox() + dz * proj.rzprojderapprox());
0715 
0716     double dphi = reco::reduceRange(phi - phiproj);
0717 
0718     double dphiapprox = reco::reduceRange(phi - (proj.phiprojapprox() + dz * proj.phiprojderapprox()));
0719 
0720     double drphi = dphi * stub->r();
0721     double drphiapprox = dphiapprox * stub->r();
0722 
0723     if (!stub->isPSmodule()) {
0724       double alphanorm = stub->alphanorm();
0725       dphi += dr * alphanorm * settings_.half2SmoduleWidth() / stub->r2();
0726       ;
0727       dphiapprox += drapprox * alphanorm * settings_.half2SmoduleWidth() / stub->r2();
0728 
0729       drphi += dr * alphanorm * settings_.half2SmoduleWidth() / stub->r();
0730       drphiapprox += dr * alphanorm * settings_.half2SmoduleWidth() / stub->r();
0731     }
0732 
0733     int seedindex = tracklet->getISeed();
0734 
0735     int idrphicut = rphicutPStable_.lookup(seedindex);
0736     int idrcut = rcutPStable_.lookup(seedindex);
0737     if (!stub->isPSmodule()) {
0738       idrphicut = rphicut2Stable_.lookup(seedindex);
0739       idrcut = rcut2Stable_.lookup(seedindex);
0740     }
0741 
0742     curr_tracklet = next_tracklet;
0743     next_tracklet = tracklet;
0744     // Do we have a new tracklet?
0745     bool newtracklet = (istep == 0 || tracklet != curr_tracklet);
0746     // If so, replace the "best" values with the cut tables
0747     if (newtracklet) {
0748       best_ideltaphi_disk = idrphicut;
0749       best_ideltar_disk = idrcut;
0750     }
0751 
0752     double drphicut = idrphicut * settings_.kphi() * settings_.kr();
0753     double drcut = idrcut * settings_.krprojshiftdisk();
0754 
0755     if (settings_.writeMonitorData("Residuals")) {
0756       double pt = 0.01 * settings_.c() * settings_.bfield() / std::abs(tracklet->rinv());
0757 
0758       globals_->ofstream("diskresiduals.txt")
0759           << layerdisk_ - N_LAYER + 1 << " " << stub->isPSmodule() << " " << tracklet->layer() << " "
0760           << abs(tracklet->disk()) << " " << pt << " " << ideltaphi * settings_.kphi() * stub->r() << " " << drphiapprox
0761           << " " << drphicut << " " << ideltar * settings_.krprojshiftdisk() << " " << deltar << " " << drcut << " "
0762           << endl;
0763     }
0764 
0765     bool match = (std::abs(drphi) < drphicut) && (std::abs(deltar) < drcut);
0766     bool imatch = (std::abs(ideltaphi * irstub) < best_ideltaphi_disk) && (std::abs(ideltar) < best_ideltar_disk);
0767     // Update the "best" values
0768     if (imatch) {
0769       best_ideltaphi_disk = std::abs(ideltaphi) * irstub;
0770       best_ideltar_disk = std::abs(ideltar);
0771     }
0772 
0773     if (settings_.debugTracklet()) {
0774       edm::LogVerbatim("Tracklet") << "imatch match disk: " << imatch << " " << match << " " << std::abs(ideltaphi)
0775                                    << " " << drphicut / (settings_.kphi() * stub->r()) << " " << std::abs(ideltar)
0776                                    << " " << drcut / settings_.krprojshiftdisk() << " r = " << stub->r();
0777     }
0778 
0779     bool keep = true;
0780     if (!settings_.doKF() || !settings_.doMultipleMatches()) {
0781       // Case of allowing only one stub per track per layer (or no KF which implies the same).
0782       if (imatch && tracklet->match(layerdisk_)) {
0783         // Veto match if is not the best one for this tracklet (in given layer)
0784         auto res = tracklet->resid(layerdisk_);
0785         keep = abs(ideltaphi) < abs(res.fpgaphiresid().value());
0786         imatch = keep;
0787       }
0788     }
0789 
0790     if (imatch) {
0791       if (settings_.debugTracklet()) {
0792         edm::LogVerbatim("Tracklet") << "MatchCalculator found match in disk " << getName();
0793       }
0794 
0795       if (std::abs(dphi) >= third * settings_.dphisectorHG()) {
0796         edm::LogPrint("Tracklet") << "dphi " << dphi << " ISeed " << tracklet->getISeed();
0797       }
0798       assert(std::abs(dphi) < third * settings_.dphisectorHG());
0799       assert(std::abs(dphiapprox) < third * settings_.dphisectorHG());
0800 
0801       tracklet->addMatch(layerdisk_,
0802                          ideltaphi,
0803                          ideltar,
0804                          drphi / stub->r(),
0805                          dr,
0806                          drphiapprox / stub->r(),
0807                          drapprox,
0808                          (phiregion_ << N_BITSMEMADDRESS) + fpgastub->stubindex().value(),
0809                          fpgastub);
0810 
0811       if (settings_.debugTracklet()) {
0812         edm::LogVerbatim("Tracklet") << "Accepted full match in disk " << getName() << " " << tracklet;
0813       }
0814 
0815       int iSeed = tracklet->getISeed();
0816       assert(fullmatches_[iSeed] != nullptr);
0817       fullmatches_[iSeed]->addMatch(tracklet, fpgastub);
0818 
0819       return true;
0820     } else {
0821       return false;
0822     }
0823   }
0824 }