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