Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:22:00

0001 //////////////////////////////////////////////////////////////////
0002 // MatchCalculator
0003 // This class loads pairs for tracklet/stubs and looks for the
0004 // best possible match.
0005 // Variables such as `best_ideltaphi_barrel` store the "global"
0006 // best value for delta phi, r, z, and r*phi, for instances
0007 // where the same tracklet has multiple stub pairs. This allows
0008 // us to find the truly best match
0009 //////////////////////////////////////////////////////////////////
0010 
0011 #include "L1Trigger/TrackFindingTracklet/interface/MatchCalculator.h"
0012 #include "L1Trigger/TrackFindingTracklet/interface/Globals.h"
0013 #include "L1Trigger/TrackFindingTracklet/interface/Util.h"
0014 #include "L1Trigger/TrackFindingTracklet/interface/CandidateMatchMemory.h"
0015 #include "L1Trigger/TrackFindingTracklet/interface/FullMatchMemory.h"
0016 #include "L1Trigger/TrackFindingTracklet/interface/AllStubsMemory.h"
0017 #include "L1Trigger/TrackFindingTracklet/interface/AllProjectionsMemory.h"
0018 #include "L1Trigger/TrackFindingTracklet/interface/Tracklet.h"
0019 #include "L1Trigger/TrackFindingTracklet/interface/Stub.h"
0020 #include "L1Trigger/TrackFindingTracklet/interface/HistBase.h"
0021 
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023 #include "FWCore/Utilities/interface/Exception.h"
0024 #include "DataFormats/Math/interface/deltaPhi.h"
0025 
0026 #include <filesystem>
0027 #include <algorithm>
0028 
0029 using namespace std;
0030 using namespace trklet;
0031 
0032 MatchCalculator::MatchCalculator(string name, Settings const& settings, Globals* global)
0033     : ProcessBase(name, settings, global),
0034       phimatchcuttable_(settings),
0035       zmatchcuttable_(settings),
0036       rphicutPStable_(settings),
0037       rphicut2Stable_(settings),
0038       rcutPStable_(settings),
0039       rcut2Stable_(settings),
0040       alphainner_(settings),
0041       alphaouter_(settings),
0042       rSSinner_(settings),
0043       rSSouter_(settings) {
0044   phiregion_ = name[8] - 'A';
0045   layerdisk_ = initLayerDisk(3);
0046 
0047   fullMatches_.resize(12, nullptr);
0048 
0049   //TODO - need to sort out constants here
0050   icorrshift_ = 7;
0051 
0052   if (layerdisk_ < N_PSLAYER) {
0053     icorzshift_ = -1 - settings_.PS_zderL_shift();
0054   } else {
0055     icorzshift_ = -1 - settings_.SS_zderL_shift();
0056   }
0057   phi0shift_ = 3;
0058   fact_ = 1;
0059   if (layerdisk_ >= N_PSLAYER && layerdisk_ < N_LAYER) {
0060     fact_ = (1 << (settings_.nzbitsstub(0) - settings_.nzbitsstub(5)));
0061     icorrshift_ -= (10 - settings_.nrbitsstub(layerdisk_));
0062     icorzshift_ += (settings_.nzbitsstub(0) - settings_.nzbitsstub(5) + settings_.nrbitsstub(layerdisk_) -
0063                     settings_.nrbitsstub(0));
0064     phi0shift_ = 0;
0065   }
0066 
0067   unsigned int region = getName()[8] - 'A';
0068   assert(region < settings_.nallstubs(layerdisk_));
0069 
0070   if (layerdisk_ < N_LAYER) {
0071     phimatchcuttable_.initmatchcut(layerdisk_, TrackletLUT::MatchType::barrelphi, region);
0072     zmatchcuttable_.initmatchcut(layerdisk_, TrackletLUT::MatchType::barrelz, region);
0073   } else {
0074     rphicutPStable_.initmatchcut(layerdisk_, TrackletLUT::MatchType::diskPSphi, region);
0075     rphicut2Stable_.initmatchcut(layerdisk_, TrackletLUT::MatchType::disk2Sphi, region);
0076     rcutPStable_.initmatchcut(layerdisk_, TrackletLUT::MatchType::diskPSr, region);
0077     rcut2Stable_.initmatchcut(layerdisk_, TrackletLUT::MatchType::disk2Sr, region);
0078     alphainner_.initmatchcut(layerdisk_, TrackletLUT::MatchType::alphainner, region);
0079     alphaouter_.initmatchcut(layerdisk_, TrackletLUT::MatchType::alphaouter, region);
0080     rSSinner_.initmatchcut(layerdisk_, TrackletLUT::MatchType::rSSinner, region);
0081     rSSouter_.initmatchcut(layerdisk_, TrackletLUT::MatchType::rSSouter, region);
0082   }
0083 
0084   for (unsigned int i = 0; i < N_DSS_MOD * 2; i++) {
0085     ialphafactinner_[i] = (1 << settings_.alphashift()) * settings_.krprojshiftdisk() * settings_.half2SmoduleWidth() /
0086                           (1 << (settings_.nbitsalpha() - 1)) / (settings_.rDSSinner(i) * settings_.rDSSinner(i)) /
0087                           settings_.kphi();
0088     ialphafactouter_[i] = (1 << settings_.alphashift()) * settings_.krprojshiftdisk() * settings_.half2SmoduleWidth() /
0089                           (1 << (settings_.nbitsalpha() - 1)) / (settings_.rDSSouter(i) * settings_.rDSSouter(i)) /
0090                           settings_.kphi();
0091   }
0092 }
0093 
0094 void MatchCalculator::addOutput(MemoryBase* memory, string output) {
0095   if (settings_.writetrace()) {
0096     edm::LogVerbatim("Tracklet") << "In " << name_ << " adding output to " << memory->getName() << " to output "
0097                                  << output;
0098   }
0099   if (output.substr(0, 8) == "matchout") {
0100     auto* tmp = dynamic_cast<FullMatchMemory*>(memory);
0101     assert(tmp != nullptr);
0102     unsigned int iSeed = getISeed(memory->getName());
0103     fullMatches_[iSeed] = tmp;
0104     return;
0105   }
0106   throw cms::Exception("BadConfig") << __FILE__ << " " << __LINE__ << " could not find output " << output;
0107 }
0108 
0109 void MatchCalculator::addInput(MemoryBase* memory, string input) {
0110   if (settings_.writetrace()) {
0111     edm::LogVerbatim("Tracklet") << "In " << name_ << " adding input from " << memory->getName() << " to input "
0112                                  << input;
0113   }
0114   if (input == "allstubin") {
0115     auto* tmp = dynamic_cast<AllStubsMemory*>(memory);
0116     assert(tmp != nullptr);
0117     allstubs_ = tmp;
0118     return;
0119   }
0120   if (input == "allprojin") {
0121     auto* tmp = dynamic_cast<AllProjectionsMemory*>(memory);
0122     assert(tmp != nullptr);
0123     allprojs_ = tmp;
0124     return;
0125   }
0126   if (input.substr(0, 5) == "match" && input.substr(input.size() - 2, 2) == "in") {
0127     auto* tmp = dynamic_cast<CandidateMatchMemory*>(memory);
0128     assert(tmp != nullptr);
0129     matches_.push_back(tmp);
0130     return;
0131   }
0132   throw cms::Exception("BadConfig") << __FILE__ << " " << __LINE__ << " could not find input " << input;
0133 }
0134 
0135 void MatchCalculator::execute(unsigned int iSector, double phioffset) {
0136   unsigned int countall = 0;
0137   unsigned int countsel = 0;
0138 
0139   //bool print = getName() == "MC_L4PHIC" && iSector == 3;
0140 
0141   Tracklet* oldTracklet = nullptr;
0142 
0143   // Get all tracklet/stub pairs
0144   std::vector<std::pair<std::pair<Tracklet*, int>, const Stub*> > mergedMatches = mergeMatches(matches_);
0145 
0146   // Number of clock cycles the pipeline in HLS takes to process the projection merging to
0147   // produce the first projection
0148   unsigned int mergedepth = 3;
0149 
0150   unsigned int maxProc = std::min(settings_.maxStep("MC") - mergedepth, (unsigned int)mergedMatches.size());
0151 
0152   // Pick some initial large values
0153   int best_ideltaphi_barrel = 0xFFFF;
0154   int best_ideltaz_barrel = 0xFFFF;
0155   int best_ideltaphi_disk = 0xFFFF;
0156   int best_ideltar_disk = 0xFFFF;
0157   unsigned int curr_projid = -1;
0158   unsigned int next_projid = -1;
0159 
0160   for (unsigned int j = 0; j < maxProc; j++) {
0161     if (settings_.debugTracklet() && j == 0) {
0162       edm::LogVerbatim("Tracklet") << getName() << " has " << mergedMatches.size() << " candidate matches";
0163     }
0164 
0165     countall++;
0166 
0167     const Stub* fpgastub = mergedMatches[j].second;
0168     Tracklet* tracklet = mergedMatches[j].first.first;
0169     const L1TStub* stub = fpgastub->l1tstub();
0170 
0171     //check that the matches are ordered correctly
0172     //allow equal here since we can have more than one candidate match per tracklet projection
0173     if (oldTracklet != nullptr) {
0174       assert(oldTracklet->TCID() <= tracklet->TCID());
0175     }
0176     oldTracklet = tracklet;
0177 
0178     if (layerdisk_ < N_LAYER) {
0179       //Integer calculation
0180 
0181       const Projection& proj = tracklet->proj(layerdisk_);
0182 
0183       int ir = fpgastub->r().value();
0184       int iphi = proj.fpgaphiproj().value();
0185       int icorr = (ir * proj.fpgaphiprojder().value()) >> icorrshift_;
0186       iphi += icorr;
0187 
0188       int iz = proj.fpgarzproj().value();
0189       int izcor = (ir * proj.fpgarzprojder().value() + (1 << (icorzshift_ - 1))) >> icorzshift_;
0190       iz += izcor;
0191 
0192       int ideltaz = fpgastub->z().value() - iz;
0193       int ideltaphi = (fpgastub->phi().value() << phi0shift_) - (iphi << (settings_.phi0bitshift() - 1 + phi0shift_));
0194 
0195       //Floating point calculations
0196 
0197       double phi = stub->phi() - phioffset;
0198       double r = stub->r();
0199       double z = stub->z();
0200 
0201       if (settings_.useapprox()) {
0202         double dphi = reco::reduceRange(phi - fpgastub->phiapprox(0.0, 0.0));
0203         assert(std::abs(dphi) < 0.001);
0204         phi = fpgastub->phiapprox(0.0, 0.0);
0205         z = fpgastub->zapprox();
0206         r = fpgastub->rapprox();
0207       }
0208 
0209       if (phi < 0)
0210         phi += 2 * M_PI;
0211 
0212       double dr = r - settings_.rmean(layerdisk_);
0213       assert(std::abs(dr) < settings_.drmax());
0214 
0215       double dphi = reco::reduceRange(phi - (proj.phiproj() + dr * proj.phiprojder()));
0216 
0217       double dz = z - (proj.rzproj() + dr * proj.rzprojder());
0218 
0219       double dphiapprox = reco::reduceRange(phi - (proj.phiprojapprox() + dr * proj.phiprojderapprox()));
0220 
0221       double dzapprox = z - (proj.rzprojapprox() + dr * proj.rzprojderapprox());
0222 
0223       int seedindex = tracklet->getISeed();
0224       unsigned int projindex = mergedMatches[j].first.second;  // Allproj index
0225       curr_projid = next_projid;
0226       next_projid = projindex;
0227 
0228       // Do we have a new tracklet?
0229       bool newtracklet = (j == 0 || projindex != curr_projid);
0230       if (j == 0)
0231         best_ideltar_disk = (1 << (fpgastub->r().nbits() - 1));  // Set to the maximum possible
0232       // If so, replace the "best" values with the cut tables
0233       if (newtracklet) {
0234         best_ideltaphi_barrel = (int)phimatchcuttable_.lookup(seedindex);
0235         best_ideltaz_barrel = (int)zmatchcuttable_.lookup(seedindex);
0236       }
0237 
0238       assert(phimatchcuttable_.lookup(seedindex) > 0);
0239       assert(zmatchcuttable_.lookup(seedindex) > 0);
0240 
0241       if (settings_.bookHistos()) {
0242         bool truthmatch = tracklet->stubtruthmatch(stub);
0243 
0244         HistBase* hists = globals_->histograms();
0245         hists->FillLayerResidual(layerdisk_ + 1,
0246                                  seedindex,
0247                                  dphiapprox * settings_.rmean(layerdisk_),
0248                                  ideltaphi * settings_.kphi1() * settings_.rmean(layerdisk_),
0249                                  ideltaz * fact_ * settings_.kz(),
0250                                  dz,
0251                                  truthmatch);
0252       }
0253 
0254       //This would catch significant consistency problems in the configuration - helps to debug if there are problems.
0255       if (std::abs(dphi) > 0.5 * settings_.dphisectorHG() || std::abs(dphiapprox) > 0.5 * settings_.dphisectorHG()) {
0256         throw cms::Exception("LogicError")
0257             << "WARNING dphi and/or dphiapprox too large : " << dphi << " " << dphiapprox << endl;
0258       }
0259 
0260       if (settings_.writeMonitorData("Residuals")) {
0261         double pt = 0.01 * settings_.c() * settings_.bfield() / std::abs(tracklet->rinv());
0262 
0263         globals_->ofstream("layerresiduals.txt")
0264             << layerdisk_ + 1 << " " << seedindex << " " << pt << " "
0265             << ideltaphi * settings_.kphi1() * settings_.rmean(layerdisk_) << " "
0266             << dphiapprox * settings_.rmean(layerdisk_) << " "
0267             << phimatchcuttable_.lookup(seedindex) * settings_.kphi1() * settings_.rmean(layerdisk_) << "   "
0268             << ideltaz * fact_ * settings_.kz() << " " << dz << " "
0269             << zmatchcuttable_.lookup(seedindex) * settings_.kz() << endl;
0270       }
0271 
0272       // integer match
0273       bool imatch = (std::abs(ideltaphi) <= best_ideltaphi_barrel) && (ideltaz * fact_ < best_ideltaz_barrel) &&
0274                     (ideltaz * fact_ >= -best_ideltaz_barrel);
0275       // Update the "best" values
0276       if (imatch) {
0277         best_ideltaphi_barrel = std::abs(ideltaphi);
0278         best_ideltaz_barrel = std::abs(ideltaz * fact_);
0279       }
0280 
0281       if (settings_.debugTracklet()) {
0282         edm::LogVerbatim("Tracklet") << getName() << " imatch = " << imatch << " ideltaphi cut " << ideltaphi << " "
0283                                      << phimatchcuttable_.lookup(seedindex) << " ideltaz*fact cut " << ideltaz * fact_
0284                                      << " " << zmatchcuttable_.lookup(seedindex);
0285       }
0286 
0287       if (imatch) {
0288         countsel++;
0289 
0290         tracklet->addMatch(layerdisk_,
0291                            ideltaphi,
0292                            ideltaz,
0293                            dphi,
0294                            dz,
0295                            dphiapprox,
0296                            dzapprox,
0297                            (phiregion_ << 7) + fpgastub->stubindex().value(),
0298                            mergedMatches[j].second);
0299 
0300         if (settings_.debugTracklet()) {
0301           edm::LogVerbatim("Tracklet") << "Accepted full match in layer " << getName() << " " << tracklet;
0302         }
0303 
0304         fullMatches_[seedindex]->addMatch(tracklet, mergedMatches[j].second);
0305       }
0306     } else {  //disk matches
0307       //check that stubs and projections in same half of detector
0308       assert(stub->z() * tracklet->t() > 0.0);
0309 
0310       int sign = (tracklet->t() > 0.0) ? 1 : -1;
0311       int disk = sign * (layerdisk_ - (N_LAYER - 1));
0312       assert(disk != 0);
0313 
0314       //Perform integer calculations here
0315 
0316       const Projection& proj = tracklet->proj(layerdisk_);
0317 
0318       int iz = fpgastub->z().value();
0319       int iphi = proj.fpgaphiproj().value();
0320 
0321       //TODO - need to express in terms of constants
0322       int shifttmp = 6;
0323       int iphicorr = (iz * proj.fpgaphiprojder().value()) >> shifttmp;
0324 
0325       iphi += iphicorr;
0326 
0327       int ir = proj.fpgarzproj().value();
0328 
0329       //TODO - need to express in terms of constants
0330       int shifttmp2 = 7;
0331       int ircorr = (iz * proj.fpgarzprojder().value()) >> shifttmp2;
0332 
0333       ir += ircorr;
0334 
0335       int ideltaphi = fpgastub->phi().value() * settings_.kphi() / settings_.kphi() - iphi;
0336 
0337       int irstub = fpgastub->r().value();
0338       int ialphafact = 0;
0339       if (!stub->isPSmodule()) {
0340         assert(irstub < (int)N_DSS_MOD * 2);
0341         if (abs(disk) <= 2) {
0342           ialphafact = ialphafactinner_[irstub];
0343           irstub = settings_.rDSSinner(irstub) / settings_.kr();
0344         } else {
0345           ialphafact = ialphafactouter_[irstub];
0346           irstub = settings_.rDSSouter(irstub) / settings_.kr();
0347         }
0348       }
0349 
0350       //TODO stub and projection r should not use different # bits...
0351       int ideltar = (irstub >> 1) - ir;
0352 
0353       if (!stub->isPSmodule()) {
0354         int ialpha = fpgastub->alpha().value();
0355         int iphialphacor = ((ideltar * ialpha * ialphafact) >> settings_.alphashift());
0356         ideltaphi += iphialphacor;
0357       }
0358 
0359       //Perform floating point calculations here
0360 
0361       double phi = stub->phi() - phioffset;
0362       double z = stub->z();
0363       double r = stub->r();
0364 
0365       if (settings_.useapprox()) {
0366         double dphi = reco::reduceRange(phi - fpgastub->phiapprox(0.0, 0.0));
0367         assert(std::abs(dphi) < 0.001);
0368         phi = fpgastub->phiapprox(0.0, 0.0);
0369         z = fpgastub->zapprox();
0370         r = fpgastub->rapprox();
0371       }
0372 
0373       if (phi < 0)
0374         phi += 2 * M_PI;
0375 
0376       double dz = z - sign * settings_.zmean(layerdisk_ - N_LAYER);
0377 
0378       if (std::abs(dz) > settings_.dzmax()) {
0379         throw cms::Exception("LogicError") << __FILE__ << " " << __LINE__ << " " << name_ << " " << tracklet->getISeed()
0380                                            << "\n stub " << stub->z() << " disk " << disk << " " << dz;
0381       }
0382 
0383       double phiproj = proj.phiproj() + dz * proj.phiprojder();
0384 
0385       double rproj = proj.rzproj() + dz * proj.rzprojder();
0386 
0387       double deltar = r - rproj;
0388 
0389       double dr = stub->r() - rproj;
0390 
0391       double dphi = reco::reduceRange(phi - phiproj);
0392 
0393       double dphiapprox = reco::reduceRange(phi - (proj.phiprojapprox() + dz * proj.phiprojderapprox()));
0394 
0395       double drapprox = stub->r() - (proj.rzprojapprox() + dz * proj.rzprojderapprox());
0396 
0397       double drphi = dphi * stub->r();
0398       double drphiapprox = dphiapprox * stub->r();
0399 
0400       if (!stub->isPSmodule()) {
0401         double alphanorm = stub->alphanorm();
0402         dphi += dr * alphanorm * settings_.half2SmoduleWidth() / stub->r2();
0403         dphiapprox += drapprox * alphanorm * settings_.half2SmoduleWidth() / stub->r2();
0404 
0405         drphi += dr * alphanorm * settings_.half2SmoduleWidth() / stub->r();
0406         drphiapprox += dr * alphanorm * settings_.half2SmoduleWidth() / stub->r();
0407       }
0408 
0409       int seedindex = tracklet->getISeed();
0410 
0411       int idrphicut = rphicutPStable_.lookup(seedindex);
0412       int idrcut = rcutPStable_.lookup(seedindex);
0413       if (!stub->isPSmodule()) {
0414         idrphicut = rphicut2Stable_.lookup(seedindex);
0415         idrcut = rcut2Stable_.lookup(seedindex);
0416       }
0417 
0418       unsigned int projindex = mergedMatches[j].first.second;  // Allproj index
0419       curr_projid = next_projid;
0420       next_projid = projindex;
0421       // Do we have a new tracklet?
0422       bool newtracklet = (j == 0 || projindex != curr_projid);
0423       // If so, replace the "best" values with the cut tables
0424       if (newtracklet) {
0425         best_ideltaphi_disk = idrphicut;
0426         best_ideltar_disk = idrcut;
0427       }
0428 
0429       double drphicut = idrphicut * settings_.kphi() * settings_.kr();
0430       double drcut = idrcut * settings_.krprojshiftdisk();
0431 
0432       bool match, imatch;
0433       if (std::abs(dphi) < third * settings_.dphisectorHG() &&
0434           std::abs(dphiapprox) < third * settings_.dphisectorHG()) {  //1/3 of sector size to catch errors
0435         if (settings_.writeMonitorData("Residuals")) {
0436           double pt = 0.01 * settings_.c() * settings_.bfield() / std::abs(tracklet->rinv());
0437 
0438           globals_->ofstream("diskresiduals.txt")
0439               << disk << " " << stub->isPSmodule() << " " << tracklet->layer() << " " << abs(tracklet->disk()) << " "
0440               << pt << " " << ideltaphi * settings_.kphi() * stub->r() << " " << drphiapprox << " " << drphicut << " "
0441               << ideltar * settings_.krprojshiftdisk() << " " << deltar << " " << drcut << " " << endl;
0442         }
0443 
0444         // floating point match
0445         match = (std::abs(drphi) < drphicut) && (std::abs(deltar) < drcut);
0446 
0447         // integer match
0448         imatch = (std::abs(ideltaphi) * irstub < best_ideltaphi_disk) && (std::abs(ideltar) < best_ideltar_disk);
0449         // Update the "best" values
0450         if (imatch) {
0451           best_ideltaphi_disk = std::abs(ideltaphi) * irstub;
0452           best_ideltar_disk = std::abs(ideltar);
0453         }
0454       } else {
0455         edm::LogProblem("Tracklet") << "WARNING dphi and/or dphiapprox too large : " << dphi << " " << dphiapprox
0456                                     << "dphi " << dphi << " Seed / ISeed " << tracklet->getISeed() << endl;
0457         match = false;
0458         imatch = false;
0459       }
0460       if (settings_.debugTracklet()) {
0461         edm::LogVerbatim("Tracklet") << "imatch match disk: " << imatch << " " << match << " " << std::abs(ideltaphi)
0462                                      << " " << drphicut / (settings_.kphi() * stub->r()) << " " << std::abs(ideltar)
0463                                      << " " << drcut / settings_.krprojshiftdisk() << " r = " << stub->r();
0464       }
0465 
0466       if (imatch) {
0467         countsel++;
0468 
0469         if (settings_.debugTracklet()) {
0470           edm::LogVerbatim("Tracklet") << "MatchCalculator found match in disk " << getName();
0471         }
0472 
0473         tracklet->addMatch(layerdisk_,
0474                            ideltaphi,
0475                            ideltar,
0476                            drphi / stub->r(),
0477                            dr,
0478                            drphiapprox / stub->r(),
0479                            drapprox,
0480                            (phiregion_ << 7) + fpgastub->stubindex().value(),
0481                            fpgastub);
0482 
0483         if (settings_.debugTracklet()) {
0484           edm::LogVerbatim("Tracklet") << "Accepted full match in disk " << getName() << " " << tracklet;
0485         }
0486 
0487         fullMatches_[seedindex]->addMatch(tracklet, mergedMatches[j].second);
0488       }
0489     }
0490     if (countall >= settings_.maxStep("MC"))
0491       break;
0492   }
0493 
0494   if (settings_.writeMonitorData("MC")) {
0495     globals_->ofstream("matchcalculator.txt") << getName() << " " << countall << " " << countsel << endl;
0496   }
0497 }
0498 
0499 // Combines all tracklet/stub pairs into a vector
0500 std::vector<std::pair<std::pair<Tracklet*, int>, const Stub*> > MatchCalculator::mergeMatches(
0501     vector<CandidateMatchMemory*>& candmatch) {
0502   std::vector<std::pair<std::pair<Tracklet*, int>, const Stub*> > tmp;
0503 
0504   std::vector<unsigned int> indexArray;
0505   indexArray.reserve(candmatch.size());
0506   for (unsigned int i = 0; i < candmatch.size(); i++) {
0507     indexArray.push_back(0);
0508   }
0509 
0510   int bestIndex = -1;
0511   do {
0512     int bestSector = 100;
0513     int bestTCID = -1;
0514     bestIndex = -1;
0515     for (unsigned int i = 0; i < candmatch.size(); i++) {
0516       if (indexArray[i] >= candmatch[i]->nMatches()) {
0517         // skip as we were at the end
0518         continue;
0519       }
0520       int TCID = candmatch[i]->getMatch(indexArray[i]).first.first->TCID();
0521       int dSector = 0;
0522       if (dSector > 2)
0523         dSector -= N_SECTOR;
0524       if (dSector < -2)
0525         dSector += N_SECTOR;
0526       assert(abs(dSector) < 2);
0527       if (dSector == -1)
0528         dSector = 2;
0529       if (dSector < bestSector) {
0530         bestSector = dSector;
0531         bestTCID = TCID;
0532         bestIndex = i;
0533       }
0534       if (dSector == bestSector) {
0535         if (TCID < bestTCID || bestTCID < 0) {
0536           bestTCID = TCID;
0537           bestIndex = i;
0538         }
0539       }
0540     }
0541     if (bestIndex != -1) {
0542       tmp.push_back(candmatch[bestIndex]->getMatch(indexArray[bestIndex]));
0543       indexArray[bestIndex]++;
0544     }
0545   } while (bestIndex != -1);
0546 
0547   if (layerdisk_ < N_LAYER) {
0548     int lastTCID = -1;
0549     bool error = false;
0550 
0551     //Allow equal TCIDs since we can have multiple candidate matches
0552     for (unsigned int i = 1; i < tmp.size(); i++) {
0553       if (lastTCID > tmp[i].first.first->TCID()) {
0554         edm::LogProblem("Tracklet") << "Wrong TCID ordering for projections in " << getName() << " last " << lastTCID
0555                                     << " " << tmp[i].first.first->TCID();
0556         error = true;
0557       } else {
0558         lastTCID = tmp[i].first.first->TCID();
0559       }
0560     }
0561 
0562     if (error) {
0563       for (unsigned int i = 1; i < tmp.size(); i++) {
0564         edm::LogProblem("Tracklet") << "Wrong order for in " << getName() << " " << i << " " << tmp[i].first.first
0565                                     << " " << tmp[i].first.first->TCID();
0566       }
0567     }
0568   }
0569 
0570   for (unsigned int i = 0; i < tmp.size(); i++) {
0571     if (i > 0) {
0572       //This allows for equal TCIDs. This means that we can e.g. have a track seeded
0573       //in L1L2 that projects to both L3 and D4. The algorithm will pick up the first hit and
0574       //drop the second
0575 
0576       assert(tmp[i - 1].first.first->TCID() <= tmp[i].first.first->TCID());
0577     }
0578   }
0579 
0580   return tmp;
0581 }