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