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
0079
0080
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
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
0120 if (settings_.removalType() == "merge") {
0121 for (const auto& it : trackstublist) {
0122 int layer = it->layer().value() + 1;
0123 if (it->layer().value() < 0) {
0124 layer = it->disk().value() + 10 * it->disk().value() / abs(it->disk().value());
0125 }
0126 stubidslist.push_back(std::make_pair(layer, it->phiregionaddress()));
0127 }
0128
0129
0130
0131 } else {
0132
0133
0134 HybridFit hybridFitter(iSector_, settings_, globals_);
0135 hybridFitter.Fit(tracklet, trackstublist);
0136 }
0137 }
0138 }
0139
0140 #else
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
0159 int layers[N_LAYER];
0160 double r[N_LAYER];
0161 unsigned int nlayers = 0;
0162 int disks[N_DISK];
0163 double z[N_DISK];
0164 unsigned int ndisks = 0;
0165
0166
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;
0187 std::bitset<N_DISK * 2> dmatches;
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
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()) {
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
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;
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++) {
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
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
0751 ichisqfit = ichisqfit >> 7;
0752
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
0764
0765
0766
0767
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
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
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
0863
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
0875
0876
0877
0878
0879 void FitTrack::execute(deque<string>& streamTrackRaw,
0880 vector<deque<StubStreamData>>& streamsStubRaw,
0881 unsigned int iSector) {
0882
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
0951 int nMatches = 0;
0952
0953
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
1004 if ((bestTracklet->getISeed() >= (int)N_SEED_PROMPT && nMatchesUniq >= 1) ||
1005 nMatchesUniq >= 2) {
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
1041 if (settings_.storeTrackBuilderOutput() && bestTracklet) {
1042
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
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
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
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
1076 streamsStubRaw[ihit++].emplace_back(seedType, *stub, valid + stubId + r + phi + rz);
1077 }
1078 }
1079
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
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 }