TrackFindingTracklet/src/FitTrack.cc

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