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