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