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