File indexing completed on 2022-10-14 01:44:14
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 #include "L1Trigger/TrackFindingTracklet/interface/Settings.h"
0020
0021 #include <TProfile.h>
0022 #include <TProfile2D.h>
0023 #include <TH1F.h>
0024 #include <TH2F.h>
0025
0026 #include <vector>
0027 #include <deque>
0028 #include <set>
0029 #include <cmath>
0030 #include <numeric>
0031 #include <sstream>
0032
0033 using namespace std;
0034 using namespace edm;
0035 using namespace trackerTFP;
0036 using namespace tt;
0037
0038 namespace trklet {
0039
0040
0041 enum Resolution { Phi, Z, NumResolution };
0042 constexpr initializer_list<Resolution> AllResolution = {Phi, Z};
0043 constexpr auto NameResolution = {"Phi", "Z"};
0044 inline string name(Resolution r) { return string(*(NameResolution.begin() + r)); }
0045
0046
0047
0048
0049
0050
0051 class AnalyzerTBout : public one::EDAnalyzer<one::WatchRuns, one::SharedResources> {
0052 public:
0053 AnalyzerTBout(const ParameterSet& iConfig);
0054 void beginJob() override {}
0055 void beginRun(const Run& iEvent, const EventSetup& iSetup) override;
0056 void analyze(const Event& iEvent, const EventSetup& iSetup) override;
0057 void endRun(const Run& iEvent, const EventSetup& iSetup) override {}
0058 void endJob() override;
0059
0060 private:
0061
0062 void formTracks(const StreamsTrack& streamsTrack,
0063 const StreamsStub& streamsStubs,
0064 vector<vector<TTStubRef>>& tracks,
0065 int channel);
0066
0067 void associate(const vector<vector<TTStubRef>>& tracks, const StubAssociation* ass, set<TPPtr>& tps, int& sum) const;
0068
0069 void fill(const FrameTrack& frameTrack, const FrameStub& frameStub);
0070
0071
0072 EDGetTokenT<TTDTC> edGetTokenTTDTC_;
0073
0074 EDGetTokenT<StreamsStub> edGetTokenAcceptedStubs_;
0075
0076 EDGetTokenT<StreamsTrack> edGetTokenAcceptedTracks_;
0077
0078 EDGetTokenT<StreamsStub> edGetTokenLostStubs_;
0079
0080 EDGetTokenT<StreamsTrack> edGetTokenLostTracks_;
0081
0082 EDGetTokenT<StubAssociation> edGetTokenSelection_;
0083
0084 EDGetTokenT<StubAssociation> edGetTokenReconstructable_;
0085
0086 ESGetToken<Setup, SetupRcd> esGetTokenSetup_;
0087
0088 ESGetToken<DataFormats, DataFormatsRcd> esGetTokenDataFormats_;
0089
0090 ESGetToken<ChannelAssignment, ChannelAssignmentRcd> esGetTokenChannelAssignment_;
0091
0092 const Setup* setup_ = nullptr;
0093
0094 const Settings settings_;
0095
0096 const DataFormats* dataFormats_ = nullptr;
0097
0098 const ChannelAssignment* channelAssignment_ = nullptr;
0099
0100 bool useMCTruth_;
0101
0102 int nEvents_ = 0;
0103
0104 vector<deque<FrameStub>> regionStubs_;
0105
0106 int region_;
0107
0108 vector<int> nOverflows_;
0109
0110
0111
0112 TProfile* prof_;
0113 TProfile* profChannel_;
0114 TH1F* hisChannel_;
0115 vector<TH1F*> hisResolution_;
0116 vector<TProfile2D*> profResolution_;
0117 vector<TH1F*> hisResolutionMe_;
0118 vector<TH1F*> hisResolutionThey_;
0119 vector<TH2F*> his2Resolution_;
0120
0121
0122 stringstream log_;
0123 };
0124
0125 AnalyzerTBout::AnalyzerTBout(const ParameterSet& iConfig) : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")) {
0126 usesResource("TFileService");
0127
0128 const InputTag& inputTag = iConfig.getParameter<InputTag>("InputTagDTC");
0129 const string& label = iConfig.getParameter<string>("LabelTBout");
0130 const string& branchAcceptedStubs = iConfig.getParameter<string>("BranchAcceptedStubs");
0131 const string& branchAcceptedTracks = iConfig.getParameter<string>("BranchAcceptedTracks");
0132 const string& branchLostStubs = iConfig.getParameter<string>("BranchLostStubs");
0133 const string& branchLostTracks = iConfig.getParameter<string>("BranchLostTracks");
0134 edGetTokenTTDTC_ = consumes<TTDTC>(inputTag);
0135 edGetTokenAcceptedStubs_ = consumes<StreamsStub>(InputTag(label, branchAcceptedStubs));
0136 edGetTokenAcceptedTracks_ = consumes<StreamsTrack>(InputTag(label, branchAcceptedTracks));
0137 edGetTokenLostStubs_ = consumes<StreamsStub>(InputTag(label, branchLostStubs));
0138 edGetTokenLostTracks_ = consumes<StreamsTrack>(InputTag(label, branchLostTracks));
0139 if (useMCTruth_) {
0140 const auto& inputTagSelecttion = iConfig.getParameter<InputTag>("InputTagSelection");
0141 const auto& inputTagReconstructable = iConfig.getParameter<InputTag>("InputTagReconstructable");
0142 edGetTokenSelection_ = consumes<StubAssociation>(inputTagSelecttion);
0143 edGetTokenReconstructable_ = consumes<StubAssociation>(inputTagReconstructable);
0144 }
0145
0146 esGetTokenSetup_ = esConsumes<Setup, SetupRcd, Transition::BeginRun>();
0147 esGetTokenDataFormats_ = esConsumes<DataFormats, DataFormatsRcd, Transition::BeginRun>();
0148 esGetTokenChannelAssignment_ = esConsumes<ChannelAssignment, ChannelAssignmentRcd, Transition::BeginRun>();
0149
0150 nOverflows_ = vector<int>(2, 0);
0151
0152 log_.setf(ios::fixed, ios::floatfield);
0153 log_.precision(4);
0154 }
0155
0156 void AnalyzerTBout::beginRun(const Run& iEvent, const EventSetup& iSetup) {
0157
0158 setup_ = &iSetup.getData(esGetTokenSetup_);
0159
0160 dataFormats_ = &iSetup.getData(esGetTokenDataFormats_);
0161
0162 channelAssignment_ = &iSetup.getData(esGetTokenChannelAssignment_);
0163
0164 Service<TFileService> fs;
0165 TFileDirectory dir;
0166 dir = fs->mkdir("TBout");
0167 prof_ = dir.make<TProfile>("Counts", ";", 9, 0.5, 9.5);
0168 prof_->GetXaxis()->SetBinLabel(1, "Stubs");
0169 prof_->GetXaxis()->SetBinLabel(2, "Tracks");
0170 prof_->GetXaxis()->SetBinLabel(3, "Lost Tracks");
0171 prof_->GetXaxis()->SetBinLabel(4, "Matched Tracks");
0172 prof_->GetXaxis()->SetBinLabel(5, "All Tracks");
0173 prof_->GetXaxis()->SetBinLabel(6, "Found TPs");
0174 prof_->GetXaxis()->SetBinLabel(7, "Found selected TPs");
0175 prof_->GetXaxis()->SetBinLabel(8, "Lost TPs");
0176 prof_->GetXaxis()->SetBinLabel(9, "All TPs");
0177
0178 constexpr int maxOcc = 180;
0179 const int numChannels = channelAssignment_->numChannelsStub() * setup_->numRegions();
0180 hisChannel_ = dir.make<TH1F>("His Channel Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0181 profChannel_ = dir.make<TProfile>("Prof Channel Occupancy", ";", numChannels, -.5, numChannels - .5);
0182
0183 constexpr int bins = 400;
0184 constexpr int binsHis = 100;
0185 constexpr double maxZ = 300.;
0186 constexpr double maxR = 120.;
0187 constexpr array<double, NumResolution> ranges{{.01, 5.}};
0188 hisResolution_.reserve(NumResolution);
0189 profResolution_.reserve(NumResolution);
0190 for (Resolution r : AllResolution) {
0191 hisResolution_.emplace_back(dir.make<TH1F>(("HisRes" + name(r)).c_str(), ";", binsHis, -ranges[r], ranges[r]));
0192 profResolution_.emplace_back(
0193 dir.make<TProfile2D>(("ProfRes" + name(r)).c_str(), ";;", bins, -maxZ, maxZ, bins, 0., maxR));
0194 hisResolutionMe_.emplace_back(
0195 dir.make<TH1F>(("HisResMe" + name(r)).c_str(), ";", binsHis, -ranges[r], ranges[r]));
0196 hisResolutionThey_.emplace_back(
0197 dir.make<TH1F>(("HisResThey" + name(r)).c_str(), ";", binsHis, -ranges[r], ranges[r]));
0198 his2Resolution_.emplace_back(
0199 dir.make<TH2F>(("His2" + name(r)).c_str(), ";;", bins, -ranges[r], ranges[r], bins, -ranges[r], ranges[r]));
0200 }
0201 regionStubs_ = vector<deque<FrameStub>>(setup_->numRegions());
0202 }
0203
0204 void AnalyzerTBout::analyze(const Event& iEvent, const EventSetup& iSetup) {
0205
0206 Handle<TTDTC> handleTTDTC;
0207 iEvent.getByToken<TTDTC>(edGetTokenTTDTC_, handleTTDTC);
0208 const TTDTC& ttDTC = *handleTTDTC;
0209 for (deque<FrameStub>& region : regionStubs_)
0210 region.clear();
0211 for (int r : ttDTC.tfpRegions()) {
0212 for (int c : ttDTC.tfpChannels()) {
0213 const StreamStub& s = ttDTC.stream(r, c);
0214 copy_if(
0215 s.begin(), s.end(), back_inserter(regionStubs_[r]), [](const FrameStub& f) { return f.first.isNonnull(); });
0216 }
0217 }
0218
0219 Handle<StreamsStub> handleAcceptedStubs;
0220 iEvent.getByToken<StreamsStub>(edGetTokenAcceptedStubs_, handleAcceptedStubs);
0221 const StreamsStub& acceptedStubs = *handleAcceptedStubs;
0222 Handle<StreamsTrack> handleAcceptedTracks;
0223 iEvent.getByToken<StreamsTrack>(edGetTokenAcceptedTracks_, handleAcceptedTracks);
0224 const StreamsTrack& acceptedTracks = *handleAcceptedTracks;
0225 Handle<StreamsStub> handleLostStubs;
0226 iEvent.getByToken<StreamsStub>(edGetTokenLostStubs_, handleLostStubs);
0227 const StreamsStub& lostStubs = *handleLostStubs;
0228 Handle<StreamsTrack> handleLostTracks;
0229 iEvent.getByToken<StreamsTrack>(edGetTokenLostTracks_, handleLostTracks);
0230 const StreamsTrack& lostTracks = *handleLostTracks;
0231
0232 const StubAssociation* selection = nullptr;
0233 const StubAssociation* reconstructable = nullptr;
0234 if (useMCTruth_) {
0235 Handle<StubAssociation> handleSelection;
0236 iEvent.getByToken<StubAssociation>(edGetTokenSelection_, handleSelection);
0237 selection = handleSelection.product();
0238 prof_->Fill(9, selection->numTPs());
0239 Handle<StubAssociation> handleReconstructable;
0240 iEvent.getByToken<StubAssociation>(edGetTokenReconstructable_, handleReconstructable);
0241 reconstructable = handleReconstructable.product();
0242 }
0243
0244 set<TPPtr> tpPtrs;
0245 set<TPPtr> tpPtrsSelection;
0246 set<TPPtr> tpPtrsLost;
0247 int allMatched(0);
0248 int allTracks(0);
0249 for (region_ = 0; region_ < setup_->numRegions(); region_++) {
0250 const int offset = region_ * channelAssignment_->numChannelsTrack();
0251 int nStubs(0);
0252 int nTracks(0);
0253 int nLost(0);
0254 for (int channel = 0; channel < channelAssignment_->numChannelsTrack(); channel++) {
0255 vector<vector<TTStubRef>> tracks;
0256 formTracks(acceptedTracks, acceptedStubs, tracks, offset + channel);
0257 vector<vector<TTStubRef>> lost;
0258 formTracks(lostTracks, lostStubs, lost, offset + channel);
0259 nTracks += tracks.size();
0260 nStubs +=
0261 accumulate(tracks.begin(), tracks.end(), 0, [](int& sum, const auto& v) { return sum += (int)v.size(); });
0262 nLost += lost.size();
0263 allTracks += tracks.size();
0264 if (!useMCTruth_)
0265 continue;
0266 int tmp(0);
0267 associate(tracks, selection, tpPtrsSelection, tmp);
0268 associate(lost, selection, tpPtrsLost, tmp);
0269 associate(tracks, reconstructable, tpPtrs, allMatched);
0270 }
0271 prof_->Fill(1, nStubs);
0272 prof_->Fill(2, nTracks);
0273 prof_->Fill(3, nLost);
0274 }
0275 vector<TPPtr> recovered;
0276 recovered.reserve(tpPtrsLost.size());
0277 set_intersection(tpPtrsLost.begin(), tpPtrsLost.end(), tpPtrs.begin(), tpPtrs.end(), back_inserter(recovered));
0278 for (const TPPtr& tpPtr : recovered)
0279 tpPtrsLost.erase(tpPtr);
0280 prof_->Fill(4, allMatched);
0281 prof_->Fill(5, allTracks);
0282 prof_->Fill(6, tpPtrs.size());
0283 prof_->Fill(7, tpPtrsSelection.size());
0284 prof_->Fill(8, tpPtrsLost.size());
0285 nEvents_++;
0286 }
0287
0288 void AnalyzerTBout::endJob() {
0289 if (nEvents_ == 0)
0290 return;
0291
0292 const double totalTPs = prof_->GetBinContent(9);
0293 const double numStubs = prof_->GetBinContent(1);
0294 const double numTracks = prof_->GetBinContent(2);
0295 const double numTracksLost = prof_->GetBinContent(3);
0296 const double totalTracks = prof_->GetBinContent(5);
0297 const double numTracksMatched = prof_->GetBinContent(4);
0298 const double numTPsAll = prof_->GetBinContent(6);
0299 const double numTPsEff = prof_->GetBinContent(7);
0300 const double numTPsLost = prof_->GetBinContent(8);
0301 const double errStubs = prof_->GetBinError(1);
0302 const double errTracks = prof_->GetBinError(2);
0303 const double errTracksLost = prof_->GetBinError(3);
0304 const double fracFake = (totalTracks - numTracksMatched) / totalTracks;
0305 const double fracDup = (numTracksMatched - numTPsAll) / totalTracks;
0306 const double eff = numTPsEff / totalTPs;
0307 const double errEff = sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0308 const double effLoss = numTPsLost / totalTPs;
0309 const double errEffLoss = sqrt(effLoss * (1. - effLoss) / totalTPs / nEvents_);
0310 const vector<double> nums = {numStubs, numTracks, numTracksLost};
0311 const vector<double> errs = {errStubs, errTracks, errTracksLost};
0312 const int wNums = ceil(log10(*max_element(nums.begin(), nums.end()))) + 5;
0313 const int wErrs = ceil(log10(*max_element(errs.begin(), errs.end()))) + 5;
0314 log_ << " TBout SUMMARY " << endl;
0315 log_ << "number of stubs per TFP = " << setw(wNums) << numStubs << " +- " << setw(wErrs) << errStubs << endl;
0316 log_ << "number of tracks per TFP = " << setw(wNums) << numTracks << " +- " << setw(wErrs) << errTracks
0317 << endl;
0318 log_ << "number of lost tracks per TFP = " << setw(wNums) << numTracksLost << " +- " << setw(wErrs) << errTracksLost
0319 << endl;
0320 log_ << " max tracking efficiency = " << setw(wNums) << eff << " +- " << setw(wErrs) << errEff << endl;
0321 log_ << " lost tracking efficiency = " << setw(wNums) << effLoss << " +- " << setw(wErrs) << errEffLoss << endl;
0322 log_ << " fake rate = " << setw(wNums) << fracFake << endl;
0323 log_ << " duplicate rate = " << setw(wNums) << fracDup << endl;
0324 log_ << "number of overflowed phi residuals: " << nOverflows_[0] << " and z: " << nOverflows_[1] << endl;
0325 log_ << "=============================================================";
0326 LogPrint("L1Trigger/TrackerTFP") << log_.str();
0327 }
0328
0329
0330 void AnalyzerTBout::formTracks(const StreamsTrack& streamsTrack,
0331 const StreamsStub& streamsStubs,
0332 vector<vector<TTStubRef>>& tracks,
0333 int channel) {
0334 const int seedType = channel % channelAssignment_->numChannelsTrack();
0335 const int offset = channelAssignment_->offsetStub(channel);
0336 const StreamTrack& streamTrack = streamsTrack[channel];
0337 const int numTracks = accumulate(streamTrack.begin(), streamTrack.end(), 0, [](int& sum, const FrameTrack& frame) {
0338 return sum += (frame.first.isNonnull() ? 1 : 0);
0339 });
0340 tracks.reserve(numTracks);
0341 for (int frame = 0; frame < (int)streamTrack.size(); frame++) {
0342 const FrameTrack& frameTrack = streamTrack[frame];
0343 if (frameTrack.first.isNull())
0344 continue;
0345 vector<TTStubRef> ttStubRefs;
0346 ttStubRefs.reserve(channelAssignment_->numProjectionLayers(seedType) +
0347 channelAssignment_->seedingLayers(seedType).size());
0348 for (int layer = 0; layer < channelAssignment_->numProjectionLayers(seedType); layer++) {
0349 const FrameStub& stub = streamsStubs[offset + layer][frame];
0350 if (stub.first.isNonnull()) {
0351 this->fill(frameTrack, stub);
0352 ttStubRefs.push_back(stub.first);
0353 }
0354 }
0355 for (const TTStubRef& ttStubRef : frameTrack.first->getStubRefs()) {
0356 int layer(0);
0357 if (!channelAssignment_->layerId(seedType, ttStubRef, layer))
0358 ttStubRefs.push_back(ttStubRef);
0359 }
0360 tracks.push_back(ttStubRefs);
0361 }
0362 }
0363
0364
0365 void AnalyzerTBout::associate(const vector<vector<TTStubRef>>& tracks,
0366 const StubAssociation* ass,
0367 set<TPPtr>& tps,
0368 int& sum) const {
0369 for (const vector<TTStubRef>& ttStubRefs : tracks) {
0370 const vector<TPPtr>& tpPtrs = ass->associate(ttStubRefs);
0371 if (tpPtrs.empty())
0372 continue;
0373 sum++;
0374 copy(tpPtrs.begin(), tpPtrs.end(), inserter(tps, tps.begin()));
0375 }
0376 }
0377
0378
0379 void AnalyzerTBout::fill(const FrameTrack& frameTrack, const FrameStub& frameStub) {
0380
0381 const deque<FrameStub>& region = regionStubs_[region_];
0382 const auto it =
0383 find_if(region.begin(), region.end(), [&frameStub](const FrameStub& f) { return f.first == frameStub.first; });
0384 if (it == region.end())
0385 throw cms::Exception("LgociError.") << "Stub on track was not in DTC collection.";
0386 const GlobalPoint ttPos = setup_->stubPos(frameStub.first);
0387 const GlobalPoint pos = setup_->stubPos(true, *it, region_);
0388 static constexpr int widthPhi = 12;
0389 static constexpr int widthZ = 9;
0390 static constexpr int widthR = 7;
0391 const bool barrel = setup_->barrel(frameStub.first);
0392 const int layerIdTracklet = setup_->trackletLayerId(frameStub.first);
0393 static const double baseR = settings_.kz();
0394 const double basePhi = barrel ? settings_.kphi1() : settings_.kphi(layerIdTracklet);
0395 const double baseZ = settings_.kz(layerIdTracklet);
0396 static const double baseInvR = settings_.kphi1() / settings_.kr() * pow(2, settings_.rinv_shift());
0397 static const double basePhi0 = settings_.kphi1() * pow(2, settings_.phi0_shift());
0398 static const double baseZ0 = settings_.kz() * pow(2, settings_.z0_shift());
0399 static const double baseTanL = settings_.kz() / settings_.kr() * pow(2, settings_.t_shift());
0400 const int widthRZ = barrel ? widthZ : widthR;
0401 const double baseRZ = barrel ? baseZ : baseR;
0402
0403 const double rInv = (frameTrack.first->rInv() / baseInvR + .5) * baseInvR;
0404 const double phi0 = (frameTrack.first->phi() / basePhi0 + .5) * basePhi0;
0405 const double z0 = (frameTrack.first->z0() / baseZ0 + .5) * baseZ0;
0406 const double tanL = (frameTrack.first->tanL() / baseTanL + .5) * baseTanL;
0407 const double phi = deltaPhi(pos.phi() - (phi0 - rInv * pos.perp() / 2.));
0408 const double r = pos.perp() - (pos.z() - z0) / tanL;
0409 const double z = pos.z() - (z0 + tanL * pos.perp());
0410 const double rz = barrel ? z : r;
0411 const int phii = floor(phi / basePhi);
0412 const int rzi = floor(rz / baseRZ);
0413 const double phid = (phii + .5) * basePhi;
0414 const double rzd = (rzi + .5) * baseRZ;
0415
0416 TTBV hw(frameStub.second);
0417 const TTBV hwRZ(hw, widthRZ, 0, true);
0418 hw >>= widthRZ;
0419 const TTBV hwPhi(hw, widthPhi, 0, true);
0420 bool overflowPhi = abs(phii - hwPhi.val()) > pow(2, widthPhi) * 7. / 8.;
0421 bool overflowZ = abs(rzi - hwRZ.val()) > pow(2, widthRZ) * 7. / 8.;
0422 const double hwPhid = hwPhi.val(basePhi);
0423 const double hwRZd = hwRZ.val(baseRZ);
0424 const vector<double> resolutions = {phid - hwPhid, rzd - hwRZd};
0425 const vector<bool> overflows = {overflowPhi, overflowZ};
0426 for (Resolution r : AllResolution) {
0427 if (overflows[r]) {
0428 cout << rzi << " " << hwRZ.val() << " " << barrel << " " << setup_->psModule(frameStub.first) << " "
0429 << frameStub.second << endl;
0430 nOverflows_[r]++;
0431 continue;
0432 }
0433 hisResolution_[r]->Fill(resolutions[r]);
0434 profResolution_[r]->Fill(ttPos.z(), ttPos.perp(), abs(resolutions[r]));
0435 }
0436 hisResolutionMe_[0]->Fill(phid);
0437 hisResolutionMe_[1]->Fill(rzd);
0438 hisResolutionThey_[0]->Fill(hwPhid);
0439 hisResolutionThey_[1]->Fill(hwRZd);
0440 his2Resolution_[0]->Fill(phid, hwPhid);
0441 his2Resolution_[1]->Fill(rzd, hwRZd);
0442 }
0443
0444 }
0445
0446 DEFINE_FWK_MODULE(trklet::AnalyzerTBout);