TrackFindingTracklet/test/AnalyzerTracklet.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 #include <TEfficiency.h>
0023
0024 #include <vector>
0025 #include <deque>
0026 #include <set>
0027 #include <cmath>
0028 #include <numeric>
0029 #include <sstream>
0030
0031 namespace trklet {
0032
0033   /*! \class  trklet::AnalyzerTracklet
0034    *  \brief  Class to analyze TTTracks found by tracklet pattern recognition
0035    *  \author Thomas Schuh
0036    *  \date   2020, Oct
0037    */
0038   class AnalyzerTracklet : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0039   public:
0040     AnalyzerTracklet(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     // gets all TPs associated too any of the tracks & number of tracks matching at least one TP
0049     void associate(const std::vector<std::vector<TTStubRef>>& tracks,
0050                    const tt::StubAssociation* ass,
0051                    std::set<TPPtr>& tps,
0052                    int& nMatchTrk,
0053                    bool perfect = false) const;
0054
0055     // ED input token of tracks
0056     edm::EDGetTokenT<tt::TTTracks> edGetToken_;
0057     // ED input token of TTStubRef to TPPtr association for tracking efficiency
0058     edm::EDGetTokenT<tt::StubAssociation> edGetTokenSelection_;
0059     // ED input token of TTStubRef to recontructable TPPtr association
0060     edm::EDGetTokenT<tt::StubAssociation> edGetTokenReconstructable_;
0061     // Setup token
0062     edm::ESGetToken<tt::Setup, tt::SetupRcd> esGetTokenSetup_;
0063     // DataFormats token
0064     edm::ESGetToken<DataFormats, ChannelAssignmentRcd> esGetTokenDataFormats_;
0065     // ChannelAssignment token
0066     edm::ESGetToken<ChannelAssignment, ChannelAssignmentRcd> esGetTokenChannelAssignment_;
0067     // stores, calculates and provides run-time constants
0068     const tt::Setup* setup_ = nullptr;
0069     // helper class to extract structured data from tt::Frames
0070     const DataFormats* dataFormats_ = nullptr;
0071     // helper class to assign tracks to channel
0072     const ChannelAssignment* channelAssignment_ = nullptr;
0073     // enables analyze of TPs
0074     bool useMCTruth_;
0075     //
0076     int nEvents_ = 0;
0077
0078     // Histograms
0079
0080     // counts per TFP (processing nonant and event)
0081     TProfile* prof_;
0082     // no. of tracks per nonant
0083     TProfile* profChannel_;
0084     TH1F* hisChannel_;
0085     TH1F* hisEff_;
0086     TH1F* hisEffTotal_;
0087     TEfficiency* eff_;
0088
0089     // printout
0090     std::stringstream log_;
0091   };
0092
0093   AnalyzerTracklet::AnalyzerTracklet(const edm::ParameterSet& iConfig)
0094       : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")) {
0095     usesResource("TFileService");
0096     // book in- and output ED products
0097     const edm::InputTag& inputTag = iConfig.getParameter<edm::InputTag>("InputTag");
0098     edGetToken_ = consumes<tt::TTTracks>(inputTag);
0099     if (useMCTruth_) {
0100       const auto& inputTagSelecttion = iConfig.getParameter<edm::InputTag>("InputTagSelection");
0101       const auto& inputTagReconstructable = iConfig.getParameter<edm::InputTag>("InputTagReconstructable");
0102       edGetTokenSelection_ = consumes<tt::StubAssociation>(inputTagSelecttion);
0103       edGetTokenReconstructable_ = consumes<tt::StubAssociation>(inputTagReconstructable);
0104     }
0105     // book ES products
0106     esGetTokenSetup_ = esConsumes<edm::Transition::BeginRun>();
0107     esGetTokenDataFormats_ = esConsumes<edm::Transition::BeginRun>();
0108     esGetTokenChannelAssignment_ = esConsumes<edm::Transition::BeginRun>();
0109     // log config
0110     log_.setf(std::ios::fixed, std::ios::floatfield);
0111     log_.precision(4);
0112   }
0113
0114   void AnalyzerTracklet::beginRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) {
0115     // helper class to store configurations
0116     setup_ = &iSetup.getData(esGetTokenSetup_);
0117     // helper class to extract structured data from tt::Frames
0118     dataFormats_ = &iSetup.getData(esGetTokenDataFormats_);
0119     // helper class to assign tracks to channel
0120     channelAssignment_ = &iSetup.getData(esGetTokenChannelAssignment_);
0121     // book histograms
0122     edm::Service<TFileService> fs;
0123     TFileDirectory dir;
0124     dir = fs->mkdir("Tracklet");
0125     prof_ = dir.make<TProfile>("Counts", ";", 10, 0.5, 10.5);
0126     prof_->GetXaxis()->SetBinLabel(1, "Stubs");
0127     prof_->GetXaxis()->SetBinLabel(2, "Tracks");
0128     prof_->GetXaxis()->SetBinLabel(3, "Lost Tracks");
0129     prof_->GetXaxis()->SetBinLabel(4, "Matched Tracks");
0130     prof_->GetXaxis()->SetBinLabel(5, "All Tracks");
0131     prof_->GetXaxis()->SetBinLabel(6, "Found TPs");
0132     prof_->GetXaxis()->SetBinLabel(7, "Found selected TPs");
0133     prof_->GetXaxis()->SetBinLabel(8, "Lost TPs");
0134     prof_->GetXaxis()->SetBinLabel(9, "All TPs");
0135     prof_->GetXaxis()->SetBinLabel(10, "Perfectly Found selected TPs");
0136     // channel occupancy
0137     constexpr int maxOcc = 180;
0138     const int numChannels = channelAssignment_->numSeedTypes();
0139     hisChannel_ = dir.make<TH1F>("His Channel Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0140     profChannel_ = dir.make<TProfile>("Prof Channel Occupancy", ";", numChannels, -.5, numChannels - .5);
0141     // Efficiencies
0142     hisEffTotal_ = dir.make<TH1F>("HisTPEtaTotal", ";", 128, -2.5, 2.5);
0143     hisEff_ = dir.make<TH1F>("HisTPEta", ";", 128, -2.5, 2.5);
0144     eff_ = dir.make<TEfficiency>("EffEta", ";", 128, -2.5, 2.5);
0145   }
0146
0147   void AnalyzerTracklet::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0148     auto fill = [](const TPPtr& tpPtr, TH1F* his) { his->Fill(tpPtr->eta()); };
0149     // read in tracklet products
0150     edm::Handle<tt::TTTracks> handle;
0151     iEvent.getByToken<tt::TTTracks>(edGetToken_, handle);
0152     // read in MCTruth
0153     const tt::StubAssociation* selection = nullptr;
0154     const tt::StubAssociation* reconstructable = nullptr;
0155     if (useMCTruth_) {
0156       edm::Handle<tt::StubAssociation> handleSelection;
0157       iEvent.getByToken<tt::StubAssociation>(edGetTokenSelection_, handleSelection);
0158       selection = handleSelection.product();
0159       prof_->Fill(9, selection->numTPs());
0160       edm::Handle<tt::StubAssociation> handleReconstructable;
0161       iEvent.getByToken<tt::StubAssociation>(edGetTokenReconstructable_, handleReconstructable);
0162       reconstructable = handleReconstructable.product();
0163       for (const auto& p : selection->getTrackingParticleToTTStubsMap())
0164         fill(p.first, hisEffTotal_);
0165     }
0166     //
0167     const tt::TTTracks& ttTracks = *handle.product();
0168     std::vector<std::vector<TTTrackRef>> ttTrackRefsRegions(setup_->numRegions());
0169     std::vector<int> nTTTracksRegions(setup_->numRegions(), 0);
0170     for (const TTTrack<Ref_Phase2TrackerDigi_>& ttTrack : ttTracks)
0171       nTTTracksRegions[ttTrack.phiSector()]++;
0172     for (int region = 0; region < setup_->numRegions(); region++)
0173       ttTrackRefsRegions[region].reserve(nTTTracksRegions[region]);
0174     int i(0);
0175     for (const TTTrack<Ref_Phase2TrackerDigi_>& ttTrack : ttTracks)
0176       ttTrackRefsRegions[ttTrack.phiSector()].emplace_back(TTTrackRef(handle, i++));
0177     for (int region = 0; region < setup_->numRegions(); region++) {
0178       const std::vector<TTTrackRef>& ttTrackRefs = ttTrackRefsRegions[region];
0179       const int nStubs =
0180           std::accumulate(ttTrackRefs.begin(), ttTrackRefs.end(), 0, [](int sum, const TTTrackRef& ttTrackRef) {
0181             return sum + ttTrackRef->getStubRefs().size();
0182           });
0183       const int nTracks = ttTrackRefs.size();
0184       prof_->Fill(1, nStubs);
0185       prof_->Fill(2, nTracks);
0186       // no access to lost tracks
0187       prof_->Fill(3, 0);
0188       for (int seedType = 0; seedType < channelAssignment_->numSeedTypes(); seedType++) {
0189         const int nTracks = std::accumulate(
0190             ttTrackRefs.begin(), ttTrackRefs.end(), 0, [seedType](int sum, const TTTrackRef& ttTrackRef) {
0191               return sum + (static_cast<int>(ttTrackRef->trackSeedType()) == seedType ? 1 : 0);
0192             });
0193         hisChannel_->Fill(nTracks);
0194         profChannel_->Fill(seedType, nTracks);
0195       }
0196     }
0197     // analyze tracklet products and associate found tracks with reconstrucable TrackingParticles
0198     std::set<TPPtr> tpPtrs;
0199     std::set<TPPtr> tpPtrsSelection;
0200     std::set<TPPtr> tpPtrsPerfect;
0201     int nAllMatched(0);
0202     // convert vector of tracks to vector of vector of associated stubs
0203     std::vector<std::vector<TTStubRef>> tracks;
0204     tracks.reserve(ttTracks.size());
0205     std::transform(ttTracks.begin(),
0206                    ttTracks.end(),
0207                    std::back_inserter(tracks),
0208                    [](const TTTrack<Ref_Phase2TrackerDigi_>& ttTrack) { return ttTrack.getStubRefs(); });
0209     if (useMCTruth_) {
0210       int tmp(0);
0211       associate(tracks, selection, tpPtrsSelection, tmp);
0212       associate(tracks, selection, tpPtrsPerfect, tmp, true);
0213       associate(tracks, reconstructable, tpPtrs, nAllMatched);
0214     }
0215     for (const TPPtr& tpPtr : tpPtrsSelection)
0216       fill(tpPtr, hisEff_);
0217     prof_->Fill(4, nAllMatched);
0218     prof_->Fill(5, ttTracks.size());
0219     prof_->Fill(6, tpPtrs.size());
0220     prof_->Fill(7, tpPtrsSelection.size());
0221     // no access to lost tp
0222     prof_->Fill(8, 0);
0223     prof_->Fill(10, tpPtrsPerfect.size());
0224     nEvents_++;
0225   }
0226
0227   void AnalyzerTracklet::endJob() {
0228     if (nEvents_ == 0)
0229       return;
0230     // effi
0231     eff_->SetPassedHistogram(*hisEff_, "f");
0232     eff_->SetTotalHistogram(*hisEffTotal_, "f");
0233     // printout SF summary
0234     const double totalTPs = prof_->GetBinContent(9);
0235     const double numStubs = prof_->GetBinContent(1);
0236     const double numTracks = prof_->GetBinContent(2);
0237     const double totalTracks = prof_->GetBinContent(5);
0238     const double numTracksMatched = prof_->GetBinContent(4);
0239     const double numTPsAll = prof_->GetBinContent(6);
0240     const double numTPsEff = prof_->GetBinContent(7);
0241     const double numTPsEffPerfect = prof_->GetBinContent(10);
0242     const double errStubs = prof_->GetBinError(1);
0243     const double errTracks = prof_->GetBinError(2);
0244     const double fracFake = (totalTracks - numTracksMatched) / totalTracks;
0245     const double fracDup = (numTracksMatched - numTPsAll) / totalTracks;
0246     const double eff = numTPsEff / totalTPs;
0247     const double errEff = std::sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0248     const double effPerfect = numTPsEffPerfect / totalTPs;
0249     const double errEffPerfect = std::sqrt(effPerfect * (1. - effPerfect) / totalTPs / nEvents_);
0250     const std::vector<double> nums = {numStubs, numTracks};
0251     const std::vector<double> errs = {errStubs, errTracks};
0252     const int wNums = std::ceil(std::log10(*std::max_element(nums.begin(), nums.end()))) + 5;
0253     const int wErrs = std::ceil(std::log10(*std::max_element(errs.begin(), errs.end()))) + 5;
0254     log_ << "                      Tracklet  SUMMARY                      " << std::endl;
0255     log_ << "number of stubs     per TFP = " << std::setw(wNums) << numStubs << " +- " << std::setw(wErrs) << errStubs
0256          << std::endl;
0257     log_ << "number of tracks    per TFP = " << std::setw(wNums) << numTracks << " +- " << std::setw(wErrs) << errTracks
0258          << std::endl;
0259     log_ << "current tracking efficiency = " << std::setw(wNums) << effPerfect << " +- " << std::setw(wErrs)
0260          << errEffPerfect << std::endl;
0261     log_ << "max     tracking efficiency = " << std::setw(wNums) << eff << " +- " << std::setw(wErrs) << errEff
0262          << std::endl;
0263     log_ << "                  fake rate = " << std::setw(wNums) << fracFake << std::endl;
0264     log_ << "             duplicate rate = " << std::setw(wNums) << fracDup << std::endl;
0265     log_ << "=============================================================";
0266     edm::LogPrint(moduleDescription().moduleName()) << log_.str();
0267   }
0268
0269   // gets all TPs associated too any of the tracks & number of tracks matching at least one TP
0270   void AnalyzerTracklet::associate(const std::vector<std::vector<TTStubRef>>& tracks,
0271                                    const tt::StubAssociation* ass,
0272                                    std::set<TPPtr>& tps,
0273                                    int& nMatchTrk,
0274                                    bool perfect) const {
0275     for (const std::vector<TTStubRef>& ttStubRefs : tracks) {
0276       const std::vector<TPPtr>& tpPtrs = perfect ? ass->associateFinal(ttStubRefs) : ass->associate(ttStubRefs);
0277       if (tpPtrs.empty())
0278         continue;
0279       nMatchTrk++;
0280       std::copy(tpPtrs.begin(), tpPtrs.end(), std::inserter(tps, tps.begin()));
0281     }
0282   }
0283
0284 }  // namespace trklet
0285
0286 DEFINE_FWK_MODULE(trklet::AnalyzerTracklet);