Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-11-15 23:40:45

0001 #include "L1Trigger/TrackFindingTracklet/interface/FitTrack.h"
0002 #include "L1Trigger/TrackFindingTracklet/interface/Globals.h"
0003 #include "L1Trigger/TrackFindingTracklet/interface/TrackDerTable.h"
0004 #include "L1Trigger/TrackFindingTracklet/interface/HybridFit.h"
0005 #include "L1Trigger/TrackFindingTracklet/interface/Tracklet.h"
0006 #include "L1Trigger/TrackFindingTracklet/interface/Stub.h"
0007 #include "L1Trigger/TrackFindingTracklet/interface/StubStreamData.h"
0008 #include "DataFormats/L1TrackTrigger/interface/TTBV.h"
0009 
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "FWCore/Utilities/interface/Exception.h"
0012 
0013 using namespace std;
0014 using namespace trklet;
0015 
0016 FitTrack::FitTrack(string name, Settings const& settings, Globals* global)
0017     : ProcessBase(name, settings, global), trackfit_(nullptr) {}
0018 
0019 void FitTrack::addOutput(MemoryBase* memory, string output) {
0020   if (settings_.writetrace()) {
0021     edm::LogVerbatim("Tracklet") << "In " << name_ << " adding output to " << memory->getName() << " to output "
0022                                  << output;
0023   }
0024   if (output == "trackout") {
0025     TrackFitMemory* tmp = dynamic_cast<TrackFitMemory*>(memory);
0026     assert(tmp != nullptr);
0027     trackfit_ = tmp;
0028     return;
0029   }
0030 
0031   throw cms::Exception("BadConfig") << __FILE__ << " " << __LINE__ << " addOutput, output = " << output << " not known";
0032 }
0033 
0034 void FitTrack::addInput(MemoryBase* memory, string input) {
0035   if (settings_.writetrace()) {
0036     edm::LogVerbatim("Tracklet") << "In " << name_ << " adding input from " << memory->getName() << " to input "
0037                                  << input;
0038   }
0039   if (input.substr(0, 4) == "tpar") {
0040     auto* tmp = dynamic_cast<TrackletParametersMemory*>(memory);
0041     assert(tmp != nullptr);
0042     seedtracklet_.push_back(tmp);
0043     return;
0044   }
0045   if (input.substr(0, 10) == "fullmatch1") {
0046     auto* tmp = dynamic_cast<FullMatchMemory*>(memory);
0047     assert(tmp != nullptr);
0048     fullmatch1_.push_back(tmp);
0049     return;
0050   }
0051   if (input.substr(0, 10) == "fullmatch2") {
0052     auto* tmp = dynamic_cast<FullMatchMemory*>(memory);
0053     assert(tmp != nullptr);
0054     fullmatch2_.push_back(tmp);
0055     return;
0056   }
0057   if (input.substr(0, 10) == "fullmatch3") {
0058     auto* tmp = dynamic_cast<FullMatchMemory*>(memory);
0059     assert(tmp != nullptr);
0060     fullmatch3_.push_back(tmp);
0061     return;
0062   }
0063   if (input.substr(0, 10) == "fullmatch4") {
0064     auto* tmp = dynamic_cast<FullMatchMemory*>(memory);
0065     assert(tmp != nullptr);
0066     fullmatch4_.push_back(tmp);
0067     return;
0068   }
0069 
0070   throw cms::Exception("BadConfig") << __FILE__ << " " << __LINE__ << " input = " << input << " not found";
0071 }
0072 
0073 #ifdef USEHYBRID
0074 void FitTrack::trackFitKF(Tracklet* tracklet,
0075                           std::vector<const Stub*>& trackstublist,
0076                           std::vector<std::pair<int, int>>& stubidslist) {
0077   if (settings_.doKF()) {
0078     // From full match lists, collect all the stubs associated with the tracklet seed
0079 
0080     // Get seed stubs first
0081     trackstublist.emplace_back(tracklet->innerFPGAStub());
0082     if (tracklet->getISeed() >= (int)N_TRKLSEED + 1)
0083       trackstublist.emplace_back(tracklet->middleFPGAStub());
0084     trackstublist.emplace_back(tracklet->outerFPGAStub());
0085 
0086     // Now get ALL matches (can have multiple per layer)
0087     for (const auto& i : fullmatch1_) {
0088       for (unsigned int j = 0; j < i->nMatches(); j++) {
0089         if (i->getTracklet(j)->TCID() == tracklet->TCID()) {
0090           trackstublist.push_back(i->getMatch(j).second);
0091         }
0092       }
0093     }
0094 
0095     for (const auto& i : fullmatch2_) {
0096       for (unsigned int j = 0; j < i->nMatches(); j++) {
0097         if (i->getTracklet(j)->TCID() == tracklet->TCID()) {
0098           trackstublist.push_back(i->getMatch(j).second);
0099         }
0100       }
0101     }
0102 
0103     for (const auto& i : fullmatch3_) {
0104       for (unsigned int j = 0; j < i->nMatches(); j++) {
0105         if (i->getTracklet(j)->TCID() == tracklet->TCID()) {
0106           trackstublist.push_back(i->getMatch(j).second);
0107         }
0108       }
0109     }
0110 
0111     for (const auto& i : fullmatch4_) {
0112       for (unsigned int j = 0; j < i->nMatches(); j++) {
0113         if (i->getTracklet(j)->TCID() == tracklet->TCID()) {
0114           trackstublist.push_back(i->getMatch(j).second);
0115         }
0116       }
0117     }
0118 
0119     // For merge removal, loop through the resulting list of stubs to calculate their stubids
0120     if (settings_.removalType() == "merge") {
0121       for (const auto& it : trackstublist) {
0122         int layer = it->layer().value() + 1;  // Assume layer (1-6) stub first
0123         if (it->layer().value() < 0) {        // if disk stub, though...
0124           layer = it->disk().value() + 10 * it->disk().value() / abs(it->disk().value());  //disk = +/- 11-15
0125         }
0126         stubidslist.push_back(std::make_pair(layer, it->phiregionaddress()));
0127       }
0128 
0129       // And that's all we need! The rest is just for track fit (done in PurgeDuplicate)
0130 
0131     } else {
0132       // Track fit only called here if not running duplicate removal
0133       // before fit. (e.g. If skipping duplicate removal).
0134       HybridFit hybridFitter(iSector_, settings_, globals_);
0135       hybridFitter.Fit(tracklet, trackstublist);
0136     }
0137   }
0138 }
0139 
0140 #else  // Code for pure Tracklet algo.
0141 
0142 void FitTrack::trackFitChisq(Tracklet* tracklet, std::vector<const Stub*>&, std::vector<std::pair<int, int>>&) {
0143   if (globals_->trackDerTable() == nullptr) {
0144     TrackDerTable* derTablePtr = new TrackDerTable(settings_);
0145 
0146     derTablePtr->readPatternFile(settings_.fitPatternFile());
0147     derTablePtr->fillTable();
0148     if (settings_.debugTracklet()) {
0149       edm::LogVerbatim("Tracklet") << "Number of entries in derivative table: " << derTablePtr->getEntries();
0150     }
0151     assert(derTablePtr->getEntries() != 0);
0152 
0153     globals_->trackDerTable() = derTablePtr;
0154   }
0155 
0156   const TrackDerTable& derTable = *globals_->trackDerTable();
0157 
0158   //First step is to build list of layers and disks.
0159   int layers[N_LAYER];
0160   double r[N_LAYER];
0161   unsigned int nlayers = 0;  // layers with found stub-projections
0162   int disks[N_DISK];
0163   double z[N_DISK];
0164   unsigned int ndisks = 0;  // disks with found stub-projections
0165 
0166   // residuals for each stub
0167   double phiresid[N_FITSTUB];
0168   double zresid[N_FITSTUB];
0169   double phiresidexact[N_FITSTUB];
0170   double zresidexact[N_FITSTUB];
0171   int iphiresid[N_FITSTUB];
0172   int izresid[N_FITSTUB];
0173   double alpha[N_FITSTUB];
0174 
0175   for (unsigned int i = 0; i < N_FITSTUB; i++) {
0176     iphiresid[i] = 0;
0177     izresid[i] = 0;
0178     alpha[i] = 0.0;
0179 
0180     phiresid[i] = 0.0;
0181     zresid[i] = 0.0;
0182     phiresidexact[i] = 0.0;
0183     zresidexact[i] = 0.0;
0184   }
0185 
0186   std::bitset<N_LAYER> lmatches;     //layer matches
0187   std::bitset<N_DISK * 2> dmatches;  //disk matches (2 per disk to separate 2S from PS)
0188 
0189   int mult = 1;
0190 
0191   unsigned int layermask = 0;
0192   unsigned int diskmask = 0;
0193   unsigned int alphaindex = 0;
0194   unsigned int power = 1;
0195 
0196   double t = tracklet->t();
0197   double rinv = tracklet->rinv();
0198 
0199   if (tracklet->isBarrel()) {
0200     for (unsigned int l = 1; l <= N_LAYER; l++) {
0201       if (l == (unsigned int)tracklet->layer() || l == (unsigned int)tracklet->layer() + 1) {
0202         lmatches.set(N_LAYER - l);
0203         layermask |= (1 << (N_LAYER - l));
0204         layers[nlayers++] = l;
0205         continue;
0206       }
0207       if (tracklet->match(l - 1)) {
0208         const Residual& resid = tracklet->resid(l - 1);
0209         lmatches.set(N_LAYER - l);
0210         layermask |= (1 << (N_LAYER - l));
0211         phiresid[nlayers] = resid.phiresidapprox();
0212         zresid[nlayers] = resid.rzresidapprox();
0213         phiresidexact[nlayers] = resid.phiresid();
0214         zresidexact[nlayers] = resid.rzresid();
0215         iphiresid[nlayers] = resid.fpgaphiresid().value();
0216         izresid[nlayers] = resid.fpgarzresid().value();
0217 
0218         layers[nlayers++] = l;
0219       }
0220     }
0221 
0222     for (unsigned int d = 1; d <= N_DISK; d++) {
0223       if (layermask & (1 << (d - 1)))
0224         continue;
0225 
0226       if (mult == 1 << (3 * settings_.alphaBitsTable()))
0227         continue;
0228 
0229       if (ndisks + nlayers >= N_FITSTUB)
0230         continue;
0231       if (tracklet->match(N_LAYER + d - 1)) {
0232         const Residual& resid = tracklet->resid(N_LAYER + d - 1);
0233         double pitch = settings_.stripPitch(resid.stubptr()->l1tstub()->isPSmodule());
0234         if (std::abs(resid.stubptr()->l1tstub()->alpha(pitch)) < 1e-20) {
0235           dmatches.set(2 * d - 1);
0236           diskmask |= (1 << (2 * (N_DISK - d) + 1));
0237         } else {
0238           int ialpha = resid.stubptr()->alpha().value();
0239           int nalpha = resid.stubptr()->alpha().nbits();
0240           nalpha = nalpha - settings_.alphaBitsTable();
0241           ialpha = (1 << (settings_.alphaBitsTable() - 1)) + (ialpha >> nalpha);
0242 
0243           alphaindex += ialpha * power;
0244           power = power << settings_.alphaBitsTable();
0245           dmatches.set(2 * (N_DISK - d));
0246           diskmask |= (1 << (2 * (N_DISK - d)));
0247           mult = mult << settings_.alphaBitsTable();
0248         }
0249         alpha[ndisks] = resid.stubptr()->l1tstub()->alpha(pitch);
0250         phiresid[nlayers + ndisks] = resid.phiresidapprox();
0251         zresid[nlayers + ndisks] = resid.rzresidapprox();
0252         phiresidexact[nlayers + ndisks] = resid.phiresid();
0253         zresidexact[nlayers + ndisks] = resid.rzresid();
0254         iphiresid[nlayers + ndisks] = resid.fpgaphiresid().value();
0255         izresid[nlayers + ndisks] = resid.fpgarzresid().value();
0256 
0257         disks[ndisks++] = d;
0258       }
0259     }
0260 
0261     if (settings_.writeMonitorData("HitPattern")) {
0262       if (mult <= 1 << (3 * settings_.alphaBitsTable())) {
0263         globals_->ofstream("hitpattern.txt")
0264             << lmatches.to_string() << " " << dmatches.to_string() << " " << mult << endl;
0265       }
0266     }
0267   }
0268 
0269   if (tracklet->isDisk()) {
0270     for (unsigned int l = 1; l <= 2; l++) {
0271       if (tracklet->match(l - 1)) {
0272         lmatches.set(N_LAYER - l);
0273 
0274         layermask |= (1 << (N_LAYER - l));
0275         const Residual& resid = tracklet->resid(l - 1);
0276         phiresid[nlayers] = resid.phiresidapprox();
0277         zresid[nlayers] = resid.rzresidapprox();
0278         phiresidexact[nlayers] = resid.phiresid();
0279         zresidexact[nlayers] = resid.rzresid();
0280         iphiresid[nlayers] = resid.fpgaphiresid().value();
0281         izresid[nlayers] = resid.fpgarzresid().value();
0282 
0283         layers[nlayers++] = l;
0284       }
0285     }
0286 
0287     for (int d1 = 1; d1 <= N_DISK; d1++) {
0288       int d = d1;
0289 
0290       // skip F/B5 if there's already a L2 match
0291       if (d == 5 and layermask & (1 << 4))
0292         continue;
0293 
0294       if (tracklet->fpgat().value() < 0.0)
0295         d = -d1;
0296       if (d1 == abs(tracklet->disk()) || d1 == abs(tracklet->disk()) + 1) {
0297         dmatches.set(2 * d1 - 1);
0298         diskmask |= (1 << (2 * (N_DISK - d1) + 1));
0299         alpha[ndisks] = 0.0;
0300         disks[ndisks++] = d;
0301         continue;
0302       }
0303 
0304       if (ndisks + nlayers >= N_FITSTUB)
0305         continue;
0306       if (tracklet->match(N_LAYER + abs(d) - 1)) {
0307         const Residual& resid = tracklet->resid(N_LAYER + abs(d) - 1);
0308         double pitch = settings_.stripPitch(resid.stubptr()->l1tstub()->isPSmodule());
0309         if (std::abs(resid.stubptr()->l1tstub()->alpha(pitch)) < 1e-20) {
0310           dmatches.set(2 * d1 - 1);
0311           diskmask |= (1 << (2 * (N_DISK - d1) + 1));
0312         } else {
0313           int ialpha = resid.stubptr()->alpha().value();
0314           int nalpha = resid.stubptr()->alpha().nbits();
0315           nalpha = nalpha - settings_.alphaBitsTable();
0316           ialpha = (1 << (settings_.alphaBitsTable() - 1)) + (ialpha >> nalpha);
0317 
0318           alphaindex += ialpha * power;
0319           power = power << settings_.alphaBitsTable();
0320           dmatches.set(2 * (N_DISK - d1));
0321           diskmask |= (1 << (2 * (N_DISK - d1)));
0322           mult = mult << settings_.alphaBitsTable();
0323         }
0324 
0325         alpha[ndisks] = resid.stubptr()->l1tstub()->alpha(pitch);
0326         assert(std::abs(resid.phiresidapprox()) < 0.2);
0327         phiresid[nlayers + ndisks] = resid.phiresidapprox();
0328         zresid[nlayers + ndisks] = resid.rzresidapprox();
0329         assert(std::abs(resid.phiresid()) < 0.2);
0330         phiresidexact[nlayers + ndisks] = resid.phiresid();
0331         zresidexact[nlayers + ndisks] = resid.rzresid();
0332         iphiresid[nlayers + ndisks] = resid.fpgaphiresid().value();
0333         izresid[nlayers + ndisks] = resid.fpgarzresid().value();
0334 
0335         disks[ndisks++] = d;
0336       }
0337     }
0338   }
0339 
0340   if (tracklet->isOverlap()) {
0341     for (unsigned int l = 1; l <= 2; l++) {
0342       if (l == (unsigned int)tracklet->layer()) {
0343         lmatches.set(N_LAYER - l);
0344         layermask |= (1 << (N_LAYER - l));
0345         layers[nlayers++] = l;
0346         continue;
0347       }
0348       if (tracklet->match(l - 1)) {
0349         lmatches.set(N_LAYER - l);
0350         layermask |= (1 << (N_LAYER - l));
0351         const Residual& resid = tracklet->resid(l - 1);
0352         assert(std::abs(resid.phiresidapprox()) < 0.2);
0353         phiresid[nlayers] = resid.phiresidapprox();
0354         zresid[nlayers] = resid.rzresidapprox();
0355         assert(std::abs(resid.phiresid()) < 0.2);
0356         phiresidexact[nlayers] = resid.phiresid();
0357         zresidexact[nlayers] = resid.rzresid();
0358         iphiresid[nlayers] = resid.fpgaphiresid().value();
0359         izresid[nlayers] = resid.fpgarzresid().value();
0360 
0361         layers[nlayers++] = l;
0362       }
0363     }
0364 
0365     for (unsigned int d1 = 1; d1 <= N_DISK; d1++) {
0366       if (mult == 1 << (3 * settings_.alphaBitsTable()))
0367         continue;
0368       int d = d1;
0369       if (tracklet->fpgat().value() < 0.0)
0370         d = -d1;
0371       if (d == tracklet->disk()) {  //All seeds in PS modules
0372         disks[ndisks] = tracklet->disk();
0373         dmatches.set(2 * d1 - 1);
0374         diskmask |= (1 << (2 * (N_DISK - d1) + 1));
0375         ndisks++;
0376         continue;
0377       }
0378 
0379       if (ndisks + nlayers >= N_FITSTUB)
0380         continue;
0381       if (tracklet->match(N_LAYER + abs(d) - 1)) {
0382         const Residual& resid = tracklet->resid(N_LAYER + abs(d) - 1);
0383         double pitch = settings_.stripPitch(resid.stubptr()->l1tstub()->isPSmodule());
0384         if (std::abs(resid.stubptr()->l1tstub()->alpha(pitch)) < 1e-20) {
0385           dmatches.set(2 * (N_DISK - d1));
0386           diskmask |= (1 << (2 * (N_DISK - d1) + 1));
0387           FPGAWord tmp;
0388           tmp.set(diskmask, 10);
0389         } else {
0390           int ialpha = resid.stubptr()->alpha().value();
0391           int nalpha = resid.stubptr()->alpha().nbits();
0392           nalpha = nalpha - settings_.alphaBitsTable();
0393           ialpha = (1 << (settings_.alphaBitsTable() - 1)) + (ialpha >> nalpha);
0394 
0395           alphaindex += ialpha * power;
0396           power = power << settings_.alphaBitsTable();
0397           dmatches.set(2 * (N_DISK - d1));
0398           diskmask |= (1 << (2 * (N_DISK - d1)));
0399           FPGAWord tmp;
0400           tmp.set(diskmask, 10);
0401           mult = mult << settings_.alphaBitsTable();
0402         }
0403 
0404         alpha[ndisks] = resid.stubptr()->l1tstub()->alpha(pitch);
0405         assert(std::abs(resid.phiresidapprox()) < 0.2);
0406         phiresid[nlayers + ndisks] = resid.phiresidapprox();
0407         zresid[nlayers + ndisks] = resid.rzresidapprox();
0408         assert(std::abs(resid.phiresid()) < 0.2);
0409         phiresidexact[nlayers + ndisks] = resid.phiresid();
0410         zresidexact[nlayers + ndisks] = resid.rzresid();
0411         iphiresid[nlayers + ndisks] = resid.fpgaphiresid().value();
0412         izresid[nlayers + ndisks] = resid.fpgarzresid().value();
0413 
0414         disks[ndisks++] = d;
0415       }
0416     }
0417   }
0418 
0419   int rinvindex =
0420       (1 << (settings_.nrinvBitsTable() - 1)) * rinv / settings_.rinvmax() + (1 << (settings_.nrinvBitsTable() - 1));
0421   if (rinvindex < 0)
0422     rinvindex = 0;
0423   if (rinvindex >= (1 << settings_.nrinvBitsTable()))
0424     rinvindex = (1 << settings_.nrinvBitsTable()) - 1;
0425 
0426   const TrackDer* derivatives = derTable.getDerivatives(layermask, diskmask, alphaindex, rinvindex);
0427 
0428   if (derivatives == nullptr) {
0429     if (settings_.warnNoDer()) {
0430       FPGAWord tmpl, tmpd;
0431       tmpl.set(layermask, 6);
0432       tmpd.set(diskmask, 10);
0433       edm::LogVerbatim("Tracklet") << "No derivative for layermask, diskmask : " << layermask << " " << tmpl.str()
0434                                    << " " << diskmask << " " << tmpd.str() << " eta = " << asinh(t);
0435     }
0436     return;
0437   }
0438 
0439   double ttabi = TrackDerTable::tpar(settings_, diskmask, layermask);
0440   if (t < 0.0)
0441     ttabi = -ttabi;
0442   double ttab = ttabi;
0443 
0444   if (settings_.debugTracklet()) {
0445     edm::LogVerbatim("Tracklet") << "Doing trackfit in  " << getName();
0446   }
0447 
0448   int sign = 1;
0449   if (t < 0.0)
0450     sign = -1;
0451 
0452   double rstub[6];
0453 
0454   double realrstub[3];
0455   realrstub[0] = -1.0;
0456   realrstub[1] = -1.0;
0457   realrstub[2] = -1.0;
0458 
0459   for (unsigned i = 0; i < nlayers; i++) {
0460     r[i] = settings_.rmean(layers[i] - 1);
0461     if (layers[i] == tracklet->layer()) {
0462       if (tracklet->isOverlap()) {
0463         realrstub[i] = tracklet->outerFPGAStub()->l1tstub()->r();
0464       } else {
0465         realrstub[i] = tracklet->innerFPGAStub()->l1tstub()->r();
0466       }
0467     }
0468     if (layers[i] == tracklet->layer() + 1) {
0469       realrstub[i] = tracklet->outerFPGAStub()->l1tstub()->r();
0470     }
0471     if (tracklet->match(layers[i] - 1) && layers[i] < 4) {
0472       const Stub* stubptr = tracklet->resid(layers[i] - 1).stubptr();
0473       realrstub[i] = stubptr->l1tstub()->r();
0474       assert(std::abs(realrstub[i] - r[i]) < 5.0);
0475     }
0476     rstub[i] = r[i];
0477   }
0478   for (unsigned i = 0; i < ndisks; i++) {
0479     z[i] = sign * settings_.zmean(abs(disks[i]) - 1);
0480     rstub[i + nlayers] = z[i] / ttabi;
0481   }
0482 
0483   double D[N_FITPARAM][N_FITSTUB * 2];
0484   double MinvDt[N_FITPARAM][N_FITSTUB * 2];
0485   int iD[N_FITPARAM][N_FITSTUB * 2];
0486   int iMinvDt[N_FITPARAM][N_FITSTUB * 2];
0487   double sigma[N_FITSTUB * 2];
0488   double kfactor[N_FITSTUB * 2];
0489 
0490   unsigned int n = nlayers + ndisks;
0491 
0492   if (settings_.exactderivatives()) {
0493     TrackDerTable::calculateDerivatives(
0494         settings_, nlayers, r, ndisks, z, alpha, t, rinv, D, iD, MinvDt, iMinvDt, sigma, kfactor);
0495     ttabi = t;
0496     ttab = t;
0497   } else {
0498     if (settings_.exactderivativesforfloating()) {
0499       TrackDerTable::calculateDerivatives(
0500           settings_, nlayers, r, ndisks, z, alpha, t, rinv, D, iD, MinvDt, iMinvDt, sigma, kfactor);
0501 
0502       double MinvDtDummy[N_FITPARAM][N_FITSTUB * 2];
0503       derivatives->fill(tracklet->fpgat().value(), MinvDtDummy, iMinvDt);
0504       ttab = t;
0505     } else {
0506       derivatives->fill(tracklet->fpgat().value(), MinvDt, iMinvDt);
0507     }
0508   }
0509 
0510   if (!settings_.exactderivatives()) {
0511     for (unsigned int i = 0; i < nlayers; i++) {
0512       if (r[i] > settings_.rPS2S())
0513         continue;
0514       for (unsigned int ii = 0; ii < nlayers; ii++) {
0515         if (r[ii] > settings_.rPS2S())
0516           continue;
0517 
0518         double tder = derivatives->tdzcorr(i, ii);
0519         double zder = derivatives->z0dzcorr(i, ii);
0520 
0521         double dr = realrstub[i] - r[i];
0522 
0523         MinvDt[2][2 * ii + 1] += dr * tder;
0524         MinvDt[3][2 * ii + 1] += dr * zder;
0525 
0526         int itder = derivatives->itdzcorr(i, ii);
0527         int izder = derivatives->iz0dzcorr(i, ii);
0528 
0529         int idr = dr / settings_.kr();
0530 
0531         iMinvDt[2][2 * ii + 1] += ((idr * itder) >> settings_.rcorrbits());
0532         iMinvDt[3][2 * ii + 1] += ((idr * izder) >> settings_.rcorrbits());
0533       }
0534     }
0535   }
0536 
0537   double rinvseed = tracklet->rinvapprox();
0538   double phi0seed = tracklet->phi0approx();
0539   double tseed = tracklet->tapprox();
0540   double z0seed = tracklet->z0approx();
0541 
0542   double rinvseedexact = tracklet->rinv();
0543   double phi0seedexact = tracklet->phi0();
0544   double tseedexact = tracklet->t();
0545   double z0seedexact = tracklet->z0();
0546 
0547   double chisqseedexact = 0.0;
0548 
0549   double delta[2 * N_FITSTUB];
0550   double deltaexact[2 * N_FITSTUB];
0551   int idelta[2 * N_FITSTUB];
0552 
0553   for (unsigned int i = 0; i < 2 * N_FITSTUB; i++) {
0554     delta[i] = 0.0;
0555     deltaexact[i] = 0.0;
0556     idelta[i] = 0;
0557   }
0558 
0559   int j = 0;
0560 
0561   for (unsigned int i = 0; i < n; i++) {
0562     if (i >= nlayers) {
0563       iphiresid[i] *= (t / ttabi);
0564       phiresid[i] *= (t / ttab);
0565       phiresidexact[i] *= (t / ttab);
0566     }
0567 
0568     idelta[j] = iphiresid[i];
0569     delta[j] = phiresid[i];
0570     if (std::abs(phiresid[i]) > 0.2) {
0571       edm::LogWarning("Tracklet") << getName() << " WARNING too large phiresid: " << phiresid[i] << " "
0572                                   << phiresidexact[i];
0573     }
0574     assert(std::abs(phiresid[i]) < 1.0);
0575     assert(std::abs(phiresidexact[i]) < 1.0);
0576     deltaexact[j++] = phiresidexact[i];
0577 
0578     idelta[j] = izresid[i];
0579     delta[j] = zresid[i];
0580     deltaexact[j++] = zresidexact[i];
0581 
0582     chisqseedexact += (deltaexact[j - 2] * deltaexact[j - 2] + deltaexact[j - 1] * deltaexact[j - 1]);
0583   }
0584   assert(j <= 12);
0585 
0586   double drinv = 0.0;
0587   double dphi0 = 0.0;
0588   double dt = 0.0;
0589   double dz0 = 0.0;
0590 
0591   double drinvexact = 0.0;
0592   double dphi0exact = 0.0;
0593   double dtexact = 0.0;
0594   double dz0exact = 0.0;
0595 
0596   int idrinv = 0;
0597   int idphi0 = 0;
0598   int idt = 0;
0599   int idz0 = 0;
0600 
0601   double drinv_cov = 0.0;
0602   double dphi0_cov = 0.0;
0603   double dt_cov = 0.0;
0604   double dz0_cov = 0.0;
0605 
0606   double drinv_covexact = 0.0;
0607   double dphi0_covexact = 0.0;
0608   double dt_covexact = 0.0;
0609   double dz0_covexact = 0.0;
0610 
0611   for (unsigned int j = 0; j < 2 * n; j++) {
0612     drinv -= MinvDt[0][j] * delta[j];
0613     dphi0 -= MinvDt[1][j] * delta[j];
0614     dt -= MinvDt[2][j] * delta[j];
0615     dz0 -= MinvDt[3][j] * delta[j];
0616 
0617     drinv_cov += D[0][j] * delta[j];
0618     dphi0_cov += D[1][j] * delta[j];
0619     dt_cov += D[2][j] * delta[j];
0620     dz0_cov += D[3][j] * delta[j];
0621 
0622     drinvexact -= MinvDt[0][j] * deltaexact[j];
0623     dphi0exact -= MinvDt[1][j] * deltaexact[j];
0624     dtexact -= MinvDt[2][j] * deltaexact[j];
0625     dz0exact -= MinvDt[3][j] * deltaexact[j];
0626 
0627     drinv_covexact += D[0][j] * deltaexact[j];
0628     dphi0_covexact += D[1][j] * deltaexact[j];
0629     dt_covexact += D[2][j] * deltaexact[j];
0630     dz0_covexact += D[3][j] * deltaexact[j];
0631 
0632     idrinv += ((iMinvDt[0][j] * idelta[j]));
0633     idphi0 += ((iMinvDt[1][j] * idelta[j]));
0634     idt += ((iMinvDt[2][j] * idelta[j]));
0635     idz0 += ((iMinvDt[3][j] * idelta[j]));
0636 
0637     if (false && j % 2 == 0) {
0638       edm::LogVerbatim("Tracklet") << "DEBUG CHI2FIT " << j << " " << rinvseed << " + " << MinvDt[0][j] * delta[j]
0639                                    << " " << MinvDt[0][j] << " " << delta[j] * rstub[j / 2] * 10000 << " \n"
0640                                    << j << " " << tracklet->fpgarinv().value() * settings_.krinvpars() << " + "
0641                                    << ((iMinvDt[0][j] * idelta[j])) * settings_.krinvpars() / 1024.0 << " "
0642                                    << iMinvDt[0][j] * settings_.krinvpars() / settings_.kphi() / 1024.0 << " "
0643                                    << idelta[j] * settings_.kphi() * rstub[j / 2] * 10000 << " " << idelta[j];
0644     }
0645   }
0646 
0647   double deltaChisqexact =
0648       drinvexact * drinv_covexact + dphi0exact * dphi0_covexact + dtexact * dt_covexact + dz0exact * dz0_covexact;
0649 
0650   int irinvseed = tracklet->fpgarinv().value();
0651   int iphi0seed = tracklet->fpgaphi0().value();
0652 
0653   int itseed = tracklet->fpgat().value();
0654   int iz0seed = tracklet->fpgaz0().value();
0655 
0656   int irinvfit = irinvseed + ((idrinv + (1 << settings_.fitrinvbitshift())) >> settings_.fitrinvbitshift());
0657 
0658   int iphi0fit = iphi0seed + (idphi0 >> settings_.fitphi0bitshift());
0659   int itfit = itseed + (idt >> settings_.fittbitshift());
0660 
0661   int iz0fit = iz0seed + (idz0 >> settings_.fitz0bitshift());
0662 
0663   double rinvfit = rinvseed - drinv;
0664   double phi0fit = phi0seed - dphi0;
0665 
0666   double tfit = tseed - dt;
0667   double z0fit = z0seed - dz0;
0668 
0669   double rinvfitexact = rinvseedexact - drinvexact;
0670   double phi0fitexact = phi0seedexact - dphi0exact;
0671 
0672   double tfitexact = tseedexact - dtexact;
0673   double z0fitexact = z0seedexact - dz0exact;
0674 
0675   double chisqfitexact = chisqseedexact + deltaChisqexact;
0676 
0677   ////////////// NEW CHISQ /////////////////////
0678   bool NewChisqDebug = false;
0679   double chisqfit = 0.0;
0680   uint ichisqfit = 0;
0681 
0682   double phifactor;
0683   double rzfactor;
0684   double iphifactor;
0685   double irzfactor;
0686   int k = 0;  // column index of D matrix
0687 
0688   if (NewChisqDebug) {
0689     edm::LogVerbatim("Tracklet") << "OG chisq: \n"
0690                                  << "drinv/cov = " << drinv << "/" << drinv_cov << " \n"
0691                                  << "dphi0/cov = " << drinv << "/" << dphi0_cov << " \n"
0692                                  << "dt/cov = " << drinv << "/" << dt_cov << " \n"
0693                                  << "dz0/cov = " << drinv << "/" << dz0_cov << "\n";
0694     std::string myout = "D[0][k]= ";
0695     for (unsigned int i = 0; i < 2 * n; i++) {
0696       myout += std::to_string(D[0][i]);
0697       myout += ", ";
0698     }
0699     edm::LogVerbatim("Tracklet") << myout;
0700   }
0701 
0702   for (unsigned int i = 0; i < n; i++) {  // loop over stubs
0703 
0704     phifactor = rstub[k / 2] * delta[k] / sigma[k] + D[0][k] * drinv + D[1][k] * dphi0 + D[2][k] * dt + D[3][k] * dz0;
0705     iphifactor = kfactor[k] * rstub[k / 2] * idelta[k] * (1 << settings_.chisqphifactbits()) / sigma[k] -
0706                  iD[0][k] * idrinv - iD[1][k] * idphi0 - iD[2][k] * idt - iD[3][k] * idz0;
0707 
0708     if (NewChisqDebug) {
0709       edm::LogVerbatim("Tracklet") << "delta[k]/sigma = " << delta[k] / sigma[k] << "  delta[k] = " << delta[k] << "\n"
0710                                    << "sum = " << phifactor - delta[k] / sigma[k] << "  drinvterm = " << D[0][k] * drinv
0711                                    << "  dphi0term = " << D[1][k] * dphi0 << "  dtterm = " << D[2][k] * dt
0712                                    << "  dz0term = " << D[3][k] * dz0 << "\n  phifactor = " << phifactor;
0713     }
0714 
0715     chisqfit += phifactor * phifactor;
0716     ichisqfit += iphifactor * iphifactor / (1 << (2 * settings_.chisqphifactbits() - 4));
0717 
0718     k++;
0719 
0720     rzfactor = delta[k] / sigma[k] + D[0][k] * drinv + D[1][k] * dphi0 + D[2][k] * dt + D[3][k] * dz0;
0721     irzfactor = kfactor[k] * idelta[k] * (1 << settings_.chisqzfactbits()) / sigma[k] - iD[0][k] * idrinv -
0722                 iD[1][k] * idphi0 - iD[2][k] * idt - iD[3][k] * idz0;
0723 
0724     if (NewChisqDebug) {
0725       edm::LogVerbatim("Tracklet") << "delta[k]/sigma = " << delta[k] / sigma[k] << "  delta[k] = " << delta[k] << "\n"
0726                                    << "sum = " << rzfactor - delta[k] / sigma[k] << "  drinvterm = " << D[0][k] * drinv
0727                                    << "  dphi0term = " << D[1][k] * dphi0 << "  dtterm = " << D[2][k] * dt
0728                                    << "  dz0term = " << D[3][k] * dz0 << "\n  rzfactor = " << rzfactor;
0729     }
0730 
0731     chisqfit += rzfactor * rzfactor;
0732     ichisqfit += irzfactor * irzfactor / (1 << (2 * settings_.chisqzfactbits() - 4));
0733 
0734     k++;
0735   }
0736 
0737   if (settings_.writeMonitorData("ChiSq")) {
0738     globals_->ofstream("chisq.txt") << asinh(itfit * settings_.ktpars()) << " " << chisqfit << " " << ichisqfit / 16.0
0739                                     << endl;
0740   }
0741 
0742   // Chisquare per DOF capped out at 11 bits, so 15 is an educated guess
0743   if (ichisqfit >= (1 << 15)) {
0744     if (NewChisqDebug) {
0745       edm::LogVerbatim("Tracklet") << "CHISQUARE (" << ichisqfit << ") LARGER THAN 11 BITS!";
0746     }
0747     ichisqfit = (1 << 15) - 1;
0748   }
0749 
0750   // Eliminate lower bits to fit in 8 bits
0751   ichisqfit = ichisqfit >> 7;
0752   // Probably redundant... enforce 8 bit cap
0753   if (ichisqfit >= (1 << 8))
0754     ichisqfit = (1 << 8) - 1;
0755 
0756   double phicrit = phi0fit - asin(0.5 * settings_.rcrit() * rinvfit);
0757   bool keep = (phicrit > settings_.phicritmin()) && (phicrit < settings_.phicritmax());
0758 
0759   if (!keep) {
0760     return;
0761   }
0762 
0763   // NOTE: setFitPars in Tracklet.h now accepts chi2 r-phi and chi2 r-z values. This class only has access
0764   // to the composite chi2. When setting fit parameters on a tracklet, this places all of the chi2 into the
0765   // r-phi fit, and sets the r-z fit value to zero.
0766   //
0767   // This is also true for the call to setFitPars in trackFitFake.
0768   tracklet->setFitPars(rinvfit,
0769                        phi0fit,
0770                        0.0,
0771                        tfit,
0772                        z0fit,
0773                        chisqfit,
0774                        0.0,
0775                        rinvfitexact,
0776                        phi0fitexact,
0777                        0.0,
0778                        tfitexact,
0779                        z0fitexact,
0780                        chisqfitexact,
0781                        0.0,
0782                        irinvfit,
0783                        iphi0fit,
0784                        0,
0785                        itfit,
0786                        iz0fit,
0787                        ichisqfit,
0788                        0,
0789                        0);
0790 }
0791 
0792 #endif
0793 
0794 void FitTrack::trackFitFake(Tracklet* tracklet, std::vector<const Stub*>&, std::vector<std::pair<int, int>>&) {
0795   tracklet->setFitPars(tracklet->rinvapprox(),
0796                        tracklet->phi0approx(),
0797                        tracklet->d0approx(),
0798                        tracklet->tapprox(),
0799                        tracklet->z0approx(),
0800                        0.0,
0801                        0.0,
0802                        tracklet->rinv(),
0803                        tracklet->phi0(),
0804                        tracklet->d0(),
0805                        tracklet->t(),
0806                        tracklet->z0(),
0807                        0.0,
0808                        0.0,
0809                        tracklet->fpgarinv().value(),
0810                        tracklet->fpgaphi0().value(),
0811                        tracklet->fpgad0().value(),
0812                        tracklet->fpgat().value(),
0813                        tracklet->fpgaz0().value(),
0814                        0,
0815                        0,
0816                        0);
0817   return;
0818 }
0819 
0820 std::vector<Tracklet*> FitTrack::orderedMatches(vector<FullMatchMemory*>& fullmatch) {
0821   std::vector<Tracklet*> tmp;
0822 
0823   std::vector<unsigned int> indexArray;
0824   for (auto& imatch : fullmatch) {
0825     //check that we have correct order
0826     if (imatch->nMatches() > 1) {
0827       for (unsigned int j = 0; j < imatch->nMatches() - 1; j++) {
0828         assert(imatch->getTracklet(j)->TCID() <= imatch->getTracklet(j + 1)->TCID());
0829       }
0830     }
0831 
0832     if (settings_.debugTracklet() && imatch->nMatches() != 0) {
0833       edm::LogVerbatim("Tracklet") << "orderedMatches: " << imatch->getName() << " " << imatch->nMatches();
0834     }
0835 
0836     indexArray.push_back(0);
0837   }
0838 
0839   int bestIndex = -1;
0840   do {
0841     int bestTCID = -1;
0842     bestIndex = -1;
0843     for (unsigned int i = 0; i < fullmatch.size(); i++) {
0844       if (indexArray[i] >= fullmatch[i]->nMatches()) {
0845         //skip as we were at the end
0846         continue;
0847       }
0848       int TCID = fullmatch[i]->getTracklet(indexArray[i])->TCID();
0849       if (TCID < bestTCID || bestTCID < 0) {
0850         bestTCID = TCID;
0851         bestIndex = i;
0852       }
0853     }
0854     if (bestIndex != -1) {
0855       tmp.push_back(fullmatch[bestIndex]->getTracklet(indexArray[bestIndex]));
0856       indexArray[bestIndex]++;
0857     }
0858   } while (bestIndex != -1);
0859 
0860   for (unsigned int i = 0; i < tmp.size(); i++) {
0861     if (i > 0) {
0862       //This allows for equal TCIDs. This means that we can e.g. have a track seeded in L1L2 that projects to both L3 and D4.
0863       //The algorithm will pick up the first hit and drop the second.
0864       if (tmp[i - 1]->TCID() > tmp[i]->TCID()) {
0865         edm::LogVerbatim("Tracklet") << "Wrong TCID ordering in " << getName() << " : " << tmp[i - 1]->TCID() << " "
0866                                      << tmp[i]->TCID();
0867       }
0868     }
0869   }
0870 
0871   return tmp;
0872 }
0873 
0874 // Adds the fitted track to the output memories to be used by pure Tracklet algo.
0875 // (Also used by Hybrid algo with non-exact Old KF emulation)
0876 // Also create output streams, that bypass these memories, (so can include gaps in time),
0877 // to be used by Hybrid case with exact New KF emulation.
0878 
0879 void FitTrack::execute(deque<string>& streamTrackRaw,
0880                        vector<deque<StubStreamData>>& streamsStubRaw,
0881                        unsigned int iSector) {
0882   // merge
0883   const std::vector<Tracklet*>& matches1 = orderedMatches(fullmatch1_);
0884   const std::vector<Tracklet*>& matches2 = orderedMatches(fullmatch2_);
0885   const std::vector<Tracklet*>& matches3 = orderedMatches(fullmatch3_);
0886   const std::vector<Tracklet*>& matches4 = orderedMatches(fullmatch4_);
0887 
0888   iSector_ = iSector;
0889 
0890   if (settings_.debugTracklet() && (matches1.size() + matches2.size() + matches3.size() + matches4.size()) > 0) {
0891     for (auto& imatch : fullmatch1_) {
0892       edm::LogVerbatim("Tracklet") << imatch->getName() << " " << imatch->nMatches();
0893     }
0894     edm::LogVerbatim("Tracklet") << getName() << " matches : " << matches1.size() << " " << matches2.size() << " "
0895                                  << matches3.size() << " " << matches4.size();
0896   }
0897 
0898   unsigned int indexArray[4];
0899   for (unsigned int i = 0; i < 4; i++) {
0900     indexArray[i] = 0;
0901   }
0902 
0903   unsigned int countAll = 0;
0904   unsigned int countFit = 0;
0905 
0906   Tracklet* bestTracklet = nullptr;
0907   do {
0908     countAll++;
0909     bestTracklet = nullptr;
0910 
0911     if (indexArray[0] < matches1.size()) {
0912       if (bestTracklet == nullptr) {
0913         bestTracklet = matches1[indexArray[0]];
0914       } else {
0915         if (matches1[indexArray[0]]->TCID() < bestTracklet->TCID())
0916           bestTracklet = matches1[indexArray[0]];
0917       }
0918     }
0919 
0920     if (indexArray[1] < matches2.size()) {
0921       if (bestTracklet == nullptr) {
0922         bestTracklet = matches2[indexArray[1]];
0923       } else {
0924         if (matches2[indexArray[1]]->TCID() < bestTracklet->TCID())
0925           bestTracklet = matches2[indexArray[1]];
0926       }
0927     }
0928 
0929     if (indexArray[2] < matches3.size()) {
0930       if (bestTracklet == nullptr) {
0931         bestTracklet = matches3[indexArray[2]];
0932       } else {
0933         if (matches3[indexArray[2]]->TCID() < bestTracklet->TCID())
0934           bestTracklet = matches3[indexArray[2]];
0935       }
0936     }
0937 
0938     if (indexArray[3] < matches4.size()) {
0939       if (bestTracklet == nullptr) {
0940         bestTracklet = matches4[indexArray[3]];
0941       } else {
0942         if (matches4[indexArray[3]]->TCID() < bestTracklet->TCID())
0943           bestTracklet = matches4[indexArray[3]];
0944       }
0945     }
0946 
0947     if (bestTracklet == nullptr)
0948       break;
0949 
0950     //Counts total number of matched hits
0951     int nMatches = 0;
0952 
0953     //Counts unique hits in each layer
0954     int nMatchesUniq = 0;
0955     bool match = false;
0956 
0957     while (indexArray[0] < matches1.size() && matches1[indexArray[0]] == bestTracklet) {
0958       indexArray[0]++;
0959       nMatches++;
0960       match = true;
0961     }
0962 
0963     if (match)
0964       nMatchesUniq++;
0965     match = false;
0966 
0967     while (indexArray[1] < matches2.size() && matches2[indexArray[1]] == bestTracklet) {
0968       indexArray[1]++;
0969       nMatches++;
0970       match = true;
0971     }
0972 
0973     if (match)
0974       nMatchesUniq++;
0975     match = false;
0976 
0977     while (indexArray[2] < matches3.size() && matches3[indexArray[2]] == bestTracklet) {
0978       indexArray[2]++;
0979       nMatches++;
0980       match = true;
0981     }
0982 
0983     if (match)
0984       nMatchesUniq++;
0985     match = false;
0986 
0987     while (indexArray[3] < matches4.size() && matches4[indexArray[3]] == bestTracklet) {
0988       indexArray[3]++;
0989       nMatches++;
0990       match = true;
0991     }
0992 
0993     if (match)
0994       nMatchesUniq++;
0995 
0996     if (settings_.debugTracklet()) {
0997       edm::LogVerbatim("Tracklet") << getName() << " : nMatches = " << nMatches << " nMatchesUniq = " << nMatchesUniq
0998                                    << " " << asinh(bestTracklet->t());
0999     }
1000 
1001     std::vector<const Stub*> trackstublist;
1002     std::vector<std::pair<int, int>> stubidslist;
1003     // Track Builder cut of >= 4 layers with stubs.
1004     if ((bestTracklet->getISeed() >= (int)N_SEED_PROMPT && nMatchesUniq >= 1) ||
1005         nMatchesUniq >= 2) {  //For seeds index >=8 (triplet seeds), there are three stubs associated from start.
1006       countFit++;
1007 
1008 #ifdef USEHYBRID
1009       if (settings_.fakefit()) {
1010         trackFitFake(bestTracklet, trackstublist, stubidslist);
1011       } else {
1012         trackFitKF(bestTracklet, trackstublist, stubidslist);
1013       }
1014 #else
1015       if (settings_.fakefit()) {
1016         trackFitFake(bestTracklet, trackstublist, stubidslist);
1017       } else {
1018         trackFitChisq(bestTracklet, trackstublist, stubidslist);
1019       }
1020 #endif
1021 
1022       if (settings_.removalType() == "merge") {
1023         trackfit_->addStubList(trackstublist);
1024         trackfit_->addStubidsList(stubidslist);
1025         bestTracklet->setTrackIndex(trackfit_->nTracks());
1026         trackfit_->addTrack(bestTracklet);
1027       } else if (bestTracklet->fit()) {
1028         assert(trackfit_ != nullptr);
1029         if (settings_.writeMonitorData("Seeds")) {
1030           ofstream fout("seeds.txt", ofstream::app);
1031           fout << __FILE__ << ":" << __LINE__ << " " << name_ << "_"
1032                << " " << bestTracklet->getISeed() << endl;
1033           fout.close();
1034         }
1035         bestTracklet->setTrackIndex(trackfit_->nTracks());
1036         trackfit_->addTrack(bestTracklet);
1037       }
1038     }
1039 
1040     // store bit and clock accurate TB output
1041     if (settings_.storeTrackBuilderOutput() && bestTracklet) {
1042       // add gap if TrackBuilder rejected track (due to too few stub layers).
1043       if (!bestTracklet->fit()) {
1044         static const string invalid = "0";
1045         streamTrackRaw.emplace_back(invalid);
1046         for (auto& stream : streamsStubRaw)
1047           stream.emplace_back(StubStreamData());
1048         continue;
1049       }
1050       // convert Track word
1051       const string rinv = bestTracklet->fpgarinv().str();
1052       const string phi0 = bestTracklet->fpgaphi0().str();
1053       const string z0 = bestTracklet->fpgaz0().str();
1054       const string t = bestTracklet->fpgat().str();
1055       const int seedType = bestTracklet->getISeed();
1056       const string seed = TTBV(seedType, settings_.nbitsseed()).str();
1057       const string valid("1");
1058       streamTrackRaw.emplace_back(valid + seed + rinv + phi0 + z0 + t);
1059 
1060       // convert projected stubs
1061       unsigned int ihit(0);
1062       for (unsigned int ilayer = 0; ilayer < N_LAYER + N_DISK; ilayer++) {
1063         if (bestTracklet->match(ilayer)) {
1064           const Residual& resid = bestTracklet->resid(ilayer);
1065           // create bit accurate 64 bit word
1066           string r = resid.stubptr()->r().str();
1067           const string& phi = resid.fpgaphiresid().str();
1068           const string& rz = resid.fpgarzresid().str();
1069           const L1TStub* stub = resid.stubptr()->l1tstub();
1070           static constexpr int widthDisk2Sidentifier = 8;
1071           bool disk2S = (stub->disk() != 0) && (stub->isPSmodule() == 0);
1072           if (disk2S)
1073             r = string(widthDisk2Sidentifier, '0') + r;
1074           const string& stubId = resid.fpgastubid().str();
1075           // store seed, L1TStub, and bit accurate 64 bit word in clock accurate output
1076           streamsStubRaw[ihit++].emplace_back(seedType, *stub, valid + stubId + r + phi + rz);
1077         }
1078       }
1079       // convert seed stubs
1080       const string& stubId0 = bestTracklet->innerFPGAStub()->stubindex().str();
1081       const L1TStub* stub0 = bestTracklet->innerFPGAStub()->l1tstub();
1082       streamsStubRaw[ihit++].emplace_back(seedType, *stub0, valid + stubId0);
1083       const string& stubId1 = bestTracklet->outerFPGAStub()->stubindex().str();
1084       const L1TStub* stub1 = bestTracklet->outerFPGAStub()->l1tstub();
1085       streamsStubRaw[ihit++].emplace_back(seedType, *stub1, valid + stubId1);
1086       // fill all layers that have no stubs with gaps
1087       while (ihit < streamsStubRaw.size()) {
1088         streamsStubRaw[ihit++].emplace_back();
1089       }
1090     }
1091 
1092   } while (bestTracklet != nullptr && countAll < settings_.maxStep("TB"));
1093 
1094   if (settings_.writeMonitorData("FT")) {
1095     globals_->ofstream("fittrack.txt") << getName() << " " << countAll << " " << countFit << endl;
1096   }
1097 }