File indexing completed on 2025-06-03 00:12:26
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
0034
0035
0036
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
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
0056 edm::EDGetTokenT<tt::TTTracks> edGetToken_;
0057
0058 edm::EDGetTokenT<tt::StubAssociation> edGetTokenSelection_;
0059
0060 edm::EDGetTokenT<tt::StubAssociation> edGetTokenReconstructable_;
0061
0062 edm::ESGetToken<tt::Setup, tt::SetupRcd> esGetTokenSetup_;
0063
0064 edm::ESGetToken<DataFormats, ChannelAssignmentRcd> esGetTokenDataFormats_;
0065
0066 edm::ESGetToken<ChannelAssignment, ChannelAssignmentRcd> esGetTokenChannelAssignment_;
0067
0068 const tt::Setup* setup_ = nullptr;
0069
0070 const DataFormats* dataFormats_ = nullptr;
0071
0072 const ChannelAssignment* channelAssignment_ = nullptr;
0073
0074 bool useMCTruth_;
0075
0076 int nEvents_ = 0;
0077
0078
0079
0080
0081 TProfile* prof_;
0082
0083 TProfile* profChannel_;
0084 TH1F* hisChannel_;
0085 TH1F* hisEff_;
0086 TH1F* hisEffTotal_;
0087 TEfficiency* eff_;
0088
0089
0090 std::stringstream log_;
0091 };
0092
0093 AnalyzerTracklet::AnalyzerTracklet(const edm::ParameterSet& iConfig)
0094 : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")) {
0095 usesResource("TFileService");
0096
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
0106 esGetTokenSetup_ = esConsumes<edm::Transition::BeginRun>();
0107 esGetTokenDataFormats_ = esConsumes<edm::Transition::BeginRun>();
0108 esGetTokenChannelAssignment_ = esConsumes<edm::Transition::BeginRun>();
0109
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
0116 setup_ = &iSetup.getData(esGetTokenSetup_);
0117
0118 dataFormats_ = &iSetup.getData(esGetTokenDataFormats_);
0119
0120 channelAssignment_ = &iSetup.getData(esGetTokenChannelAssignment_);
0121
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
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
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
0150 edm::Handle<tt::TTTracks> handle;
0151 iEvent.getByToken<tt::TTTracks>(edGetToken_, handle);
0152
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
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
0198 std::set<TPPtr> tpPtrs;
0199 std::set<TPPtr> tpPtrsSelection;
0200 std::set<TPPtr> tpPtrsPerfect;
0201 int nAllMatched(0);
0202
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
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
0231 eff_->SetPassedHistogram(*hisEff_, "f");
0232 eff_->SetTotalHistogram(*hisEffTotal_, "f");
0233
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
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 }
0285
0286 DEFINE_FWK_MODULE(trklet::AnalyzerTracklet);