TrackFindingTMTT/plugins/TMTrackProducer.cc

0001 #include "L1Trigger/TrackFindingTMTT/plugins/TMTrackProducer.h"
0002 #include "L1Trigger/TrackFindingTMTT/interface/InputData.h"
0003 #include "L1Trigger/TrackFindingTMTT/interface/Sector.h"
0004 #include "L1Trigger/TrackFindingTMTT/interface/HTrphi.h"
0005 #include "L1Trigger/TrackFindingTMTT/interface/Make3Dtracks.h"
0006 #include "L1Trigger/TrackFindingTMTT/interface/DupFitTrkKiller.h"
0007 #include "L1Trigger/TrackFindingTMTT/interface/TrackFitFactory.h"
0008 #include "L1Trigger/TrackFindingTMTT/interface/L1fittedTrack.h"
0009 #include "L1Trigger/TrackFindingTMTT/interface/ConverterToTTTrack.h"
0010 #include "L1Trigger/TrackFindingTMTT/interface/HTcell.h"
0011 #include "L1Trigger/TrackFindingTMTT/interface/MuxHToutputs.h"
0012 #include "L1Trigger/TrackFindingTMTT/interface/MiniHTstage.h"
0013 #include "L1Trigger/TrackFindingTMTT/interface/Array2D.h"
0014 #include "L1Trigger/TrackFindingTMTT/interface/PrintL1trk.h"
0015
0016 #include "FWCore/Framework/interface/ESHandle.h"
0017
0018 #include <iostream>
0019 #include <vector>
0020 #include <list>
0021 #include <set>
0022 #include <sstream>
0023 #include <mutex>
0024
0025 using namespace std;
0026
0027 namespace tmtt {
0028
0029   namespace {
0030     std::once_flag printOnce;
0031   }  // namespace
0032
0033   std::unique_ptr<GlobalCacheTMTT> TMTrackProducer::initializeGlobalCache(edm::ParameterSet const& iConfig) {
0034     return std::make_unique<GlobalCacheTMTT>(iConfig);
0035   }
0036
0037   TMTrackProducer::TMTrackProducer(const edm::ParameterSet& iConfig, GlobalCacheTMTT const* globalCacheTMTT)
0038       : settings_(iConfig),                                        // Set configuration parameters
0039         stubWindowSuggest_(globalCacheTMTT->stubWindowSuggest()),  // For tuning FE stub window sizes
0040         hists_(globalCacheTMTT->hists()),                          // Initialize histograms
0041         htRphiErrMon_(globalCacheTMTT->htRphiErrMon()),            // rphi HT error monitoring
0042         debug_(true)                                               // Debug printout
0043   {
0044     using namespace edm;
0045
0046     // Get tokens for ES data access.
0047     magneticFieldToken_ =
0048         esConsumes<MagneticField, IdealMagneticFieldRecord, Transition::BeginRun>(settings_.magneticFieldInputTag());
0049     trackerGeometryToken_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord, Transition::BeginRun>(
0050         settings_.trackerGeometryInputTag());
0051     trackerTopologyToken_ =
0052         esConsumes<TrackerTopology, TrackerTopologyRcd, Transition::BeginRun>(settings_.trackerTopologyInputTag());
0053     ttStubAlgoToken_ =
0054         esConsumes<StubAlgorithm, TTStubAlgorithmRecord, Transition::BeginRun>(settings_.ttStubAlgoInputTag());
0055
0056     // Get tokens for ED data access.
0057     stubToken_ = consumes<TTStubDetSetVec>(settings_.stubInputTag());
0058     if (settings_.enableMCtruth()) {
0059       // These lines use lots of CPU, even if no use of truth info is made later.
0060       tpToken_ = consumes<TrackingParticleCollection>(settings_.tpInputTag());
0061       stubTruthToken_ = consumes<TTStubAssMap>(settings_.stubTruthInputTag());
0062       clusterTruthToken_ = consumes<TTClusterAssMap>(settings_.clusterTruthInputTag());
0063       genJetToken_ = consumes<reco::GenJetCollection>(settings_.genJetInputTag());
0064     }
0065
0066     trackFitters_ = settings_.trackFitters();
0067     useRZfilter_ = settings_.useRZfilter();
0068     runRZfilter_ = (not useRZfilter_.empty());  // Do any fitters require an r-z track filter to be run?
0069
0070     // Book histograms.
0071     //hists_.book();
0072
0073     // Create track fitting algorithm
0074     for (const string& fitterName : trackFitters_) {
0075       fitterWorkerMap_[fitterName] = trackFitFactory::create(fitterName, &settings_);
0076     }
0077
0078     //--- Define EDM output to be written to file (if required)
0079
0080     if (settings_.enableOutputIntermediateTTTracks()) {
0081       // L1 tracks found by Hough Transform
0082       produces<TTTrackCollection>("TML1TracksHT").setBranchAlias("TML1TracksHT");
0083       // L1 tracks found by r-z track filter.
0084       if (runRZfilter_)
0085         produces<TTTrackCollection>("TML1TracksRZ").setBranchAlias("TML1TracksRZ");
0086     }
0087     // L1 tracks after track fit by each of the fitting algorithms under study
0088     for (const string& fitterName : trackFitters_) {
0089       string edmName = string("TML1Tracks") + fitterName;
0090       produces<TTTrackCollection>(edmName).setBranchAlias(edmName);
0091     }
0092   }
0093
0094   //=== Run every run
0095
0096   void TMTrackProducer::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0097     // Get the B-field and store its value in the Settings class.
0098     const MagneticField* theMagneticField = &(iSetup.getData(magneticFieldToken_));
0099     float bField = theMagneticField->inTesla(GlobalPoint(0, 0, 0)).z();  // B field in Tesla.
0100     settings_.setMagneticField(bField);
0101
0102     // Set also B field in GlobalCacheTMTT (used only for Histogramming)
0103     globalCache()->settings().setMagneticField(bField);
0104
0105     std::stringstream text;
0106     text << "\n--- B field = " << bField << " Tesla ---\n";
0107     std::call_once(printOnce, [](string t) { PrintL1trk() << t; }, text.str());
0108
0109     // Get tracker geometry
0110     trackerGeometry_ = &(iSetup.getData(trackerGeometryToken_));
0111     trackerTopology_ = &(iSetup.getData(trackerTopologyToken_));
0112
0113     // Loop over tracker modules to get module info.
0114
0115     // Identifies tracker module type for firmware.
0116     TrackerModule::ModuleTypeCfg moduleTypeCfg;
0117     moduleTypeCfg.pitchVsType = settings_.pitchVsType();
0118     moduleTypeCfg.spaceVsType = settings_.spaceVsType();
0119     moduleTypeCfg.barrelVsType = settings_.barrelVsType();
0120     moduleTypeCfg.psVsType = settings_.psVsType();
0121     moduleTypeCfg.tiltedVsType = settings_.tiltedVsType();
0122
0123     listTrackerModule_.clear();
0124     for (const GeomDet* gd : trackerGeometry_->dets()) {
0125       DetId detId = gd->geographicalId();
0126       // Phase 2 Outer Tracker uses TOB for entire barrel & TID for entire endcap.
0127       if (detId.subdetId() != StripSubdetector::TOB && detId.subdetId() != StripSubdetector::TID)
0128         continue;
0129       if (trackerTopology_->isLower(detId)) {  // Select only lower of the two sensors in a module.
0130         // Store info about this tracker module.
0131         listTrackerModule_.emplace_back(trackerGeometry_, trackerTopology_, moduleTypeCfg, detId);
0132       }
0133     }
0134
0135     // Takes one copy of this to GlobalCacheTMTT for later histogramming.
0136     globalCache()->setListTrackerModule(listTrackerModule_);
0137
0138     // Get TTStubProducerAlgorithm algorithm, to adjust stub bend FE encoding.
0139     stubAlgo_ = dynamic_cast<const StubAlgorithmOfficial*>(&iSetup.getData(ttStubAlgoToken_));
0140     // Get FE stub window size from TTStub producer configuration
0141     const edm::ESHandle<StubAlgorithm> stubAlgoHandle = iSetup.getHandle(ttStubAlgoToken_);
0142     const edm::ParameterSet& pSetStubAlgo = getParameterSet(stubAlgoHandle.description()->pid_);
0143     stubFEWindows_ = std::make_unique<StubFEWindows>(pSetStubAlgo);
0144     // Initialize utilities needing FE window size.
0145     stubWindowSuggest_.setFEWindows(stubFEWindows_.get());
0146     degradeBend_ = std::make_unique<DegradeBend>(trackerTopology_, stubFEWindows_.get(), stubAlgo_);
0147   }
0148
0149   //=== Run every event
0150
0151   void TMTrackProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0152     // Note useful info about MC truth particles and about reconstructed stubs .
0153     InputData inputData(iEvent,
0154                         iSetup,
0155                         &settings_,
0156                         &stubWindowSuggest_,
0157                         degradeBend_.get(),
0158                         trackerGeometry_,
0159                         trackerTopology_,
0160                         listTrackerModule_,
0161                         tpToken_,
0162                         stubToken_,
0163                         stubTruthToken_,
0164                         clusterTruthToken_,
0165                         genJetToken_);
0166
0167     const list<TP>& vTPs = inputData.getTPs();
0168     const list<Stub*>& vStubs = inputData.stubs();
0169
0170     // Creates matrix of Sector objects, which decide which stubs are in which (eta,phi) sector
0171     Array2D<unique_ptr<Sector>> mSectors(settings_.numPhiSectors(), settings_.numEtaRegions());
0172     // Create matrix of r-phi Hough-Transform arrays, with one-to-one correspondence to sectors.
0173     Array2D<unique_ptr<HTrphi>> mHtRphis(settings_.numPhiSectors(), settings_.numEtaRegions());
0174     // Create matrix of Make3Dtracks objects, to run optional r-z track filter, with one-to-one correspondence to sectors.
0175     Array2D<unique_ptr<Make3Dtracks>> mMake3Dtrks(settings_.numPhiSectors(), settings_.numEtaRegions());
0176     // Create matrix of tracks from each fitter in each sector
0177     Array2D<map<string, std::list<L1fittedTrack>>> mapmFitTrks(settings_.numPhiSectors(), settings_.numEtaRegions());
0178     // Final tracks after duplicate removal from each track fitter in entire tracker.
0179     map<string, list<const L1fittedTrack*>> mapFinalTracks;
0180
0181     //=== Initialization
0182     // Create utility for converting L1 tracks from our private format to official CMSSW EDM format.
0183     const ConverterToTTTrack converter(&settings_);
0184
0185     // Pointers to TTTrack collections for ED output.
0186     auto htTTTracksForOutput = std::make_unique<TTTrackCollection>();
0187     auto rzTTTracksForOutput = std::make_unique<TTTrackCollection>();
0188     map<string, unique_ptr<TTTrackCollection>> allFitTTTracksForOutput;
0189     for (const string& fitterName : trackFitters_) {
0190       auto fitTTTracksForOutput = std::make_unique<TTTrackCollection>();
0191       allFitTTTracksForOutput[fitterName] = std::move(fitTTTracksForOutput);
0192     }
0193
0194     //=== Do tracking in the r-phi Hough transform within each sector.
0195
0196     // Fill Hough-Transform arrays with stubs.
0197     for (unsigned int iPhiSec = 0; iPhiSec < settings_.numPhiSectors(); iPhiSec++) {
0198       for (unsigned int iEtaReg = 0; iEtaReg < settings_.numEtaRegions(); iEtaReg++) {
0199         // Initialize constants for this sector.
0200         mSectors(iPhiSec, iEtaReg) = std::make_unique<Sector>(&settings_, iPhiSec, iEtaReg);
0201         Sector* sector = mSectors(iPhiSec, iEtaReg).get();
0202
0203         mHtRphis(iPhiSec, iEtaReg) = std::make_unique<HTrphi>(
0204             &settings_, iPhiSec, iEtaReg, sector->etaMin(), sector->etaMax(), sector->phiCentre(), &htRphiErrMon_);
0205         HTrphi* htRphi = mHtRphis(iPhiSec, iEtaReg).get();
0206
0207         // Check sector is enabled (always true, except if user disabled some for special studies).
0208         if (settings_.isHTRPhiEtaRegWhitelisted(iEtaReg)) {
0209           for (Stub* stub : vStubs) {
0210             // Digitize stub as would be at input to GP. This doesn't need the nonant number, since we assumed an integer number of
0211             // phi digitisation  bins inside an nonant. N.B. This changes the coordinates & bend stored in the stub.
0212
0213             if (settings_.enableDigitize())
0214               stub->digitize(iPhiSec, Stub::DigiStage::GP);
0215
0216             // Check if stub is inside this sector
0217             bool inside = sector->inside(stub);
0218
0219             if (inside) {
0220               // Check which eta subsectors within the sector the stub is compatible with (if subsectors being used).
0221               const vector<bool> inEtaSubSecs = sector->insideEtaSubSecs(stub);
0222
0223               // Digitize stub if as would be at input to HT, which slightly degrades its coord. & bend resolution, affecting the HT performance.
0224               if (settings_.enableDigitize())
0225                 stub->digitize(iPhiSec, Stub::DigiStage::HT);
0226
0227               // Store stub in Hough transform array for this sector, indicating its compatibility with eta subsectors with sector.
0228               htRphi->store(stub, inEtaSubSecs);
0229             }
0230           }
0231         }
0232
0233         // Find tracks in r-phi HT array.
0234         htRphi->end();  // Calls htArrayRphi_.end() -> HTBase::end()
0235       }
0236     }
0237
0238     if (settings_.muxOutputsHT() > 0) {
0239       // Multiplex outputs of several HT onto one pair of output opto-links.
0240       // This only affects tracking performance if option busySectorKill is enabled, so that tracks that
0241       // can't be sent down the link within the time-multiplexed period are killed.
0242       MuxHToutputs muxHT(&settings_);
0243       muxHT.exec(mHtRphis);
0244     }
0245
0246     // Optionally, run 2nd stage mini HT -- WITHOUT TRUNCATION ???
0247     if (settings_.miniHTstage()) {
0248       MiniHTstage miniHTstage(&settings_);
0249       miniHTstage.exec(mHtRphis);
0250     }
0251
0252     //=== Make 3D tracks, optionally running r-z track filters (such as Seed Filter) & duplicate track removal.
0253
0254     for (unsigned int iPhiSec = 0; iPhiSec < settings_.numPhiSectors(); iPhiSec++) {
0255       for (unsigned int iEtaReg = 0; iEtaReg < settings_.numEtaRegions(); iEtaReg++) {
0256         const Sector* sector = mSectors(iPhiSec, iEtaReg).get();
0257
0258         // Get tracks found by r-phi HT.
0259         const HTrphi* htRphi = mHtRphis(iPhiSec, iEtaReg).get();
0260         const list<L1track2D>& vecTracksRphi = htRphi->trackCands2D();
0261
0262         // Initialize utility for making 3D tracks from 2D ones.
0263         mMake3Dtrks(iPhiSec, iEtaReg) = std::make_unique<Make3Dtracks>(
0264             &settings_, iPhiSec, iEtaReg, sector->etaMin(), sector->etaMax(), sector->phiCentre());
0265         Make3Dtracks* make3Dtrk = mMake3Dtrks(iPhiSec, iEtaReg).get();
0266
0267         // Convert 2D tracks found by HT to 3D tracks (optionally by running r-z filters & duplicate track removal)
0268         make3Dtrk->run(vecTracksRphi);
0269
0270         if (settings_.enableOutputIntermediateTTTracks()) {
0271           // Convert these tracks to EDM format for output (used for collaborative work outside TMTT group).
0272           // Do this for tracks output by HT & optionally also for those output by r-z track filter.
0273           const list<L1track3D>& vecTrk3D_ht = make3Dtrk->trackCands3D(false);
0274           for (const L1track3D& trk : vecTrk3D_ht) {
0275             TTTrack<Ref_Phase2TrackerDigi_> htTTTrack = converter.makeTTTrack(&trk, iPhiSec, iEtaReg);
0276             htTTTracksForOutput->push_back(htTTTrack);
0277           }
0278
0279           if (runRZfilter_) {
0280             const list<L1track3D>& vecTrk3D_rz = make3Dtrk->trackCands3D(true);
0281             for (const L1track3D& trk : vecTrk3D_rz) {
0282               TTTrack<Ref_Phase2TrackerDigi_> rzTTTrack = converter.makeTTTrack(&trk, iPhiSec, iEtaReg);
0283               rzTTTracksForOutput->push_back(rzTTTrack);
0284             }
0285           }
0286         }
0287       }
0288     }
0289
0290     //=== Do a helix fit to all the track candidates.
0291
0292     // Loop over all the fitting algorithms we are trying.
0293     for (const string& fitterName : trackFitters_) {
0294       for (unsigned int iPhiSec = 0; iPhiSec < settings_.numPhiSectors(); iPhiSec++) {
0295         for (unsigned int iEtaReg = 0; iEtaReg < settings_.numEtaRegions(); iEtaReg++) {
0296           const Make3Dtracks* make3Dtrk = mMake3Dtrks(iPhiSec, iEtaReg).get();
0297
0298           // Does this fitter require r-z track filter to be run before it?
0299           bool useRZfilt = (std::count(useRZfilter_.begin(), useRZfilter_.end(), fitterName) > 0);
0300
0301           // Get 3D track candidates found by Hough transform (plus optional r-z filters/duplicate removal) in this sector.
0302           const list<L1track3D>& vecTrk3D = make3Dtrk->trackCands3D(useRZfilt);
0303
0304           // Find list where fitted tracks will be stored.
0305           list<L1fittedTrack>& fitTrksInSec = mapmFitTrks(iPhiSec, iEtaReg)[fitterName];
0306
0307           // Fit all tracks in this sector
0308           for (const L1track3D& trk : vecTrk3D) {
0309             // Ensure stubs assigned to this track is digitized with respect to the phi sector the track is in.
0310             if (settings_.enableDigitize()) {
0311               const vector<Stub*>& stubsOnTrk = trk.stubs();
0312               for (Stub* s : stubsOnTrk) {
0313                 // Also digitize stub in way this specific track fitter uses it.
0314                 s->digitize(iPhiSec, Stub::DigiStage::TF);
0315               }
0316             }
0317
0318             L1fittedTrack fitTrk = fitterWorkerMap_[fitterName]->fit(trk);
0319
0320             if (fitTrk.accepted()) {  // If fitter accepted track, then store it.
0321               // Optionally digitize fitted track, degrading slightly resolution.
0322               if (settings_.enableDigitize())
0323                 fitTrk.digitizeTrack(fitterName);
0324               // Store fitted tracks, such that there is one fittedTracks corresponding to each HT tracks.
0325               fitTrksInSec.push_back(fitTrk);
0326             }
0327           }
0328         }
0329       }
0330     }
0331
0332     // Run duplicate track removal on the fitted tracks if requested.
0333
0334     // Initialize the duplicate track removal algorithm that can optionally be run after the track fit.
0335     DupFitTrkKiller killDupFitTrks(&settings_);
0336
0337     // Loop over all the fitting algorithms we used.
0338     for (const string& fitterName : trackFitters_) {
0339       for (unsigned int iPhiSec = 0; iPhiSec < settings_.numPhiSectors(); iPhiSec++) {
0340         for (unsigned int iEtaReg = 0; iEtaReg < settings_.numEtaRegions(); iEtaReg++) {
0341           // Get fitted tracks in sector
0342           const list<L1fittedTrack>& fitTrksInSec = mapmFitTrks(iPhiSec, iEtaReg)[fitterName];
0343
0344           // Run duplicate removal
0345           list<const L1fittedTrack*> filteredFitTrksInSec = killDupFitTrks.filter(fitTrksInSec);
0346
0347           // Prepare TTTrack collection.
0348           for (const L1fittedTrack* fitTrk : filteredFitTrksInSec) {
0349             // Convert these fitted tracks to EDM format for output (used for collaborative work outside TMTT group).
0350             TTTrack<Ref_Phase2TrackerDigi_> fitTTTrack = converter.makeTTTrack(fitTrk, iPhiSec, iEtaReg);
0351             allFitTTTracksForOutput[fitterName]->push_back(fitTTTrack);
0352           }
0353
0354           // Store fitted tracks from entire tracker.
0355           mapFinalTracks[fitterName].insert(
0356               mapFinalTracks[fitterName].end(), filteredFitTrksInSec.begin(), filteredFitTrksInSec.end());
0357         }
0358       }
0359     }
0360
0361     // Debug printout
0362     if (debug_) {
0363       PrintL1trk() << "INPUT #TPs = " << vTPs.size() << " #STUBs = " << vStubs.size();
0364       unsigned int numHTtracks = 0;
0365       for (unsigned int iPhiSec = 0; iPhiSec < settings_.numPhiSectors(); iPhiSec++) {
0366         for (unsigned int iEtaReg = 0; iEtaReg < settings_.numEtaRegions(); iEtaReg++) {
0367           const Make3Dtracks* make3Dtrk = mMake3Dtrks(iPhiSec, iEtaReg).get();
0368           numHTtracks += make3Dtrk->trackCands3D(false).size();
0369         }
0370       }
0371       PrintL1trk() << "Number of tracks after HT = " << numHTtracks;
0372       for (const auto& p : mapFinalTracks) {
0373         const string& fitName = p.first;
0374         const list<const L1fittedTrack*> fittedTracks = p.second;
0375         PrintL1trk() << "Number of tracks after " << fitName << " track helix fit = " << fittedTracks.size();
0376       }
0377     }
0378
0379     // Allow histogramming to plot undigitized variables.
0380     for (Stub* stub : vStubs) {
0381       if (settings_.enableDigitize())
0382         stub->setDigitizeWarningsOn(false);
0383     }
0384
0385     // Fill histograms to monitor input data & tracking performance.
0386     hists_.fill(inputData, mSectors, mHtRphis, mMake3Dtrks, mapFinalTracks);
0387
0388     //=== Store output EDM track and hardware stub collections.
0389     if (settings_.enableOutputIntermediateTTTracks()) {
0390       iEvent.put(std::move(htTTTracksForOutput), "TML1TracksHT");
0391       if (runRZfilter_)
0392         iEvent.put(std::move(rzTTTracksForOutput), "TML1TracksRZ");
0393     }
0394     for (const string& fitterName : trackFitters_) {
0395       string edmName = string("TML1Tracks") + fitterName;
0396       iEvent.put(std::move(allFitTTTracksForOutput[fitterName]), edmName);
0397     }
0398   }
0399
0400   void TMTrackProducer::globalEndJob(GlobalCacheTMTT* globalCacheTMTT) {
0401     const Settings& settings = globalCacheTMTT->settings();
0402
0403     // Print stub window sizes that TMTT recommends CMS uses in FE chips.
0404     if (settings.printStubWindows())
0405       globalCacheTMTT->stubWindowSuggest().printResults();
0406
0407     // Print (once) info about tracker geometry.
0408     globalCacheTMTT->hists().trackerGeometryAnalysis(globalCacheTMTT->listTrackerModule());
0409
0410     PrintL1trk() << "\n Number of (eta,phi) sectors used = (" << settings.numEtaRegions() << ","
0411                  << settings.numPhiSectors() << ")";
0412
0413     // Print job summary
0414     globalCacheTMTT->hists().endJobAnalysis(&(globalCacheTMTT->htRphiErrMon()));
0415   }
0416
0417 }  // namespace tmtt
0418
0419 DEFINE_FWK_MODULE(tmtt::TMTrackProducer);