File indexing completed on 2024-04-06 12:22:06
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/TrackerTFP/interface/LayerEncoding.h"
0019 #include "L1Trigger/TrackerTFP/interface/KalmanFilterFormats.h"
0020
0021 #include <TProfile.h>
0022 #include <TH1F.h>
0023
0024 #include <vector>
0025 #include <deque>
0026 #include <set>
0027 #include <cmath>
0028 #include <numeric>
0029 #include <sstream>
0030
0031 using namespace std;
0032 using namespace edm;
0033 using namespace tt;
0034 using namespace trackerTFP;
0035
0036 namespace trklet {
0037
0038
0039
0040
0041
0042
0043 class AnalyzerKFout : public one::EDAnalyzer<one::WatchRuns, one::SharedResources> {
0044 public:
0045 AnalyzerKFout(const ParameterSet& iConfig);
0046 void beginJob() override {}
0047 void beginRun(const Run& iEvent, const EventSetup& iSetup) override;
0048 void analyze(const Event& iEvent, const EventSetup& iSetup) override;
0049 void endRun(const Run& iEvent, const EventSetup& iSetup) override {}
0050 void endJob() override;
0051
0052 private:
0053
0054 void associate(const set<TTTrackRef>& ttTracks, const StubAssociation* ass, set<TPPtr>& tps, int& sum) const;
0055
0056
0057 EDGetTokenT<StreamsTrack> edGetTokenAccepted_;
0058
0059 EDGetTokenT<StreamsTrack> edGetTokenLost_;
0060
0061 EDGetTokenT<StubAssociation> edGetTokenSelection_;
0062
0063 EDGetTokenT<StubAssociation> edGetTokenReconstructable_;
0064
0065 ESGetToken<Setup, SetupRcd> esGetTokenSetup_;
0066
0067 ESGetToken<DataFormats, DataFormatsRcd> esGetTokenDataFormats_;
0068
0069 const Setup* setup_ = nullptr;
0070
0071 bool useMCTruth_;
0072
0073 int nEvents_ = 0;
0074
0075
0076
0077 TProfile* prof_;
0078 TProfile* profChannel_;
0079 TH1F* hisChannel_;
0080
0081
0082 stringstream log_;
0083 };
0084
0085 AnalyzerKFout::AnalyzerKFout(const ParameterSet& iConfig) : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")) {
0086 usesResource("TFileService");
0087
0088 const string& label = iConfig.getParameter<string>("LabelKFout");
0089 const string& branchAccepted = iConfig.getParameter<string>("BranchAcceptedTracks");
0090 const string& branchLost = iConfig.getParameter<string>("BranchLostTracks");
0091 edGetTokenAccepted_ = consumes<StreamsTrack>(InputTag(label, branchAccepted));
0092 edGetTokenLost_ = consumes<StreamsTrack>(InputTag(label, branchLost));
0093 if (useMCTruth_) {
0094 const auto& inputTagSelecttion = iConfig.getParameter<InputTag>("InputTagSelection");
0095 const auto& inputTagReconstructable = iConfig.getParameter<InputTag>("InputTagReconstructable");
0096 edGetTokenSelection_ = consumes<StubAssociation>(inputTagSelecttion);
0097 edGetTokenReconstructable_ = consumes<StubAssociation>(inputTagReconstructable);
0098 }
0099
0100 esGetTokenSetup_ = esConsumes<Setup, SetupRcd, Transition::BeginRun>();
0101
0102 log_.setf(ios::fixed, ios::floatfield);
0103 log_.precision(4);
0104 }
0105
0106 void AnalyzerKFout::beginRun(const Run& iEvent, const EventSetup& iSetup) {
0107
0108 setup_ = &iSetup.getData(esGetTokenSetup_);
0109
0110 Service<TFileService> fs;
0111 TFileDirectory dir;
0112 dir = fs->mkdir("KFout");
0113 prof_ = dir.make<TProfile>("Counts", ";", 9, 0.5, 9.5);
0114 prof_->GetXaxis()->SetBinLabel(1, "Stubs");
0115 prof_->GetXaxis()->SetBinLabel(2, "Tracks");
0116 prof_->GetXaxis()->SetBinLabel(3, "Lost Tracks");
0117 prof_->GetXaxis()->SetBinLabel(4, "Matched Tracks");
0118 prof_->GetXaxis()->SetBinLabel(5, "All Tracks");
0119 prof_->GetXaxis()->SetBinLabel(6, "Found TPs");
0120 prof_->GetXaxis()->SetBinLabel(7, "Found selected TPs");
0121 prof_->GetXaxis()->SetBinLabel(8, "Lost TPs");
0122 prof_->GetXaxis()->SetBinLabel(9, "All TPs");
0123
0124 constexpr int maxOcc = 180;
0125 const int numChannels = setup_->tfpNumChannel();
0126 hisChannel_ = dir.make<TH1F>("His Channel Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0127 profChannel_ = dir.make<TProfile>("Prof Channel Occupancy", ";", numChannels, -.5, numChannels - .5);
0128 }
0129
0130 void AnalyzerKFout::analyze(const Event& iEvent, const EventSetup& iSetup) {
0131
0132 Handle<StreamsTrack> handleAccepted;
0133 iEvent.getByToken<StreamsTrack>(edGetTokenAccepted_, handleAccepted);
0134 Handle<StreamsTrack> handleLost;
0135 iEvent.getByToken<StreamsTrack>(edGetTokenLost_, handleLost);
0136
0137 const StubAssociation* selection = nullptr;
0138 const StubAssociation* reconstructable = nullptr;
0139 if (useMCTruth_) {
0140 Handle<StubAssociation> handleSelection;
0141 iEvent.getByToken<StubAssociation>(edGetTokenSelection_, handleSelection);
0142 selection = handleSelection.product();
0143 prof_->Fill(9, selection->numTPs());
0144 Handle<StubAssociation> handleReconstructable;
0145 iEvent.getByToken<StubAssociation>(edGetTokenReconstructable_, handleReconstructable);
0146 reconstructable = handleReconstructable.product();
0147 }
0148
0149 set<TPPtr> tpPtrs;
0150 set<TPPtr> tpPtrsSelection;
0151 set<TPPtr> tpPtrsLost;
0152 int allMatched(0);
0153 int allTracks(0);
0154 for (int region = 0; region < setup_->numRegions(); region++) {
0155 int nStubsRegion(0);
0156 int nTracksRegion(0);
0157 int nLostRegion(0);
0158 for (int channel = 0; channel < setup_->tfpNumChannel(); channel++) {
0159 const int index = region * setup_->tfpNumChannel() + channel;
0160 const StreamTrack& accepted = handleAccepted->at(index);
0161 const StreamTrack& lost = handleLost->at(index);
0162 hisChannel_->Fill(accepted.size());
0163 profChannel_->Fill(channel, accepted.size());
0164 set<TTTrackRef> tracks;
0165 for (const FrameTrack& frame : accepted)
0166 if (frame.first.isNonnull())
0167 tracks.insert(frame.first);
0168 nTracksRegion += tracks.size();
0169 nStubsRegion += accumulate(tracks.begin(), tracks.end(), 0, [](int sum, const TTTrackRef& ttTrackRef) {
0170 return sum + (int)ttTrackRef->getStubRefs().size();
0171 });
0172 set<TTTrackRef> tracksLost;
0173 for (const FrameTrack& frame : lost)
0174 if (frame.first.isNonnull())
0175 tracksLost.insert(frame.first);
0176 nLostRegion += tracksLost.size();
0177 allTracks += tracks.size();
0178 int tmp(0);
0179 associate(tracks, selection, tpPtrsSelection, tmp);
0180 associate(tracksLost, selection, tpPtrsLost, tmp);
0181 associate(tracks, reconstructable, tpPtrs, allMatched);
0182 }
0183 prof_->Fill(1, nStubsRegion);
0184 prof_->Fill(2, nTracksRegion);
0185 prof_->Fill(3, nLostRegion);
0186 }
0187 deque<TPPtr> tpPtrsRealLost;
0188 set_difference(tpPtrsLost.begin(), tpPtrsLost.end(), tpPtrs.begin(), tpPtrs.end(), back_inserter(tpPtrsRealLost));
0189 prof_->Fill(4, allMatched);
0190 prof_->Fill(5, allTracks);
0191 prof_->Fill(6, tpPtrs.size());
0192 prof_->Fill(7, tpPtrsSelection.size());
0193 prof_->Fill(8, tpPtrsRealLost.size());
0194 nEvents_++;
0195 }
0196
0197 void AnalyzerKFout::endJob() {
0198 if (nEvents_ == 0)
0199 return;
0200
0201 const double totalTPs = prof_->GetBinContent(9);
0202 const double numStubs = prof_->GetBinContent(1);
0203 const double numTracks = prof_->GetBinContent(2);
0204 const double numTracksLost = prof_->GetBinContent(3);
0205 const double totalTracks = prof_->GetBinContent(5);
0206 const double numTracksMatched = prof_->GetBinContent(4);
0207 const double numTPsAll = prof_->GetBinContent(6);
0208 const double numTPsEff = prof_->GetBinContent(7);
0209 const double numTPsLost = prof_->GetBinContent(8);
0210 const double errStubs = prof_->GetBinError(1);
0211 const double errTracks = prof_->GetBinError(2);
0212 const double errTracksLost = prof_->GetBinError(3);
0213 const double fracFake = (totalTracks - numTracksMatched) / totalTracks;
0214 const double fracDup = (numTracksMatched - numTPsAll) / totalTracks;
0215 const double eff = numTPsEff / totalTPs;
0216 const double errEff = sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0217 const double effLoss = numTPsLost / totalTPs;
0218 const double errEffLoss = sqrt(effLoss * (1. - effLoss) / totalTPs / nEvents_);
0219 const vector<double> nums = {numStubs, numTracks, numTracksLost};
0220 const vector<double> errs = {errStubs, errTracks, errTracksLost};
0221 const int wNums = ceil(log10(*max_element(nums.begin(), nums.end()))) + 5;
0222 const int wErrs = ceil(log10(*max_element(errs.begin(), errs.end()))) + 5;
0223 log_ << " KFout SUMMARY " << endl;
0224
0225 log_ << "number of tracks per TFP = " << setw(wNums) << numTracks << " +- " << setw(wErrs) << errTracks
0226 << endl;
0227 log_ << "number of lost tracks per TFP = " << setw(wNums) << numTracksLost << " +- " << setw(wErrs) << errTracksLost
0228 << endl;
0229 log_ << " tracking efficiency = " << setw(wNums) << eff << " +- " << setw(wErrs) << errEff << endl;
0230 log_ << " lost tracking efficiency = " << setw(wNums) << effLoss << " +- " << setw(wErrs) << errEffLoss << endl;
0231 log_ << " fake rate = " << setw(wNums) << fracFake << endl;
0232 log_ << " duplicate rate = " << setw(wNums) << fracDup << endl;
0233 log_ << "=============================================================";
0234 LogPrint("L1Trigger/TrackerTFP") << log_.str();
0235 }
0236
0237
0238 void AnalyzerKFout::associate(const set<TTTrackRef>& ttTracks,
0239 const StubAssociation* ass,
0240 set<TPPtr>& tps,
0241 int& sum) const {
0242 for (const TTTrackRef& ttTrack : ttTracks) {
0243 const vector<TTStubRef>& ttStubRefs = ttTrack->getStubRefs();
0244 const vector<TPPtr>& tpPtrs = ass->associate(ttStubRefs);
0245 if (tpPtrs.empty())
0246 continue;
0247 sum++;
0248 copy(tpPtrs.begin(), tpPtrs.end(), inserter(tps, tps.begin()));
0249 }
0250 }
0251
0252 }
0253
0254 DEFINE_FWK_MODULE(trklet::AnalyzerKFout);