TrackerTFP/test/AnalyzerMHT.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 #include <TEfficiency.h>
0022
0023 #include <vector>
0024 #include <deque>
0025 #include <set>
0026 #include <cmath>
0027 #include <numeric>
0028 #include <sstream>
0029
0030 using namespace std;
0031 using namespace edm;
0032 using namespace tt;
0033
0034 namespace trackerTFP {
0035
0036   /*! \class  trackerTFP::AnalyzerMHT
0037    *  \brief  Class to analyze hardware like structured TTStub Collection generated by Mini Hough Transform
0038    *  \author Thomas Schuh
0039    *  \date   2020, Apr
0040    */
0041   class AnalyzerMHT : public one::EDAnalyzer<one::WatchRuns, one::SharedResources> {
0042   public:
0043     AnalyzerMHT(const ParameterSet& iConfig);
0044     void beginJob() override {}
0045     void beginRun(const Run& iEvent, const EventSetup& iSetup) override;
0046     void analyze(const Event& iEvent, const EventSetup& iSetup) override;
0047     void endRun(const Run& iEvent, const EventSetup& iSetup) override {}
0048     void endJob() override;
0049
0050   private:
0051     //
0052     void formTracks(const StreamStub& stream, vector<vector<TTStubRef>>& tracks) const;
0053     //
0054     void associate(const vector<vector<TTStubRef>>& tracks, const StubAssociation* ass, set<TPPtr>& tps, int& sum) const;
0055
0056     // ED input token of stubs
0057     EDGetTokenT<StreamsStub> edGetTokenAccepted_;
0058     // ED input token of lost stubs
0059     EDGetTokenT<StreamsStub> edGetTokenLost_;
0060     // ED input token of TTStubRef to TPPtr association for tracking efficiency
0061     EDGetTokenT<StubAssociation> edGetTokenSelection_;
0062     // ED input token of TTStubRef to recontructable TPPtr association
0063     EDGetTokenT<StubAssociation> edGetTokenReconstructable_;
0064     // Setup token
0065     ESGetToken<Setup, SetupRcd> esGetTokenSetup_;
0066     // DataFormats token
0067     ESGetToken<DataFormats, DataFormatsRcd> esGetTokenDataFormats_;
0068     // stores, calculates and provides run-time constants
0069     const Setup* setup_ = nullptr;
0070     // helper class to extract structured data from tt::Frames
0071     const DataFormats* dataFormats_ = nullptr;
0072     // enables analyze of TPs
0073     bool useMCTruth_;
0074     //
0075     int nEvents_ = 0;
0076
0077     // Histograms
0078
0079     TProfile* prof_;
0080     TProfile* profChannel_;
0081     TH1F* hisChannel_;
0082     TH1F* hisEff_;
0083     TH1F* hisEffTotal_;
0084     TEfficiency* eff_;
0085
0086     // printout
0087     stringstream log_;
0088   };
0089
0090   AnalyzerMHT::AnalyzerMHT(const ParameterSet& iConfig) : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")) {
0091     usesResource("TFileService");
0092     // book in- and output ED products
0093     const string& label = iConfig.getParameter<string>("LabelMHT");
0094     const string& branchAccepted = iConfig.getParameter<string>("BranchAcceptedStubs");
0095     const string& branchLost = iConfig.getParameter<string>("BranchLostStubs");
0096     edGetTokenAccepted_ = consumes<StreamsStub>(InputTag(label, branchAccepted));
0097     edGetTokenLost_ = consumes<StreamsStub>(InputTag(label, branchLost));
0098     if (useMCTruth_) {
0099       const auto& inputTagSelecttion = iConfig.getParameter<InputTag>("InputTagSelection");
0100       const auto& inputTagReconstructable = iConfig.getParameter<InputTag>("InputTagReconstructable");
0101       edGetTokenSelection_ = consumes<StubAssociation>(inputTagSelecttion);
0102       edGetTokenReconstructable_ = consumes<StubAssociation>(inputTagReconstructable);
0103     }
0104     // book ES products
0105     esGetTokenSetup_ = esConsumes<Setup, SetupRcd, Transition::BeginRun>();
0106     esGetTokenDataFormats_ = esConsumes<DataFormats, DataFormatsRcd, Transition::BeginRun>();
0107     // log config
0108     log_.setf(ios::fixed, ios::floatfield);
0109     log_.precision(4);
0110   }
0111
0112   void AnalyzerMHT::beginRun(const Run& iEvent, const 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     Service<TFileService> fs;
0119     TFileDirectory dir;
0120     dir = fs->mkdir("MHT");
0121     prof_ = dir.make<TProfile>("Counts", ";", 9, 0.5, 9.5);
0122     prof_->GetXaxis()->SetBinLabel(1, "Stubs");
0123     prof_->GetXaxis()->SetBinLabel(2, "Tracks");
0124     prof_->GetXaxis()->SetBinLabel(3, "Lost Tracks");
0125     prof_->GetXaxis()->SetBinLabel(4, "Matched Tracks");
0126     prof_->GetXaxis()->SetBinLabel(5, "All Tracks");
0127     prof_->GetXaxis()->SetBinLabel(6, "Found TPs");
0128     prof_->GetXaxis()->SetBinLabel(7, "Found selected TPs");
0129     prof_->GetXaxis()->SetBinLabel(8, "Lost TPs");
0130     prof_->GetXaxis()->SetBinLabel(9, "All TPs");
0131     // channel occupancy
0132     constexpr int maxOcc = 180;
0133     const int numChannels = dataFormats_->numChannel(Process::mht);
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     // Efficiencies
0137     hisEffTotal_ = dir.make<TH1F>("HisTPEtaTotal", ";", 128, -2.5, 2.5);
0138     hisEff_ = dir.make<TH1F>("HisTPEta", ";", 128, -2.5, 2.5);
0139     eff_ = dir.make<TEfficiency>("EffEta", ";", 128, -2.5, 2.5);
0140   }
0141
0142   void AnalyzerMHT::analyze(const Event& iEvent, const EventSetup& iSetup) {
0143     auto fill = [](const TPPtr& tpPtr, TH1F* his) { his->Fill(tpPtr->eta()); };
0144     // read in ht products
0145     Handle<StreamsStub> handleAccepted;
0146     iEvent.getByToken<StreamsStub>(edGetTokenAccepted_, handleAccepted);
0147     Handle<StreamsStub> handleLost;
0148     iEvent.getByToken<StreamsStub>(edGetTokenLost_, handleLost);
0149     // read in MCTruth
0150     const StubAssociation* selection = nullptr;
0151     const StubAssociation* reconstructable = nullptr;
0152     if (useMCTruth_) {
0153       Handle<StubAssociation> handleSelection;
0154       iEvent.getByToken<StubAssociation>(edGetTokenSelection_, handleSelection);
0155       selection = handleSelection.product();
0156       prof_->Fill(9, selection->numTPs());
0157       Handle<StubAssociation> handleReconstructable;
0158       iEvent.getByToken<StubAssociation>(edGetTokenReconstructable_, handleReconstructable);
0159       reconstructable = handleReconstructable.product();
0160       for (const auto& p : selection->getTrackingParticleToTTStubsMap())
0161         fill(p.first, hisEffTotal_);
0162     }
0163     // analyze ht products and associate found tracks with reconstrucable TrackingParticles
0164     set<TPPtr> tpPtrs;
0165     set<TPPtr> tpPtrsSelection;
0166     set<TPPtr> tpPtrsLost;
0167     int allMatched(0);
0168     int allTracks(0);
0169     for (int region = 0; region < setup_->numRegions(); region++) {
0170       int nStubs(0);
0171       int nTracks(0);
0172       int nLost(0);
0173       for (int channel = 0; channel < dataFormats_->numChannel(Process::mht); channel++) {
0174         const int index = region * dataFormats_->numChannel(Process::mht) + channel;
0175         const StreamStub& accepted = handleAccepted->at(index);
0176         hisChannel_->Fill(accepted.size());
0177         profChannel_->Fill(channel, accepted.size());
0178         nStubs += accumulate(accepted.begin(), accepted.end(), 0, [](int sum, const FrameStub& frame) {
0179           return sum + (frame.first.isNonnull() ? 1 : 0);
0180         });
0181         vector<vector<TTStubRef>> tracks;
0182         vector<vector<TTStubRef>> lost;
0183         formTracks(accepted, tracks);
0184         formTracks(handleLost->at(index), lost);
0185         nTracks += tracks.size();
0186         allTracks += tracks.size();
0187         nLost += lost.size();
0188         if (!useMCTruth_)
0189           continue;
0190         int tmp(0);
0191         associate(tracks, selection, tpPtrsSelection, tmp);
0192         associate(lost, selection, tpPtrsLost, tmp);
0193         associate(tracks, reconstructable, tpPtrs, allMatched);
0194       }
0195       prof_->Fill(1, nStubs);
0196       prof_->Fill(2, nTracks);
0197       prof_->Fill(3, nLost);
0198     }
0199     for (const TPPtr& tpPtr : tpPtrsSelection)
0200       fill(tpPtr, hisEff_);
0201     vector<TPPtr> recovered;
0202     recovered.reserve(tpPtrsLost.size());
0203     set_intersection(tpPtrsLost.begin(), tpPtrsLost.end(), tpPtrs.begin(), tpPtrs.end(), back_inserter(recovered));
0204     for (const TPPtr& tpPtr : recovered)
0205       tpPtrsLost.erase(tpPtr);
0206     prof_->Fill(4, allMatched);
0207     prof_->Fill(5, allTracks);
0208     prof_->Fill(6, tpPtrs.size());
0209     prof_->Fill(7, tpPtrsSelection.size());
0210     prof_->Fill(8, tpPtrsLost.size());
0211     nEvents_++;
0212   }
0213
0214   void AnalyzerMHT::endJob() {
0215     if (nEvents_ == 0)
0216       return;
0217     // effi
0218     eff_->SetPassedHistogram(*hisEff_, "f");
0219     eff_->SetTotalHistogram(*hisEffTotal_, "f");
0220     // printout MHT summary
0221     const double totalTPs = prof_->GetBinContent(9);
0222     const double numStubs = prof_->GetBinContent(1);
0223     const double numTracks = prof_->GetBinContent(2);
0224     const double numTracksLost = prof_->GetBinContent(3);
0225     const double totalTracks = prof_->GetBinContent(5);
0226     const double numTracksMatched = prof_->GetBinContent(4);
0227     const double numTPsAll = prof_->GetBinContent(6);
0228     const double numTPsEff = prof_->GetBinContent(7);
0229     const double numTPsLost = prof_->GetBinContent(8);
0230     const double errStubs = prof_->GetBinError(1);
0231     const double errTracks = prof_->GetBinError(2);
0232     const double errTracksLost = prof_->GetBinError(3);
0233     const double fracFake = (totalTracks - numTracksMatched) / totalTracks;
0234     const double fracDup = (numTracksMatched - numTPsAll) / totalTracks;
0235     const double eff = numTPsEff / totalTPs;
0236     const double errEff = sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0237     const double effLoss = numTPsLost / totalTPs;
0238     const double errEffLoss = sqrt(effLoss * (1. - effLoss) / totalTPs / nEvents_);
0239     const vector<double> nums = {numStubs, numTracks, numTracksLost};
0240     const vector<double> errs = {errStubs, errTracks, errTracksLost};
0241     const int wNums = ceil(log10(*max_element(nums.begin(), nums.end()))) + 5;
0242     const int wErrs = ceil(log10(*max_element(errs.begin(), errs.end()))) + 5;
0243     log_ << "                         MHT SUMMARY                         " << endl;
0244     log_ << "number of stubs       per TFP = " << setw(wNums) << numStubs << " +- " << setw(wErrs) << errStubs << endl;
0245     log_ << "number of tracks      per TFP = " << setw(wNums) << numTracks << " +- " << setw(wErrs) << errTracks
0246          << endl;
0247     log_ << "number of lost tracks per TFP = " << setw(wNums) << numTracksLost << " +- " << setw(wErrs) << errTracksLost
0248          << endl;
0249     log_ << "     max  tracking efficiency = " << setw(wNums) << eff << " +- " << setw(wErrs) << errEff << endl;
0250     log_ << "     lost tracking efficiency = " << setw(wNums) << effLoss << " +- " << setw(wErrs) << errEffLoss << endl;
0251     log_ << "                    fake rate = " << setw(wNums) << fracFake << endl;
0252     log_ << "               duplicate rate = " << setw(wNums) << fracDup << endl;
0253     log_ << "=============================================================";
0254     LogPrint("L1Trigger/TrackerTFP") << log_.str();
0255   }
0256
0257   //
0258   void AnalyzerMHT::formTracks(const StreamStub& stream, vector<vector<TTStubRef>>& tracks) const {
0259     vector<StubMHT> stubs;
0260     stubs.reserve(stream.size());
0261     for (const FrameStub& frame : stream)
0262       if (frame.first.isNonnull())
0263         stubs.emplace_back(frame, dataFormats_);
0264     for (auto it = stubs.begin(); it != stubs.end();) {
0265       const auto start = it;
0266       const int id = it->trackId();
0267       auto different = [id](const StubMHT& stub) { return id != stub.trackId(); };
0268       it = find_if(it, stubs.end(), different);
0269       vector<TTStubRef> ttStubRefs;
0270       ttStubRefs.reserve(distance(start, it));
0271       transform(start, it, back_inserter(ttStubRefs), [](const StubMHT& stub) { return stub.ttStubRef(); });
0272       tracks.push_back(ttStubRefs);
0273     }
0274   }
0275
0276   //
0277   void AnalyzerMHT::associate(const vector<vector<TTStubRef>>& tracks,
0278                               const StubAssociation* ass,
0279                               set<TPPtr>& tps,
0280                               int& sum) const {
0281     for (const vector<TTStubRef>& ttStubRefs : tracks) {
0282       const vector<TPPtr>& tpPtrs = ass->associate(ttStubRefs);
0283       if (tpPtrs.empty())
0284         continue;
0285       sum++;
0286       copy(tpPtrs.begin(), tpPtrs.end(), inserter(tps, tps.begin()));
0287     }
0288   }
0289
0290 }  // namespace trackerTFP
0291
0292 DEFINE_FWK_MODULE(trackerTFP::AnalyzerMHT);