TrackFindingTracklet/src/MatchProcessor.cc

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