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