TrackerTFP/test/AnalyzerDR.cc

0001 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0002 #include "FWCore/Framework/interface/Run.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/ServiceRegistry/interface/Service.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/Utilities/interface/EDGetToken.h"
0010 #include "FWCore/Utilities/interface/InputTag.h"
0011 #include "FWCore/Utilities/interface/Exception.h"
0012 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0013 #include "DataFormats/Common/interface/Handle.h"
0014
0015 #include "SimTracker/TrackTriggerAssociation/interface/StubAssociation.h"
0016 #include "L1Trigger/TrackTrigger/interface/Setup.h"
0017 #include "L1Trigger/TrackerTFP/interface/DataFormats.h"
0018
0019 #include <TProfile.h>
0020 #include <TH1F.h>
0021
0022 #include <vector>
0023 #include <deque>
0024 #include <set>
0025 #include <cmath>
0026 #include <numeric>
0027 #include <sstream>
0028
0029 namespace trackerTFP {
0030
0031   /*! \class  trackerTFP::AnalyzerDR
0032    *  \brief  Class to analyze hardware like structured track Collection generated by Duplicate Removal
0033    *  \author Thomas Schuh
0034    *  \date   2023, Feb
0035    */
0036   class AnalyzerDR : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0037   public:
0038     AnalyzerDR(const edm::ParameterSet& iConfig);
0039     void beginJob() override {}
0040     void beginRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) override;
0041     void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0042     void endRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) override {}
0043     void endJob() override;
0044
0045   private:
0046     //
0047     void formTracks(const tt::StreamsTrack& streamsTrack,
0048                     const tt::StreamsStub& streamsStubs,
0049                     std::vector<std::vector<TTStubRef>>& tracks,
0050                     int channel) const;
0051     //
0052     void associate(const std::vector<std::vector<TTStubRef>>& tracks,
0053                    const tt::StubAssociation* ass,
0054                    std::set<TPPtr>& tps,
0055                    int& sum,
0056                    bool perfect = true) const;
0057     // ED input token of stubs
0058     edm::EDGetTokenT<tt::StreamsStub> edGetTokenStubs_;
0059     // ED input token of tracks
0060     edm::EDGetTokenT<tt::StreamsTrack> edGetTokenTracks_;
0061     // ED input token of TTStubRef to TPPtr association for tracking efficiency
0062     edm::EDGetTokenT<tt::StubAssociation> edGetTokenSelection_;
0063     // ED input token of TTStubRef to recontructable TPPtr association
0064     edm::EDGetTokenT<tt::StubAssociation> edGetTokenReconstructable_;
0065     // Setup token
0066     edm::ESGetToken<tt::Setup, tt::SetupRcd> esGetTokenSetup_;
0067     // DataFormats token
0068     edm::ESGetToken<DataFormats, DataFormatsRcd> esGetTokenDataFormats_;
0069     // stores, calculates and provides run-time constants
0070     const tt::Setup* setup_ = nullptr;
0071     // helper class to extract structured data from tt::Frames
0072     const DataFormats* dataFormats_ = nullptr;
0073     // enables analyze of TPs
0074     bool useMCTruth_;
0075     //
0076     int nEvents_ = 0;
0077
0078     // Histograms
0079
0080     TProfile* prof_;
0081     TProfile* profChannel_;
0082     TProfile* profTracks_;
0083     TH1F* hisChannel_;
0084     TH1F* hisTracks_;
0085
0086     // printout
0087     std::stringstream log_;
0088   };
0089
0090   AnalyzerDR::AnalyzerDR(const edm::ParameterSet& iConfig) : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")) {
0091     usesResource("TFileService");
0092     // book in- and output ED products
0093     const std::string& label = iConfig.getParameter<std::string>("OutputLabelDR");
0094     const std::string& branchStubs = iConfig.getParameter<std::string>("BranchStubs");
0095     const std::string& branchTracks = iConfig.getParameter<std::string>("BranchTracks");
0096     edGetTokenStubs_ = consumes<tt::StreamsStub>(edm::InputTag(label, branchStubs));
0097     edGetTokenTracks_ = consumes<tt::StreamsTrack>(edm::InputTag(label, branchTracks));
0098     if (useMCTruth_) {
0099       const auto& inputTagSelecttion = iConfig.getParameter<edm::InputTag>("InputTagSelection");
0100       const auto& inputTagReconstructable = iConfig.getParameter<edm::InputTag>("InputTagReconstructable");
0101       edGetTokenSelection_ = consumes<tt::StubAssociation>(inputTagSelecttion);
0102       edGetTokenReconstructable_ = consumes<tt::StubAssociation>(inputTagReconstructable);
0103     }
0104     // book ES products
0105     esGetTokenSetup_ = esConsumes<edm::Transition::BeginRun>();
0106     esGetTokenDataFormats_ = esConsumes<edm::Transition::BeginRun>();
0107     // log config
0108     log_.setf(std::ios::fixed, std::ios::floatfield);
0109     log_.precision(4);
0110   }
0111
0112   void AnalyzerDR::beginRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) {
0113     // helper class to store configurations
0114     setup_ = &iSetup.getData(esGetTokenSetup_);
0115     // helper class to extract structured data from tt::Frames
0116     dataFormats_ = &iSetup.getData(esGetTokenDataFormats_);
0117     // book histograms
0118     edm::Service<TFileService> fs;
0119     TFileDirectory dir;
0120     dir = fs->mkdir("DR");
0121     prof_ = dir.make<TProfile>("Counts", ";", 12, 0.5, 12.5);
0122     prof_->GetXaxis()->SetBinLabel(1, "Stubs");
0123     prof_->GetXaxis()->SetBinLabel(2, "Tracks");
0124     prof_->GetXaxis()->SetBinLabel(4, "Matched Tracks");
0125     prof_->GetXaxis()->SetBinLabel(5, "All Tracks");
0126     prof_->GetXaxis()->SetBinLabel(6, "Found TPs");
0127     prof_->GetXaxis()->SetBinLabel(7, "Found selected TPs");
0128     prof_->GetXaxis()->SetBinLabel(9, "All TPs");
0129     prof_->GetXaxis()->SetBinLabel(10, "states");
0130     prof_->GetXaxis()->SetBinLabel(12, "max tp");
0131     // channel occupancy
0132     constexpr int maxOcc = 180;
0133     const int numChannels = dataFormats_->numChannel(Process::dr);
0134     hisChannel_ = dir.make<TH1F>("His Channel Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0135     profChannel_ = dir.make<TProfile>("Prof Channel Occupancy", ";", numChannels, -.5, numChannels - .5);
0136     // track occupancy
0137     hisTracks_ = dir.make<TH1F>("His Track Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0138     profTracks_ = dir.make<TProfile>("Prof Track Occupancy", ";", numChannels, -.5, numChannels - .5);
0139   }
0140
0141   void AnalyzerDR::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0142     // read in ht products
0143     edm::Handle<tt::StreamsStub> handleStubs;
0144     iEvent.getByToken<tt::StreamsStub>(edGetTokenStubs_, handleStubs);
0145     const tt::StreamsStub& acceptedStubs = *handleStubs;
0146     edm::Handle<tt::StreamsTrack> handleTracks;
0147     iEvent.getByToken<tt::StreamsTrack>(edGetTokenTracks_, handleTracks);
0148     const tt::StreamsTrack& acceptedTracks = *handleTracks;
0149     // read in MCTruth
0150     const tt::StubAssociation* selection = nullptr;
0151     const tt::StubAssociation* reconstructable = nullptr;
0152     if (useMCTruth_) {
0153       edm::Handle<tt::StubAssociation> handleSelection;
0154       iEvent.getByToken<tt::StubAssociation>(edGetTokenSelection_, handleSelection);
0155       selection = handleSelection.product();
0156       prof_->Fill(9, selection->numTPs());
0157       edm::Handle<tt::StubAssociation> handleReconstructable;
0158       iEvent.getByToken<tt::StubAssociation>(edGetTokenReconstructable_, handleReconstructable);
0159       reconstructable = handleReconstructable.product();
0160     }
0161     // analyze ht products and associate found tracks with reconstrucable TrackingParticles
0162     std::set<TPPtr> tpPtrs;
0163     std::set<TPPtr> tpPtrsSelection;
0164     std::set<TPPtr> tpPtrsMax;
0165     int allMatched(0);
0166     int allTracks(0);
0167     for (int region = 0; region < setup_->numRegions(); region++) {
0168       const int offset = region * dataFormats_->numChannel(Process::dr);
0169       int nStubs(0);
0170       int nTracks(0);
0171       for (int channel = 0; channel < dataFormats_->numChannel(Process::dr); channel++) {
0172         std::vector<std::vector<TTStubRef>> tracks;
0173         formTracks(acceptedTracks, acceptedStubs, tracks, offset + channel);
0174         hisTracks_->Fill(tracks.size());
0175         profTracks_->Fill(channel, tracks.size());
0176         nTracks += tracks.size();
0177         nStubs += std::accumulate(tracks.begin(), tracks.end(), 0, [](int sum, const std::vector<TTStubRef>& track) {
0178           return sum + track.size();
0179         });
0180         allTracks += tracks.size();
0181         if (!useMCTruth_)
0182           continue;
0183         int tmp(0);
0184         associate(tracks, selection, tpPtrsSelection, tmp);
0185         associate(tracks, reconstructable, tpPtrs, allMatched, false);
0186         associate(tracks, selection, tpPtrsMax, tmp, false);
0187         const int size = acceptedTracks[offset + channel].size();
0188         hisChannel_->Fill(size);
0189         profChannel_->Fill(channel, size);
0190       }
0191       prof_->Fill(1, nStubs);
0192       prof_->Fill(2, nTracks);
0193     }
0194     prof_->Fill(4, allMatched);
0195     prof_->Fill(5, allTracks);
0196     prof_->Fill(6, tpPtrs.size());
0197     prof_->Fill(7, tpPtrsSelection.size());
0198     prof_->Fill(12, tpPtrsMax.size());
0199     nEvents_++;
0200   }
0201
0202   void AnalyzerDR::endJob() {
0203     if (nEvents_ == 0)
0204       return;
0205     // printout DR summary
0206     const double totalTPs = prof_->GetBinContent(9);
0207     const double numStubs = prof_->GetBinContent(1);
0208     const double numTracks = prof_->GetBinContent(2);
0209     const double totalTracks = prof_->GetBinContent(5);
0210     const double numTracksMatched = prof_->GetBinContent(4);
0211     const double numTPsAll = prof_->GetBinContent(6);
0212     const double numTPsEff = prof_->GetBinContent(7);
0213     const double numTPsEffMax = prof_->GetBinContent(12);
0214     const double errStubs = prof_->GetBinError(1);
0215     const double errTracks = prof_->GetBinError(2);
0216     const double fracFake = (totalTracks - numTracksMatched) / totalTracks;
0217     const double fracDup = (numTracksMatched - numTPsAll) / totalTracks;
0218     const double eff = numTPsEff / totalTPs;
0219     const double errEff = std::sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0220     const double effMax = numTPsEffMax / totalTPs;
0221     const double errEffMax = std::sqrt(effMax * (1. - effMax) / totalTPs / nEvents_);
0222     const std::vector<double> nums = {numStubs, numTracks};
0223     const std::vector<double> errs = {errStubs, errTracks};
0224     const int wNums = std::ceil(std::log10(*std::max_element(nums.begin(), nums.end()))) + 5;
0225     const int wErrs = std::ceil(std::log10(*std::max_element(errs.begin(), errs.end()))) + 5;
0226     log_ << "                         DR  SUMMARY                         " << std::endl;
0227     log_ << "number of stubs       per TFP = " << std::setw(wNums) << numStubs << " +- " << std::setw(wErrs) << errStubs
0228          << std::endl;
0229     log_ << "number of tracks      per TFP = " << std::setw(wNums) << numTracks << " +- " << std::setw(wErrs)
0230          << errTracks << std::endl;
0231     log_ << "          tracking efficiency = " << std::setw(wNums) << eff << " +- " << std::setw(wErrs) << errEff
0232          << std::endl;
0233     log_ << "      max tracking efficiency = " << std::setw(wNums) << effMax << " +- " << std::setw(wErrs) << errEffMax
0234          << std::endl;
0235     log_ << "                    fake rate = " << std::setw(wNums) << fracFake << std::endl;
0236     log_ << "               duplicate rate = " << std::setw(wNums) << fracDup << std::endl;
0237     log_ << "=============================================================";
0238     edm::LogPrint(moduleDescription().moduleName()) << log_.str();
0239   }
0240
0241   //
0242   void AnalyzerDR::formTracks(const tt::StreamsTrack& streamsTrack,
0243                               const tt::StreamsStub& streamsStubs,
0244                               std::vector<std::vector<TTStubRef>>& tracks,
0245                               int channel) const {
0246     const int offset = channel * setup_->numLayers();
0247     const tt::StreamTrack& streamTrack = streamsTrack[channel];
0248     const int numTracks =
0249         std::accumulate(streamTrack.begin(), streamTrack.end(), 0, [](int sum, const tt::FrameTrack& frame) {
0250           return sum + (frame.first.isNonnull() ? 1 : 0);
0251         });
0252     tracks.reserve(numTracks);
0253     for (int frame = 0; frame < static_cast<int>(streamTrack.size()); frame++) {
0254       const tt::FrameTrack& frameTrack = streamTrack[frame];
0255       if (frameTrack.first.isNull())
0256         continue;
0257       std::deque<TTStubRef> stubs;
0258       for (int layer = 0; layer < setup_->numLayers(); layer++) {
0259         const tt::FrameStub& stub = streamsStubs[offset + layer][frame];
0260         if (stub.first.isNonnull())
0261           stubs.push_back(stub.first);
0262       }
0263       tracks.emplace_back(stubs.begin(), stubs.end());
0264     }
0265   }
0266
0267   //
0268   void AnalyzerDR::associate(const std::vector<std::vector<TTStubRef>>& tracks,
0269                              const tt::StubAssociation* ass,
0270                              std::set<TPPtr>& tps,
0271                              int& sum,
0272                              bool perfect) const {
0273     for (const std::vector<TTStubRef>& ttStubRefs : tracks) {
0274       const std::vector<TPPtr>& tpPtrs = perfect ? ass->associateFinal(ttStubRefs) : ass->associate(ttStubRefs);
0275       if (tpPtrs.empty())
0276         continue;
0277       sum++;
0278       std::copy(tpPtrs.begin(), tpPtrs.end(), std::inserter(tps, tps.begin()));
0279     }
0280   }
0281
0282 }  // namespace trackerTFP
0283
0284 DEFINE_FWK_MODULE(trackerTFP::AnalyzerDR);