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