Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-07-07 22:33:32 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!