File indexing completed on 2022-10-14 01:44:13
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 #include "L1Trigger/TrackFindingTracklet/interface/ChannelAssignment.h"
0019
0020 #include <TProfile.h>
0021 #include <TH1F.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
0038
0039
0040
0041
0042 class AnalyzerKFin : public one::EDAnalyzer<one::WatchRuns, one::SharedResources> {
0043 public:
0044 AnalyzerKFin(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
0053 void formTracks(const StreamsTrack& streamsTrack,
0054 const StreamsStub& streamsStubs,
0055 vector<vector<TTStubRef>>& tracks,
0056 int channel) const;
0057
0058 void associate(const vector<vector<TTStubRef>>& tracks,
0059 const StubAssociation* ass,
0060 set<TPPtr>& tps,
0061 int& sum,
0062 bool perfect = false) const;
0063
0064
0065 EDGetTokenT<StreamsStub> edGetTokenAcceptedStubs_;
0066
0067 EDGetTokenT<StreamsTrack> edGetTokenAcceptedTracks_;
0068
0069 EDGetTokenT<StreamsStub> edGetTokenLostStubs_;
0070
0071 EDGetTokenT<StreamsTrack> edGetTokenLostTracks_;
0072
0073 EDGetTokenT<StubAssociation> edGetTokenSelection_;
0074
0075 EDGetTokenT<StubAssociation> edGetTokenReconstructable_;
0076
0077 ESGetToken<Setup, SetupRcd> esGetTokenSetup_;
0078
0079 ESGetToken<DataFormats, DataFormatsRcd> esGetTokenDataFormats_;
0080
0081 ESGetToken<ChannelAssignment, ChannelAssignmentRcd> esGetTokenChannelAssignment_;
0082
0083 const Setup* setup_ = nullptr;
0084
0085 const DataFormats* dataFormats_ = nullptr;
0086
0087 const ChannelAssignment* channelAssignment_ = nullptr;
0088
0089 bool useMCTruth_;
0090
0091 int nEvents_ = 0;
0092
0093
0094
0095 TProfile* prof_;
0096 TProfile* profChannel_;
0097 TH1F* hisChannel_;
0098
0099
0100 stringstream log_;
0101 };
0102
0103 AnalyzerKFin::AnalyzerKFin(const ParameterSet& iConfig) : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")) {
0104 usesResource("TFileService");
0105
0106 const string& label = iConfig.getParameter<string>("LabelKFin");
0107 const string& branchAcceptedStubs = iConfig.getParameter<string>("BranchAcceptedStubs");
0108 const string& branchAcceptedTracks = iConfig.getParameter<string>("BranchAcceptedTracks");
0109 const string& branchLostStubs = iConfig.getParameter<string>("BranchLostStubs");
0110 const string& branchLostTracks = iConfig.getParameter<string>("BranchLostTracks");
0111 edGetTokenAcceptedStubs_ = consumes<StreamsStub>(InputTag(label, branchAcceptedStubs));
0112 edGetTokenAcceptedTracks_ = consumes<StreamsTrack>(InputTag(label, branchAcceptedTracks));
0113 edGetTokenLostStubs_ = consumes<StreamsStub>(InputTag(label, branchLostStubs));
0114 edGetTokenLostTracks_ = consumes<StreamsTrack>(InputTag(label, branchLostTracks));
0115 if (useMCTruth_) {
0116 const auto& inputTagSelecttion = iConfig.getParameter<InputTag>("InputTagSelection");
0117 const auto& inputTagReconstructable = iConfig.getParameter<InputTag>("InputTagReconstructable");
0118 edGetTokenSelection_ = consumes<StubAssociation>(inputTagSelecttion);
0119 edGetTokenReconstructable_ = consumes<StubAssociation>(inputTagReconstructable);
0120 }
0121
0122 esGetTokenSetup_ = esConsumes<Setup, SetupRcd, Transition::BeginRun>();
0123 esGetTokenDataFormats_ = esConsumes<DataFormats, DataFormatsRcd, Transition::BeginRun>();
0124 esGetTokenChannelAssignment_ = esConsumes<ChannelAssignment, ChannelAssignmentRcd, Transition::BeginRun>();
0125
0126 log_.setf(ios::fixed, ios::floatfield);
0127 log_.precision(4);
0128 }
0129
0130 void AnalyzerKFin::beginRun(const Run& iEvent, const EventSetup& iSetup) {
0131
0132 setup_ = &iSetup.getData(esGetTokenSetup_);
0133
0134 dataFormats_ = &iSetup.getData(esGetTokenDataFormats_);
0135
0136 channelAssignment_ = &iSetup.getData(esGetTokenChannelAssignment_);
0137
0138 Service<TFileService> fs;
0139 TFileDirectory dir;
0140 dir = fs->mkdir("KFin");
0141 prof_ = dir.make<TProfile>("Counts", ";", 10, 0.5, 10.5);
0142 prof_->GetXaxis()->SetBinLabel(1, "Stubs");
0143 prof_->GetXaxis()->SetBinLabel(2, "Tracks");
0144 prof_->GetXaxis()->SetBinLabel(3, "Lost Tracks");
0145 prof_->GetXaxis()->SetBinLabel(4, "Matched Tracks");
0146 prof_->GetXaxis()->SetBinLabel(5, "All Tracks");
0147 prof_->GetXaxis()->SetBinLabel(6, "Found TPs");
0148 prof_->GetXaxis()->SetBinLabel(7, "Found selected TPs");
0149 prof_->GetXaxis()->SetBinLabel(8, "Lost TPs");
0150 prof_->GetXaxis()->SetBinLabel(9, "All TPs");
0151 prof_->GetXaxis()->SetBinLabel(10, "Perfect TPs");
0152
0153 constexpr int maxOcc = 180;
0154 const int numChannels = channelAssignment_->numChannelsTrack() * setup_->numLayers() * setup_->numRegions();
0155 hisChannel_ = dir.make<TH1F>("His Channel Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0156 profChannel_ = dir.make<TProfile>("Prof Channel Occupancy", ";", numChannels, -.5, numChannels - .5);
0157 }
0158
0159 void AnalyzerKFin::analyze(const Event& iEvent, const EventSetup& iSetup) {
0160
0161 Handle<StreamsStub> handleAcceptedStubs;
0162 iEvent.getByToken<StreamsStub>(edGetTokenAcceptedStubs_, handleAcceptedStubs);
0163 const StreamsStub& acceptedStubs = *handleAcceptedStubs;
0164 Handle<StreamsTrack> handleAcceptedTracks;
0165 iEvent.getByToken<StreamsTrack>(edGetTokenAcceptedTracks_, handleAcceptedTracks);
0166 const StreamsTrack& acceptedTracks = *handleAcceptedTracks;
0167 Handle<StreamsStub> handleLostStubs;
0168 iEvent.getByToken<StreamsStub>(edGetTokenLostStubs_, handleLostStubs);
0169 const StreamsStub& lostStubs = *handleLostStubs;
0170 Handle<StreamsTrack> handleLostTracks;
0171 iEvent.getByToken<StreamsTrack>(edGetTokenLostTracks_, handleLostTracks);
0172 const StreamsTrack& lostTracks = *handleLostTracks;
0173
0174 const StubAssociation* selection = nullptr;
0175 const StubAssociation* reconstructable = nullptr;
0176 if (useMCTruth_) {
0177 Handle<StubAssociation> handleSelection;
0178 iEvent.getByToken<StubAssociation>(edGetTokenSelection_, handleSelection);
0179 selection = handleSelection.product();
0180 prof_->Fill(9, selection->numTPs());
0181 Handle<StubAssociation> handleReconstructable;
0182 iEvent.getByToken<StubAssociation>(edGetTokenReconstructable_, handleReconstructable);
0183 reconstructable = handleReconstructable.product();
0184 }
0185
0186 set<TPPtr> tpPtrs;
0187 set<TPPtr> tpPtrsSelection;
0188 set<TPPtr> tpPtrsPerfect;
0189 set<TPPtr> tpPtrsLost;
0190 int allMatched(0);
0191 int allTracks(0);
0192 for (int region = 0; region < setup_->numRegions(); region++) {
0193 const int offset = region * channelAssignment_->numChannelsTrack();
0194 int nStubs(0);
0195 int nTracks(0);
0196 int nLost(0);
0197 for (int channel = 0; channel < channelAssignment_->numChannelsTrack(); channel++) {
0198 vector<vector<TTStubRef>> tracks;
0199 formTracks(acceptedTracks, acceptedStubs, tracks, offset + channel);
0200 vector<vector<TTStubRef>> lost;
0201 formTracks(lostTracks, lostStubs, lost, offset + channel);
0202 nTracks += tracks.size();
0203 nStubs += accumulate(tracks.begin(), tracks.end(), 0, [](int& sum, const vector<TTStubRef>& track) {
0204 return sum += (int)track.size();
0205 });
0206 nLost += lost.size();
0207 allTracks += tracks.size();
0208 if (!useMCTruth_)
0209 continue;
0210 int tmp(0);
0211 associate(tracks, selection, tpPtrsSelection, tmp);
0212 associate(tracks, selection, tpPtrsPerfect, tmp, true);
0213 associate(lost, selection, tpPtrsLost, tmp);
0214 associate(tracks, reconstructable, tpPtrs, allMatched);
0215 }
0216 prof_->Fill(1, nStubs);
0217 prof_->Fill(2, nTracks);
0218 prof_->Fill(3, nLost);
0219 }
0220 vector<TPPtr> recovered;
0221 recovered.reserve(tpPtrsLost.size());
0222 set_intersection(tpPtrsLost.begin(), tpPtrsLost.end(), tpPtrs.begin(), tpPtrs.end(), back_inserter(recovered));
0223 for (const TPPtr& tpPtr : recovered)
0224 tpPtrsLost.erase(tpPtr);
0225 prof_->Fill(4, allMatched);
0226 prof_->Fill(5, allTracks);
0227 prof_->Fill(6, tpPtrs.size());
0228 prof_->Fill(7, tpPtrsSelection.size());
0229 prof_->Fill(8, tpPtrsLost.size());
0230 prof_->Fill(10, tpPtrsPerfect.size());
0231 nEvents_++;
0232 }
0233
0234 void AnalyzerKFin::endJob() {
0235 if (nEvents_ == 0)
0236 return;
0237
0238 const double totalTPs = prof_->GetBinContent(9);
0239 const double numStubs = prof_->GetBinContent(1);
0240 const double numTracks = prof_->GetBinContent(2);
0241 const double numTracksLost = prof_->GetBinContent(3);
0242 const double totalTracks = prof_->GetBinContent(5);
0243 const double numTracksMatched = prof_->GetBinContent(4);
0244 const double numTPsAll = prof_->GetBinContent(6);
0245 const double numTPsEff = prof_->GetBinContent(7);
0246 const double numTPsLost = prof_->GetBinContent(8);
0247 const double numTPsEffPerfect = prof_->GetBinContent(10);
0248 const double errStubs = prof_->GetBinError(1);
0249 const double errTracks = prof_->GetBinError(2);
0250 const double errTracksLost = prof_->GetBinError(3);
0251 const double fracFake = (totalTracks - numTracksMatched) / totalTracks;
0252 const double fracDup = (numTracksMatched - numTPsAll) / totalTracks;
0253 const double eff = numTPsEff / totalTPs;
0254 const double errEff = sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0255 const double effLoss = numTPsLost / totalTPs;
0256 const double errEffLoss = sqrt(effLoss * (1. - effLoss) / totalTPs / nEvents_);
0257 const double effPerfect = numTPsEffPerfect / totalTPs;
0258 const double errEffPerfect = sqrt(effPerfect * (1. - effPerfect) / totalTPs / nEvents_);
0259 const vector<double> nums = {numStubs, numTracks, numTracksLost};
0260 const vector<double> errs = {errStubs, errTracks, errTracksLost};
0261 const int wNums = ceil(log10(*max_element(nums.begin(), nums.end()))) + 5;
0262 const int wErrs = ceil(log10(*max_element(errs.begin(), errs.end()))) + 5;
0263 log_ << " KFin SUMMARY " << endl;
0264 log_ << "number of stubs per TFP = " << setw(wNums) << numStubs << " +- " << setw(wErrs) << errStubs << endl;
0265 log_ << "number of tracks per TFP = " << setw(wNums) << numTracks << " +- " << setw(wErrs) << errTracks
0266 << endl;
0267 log_ << "number of lost tracks per TFP = " << setw(wNums) << numTracksLost << " +- " << setw(wErrs) << errTracksLost
0268 << endl;
0269 log_ << " current tracking efficiency = " << setw(wNums) << effPerfect << " +- " << setw(wErrs) << errEffPerfect
0270 << endl;
0271 log_ << " max tracking efficiency = " << setw(wNums) << eff << " +- " << setw(wErrs) << errEff << endl;
0272 log_ << " lost tracking efficiency = " << setw(wNums) << effLoss << " +- " << setw(wErrs) << errEffLoss << endl;
0273 log_ << " fake rate = " << setw(wNums) << fracFake << endl;
0274 log_ << " duplicate rate = " << setw(wNums) << fracDup << endl;
0275 log_ << "=============================================================";
0276 LogPrint("L1Trigger/TrackerTFP") << log_.str();
0277 }
0278
0279
0280 void AnalyzerKFin::formTracks(const StreamsTrack& streamsTrack,
0281 const StreamsStub& streamsStubs,
0282 vector<vector<TTStubRef>>& tracks,
0283 int channel) const {
0284 const int offset = channel * setup_->numLayers();
0285 const StreamTrack& streamTrack = streamsTrack[channel];
0286 const int numTracks = accumulate(streamTrack.begin(), streamTrack.end(), 0, [](int& sum, const FrameTrack& frame) {
0287 return sum += (frame.first.isNonnull() ? 1 : 0);
0288 });
0289 tracks.reserve(numTracks);
0290 for (int frame = 0; frame < (int)streamTrack.size(); frame++) {
0291 const FrameTrack& frameTrack = streamTrack[frame];
0292 if (frameTrack.first.isNull())
0293 continue;
0294 vector<TTStubRef> ttStubRefs;
0295 ttStubRefs.reserve(setup_->numLayers());
0296 for (int layer = 0; layer < setup_->numLayers(); layer++) {
0297 const FrameStub& stub = streamsStubs[offset + layer][frame];
0298 if (stub.first.isNonnull())
0299 ttStubRefs.push_back(stub.first);
0300 }
0301 tracks.push_back(ttStubRefs);
0302 }
0303 }
0304
0305
0306 void AnalyzerKFin::associate(const vector<vector<TTStubRef>>& tracks,
0307 const StubAssociation* ass,
0308 set<TPPtr>& tps,
0309 int& sum,
0310 bool perfect) const {
0311 for (const vector<TTStubRef>& ttStubRefs : tracks) {
0312 const vector<TPPtr>& tpPtrs = perfect ? ass->associateFinal(ttStubRefs) : ass->associate(ttStubRefs);
0313 if (tpPtrs.empty())
0314 continue;
0315 sum++;
0316 copy(tpPtrs.begin(), tpPtrs.end(), inserter(tps, tps.begin()));
0317 }
0318 }
0319
0320 }
0321
0322 DEFINE_FWK_MODULE(trklet::AnalyzerKFin);