TrackerTFP/test/AnalyzerKFin.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 using namespace std;
0030 using namespace edm;
0031 using namespace trackerTFP;
0032 using namespace tt;
0033
0034 namespace trackerTFP {
0035
0036   /*! \class  trackerTFP::AnalyzerKFin
0037    *  \brief  Class to analyze hardware like structured TTStub Collection generated by Seed Filter
0038    *  \author Thomas Schuh
0039    *  \date   2020, April
0040    */
0041   class AnalyzerKFin : public one::EDAnalyzer<one::WatchRuns, one::SharedResources> {
0042   public:
0043     AnalyzerKFin(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 StreamsTrack& streamsTrack,
0053                     const StreamsStub& streamsStubs,
0054                     vector<vector<TTStubRef>>& tracks,
0055                     int channel) const;
0056     //
0057     void associate(const vector<vector<TTStubRef>>& tracks, const StubAssociation* ass, set<TPPtr>& tps, int& sum) const;
0058
0059     // ED input token of stubs
0060     EDGetTokenT<StreamsStub> edGetTokenAcceptedStubs_;
0061     // ED input token of tracks
0062     EDGetTokenT<StreamsTrack> edGetTokenAcceptedTracks_;
0063     // ED input token of lost stubs
0064     EDGetTokenT<StreamsStub> edGetTokenLostStubs_;
0065     // ED input token of lost tracks
0066     EDGetTokenT<StreamsTrack> edGetTokenLostTracks_;
0067     // ED input token of TTStubRef to TPPtr association for tracking efficiency
0068     EDGetTokenT<StubAssociation> edGetTokenSelection_;
0069     // ED input token of TTStubRef to recontructable TPPtr association
0070     EDGetTokenT<StubAssociation> edGetTokenReconstructable_;
0071     // Setup token
0072     ESGetToken<Setup, SetupRcd> esGetTokenSetup_;
0073     // DataFormats token
0074     ESGetToken<DataFormats, DataFormatsRcd> esGetTokenDataFormats_;
0075     // stores, calculates and provides run-time constants
0076     const Setup* setup_ = nullptr;
0077     // helper class to extract structured data from tt::Frames
0078     const DataFormats* dataFormats_ = nullptr;
0079     // enables analyze of TPs
0080     bool useMCTruth_;
0081     //
0082     int nEvents_ = 0;
0083
0084     // Histograms
0085
0086     TProfile* prof_;
0087     TProfile* profChannel_;
0088     TH1F* hisChannel_;
0089
0090     // printout
0091     stringstream log_;
0092   };
0093
0094   AnalyzerKFin::AnalyzerKFin(const ParameterSet& iConfig) : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")) {
0095     usesResource("TFileService");
0096     // book in- and output ED products
0097     const string& label = iConfig.getParameter<string>("LabelKFin");
0098     const string& branchAcceptedStubs = iConfig.getParameter<string>("BranchAcceptedStubs");
0099     const string& branchAcceptedTracks = iConfig.getParameter<string>("BranchAcceptedTracks");
0100     const string& branchLostStubs = iConfig.getParameter<string>("BranchLostStubs");
0101     const string& branchLostTracks = iConfig.getParameter<string>("BranchLostTracks");
0102     edGetTokenAcceptedStubs_ = consumes<StreamsStub>(InputTag(label, branchAcceptedStubs));
0103     edGetTokenAcceptedTracks_ = consumes<StreamsTrack>(InputTag(label, branchAcceptedTracks));
0104     edGetTokenLostStubs_ = consumes<StreamsStub>(InputTag(label, branchLostStubs));
0105     edGetTokenLostTracks_ = consumes<StreamsTrack>(InputTag(label, branchLostTracks));
0106     if (useMCTruth_) {
0107       const auto& inputTagSelecttion = iConfig.getParameter<InputTag>("InputTagSelection");
0108       const auto& inputTagReconstructable = iConfig.getParameter<InputTag>("InputTagReconstructable");
0109       edGetTokenSelection_ = consumes<StubAssociation>(inputTagSelecttion);
0110       edGetTokenReconstructable_ = consumes<StubAssociation>(inputTagReconstructable);
0111     }
0112     // book ES products
0113     esGetTokenSetup_ = esConsumes<Setup, SetupRcd, Transition::BeginRun>();
0114     esGetTokenDataFormats_ = esConsumes<DataFormats, DataFormatsRcd, Transition::BeginRun>();
0115     // log config
0116     log_.setf(ios::fixed, ios::floatfield);
0117     log_.precision(4);
0118   }
0119
0120   void AnalyzerKFin::beginRun(const Run& iEvent, const EventSetup& iSetup) {
0121     // helper class to store configurations
0122     setup_ = &iSetup.getData(esGetTokenSetup_);
0123     // helper class to extract structured data from tt::Frames
0124     dataFormats_ = &iSetup.getData(esGetTokenDataFormats_);
0125     // book histograms
0126     Service<TFileService> fs;
0127     TFileDirectory dir;
0128     dir = fs->mkdir("KFin");
0129     prof_ = dir.make<TProfile>("Counts", ";", 9, 0.5, 9.5);
0130     prof_->GetXaxis()->SetBinLabel(1, "Stubs");
0131     prof_->GetXaxis()->SetBinLabel(2, "Tracks");
0132     prof_->GetXaxis()->SetBinLabel(3, "Lost Tracks");
0133     prof_->GetXaxis()->SetBinLabel(4, "Matched Tracks");
0134     prof_->GetXaxis()->SetBinLabel(5, "All Tracks");
0135     prof_->GetXaxis()->SetBinLabel(6, "Found TPs");
0136     prof_->GetXaxis()->SetBinLabel(7, "Found selected TPs");
0137     prof_->GetXaxis()->SetBinLabel(8, "Lost TPs");
0138     prof_->GetXaxis()->SetBinLabel(9, "All TPs");
0139     // channel occupancy
0140     constexpr int maxOcc = 180;
0141     const int numChannels = dataFormats_->numChannel(Process::kfin);
0142     hisChannel_ = dir.make<TH1F>("His Channel Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0143     profChannel_ = dir.make<TProfile>("Prof Channel Occupancy", ";", numChannels, -.5, numChannels - .5);
0144   }
0145
0146   void AnalyzerKFin::analyze(const Event& iEvent, const EventSetup& iSetup) {
0147     // read in ht products
0148     Handle<StreamsStub> handleAcceptedStubs;
0149     iEvent.getByToken<StreamsStub>(edGetTokenAcceptedStubs_, handleAcceptedStubs);
0150     const StreamsStub& acceptedStubs = *handleAcceptedStubs;
0151     Handle<StreamsTrack> handleAcceptedTracks;
0152     iEvent.getByToken<StreamsTrack>(edGetTokenAcceptedTracks_, handleAcceptedTracks);
0153     const StreamsTrack& acceptedTracks = *handleAcceptedTracks;
0154     Handle<StreamsStub> handleLostStubs;
0155     iEvent.getByToken<StreamsStub>(edGetTokenLostStubs_, handleLostStubs);
0156     const StreamsStub& lostStubs = *handleLostStubs;
0157     Handle<StreamsTrack> handleLostTracks;
0158     iEvent.getByToken<StreamsTrack>(edGetTokenLostTracks_, handleLostTracks);
0159     const StreamsTrack& lostTracks = *handleLostTracks;
0160     // read in MCTruth
0161     const StubAssociation* selection = nullptr;
0162     const StubAssociation* reconstructable = nullptr;
0163     if (useMCTruth_) {
0164       Handle<StubAssociation> handleSelection;
0165       iEvent.getByToken<StubAssociation>(edGetTokenSelection_, handleSelection);
0166       selection = handleSelection.product();
0167       prof_->Fill(9, selection->numTPs());
0168       Handle<StubAssociation> handleReconstructable;
0169       iEvent.getByToken<StubAssociation>(edGetTokenReconstructable_, handleReconstructable);
0170       reconstructable = handleReconstructable.product();
0171     }
0172     // analyze ht products and associate found tracks with reconstrucable TrackingParticles
0173     set<TPPtr> tpPtrs;
0174     set<TPPtr> tpPtrsSelection;
0175     set<TPPtr> tpPtrsLost;
0176     int allMatched(0);
0177     int allTracks(0);
0178     for (int region = 0; region < setup_->numRegions(); region++) {
0179       const int offset = region * dataFormats_->numChannel(Process::kf);
0180       int nStubs(0);
0181       int nTracks(0);
0182       int nLost(0);
0183       for (int channel = 0; channel < dataFormats_->numChannel(Process::kf); channel++) {
0184         vector<vector<TTStubRef>> tracks;
0185         formTracks(acceptedTracks, acceptedStubs, tracks, offset + channel);
0186         vector<vector<TTStubRef>> lost;
0187         formTracks(lostTracks, lostStubs, lost, offset + channel);
0188         nTracks += tracks.size();
0189         nStubs += accumulate(tracks.begin(), tracks.end(), 0, [](int sum, const vector<TTStubRef>& track) {
0190           return sum + (int)track.size();
0191         });
0192         nLost += lost.size();
0193         allTracks += tracks.size();
0194         if (!useMCTruth_)
0195           continue;
0196         int tmp(0);
0197         associate(tracks, selection, tpPtrsSelection, tmp);
0198         associate(lost, selection, tpPtrsLost, tmp);
0199         associate(tracks, reconstructable, tpPtrs, allMatched);
0200       }
0201       prof_->Fill(1, nStubs);
0202       prof_->Fill(2, nTracks);
0203       prof_->Fill(3, nLost);
0204     }
0205     vector<TPPtr> recovered;
0206     recovered.reserve(tpPtrsLost.size());
0207     set_intersection(tpPtrsLost.begin(), tpPtrsLost.end(), tpPtrs.begin(), tpPtrs.end(), back_inserter(recovered));
0208     for (const TPPtr& tpPtr : recovered)
0209       tpPtrsLost.erase(tpPtr);
0210     prof_->Fill(4, allMatched);
0211     prof_->Fill(5, allTracks);
0212     prof_->Fill(6, tpPtrs.size());
0213     prof_->Fill(7, tpPtrsSelection.size());
0214     prof_->Fill(8, tpPtrsLost.size());
0215     nEvents_++;
0216   }
0217
0218   void AnalyzerKFin::endJob() {
0219     if (nEvents_ == 0)
0220       return;
0221     // printout SF summary
0222     const double totalTPs = prof_->GetBinContent(9);
0223     const double numStubs = prof_->GetBinContent(1);
0224     const double numTracks = prof_->GetBinContent(2);
0225     const double numTracksLost = prof_->GetBinContent(3);
0226     const double totalTracks = prof_->GetBinContent(5);
0227     const double numTracksMatched = prof_->GetBinContent(4);
0228     const double numTPsAll = prof_->GetBinContent(6);
0229     const double numTPsEff = prof_->GetBinContent(7);
0230     const double numTPsLost = prof_->GetBinContent(8);
0231     const double errStubs = prof_->GetBinError(1);
0232     const double errTracks = prof_->GetBinError(2);
0233     const double errTracksLost = prof_->GetBinError(3);
0234     const double fracFake = (totalTracks - numTracksMatched) / totalTracks;
0235     const double fracDup = (numTracksMatched - numTPsAll) / totalTracks;
0236     const double eff = numTPsEff / totalTPs;
0237     const double errEff = sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0238     const double effLoss = numTPsLost / totalTPs;
0239     const double errEffLoss = sqrt(effLoss * (1. - effLoss) / totalTPs / nEvents_);
0240     const vector<double> nums = {numStubs, numTracks, numTracksLost};
0241     const vector<double> errs = {errStubs, errTracks, errTracksLost};
0242     const int wNums = ceil(log10(*max_element(nums.begin(), nums.end()))) + 5;
0243     const int wErrs = ceil(log10(*max_element(errs.begin(), errs.end()))) + 5;
0244     log_ << "                        KFin  SUMMARY                        " << endl;
0245     log_ << "number of stubs       per TFP = " << setw(wNums) << numStubs << " +- " << setw(wErrs) << errStubs << endl;
0246     log_ << "number of tracks      per TFP = " << setw(wNums) << numTracks << " +- " << setw(wErrs) << errTracks
0247          << endl;
0248     log_ << "number of lost tracks per TFP = " << setw(wNums) << numTracksLost << " +- " << setw(wErrs) << errTracksLost
0249          << endl;
0250     log_ << "     max  tracking efficiency = " << setw(wNums) << eff << " +- " << setw(wErrs) << errEff << endl;
0251     log_ << "     lost tracking efficiency = " << setw(wNums) << effLoss << " +- " << setw(wErrs) << errEffLoss << endl;
0252     log_ << "                    fake rate = " << setw(wNums) << fracFake << endl;
0253     log_ << "               duplicate rate = " << setw(wNums) << fracDup << endl;
0254     log_ << "=============================================================";
0255     LogPrint("L1Trigger/TrackerTFP") << log_.str();
0256   }
0257
0258   //
0259   void AnalyzerKFin::formTracks(const StreamsTrack& streamsTrack,
0260                                 const StreamsStub& streamsStubs,
0261                                 vector<vector<TTStubRef>>& tracks,
0262                                 int channel) const {
0263     const int offset = channel * setup_->numLayers();
0264     const StreamTrack& streamTrack = streamsTrack[channel];
0265     const int numTracks = accumulate(streamTrack.begin(), streamTrack.end(), 0, [](int sum, const FrameTrack& frame) {
0266       return sum + (frame.first.isNonnull() ? 1 : 0);
0267     });
0268     tracks.reserve(numTracks);
0269     for (int frame = 0; frame < (int)streamTrack.size(); frame++) {
0270       const FrameTrack& frameTrack = streamTrack[frame];
0271       if (frameTrack.first.isNull())
0272         continue;
0273       const auto end = find_if(next(streamTrack.begin(), frame + 1), streamTrack.end(), [](const FrameTrack& frame) {
0274         return frame.first.isNonnull();
0275       });
0276       const int size = distance(next(streamTrack.begin(), frame), end);
0277       int numStubs(0);
0278       for (int layer = 0; layer < setup_->numLayers(); layer++) {
0279         const StreamStub& stream = streamsStubs[offset + layer];
0280         numStubs +=
0281             accumulate(stream.begin() + frame, stream.begin() + frame + size, 0, [](int sum, const FrameStub& frame) {
0282               return sum + (frame.first.isNonnull() ? 1 : 0);
0283             });
0284       }
0285       vector<TTStubRef> stubs;
0286       stubs.reserve(numStubs);
0287       for (int layer = 0; layer < setup_->numLayers(); layer++) {
0288         for (int f = frame; f < frame + size; f++) {
0289           const FrameStub& stub = streamsStubs[offset + layer][f];
0290           if (stub.first.isNonnull())
0291             stubs.push_back(stub.first);
0292         }
0293       }
0294       tracks.push_back(stubs);
0295     }
0296   }
0297
0298   //
0299   void AnalyzerKFin::associate(const vector<vector<TTStubRef>>& tracks,
0300                                const StubAssociation* ass,
0301                                set<TPPtr>& tps,
0302                                int& sum) const {
0303     for (const vector<TTStubRef>& ttStubRefs : tracks) {
0304       const vector<TPPtr>& tpPtrs = ass->associate(ttStubRefs);
0305       if (tpPtrs.empty())
0306         continue;
0307       sum++;
0308       copy(tpPtrs.begin(), tpPtrs.end(), inserter(tps, tps.begin()));
0309     }
0310   }
0311
0312 }  // namespace trackerTFP
0313
0314 DEFINE_FWK_MODULE(trackerTFP::AnalyzerKFin);