Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-10-14 01:44:06 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!