TrackerDTC/test/Analyzer.cc

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/DetId/interface/DetId.h"
0014 #include "DataFormats/Common/interface/Ptr.h"
0015 #include "DataFormats/Common/interface/Handle.h"
0016 #include "DataFormats/L1TrackTrigger/interface/TTDTC.h"
0017 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0018 #include "DataFormats/GeometrySurface/interface/Plane.h"
0019 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0020
0021 #include "SimTracker/TrackTriggerAssociation/interface/StubAssociation.h"
0022 #include "L1Trigger/TrackTrigger/interface/Setup.h"
0023 #include "L1Trigger/TrackerDTC/interface/LayerEncoding.h"
0024
0025 #include <TProfile.h>
0026 #include <TProfile2D.h>
0027 #include <TH1F.h>
0028 #include <TH2F.h>
0029 #include <TEfficiency.h>
0030
0031 #include <vector>
0032 #include <map>
0033 #include <utility>
0034 #include <set>
0035 #include <algorithm>
0036 #include <cmath>
0037 #include <numeric>
0038 #include <array>
0039 #include <initializer_list>
0040 #include <sstream>
0041
0042 namespace trackerDTC {
0043
0044   // stub resolution plots helper
0045   enum Resolution { R, Phi, Z, NumResolution };
0046   constexpr std::initializer_list<Resolution> AllResolution = {R, Phi, Z};
0047   constexpr auto NameResolution = {"R", "Phi", "Z"};
0048   inline std::string name(Resolution r) { return std::string(*(NameResolution.begin() + r)); }
0049   // max tracking efficiency plots helper
0050   enum Efficiency { Phi0, Pt, InvPt, D0, Z0, Eta, NumEfficiency };
0051   constexpr std::initializer_list<Efficiency> AllEfficiency = {Phi0, Pt, InvPt, D0, Z0, Eta};
0052   constexpr auto NameEfficiency = {"Phi0", "Pt", "InvPt", "D0", "Z0", "Eta"};
0053   inline std::string name(Efficiency e) { return std::string(*(NameEfficiency.begin() + e)); }
0054
0055   /*! \class  trackerDTC::Analyzer
0056    *  \brief  Class to analyze hardware like structured TTStub Collection used by Track Trigger emulators, runs DTC stub emulation, plots performance & stub occupancy
0057    *  \author Thomas Schuh
0058    *  \date   2020, Apr
0059    */
0060   class Analyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0061   public:
0062     Analyzer(const edm::ParameterSet& iConfig);
0063     void beginJob() override {}
0064     void beginRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) override;
0065     void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0066     void endRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) override {}
0067     void endJob() override;
0068
0069   private:
0070     // fills kinematic tp histograms
0071     void fill(const TPPtr& tpPtr, const std::vector<TH1F*> th1fs) const;
0072     // fill stub related histograms
0073     void fill(const tt::StreamStub& stream, int region, int channel, int& sum, TH2F* th2f);
0074     // prints out MC summary
0075     void endJobMC();
0076     // prints out DTC summary
0077     void endJobDTC();
0078
0079     // ED input token of DTC stubs
0080     edm::EDGetTokenT<TTDTC> edGetTokenTTDTCAccepted_;
0081     // ED input token of lost DTC stubs
0082     edm::EDGetTokenT<TTDTC> edGetTokenTTDTCLost_;
0083     // ED input token of TTStubRef to TPPtr association for tracking efficiency
0084     edm::EDGetTokenT<tt::StubAssociation> edGetTokenSelection_;
0085     // ED input token of TTStubRef to recontructable TPPtr association
0086     edm::EDGetTokenT<tt::StubAssociation> edGetTokenReconstructable_;
0087     // Setup token
0088     edm::ESGetToken<tt::Setup, tt::SetupRcd> esGetToken_;
0089     // stores, calculates and provides run-time constants
0090     const tt::Setup* setup_;
0091     // enables analyze of TPs
0092     bool useMCTruth_;
0093     //
0094     int nEvents_ = 0;
0095
0096     // Histograms
0097
0098     TProfile* profMC_;
0099     TProfile* profDTC_;
0100     TProfile* profChannel_;
0101     TH1F* hisChannel_;
0102     TH2F* hisRZStubs_;
0103     TH2F* hisRZStubsLost_;
0104     TH2F* hisRZStubsEff_;
0105     std::vector<TH1F*> hisResolution_;
0106     std::vector<TProfile2D*> profResolution_;
0107     std::vector<TH1F*> hisEff_;
0108     std::vector<TH1F*> hisEffMC_;
0109     std::vector<TEfficiency*> eff_;
0110
0111     // printout
0112     std::stringstream log_;
0113   };
0114
0115   Analyzer::Analyzer(const edm::ParameterSet& iConfig) : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")) {
0116     usesResource("TFileService");
0117     // book in- and output ED products
0118     const auto& inputTagAccepted = iConfig.getParameter<edm::InputTag>("InputTagAccepted");
0119     const auto& inputTagLost = iConfig.getParameter<edm::InputTag>("InputTagLost");
0120     edGetTokenTTDTCAccepted_ = consumes<TTDTC>(inputTagAccepted);
0121     edGetTokenTTDTCLost_ = consumes<TTDTC>(inputTagLost);
0122     if (useMCTruth_) {
0123       const auto& inputTagSelection = iConfig.getParameter<edm::InputTag>("InputTagSelection");
0124       const auto& inputTagReconstructable = iConfig.getParameter<edm::InputTag>("InputTagReconstructable");
0125       edGetTokenSelection_ = consumes<tt::StubAssociation>(inputTagSelection);
0126       edGetTokenReconstructable_ = consumes<tt::StubAssociation>(inputTagReconstructable);
0127     }
0128     // book ES product
0129     esGetToken_ = esConsumes<edm::Transition::BeginRun>();
0130     // log config
0131     log_.setf(std::ios::fixed, std::ios::floatfield);
0132     log_.precision(4);
0133   }
0134
0135   void Analyzer::beginRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) {
0136     // helper class to store configurations
0137     setup_ = &iSetup.getData(esGetToken_);
0138     // book histograms
0139     edm::Service<TFileService> fs;
0140     TFileDirectory dir;
0141     // mc
0142     dir = fs->mkdir("MC");
0143     profMC_ = dir.make<TProfile>("Counts", ";", 6, 0.5, 6.5);
0144     profMC_->GetXaxis()->SetBinLabel(1, "Stubs");
0145     profMC_->GetXaxis()->SetBinLabel(2, "Matched Stubs");
0146     profMC_->GetXaxis()->SetBinLabel(3, "reco TPs");
0147     profMC_->GetXaxis()->SetBinLabel(4, "eff TPs");
0148     profMC_->GetXaxis()->SetBinLabel(5, "total eff TPs");
0149     profMC_->GetXaxis()->SetBinLabel(6, "Cluster");
0150     constexpr std::array<int, NumEfficiency> binsEff{{9 * 8, 10, 16, 10, 30, 24}};
0151     constexpr std::array<std::pair<double, double>, NumEfficiency> rangesEff{
0152         {{-M_PI, M_PI}, {0., 100.}, {-1. / 3., 1. / 3.}, {-5., 5.}, {-15., 15.}, {-2.4, 2.4}}};
0153     if (useMCTruth_) {
0154       hisEffMC_.reserve(NumEfficiency);
0155       for (Efficiency e : AllEfficiency)
0156         hisEffMC_.emplace_back(
0157             dir.make<TH1F>(("HisTP" + name(e)).c_str(), ";", binsEff[e], rangesEff[e].first, rangesEff[e].second));
0158     }
0159     // dtc
0160     dir = fs->mkdir("DTC");
0161     profDTC_ = dir.make<TProfile>("Counts", ";", 3, 0.5, 3.5);
0162     profDTC_->GetXaxis()->SetBinLabel(1, "Stubs");
0163     profDTC_->GetXaxis()->SetBinLabel(2, "Lost Stubs");
0164     profDTC_->GetXaxis()->SetBinLabel(3, "TPs");
0165     // channel occupancy
0166     constexpr int maxOcc = 180;
0167     const int numChannels = setup_->numDTCs() * setup_->numOverlappingRegions();
0168     hisChannel_ = dir.make<TH1F>("His Channel Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0169     profChannel_ = dir.make<TProfile>("Prof Channel Occupancy", ";", numChannels, -.5, numChannels - .5);
0170     // max tracking efficiencies
0171     if (useMCTruth_) {
0172       dir = fs->mkdir("DTC/Effi");
0173       hisEff_.reserve(NumEfficiency);
0174       for (Efficiency e : AllEfficiency)
0175         hisEff_.emplace_back(
0176             dir.make<TH1F>(("HisTP" + name(e)).c_str(), ";", binsEff[e], rangesEff[e].first, rangesEff[e].second));
0177       eff_.reserve(NumEfficiency);
0178       for (Efficiency e : AllEfficiency)
0179         eff_.emplace_back(
0180             dir.make<TEfficiency>(("Eff" + name(e)).c_str(), ";", binsEff[e], rangesEff[e].first, rangesEff[e].second));
0181     }
0182     // lost stub fraction in r-z
0183     dir = fs->mkdir("DTC/Loss");
0184     constexpr int bins = 400;
0185     constexpr double maxZ = 300.;
0186     constexpr double maxR = 120.;
0187     hisRZStubs_ = dir.make<TH2F>("RZ Stubs", ";;", bins, -maxZ, maxZ, bins, 0., maxR);
0188     hisRZStubsLost_ = dir.make<TH2F>("RZ Stubs Lost", ";;", bins, -maxZ, maxZ, bins, 0., maxR);
0189     hisRZStubsEff_ = dir.make<TH2F>("RZ Stubs Eff", ";;", bins, -maxZ, maxZ, bins, 0., maxR);
0190     // stub parameter resolutions
0191     dir = fs->mkdir("DTC/Res");
0192     constexpr std::array<double, NumResolution> ranges{{.2, .0001, .5}};
0193     constexpr int binsHis = 100;
0194     hisResolution_.reserve(NumResolution);
0195     profResolution_.reserve(NumResolution);
0196     for (Resolution r : AllResolution) {
0197       hisResolution_.emplace_back(dir.make<TH1F>(("HisRes" + name(r)).c_str(), ";", binsHis, -ranges[r], ranges[r]));
0198       profResolution_.emplace_back(
0199           dir.make<TProfile2D>(("ProfRes" + name(r)).c_str(), ";;", bins, -maxZ, maxZ, bins, 0., maxR));
0200     }
0201   }
0202
0203   void Analyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0204     // read in dtc products
0205     edm::Handle<TTDTC> handleTTDTCAccepted;
0206     iEvent.getByToken<TTDTC>(edGetTokenTTDTCAccepted_, handleTTDTCAccepted);
0207     edm::Handle<TTDTC> handleTTDTCLost;
0208     iEvent.getByToken<TTDTC>(edGetTokenTTDTCLost_, handleTTDTCLost);
0209     // read in MCTruth
0210     const tt::StubAssociation* selection = nullptr;
0211     const tt::StubAssociation* reconstructable = nullptr;
0212     if (useMCTruth_) {
0213       edm::Handle<tt::StubAssociation> handleSelection;
0214       iEvent.getByToken<tt::StubAssociation>(edGetTokenSelection_, handleSelection);
0215       selection = handleSelection.product();
0216       edm::Handle<tt::StubAssociation> handleReconstructable;
0217       iEvent.getByToken<tt::StubAssociation>(edGetTokenReconstructable_, handleReconstructable);
0218       reconstructable = handleReconstructable.product();
0219       profMC_->Fill(3, reconstructable->numTPs() / (double)setup_->numRegions());
0220       profMC_->Fill(4, selection->numTPs() / (double)setup_->numRegions());
0221       profMC_->Fill(5, selection->numTPs());
0222       for (const auto& p : selection->getTrackingParticleToTTStubsMap())
0223         fill(p.first, hisEffMC_);
0224     }
0225     // analyze dtc products and find still reconstrucable TrackingParticles
0226     std::set<TPPtr> tpPtrs;
0227     for (int region = 0; region < setup_->numRegions(); region++) {
0228       int nStubs(0);
0229       int nLost(0);
0230       std::map<TPPtr, std::vector<TTStubRef>> mapTPsTTStubs;
0231       for (int channel = 0; channel < setup_->numDTCsPerTFP(); channel++) {
0232         const tt::StreamStub& accepted = handleTTDTCAccepted->stream(region, channel);
0233         const tt::StreamStub& lost = handleTTDTCLost->stream(region, channel);
0234         hisChannel_->Fill(accepted.size());
0235         profChannel_->Fill(channel, accepted.size());
0236         fill(accepted, region, channel, nStubs, hisRZStubs_);
0237         fill(lost, region, channel, nLost, hisRZStubsLost_);
0238         if (!useMCTruth_)
0239           continue;
0240         for (const tt::FrameStub& frame : accepted) {
0241           if (frame.first.isNull())
0242             continue;
0243           for (const TPPtr& tpPtr : selection->findTrackingParticlePtrs(frame.first)) {
0244             auto it = mapTPsTTStubs.find(tpPtr);
0245             if (it == mapTPsTTStubs.end()) {
0246               it = mapTPsTTStubs.emplace(tpPtr, std::vector<TTStubRef>()).first;
0247               it->second.reserve(selection->findTTStubRefs(tpPtr).size());
0248             }
0249             it->second.push_back(frame.first);
0250           }
0251         }
0252         for (const auto& p : mapTPsTTStubs)
0253           if (setup_->reconstructable(p.second))
0254             tpPtrs.insert(p.first);
0255       }
0256       profDTC_->Fill(1, nStubs);
0257       profDTC_->Fill(2, nLost);
0258     }
0259     for (const TPPtr& tpPtr : tpPtrs)
0260       fill(tpPtr, hisEff_);
0261     profDTC_->Fill(3, tpPtrs.size());
0262     nEvents_++;
0263   }
0264
0265   void Analyzer::endJob() {
0266     if (nEvents_ == 0)
0267       return;
0268     // create r-z stub fraction plot
0269     TH2F th2f("", ";;", 400, -300, 300., 400, 0., 120.);
0270     th2f.Add(hisRZStubsLost_);
0271     th2f.Add(hisRZStubs_);
0272     hisRZStubsEff_->Add(hisRZStubsLost_);
0273     hisRZStubsEff_->Divide(&th2f);
0274     // create efficieny plots
0275     if (useMCTruth_) {
0276       for (Efficiency e : AllEfficiency) {
0277         eff_[e]->SetPassedHistogram(*hisEff_[e], "f");
0278         eff_[e]->SetTotalHistogram(*hisEffMC_[e], "f");
0279       }
0280     }
0281     log_ << "'Lost' below refers to truncation losses" << std::endl;
0282     // printout MC summary
0283     endJobMC();
0284     // printout DTC summary
0285     endJobDTC();
0286     log_ << "=============================================================";
0287     edm::LogPrint(moduleDescription().moduleName()) << log_.str();
0288   }
0289
0290   // fills kinematic tp histograms
0291   void Analyzer::fill(const TPPtr& tpPtr, const std::vector<TH1F*> th1fs) const {
0292     const double s = sin(tpPtr->phi());
0293     const double c = cos(tpPtr->phi());
0294     const TrackingParticle::Point& v = tpPtr->vertex();
0295     const std::vector<double> x = {tpPtr->phi(),
0296                                    tpPtr->pt(),
0297                                    tpPtr->charge() / tpPtr->pt(),
0298                                    v.x() * s - v.y() * c,
0299                                    v.z() - (v.x() * c + v.y() * s) * sinh(tpPtr->eta()),
0300                                    tpPtr->eta()};
0301     for (Efficiency e : AllEfficiency)
0302       th1fs[e]->Fill(x[e]);
0303   }
0304
0305   // fill stub related histograms
0306   void Analyzer::fill(const tt::StreamStub& stream, int region, int channel, int& sum, TH2F* th2f) {
0307     for (const tt::FrameStub& frame : stream) {
0308       if (frame.first.isNull())
0309         continue;
0310       sum++;
0311       const GlobalPoint& pos = setup_->stubPos(frame, region);
0312       const GlobalPoint& ttPos = setup_->stubPos(frame.first);
0313       const std::vector<double> resolutions = {
0314           ttPos.perp() - pos.perp(), tt::deltaPhi(ttPos.phi() - pos.phi()), ttPos.z() - pos.z()};
0315       for (Resolution r : AllResolution) {
0316         hisResolution_[r]->Fill(resolutions[r]);
0317         profResolution_[r]->Fill(ttPos.z(), ttPos.perp(), abs(resolutions[r]));
0318       }
0319       th2f->Fill(ttPos.z(), ttPos.perp());
0320     }
0321   }
0322
0323   // prints out Monte Carlo summary
0324   void Analyzer::endJobMC() {
0325     const double numStubs = profMC_->GetBinContent(1);
0326     const double numStubsMatched = profMC_->GetBinContent(2);
0327     const double numTPsReco = profMC_->GetBinContent(3);
0328     const double numTPsEff = profMC_->GetBinContent(4);
0329     const double errStubs = profMC_->GetBinError(1);
0330     const double errStubsMatched = profMC_->GetBinError(2);
0331     const double errTPsReco = profMC_->GetBinError(3);
0332     const double errTPsEff = profMC_->GetBinError(4);
0333     const double numCluster = profMC_->GetBinContent(6);
0334     const double errCluster = profMC_->GetBinError(6);
0335     const std::vector<double> nums = {numStubs, numStubsMatched, numTPsReco, numTPsEff, numCluster};
0336     const std::vector<double> errs = {errStubs, errStubsMatched, errTPsReco, errTPsEff, errCluster};
0337     const int wNums = std::ceil(std::log10(*std::max_element(nums.begin(), nums.end()))) + 5;
0338     const int wErrs = std::ceil(std::log10(*std::max_element(errs.begin(), errs.end()))) + 5;
0339     log_ << "=============================================================" << std::endl;
0340     log_ << "                         Monte Carlo  SUMMARY                         " << std::endl;
0341     log_ << "number of cluster       per TFP = " << std::setw(wNums) << numCluster << " +- " << std::setw(wErrs)
0342          << errCluster << std::endl;
0343     log_ << "number of stubs         per TFP = " << std::setw(wNums) << numStubs << " +- " << std::setw(wErrs)
0344          << errStubs << std::endl;
0345     log_ << "number of matched stubs per TFP = " << std::setw(wNums) << numStubsMatched << " +- " << std::setw(wErrs)
0346          << errStubsMatched << std::endl;
0347     log_ << "number of TPs           per TFP = " << std::setw(wNums) << numTPsReco << " +- " << std::setw(wErrs)
0348          << errTPsReco << std::endl;
0349     log_ << "number of TPs for eff   per TFP = " << std::setw(wNums) << numTPsEff << " +- " << std::setw(wErrs)
0350          << errTPsEff << std::endl;
0351   }
0352
0353   // prints out DTC summary
0354   void Analyzer::endJobDTC() {
0355     const double numStubs = profDTC_->GetBinContent(1);
0356     const double numStubsLost = profDTC_->GetBinContent(2);
0357     const double numTPs = profDTC_->GetBinContent(3);
0358     const double errStubs = profDTC_->GetBinError(1);
0359     const double errStubsLost = profDTC_->GetBinError(2);
0360     const double totalTPs = profMC_->GetBinContent(5);
0361     const double eff = numTPs / totalTPs;
0362     const double errEff = sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0363     const std::vector<double> nums = {numStubs, numStubsLost};
0364     const std::vector<double> errs = {errStubs, errStubsLost};
0365     const int wNums = std::ceil(std::log10(*std::max_element(nums.begin(), nums.end()))) + 5;
0366     const int wErrs = std::ceil(std::log10(*std::max_element(errs.begin(), errs.end()))) + 5;
0367     log_ << "=============================================================" << std::endl;
0368     log_ << "                         DTC SUMMARY                         " << std::endl;
0369     log_ << "number of stubs      per TFP = " << std::setw(wNums) << numStubs << " +- " << std::setw(wErrs) << errStubs
0370          << std::endl;
0371     log_ << "number of lost stubs per TFP = " << std::setw(wNums) << numStubsLost << " +- " << std::setw(wErrs)
0372          << errStubsLost << std::endl;
0373     log_ << "     max tracking efficiency = " << std::setw(wNums) << eff << " +- " << std::setw(wErrs) << errEff
0374          << std::endl;
0375   }
0376
0377 }  // namespace trackerDTC
0378
0379 DEFINE_FWK_MODULE(trackerDTC::Analyzer);