TrackFindingTracklet/test/AnalyzerKFout.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 #include "L1Trigger/TrackerTFP/interface/LayerEncoding.h"
0019 #include "L1Trigger/TrackerTFP/interface/KalmanFilterFormats.h"
0020
0021 #include <TProfile.h>
0022 #include <TH1F.h>
0023
0024 #include <vector>
0025 #include <deque>
0026 #include <set>
0027 #include <cmath>
0028 #include <numeric>
0029 #include <sstream>
0030
0031 using namespace std;
0032 using namespace edm;
0033 using namespace tt;
0034 using namespace trackerTFP;
0035
0036 namespace trklet {
0037
0038   /*! \class  trklet::AnalyzerKFout
0039    *  \brief  Class to analyze hardware like structured Track Collection generated by Kalman Filter output module
0040    *  \author Thomas Schuh
0041    *  \date   2021, Aug
0042    */
0043   class AnalyzerKFout : public one::EDAnalyzer<one::WatchRuns, one::SharedResources> {
0044   public:
0045     AnalyzerKFout(const ParameterSet& iConfig);
0046     void beginJob() override {}
0047     void beginRun(const Run& iEvent, const EventSetup& iSetup) override;
0048     void analyze(const Event& iEvent, const EventSetup& iSetup) override;
0049     void endRun(const Run& iEvent, const EventSetup& iSetup) override {}
0050     void endJob() override;
0051
0052   private:
0053     //
0054     void associate(const set<TTTrackRef>& ttTracks, const StubAssociation* ass, set<TPPtr>& tps, int& sum) const;
0055
0056     // ED input token of accepted Stubs
0057     EDGetTokenT<StreamsTrack> edGetTokenAccepted_;
0058     // ED input token of lost Tracks
0059     EDGetTokenT<StreamsTrack> 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     // enables analyze of TPs
0071     bool useMCTruth_;
0072     //
0073     int nEvents_ = 0;
0074
0075     // Histograms
0076
0077     TProfile* prof_;
0078     TProfile* profChannel_;
0079     TH1F* hisChannel_;
0080
0081     // printout
0082     stringstream log_;
0083   };
0084
0085   AnalyzerKFout::AnalyzerKFout(const ParameterSet& iConfig) : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")) {
0086     usesResource("TFileService");
0087     // book in- and output ED products
0088     const string& label = iConfig.getParameter<string>("LabelKFout");
0089     const string& branchAccepted = iConfig.getParameter<string>("BranchAcceptedTracks");
0090     const string& branchLost = iConfig.getParameter<string>("BranchLostTracks");
0091     edGetTokenAccepted_ = consumes<StreamsTrack>(InputTag(label, branchAccepted));
0092     edGetTokenLost_ = consumes<StreamsTrack>(InputTag(label, branchLost));
0093     if (useMCTruth_) {
0094       const auto& inputTagSelecttion = iConfig.getParameter<InputTag>("InputTagSelection");
0095       const auto& inputTagReconstructable = iConfig.getParameter<InputTag>("InputTagReconstructable");
0096       edGetTokenSelection_ = consumes<StubAssociation>(inputTagSelecttion);
0097       edGetTokenReconstructable_ = consumes<StubAssociation>(inputTagReconstructable);
0098     }
0099     // book ES products
0100     esGetTokenSetup_ = esConsumes<Setup, SetupRcd, Transition::BeginRun>();
0101     // log config
0102     log_.setf(ios::fixed, ios::floatfield);
0103     log_.precision(4);
0104   }
0105
0106   void AnalyzerKFout::beginRun(const Run& iEvent, const EventSetup& iSetup) {
0107     // helper class to store configurations
0108     setup_ = &iSetup.getData(esGetTokenSetup_);
0109     // book histograms
0110     Service<TFileService> fs;
0111     TFileDirectory dir;
0112     dir = fs->mkdir("KFout");
0113     prof_ = dir.make<TProfile>("Counts", ";", 9, 0.5, 9.5);
0114     prof_->GetXaxis()->SetBinLabel(1, "Stubs");
0115     prof_->GetXaxis()->SetBinLabel(2, "Tracks");
0116     prof_->GetXaxis()->SetBinLabel(3, "Lost Tracks");
0117     prof_->GetXaxis()->SetBinLabel(4, "Matched Tracks");
0118     prof_->GetXaxis()->SetBinLabel(5, "All Tracks");
0119     prof_->GetXaxis()->SetBinLabel(6, "Found TPs");
0120     prof_->GetXaxis()->SetBinLabel(7, "Found selected TPs");
0121     prof_->GetXaxis()->SetBinLabel(8, "Lost TPs");
0122     prof_->GetXaxis()->SetBinLabel(9, "All TPs");
0123     // channel occupancy
0124     constexpr int maxOcc = 180;
0125     const int numChannels = setup_->tfpNumChannel();
0126     hisChannel_ = dir.make<TH1F>("His Channel Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0127     profChannel_ = dir.make<TProfile>("Prof Channel Occupancy", ";", numChannels, -.5, numChannels - .5);
0128   }
0129
0130   void AnalyzerKFout::analyze(const Event& iEvent, const EventSetup& iSetup) {
0131     // read in kf products
0132     Handle<StreamsTrack> handleAccepted;
0133     iEvent.getByToken<StreamsTrack>(edGetTokenAccepted_, handleAccepted);
0134     Handle<StreamsTrack> handleLost;
0135     iEvent.getByToken<StreamsTrack>(edGetTokenLost_, handleLost);
0136     // read in MCTruth
0137     const StubAssociation* selection = nullptr;
0138     const StubAssociation* reconstructable = nullptr;
0139     if (useMCTruth_) {
0140       Handle<StubAssociation> handleSelection;
0141       iEvent.getByToken<StubAssociation>(edGetTokenSelection_, handleSelection);
0142       selection = handleSelection.product();
0143       prof_->Fill(9, selection->numTPs());
0144       Handle<StubAssociation> handleReconstructable;
0145       iEvent.getByToken<StubAssociation>(edGetTokenReconstructable_, handleReconstructable);
0146       reconstructable = handleReconstructable.product();
0147     }
0148     // analyze kfout products and associate found tracks with reconstrucable TrackingParticles
0149     set<TPPtr> tpPtrs;
0150     set<TPPtr> tpPtrsSelection;
0151     set<TPPtr> tpPtrsLost;
0152     int allMatched(0);
0153     int allTracks(0);
0154     for (int region = 0; region < setup_->numRegions(); region++) {
0155       int nStubsRegion(0);
0156       int nTracksRegion(0);
0157       int nLostRegion(0);
0158       for (int channel = 0; channel < setup_->tfpNumChannel(); channel++) {
0159         const int index = region * setup_->tfpNumChannel() + channel;
0160         const StreamTrack& accepted = handleAccepted->at(index);
0161         const StreamTrack& lost = handleLost->at(index);
0162         hisChannel_->Fill(accepted.size());
0163         profChannel_->Fill(channel, accepted.size());
0164         set<TTTrackRef> tracks;
0165         for (const FrameTrack& frame : accepted)
0166           if (frame.first.isNonnull())
0167             tracks.insert(frame.first);
0168         nTracksRegion += tracks.size();
0169         nStubsRegion += accumulate(tracks.begin(), tracks.end(), 0, [](int sum, const TTTrackRef& ttTrackRef) {
0170           return sum + (int)ttTrackRef->getStubRefs().size();
0171         });
0172         set<TTTrackRef> tracksLost;
0173         for (const FrameTrack& frame : lost)
0174           if (frame.first.isNonnull())
0175             tracksLost.insert(frame.first);
0176         nLostRegion += tracksLost.size();
0177         allTracks += tracks.size();
0178         int tmp(0);
0179         associate(tracks, selection, tpPtrsSelection, tmp);
0180         associate(tracksLost, selection, tpPtrsLost, tmp);
0181         associate(tracks, reconstructable, tpPtrs, allMatched);
0182       }
0183       prof_->Fill(1, nStubsRegion);
0184       prof_->Fill(2, nTracksRegion);
0185       prof_->Fill(3, nLostRegion);
0186     }
0187     deque<TPPtr> tpPtrsRealLost;
0188     set_difference(tpPtrsLost.begin(), tpPtrsLost.end(), tpPtrs.begin(), tpPtrs.end(), back_inserter(tpPtrsRealLost));
0189     prof_->Fill(4, allMatched);
0190     prof_->Fill(5, allTracks);
0191     prof_->Fill(6, tpPtrs.size());
0192     prof_->Fill(7, tpPtrsSelection.size());
0193     prof_->Fill(8, tpPtrsRealLost.size());
0194     nEvents_++;
0195   }
0196
0197   void AnalyzerKFout::endJob() {
0198     if (nEvents_ == 0)
0199       return;
0200     // printout SF summary
0201     const double totalTPs = prof_->GetBinContent(9);
0202     const double numStubs = prof_->GetBinContent(1);
0203     const double numTracks = prof_->GetBinContent(2);
0204     const double numTracksLost = prof_->GetBinContent(3);
0205     const double totalTracks = prof_->GetBinContent(5);
0206     const double numTracksMatched = prof_->GetBinContent(4);
0207     const double numTPsAll = prof_->GetBinContent(6);
0208     const double numTPsEff = prof_->GetBinContent(7);
0209     const double numTPsLost = prof_->GetBinContent(8);
0210     const double errStubs = prof_->GetBinError(1);
0211     const double errTracks = prof_->GetBinError(2);
0212     const double errTracksLost = prof_->GetBinError(3);
0213     const double fracFake = (totalTracks - numTracksMatched) / totalTracks;
0214     const double fracDup = (numTracksMatched - numTPsAll) / totalTracks;
0215     const double eff = numTPsEff / totalTPs;
0216     const double errEff = sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0217     const double effLoss = numTPsLost / totalTPs;
0218     const double errEffLoss = sqrt(effLoss * (1. - effLoss) / totalTPs / nEvents_);
0219     const vector<double> nums = {numStubs, numTracks, numTracksLost};
0220     const vector<double> errs = {errStubs, errTracks, errTracksLost};
0221     const int wNums = ceil(log10(*max_element(nums.begin(), nums.end()))) + 5;
0222     const int wErrs = ceil(log10(*max_element(errs.begin(), errs.end()))) + 5;
0223     log_ << "                        KFout SUMMARY                        " << endl;
0224     //log_ << "number of stubs       per TFP = " << setw(wNums) << numStubs << " +- " << setw(wErrs) << errStubs << endl;
0225     log_ << "number of tracks      per TFP = " << setw(wNums) << numTracks << " +- " << setw(wErrs) << errTracks
0226          << endl;
0227     log_ << "number of lost tracks per TFP = " << setw(wNums) << numTracksLost << " +- " << setw(wErrs) << errTracksLost
0228          << endl;
0229     log_ << "          tracking efficiency = " << setw(wNums) << eff << " +- " << setw(wErrs) << errEff << endl;
0230     log_ << "     lost tracking efficiency = " << setw(wNums) << effLoss << " +- " << setw(wErrs) << errEffLoss << endl;
0231     log_ << "                    fake rate = " << setw(wNums) << fracFake << endl;
0232     log_ << "               duplicate rate = " << setw(wNums) << fracDup << endl;
0233     log_ << "=============================================================";
0234     LogPrint("L1Trigger/TrackerTFP") << log_.str();
0235   }
0236
0237   //
0238   void AnalyzerKFout::associate(const set<TTTrackRef>& ttTracks,
0239                                 const StubAssociation* ass,
0240                                 set<TPPtr>& tps,
0241                                 int& sum) const {
0242     for (const TTTrackRef& ttTrack : ttTracks) {
0243       const vector<TTStubRef>& ttStubRefs = ttTrack->getStubRefs();
0244       const vector<TPPtr>& tpPtrs = ass->associate(ttStubRefs);
0245       if (tpPtrs.empty())
0246         continue;
0247       sum++;
0248       copy(tpPtrs.begin(), tpPtrs.end(), inserter(tps, tps.begin()));
0249     }
0250   }
0251
0252 }  // namespace trklet
0253
0254 DEFINE_FWK_MODULE(trklet::AnalyzerKFout);