TrackFindingTracklet/test/AnalyzerKF.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/Common/interface/Handle.h"
0014
0015 #include "SimTracker/TrackTriggerAssociation/interface/StubAssociation.h"
0016 #include "L1Trigger/TrackTrigger/interface/Setup.h"
0017 #include "L1Trigger/TrackFindingTracklet/interface/DataFormats.h"
0018
0019 #include <TProfile.h>
0020 #include <TH1F.h>
0021 #include <TEfficiency.h>
0022
0023 #include <vector>
0024 #include <utility>
0025 #include <deque>
0026 #include <set>
0027 #include <cmath>
0028 #include <numeric>
0029 #include <sstream>
0030
0031 namespace trklet {
0032
0033   /*! \class  trklet::AnalyzerKF
0034    *  \brief  Class to analyze hardware like structured TTTrack Collection generated by Kalman Filter
0035    *  \author Thomas Schuh
0036    *  \date   2020, Sep
0037    */
0038   class AnalyzerKF : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0039   public:
0040     AnalyzerKF(const edm::ParameterSet& iConfig);
0041     void beginJob() override {}
0042     void beginRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) override;
0043     void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0044     void endRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) override {}
0045     void endJob() override;
0046
0047   private:
0048     //
0049     void associate(const std::vector<TrackKF>& tracks,
0050                    const std::vector<std::vector<StubKF*>>& stubs,
0051                    int region,
0052                    const tt::StubAssociation* ass,
0053                    std::set<TPPtr>& tps,
0054                    int& sum,
0055                    const std::vector<TH1F*>& his,
0056                    const std::vector<TProfile*>& prof,
0057                    bool perfect = true) const;
0058     //
0059     void analyzeTPz0(const tt::StubAssociation* sa);
0060     // ED input token of accepted Tracks
0061     edm::EDGetTokenT<tt::StreamsStub> edGetTokenStubs_;
0062     // ED input token of accepted Stubs
0063     edm::EDGetTokenT<tt::StreamsTrack> edGetTokenTracks_;
0064     // ED input token for number of accepted States
0065     edm::EDGetTokenT<int> edGetTokenNumStatesAccepted_;
0066     // ED input token for number of lost States
0067     edm::EDGetTokenT<int> edGetTokenNumStatesTruncated_;
0068     // ED input token of TTStubRef to TPPtr association for tracking efficiency
0069     edm::EDGetTokenT<tt::StubAssociation> edGetTokenSelection_;
0070     // ED input token of TTStubRef to recontructable TPPtr association
0071     edm::EDGetTokenT<tt::StubAssociation> edGetTokenReconstructable_;
0072     // Setup token
0073     edm::ESGetToken<tt::Setup, tt::SetupRcd> esGetTokenSetup_;
0074     // DataFormats token
0075     edm::ESGetToken<DataFormats, ChannelAssignmentRcd> esGetTokenDataFormats_;
0076     // stores, calculates and provides run-time constants
0077     const tt::Setup* setup_ = nullptr;
0078     //
0079     const DataFormats* dataFormats_ = nullptr;
0080     // enables analyze of TPs
0081     bool useMCTruth_;
0082     //
0083     int nEvents_ = 0;
0084     //
0085     std::string label_;
0086
0087     // Histograms
0088
0089     TProfile* prof_;
0090     TProfile* profChannel_;
0091     TH1F* hisChannel_;
0092     std::vector<TH1F*> hisRes_;
0093     std::vector<TProfile*> profRes_;
0094     TH1F* hisEffD0_;
0095     TH1F* hisEffD0Total_;
0096     TEfficiency* effD0_;
0097     TH1F* hisEffEta_;
0098     TH1F* hisEffEtaTotal_;
0099     TEfficiency* effEta_;
0100     TH1F* hisEffInv2R_;
0101     TH1F* hisEffInv2RTotal_;
0102     TEfficiency* effInv2R_;
0103     TH1F* hisEffPT_;
0104     TH1F* hisEffPTTotal_;
0105     TEfficiency* effPT_;
0106     TH1F* hisTracks_;
0107     TH1F* hisLayers_;
0108     TH1F* hisNumLayers_;
0109     TProfile* profNumLayers_;
0110
0111     // printout
0112     std::stringstream log_;
0113   };
0114
0115   AnalyzerKF::AnalyzerKF(const edm::ParameterSet& iConfig)
0116       : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")), nEvents_(0), hisRes_(9), profRes_(5) {
0117     usesResource("TFileService");
0118     // book in- and output ED products
0119     label_ = iConfig.getParameter<std::string>("OutputLabelKF");
0120     const std::string& branchStubs = iConfig.getParameter<std::string>("BranchStubs");
0121     const std::string& branchTracks = iConfig.getParameter<std::string>("BranchTracks");
0122     const std::string& branchTruncated = iConfig.getParameter<std::string>("BranchTruncated");
0123     edGetTokenStubs_ = consumes<tt::StreamsStub>(edm::InputTag(label_, branchStubs));
0124     edGetTokenTracks_ = consumes<tt::StreamsTrack>(edm::InputTag(label_, branchTracks));
0125     edGetTokenNumStatesAccepted_ = consumes<int>(edm::InputTag(label_, branchTracks));
0126     edGetTokenNumStatesTruncated_ = consumes<int>(edm::InputTag(label_, branchTruncated));
0127     if (useMCTruth_) {
0128       const auto& inputTagSelecttion = iConfig.getParameter<edm::InputTag>("InputTagSelection");
0129       const auto& inputTagReconstructable = iConfig.getParameter<edm::InputTag>("InputTagReconstructable");
0130       edGetTokenSelection_ = consumes<tt::StubAssociation>(inputTagSelecttion);
0131       edGetTokenReconstructable_ = consumes<tt::StubAssociation>(inputTagReconstructable);
0132     }
0133     // book ES products
0134     esGetTokenSetup_ = esConsumes<edm::Transition::BeginRun>();
0135     esGetTokenDataFormats_ = esConsumes<edm::Transition::BeginRun>();
0136     // log config
0137     log_.setf(std::ios::fixed, std::ios::floatfield);
0138     log_.precision(4);
0139   }
0140
0141   void AnalyzerKF::beginRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) {
0142     // helper class to store configurations
0143     setup_ = &iSetup.getData(esGetTokenSetup_);
0144     dataFormats_ = &iSetup.getData(esGetTokenDataFormats_);
0145     // book histograms
0146     edm::Service<TFileService> fs;
0147     TFileDirectory dir;
0148     dir = fs->mkdir("KF");
0149     prof_ = dir.make<TProfile>("Counts", ";", 14, 0.5, 14.5);
0150     prof_->GetXaxis()->SetBinLabel(1, "Stubs");
0151     prof_->GetXaxis()->SetBinLabel(2, "Tracks");
0152     prof_->GetXaxis()->SetBinLabel(4, "Matched Tracks");
0153     prof_->GetXaxis()->SetBinLabel(5, "All Tracks");
0154     prof_->GetXaxis()->SetBinLabel(6, "Found TPs");
0155     prof_->GetXaxis()->SetBinLabel(7, "Found selected TPs");
0156     prof_->GetXaxis()->SetBinLabel(9, "All TPs");
0157     prof_->GetXaxis()->SetBinLabel(10, "states");
0158     prof_->GetXaxis()->SetBinLabel(12, "max tp");
0159     prof_->GetXaxis()->SetBinLabel(13, "All electron TPs");
0160     prof_->GetXaxis()->SetBinLabel(14, "max electron tp");
0161     // channel occupancy
0162     constexpr int maxOcc = 180;
0163     const int numChannels = 1;
0164     hisChannel_ = dir.make<TH1F>("His Channel Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0165     profChannel_ = dir.make<TProfile>("Prof Channel Occupancy", ";", numChannels, -.5, numChannels - .5);
0166     // resoultions
0167     static const std::vector<std::string> names = {"phi0", "inv2R", "z0", "cot", "d0"};
0168     static const std::vector<double> ranges = {.01, .004, 5., .4, 5.};
0169     for (int i = 0; i < 5; i++) {
0170       hisRes_[i] = dir.make<TH1F>(("HisRes" + names[i]).c_str(), ";", 100, -ranges[i], ranges[i]);
0171       profRes_[i] = dir.make<TProfile>(("ProfRes" + names[i]).c_str(), ";", 32, 0, 2.4);
0172     }
0173     for (int i = 5; i < 9; i++) {
0174       const std::string name = (i < 7 ? names[2] : names[3]) + (i % 2 == 1 ? "plus" : "minus");
0175       const double range = (i < 7 ? ranges[2] : ranges[3]);
0176       hisRes_[i] = dir.make<TH1F>(("HisRes" + name).c_str(), ";", 100, -range, range);
0177     }
0178     // Efficiencies
0179     const double rangeD0 = 10.;
0180     hisEffD0_ = dir.make<TH1F>("HisTPD0", ";", 32, -rangeD0 / 2., rangeD0 / 2.);
0181     hisEffD0Total_ = dir.make<TH1F>("HisTPD0Total", ";", 32, -rangeD0 / 2., rangeD0 / 2.);
0182     effD0_ = dir.make<TEfficiency>("EffD0", ";", 32, -rangeD0 / 2., rangeD0 / 2.);
0183     hisEffEtaTotal_ = dir.make<TH1F>("HisTPEtaTotal", ";", 48, -2.4, 2.4);
0184     hisEffEta_ = dir.make<TH1F>("HisTPEta", ";", 48, -2.4, 2.4);
0185     effEta_ = dir.make<TEfficiency>("EffEta", ";", 48, -2.4, 2.4);
0186     const double rangeInv2R = dataFormats_->format(Variable::inv2R, Process::dr).range();
0187     hisEffInv2R_ = dir.make<TH1F>("HisTPInv2R", ";", 32, -rangeInv2R / 2., rangeInv2R / 2.);
0188     hisEffInv2RTotal_ = dir.make<TH1F>("HisTPInv2RTotal", ";", 32, -rangeInv2R / 2., rangeInv2R / 2.);
0189     effInv2R_ = dir.make<TEfficiency>("EffInv2R", ";", 32, -rangeInv2R / 2., rangeInv2R / 2.);
0190     hisEffPT_ = dir.make<TH1F>("HisTPPT", ";", 100, 0, 100);
0191     hisEffPTTotal_ = dir.make<TH1F>("HisTPPTTotal", ";", 100, 0, 100);
0192     effPT_ = dir.make<TEfficiency>("EffPT", ";", 100, 0, 100);
0193     // tracks
0194     hisTracks_ = dir.make<TH1F>("HisTracks", ";", 40, 0., 400);
0195     // layers
0196     hisLayers_ = dir.make<TH1F>("HisLayers", ";", 8, 0, 8);
0197     hisNumLayers_ = dir.make<TH1F>("HisNumLayers", ";", 9, 0, 9);
0198     profNumLayers_ = dir.make<TProfile>("Prof NumLayers", ";", 32, 0, 2.4);
0199   }
0200
0201   void AnalyzerKF::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0202     static const int numLayers = setup_->numLayers();
0203     auto fill = [this](const TPPtr& tpPtr, TH1F* hisEta, TH1F* hisInv2R, TH1F* hisPT, TH1F* hisD0) {
0204       hisEta->Fill(tpPtr->eta());
0205       hisInv2R->Fill(tpPtr->charge() / tpPtr->pt() * setup_->invPtToDphi());
0206       hisPT->Fill(tpPtr->pt());
0207       hisD0->Fill(tpPtr->d0());
0208     };
0209     // read in kf products
0210     edm::Handle<tt::StreamsStub> handleStubs;
0211     iEvent.getByToken<tt::StreamsStub>(edGetTokenStubs_, handleStubs);
0212     const tt::StreamsStub& allStubs = *handleStubs;
0213     edm::Handle<tt::StreamsTrack> handleTracks;
0214     iEvent.getByToken<tt::StreamsTrack>(edGetTokenTracks_, handleTracks);
0215     const tt::StreamsTrack& allTracks = *handleTracks;
0216     edm::Handle<int> handleNumStatesAccepted;
0217     iEvent.getByToken<int>(edGetTokenNumStatesAccepted_, handleNumStatesAccepted);
0218     edm::Handle<int> handleNumStatesTruncated;
0219     iEvent.getByToken<int>(edGetTokenNumStatesTruncated_, handleNumStatesTruncated);
0220     // read in MCTruth
0221     const tt::StubAssociation* selection = nullptr;
0222     const tt::StubAssociation* reconstructable = nullptr;
0223     if (useMCTruth_) {
0224       edm::Handle<tt::StubAssociation> handleSelection;
0225       iEvent.getByToken<tt::StubAssociation>(edGetTokenSelection_, handleSelection);
0226       selection = handleSelection.product();
0227       prof_->Fill(9, selection->numTPs());
0228       edm::Handle<tt::StubAssociation> handleReconstructable;
0229       iEvent.getByToken<tt::StubAssociation>(edGetTokenReconstructable_, handleReconstructable);
0230       reconstructable = handleReconstructable.product();
0231       for (const auto& p : selection->getTrackingParticleToTTStubsMap())
0232         fill(p.first, hisEffEtaTotal_, hisEffInv2RTotal_, hisEffPTTotal_, hisEffD0Total_);
0233     }
0234     // analyze kf products and associate found tracks with reconstrucable TrackingParticles
0235     std::set<TPPtr> tpPtrs;
0236     std::set<TPPtr> tpPtrsSelection;
0237     std::set<TPPtr> tpPtrsMax;
0238     int numMatched(0);
0239     int numTracks(0);
0240     for (int region = 0; region < setup_->numRegions(); region++) {
0241       int nRegionStubs(0);
0242       int nRegionTracks(0);
0243       const int offset = region * numLayers;
0244       const tt::StreamTrack& channelTracks = allTracks[region];
0245       hisChannel_->Fill(channelTracks.size());
0246       profChannel_->Fill(1, channelTracks.size());
0247       std::vector<TrackKF> tracks;
0248       std::vector<StubKF> stubs;
0249       std::vector<std::vector<StubKF*>> tracksStubs(channelTracks.size(), std::vector<StubKF*>(numLayers, nullptr));
0250       tracks.reserve(channelTracks.size());
0251       stubs.reserve(channelTracks.size() * numLayers);
0252       for (int frame = 0; frame < (int)channelTracks.size(); frame++) {
0253         tracks.emplace_back(channelTracks[frame], dataFormats_);
0254         const double eta = std::abs(std::asinh(tracks.back().cot()));
0255         int nLayers(0);
0256         for (int layer = 0; layer < numLayers; layer++) {
0257           const tt::FrameStub& fs = allStubs[offset + layer][frame];
0258           if (fs.first.isNull())
0259             continue;
0260           stubs.emplace_back(fs, dataFormats_);
0261           tracksStubs[frame][layer] = &stubs.back();
0262           hisLayers_->Fill(layer);
0263           nLayers++;
0264         }
0265         hisNumLayers_->Fill(nLayers);
0266         profNumLayers_->Fill(eta, nLayers);
0267       }
0268       nRegionStubs += stubs.size();
0269       nRegionTracks += tracks.size();
0270       if (!useMCTruth_)
0271         continue;
0272       int tmp(0);
0273       associate(tracks, tracksStubs, region, selection, tpPtrsSelection, tmp, hisRes_, profRes_);
0274       associate(tracks,
0275                 tracksStubs,
0276                 region,
0277                 reconstructable,
0278                 tpPtrs,
0279                 numMatched,
0280                 std::vector<TH1F*>(),
0281                 std::vector<TProfile*>());
0282       associate(
0283           tracks, tracksStubs, region, selection, tpPtrsMax, tmp, std::vector<TH1F*>(), std::vector<TProfile*>(), false);
0284       numTracks += nRegionTracks;
0285       prof_->Fill(1, nRegionStubs);
0286       prof_->Fill(2, nRegionTracks);
0287     }
0288     for (const TPPtr& tpPtr : tpPtrsSelection)
0289       fill(tpPtr, hisEffEta_, hisEffInv2R_, hisEffPT_, hisEffD0_);
0290     prof_->Fill(4, numMatched);
0291     prof_->Fill(5, numTracks);
0292     prof_->Fill(6, tpPtrs.size());
0293     prof_->Fill(7, tpPtrsSelection.size());
0294     prof_->Fill(10, *handleNumStatesAccepted);
0295     prof_->Fill(11, *handleNumStatesTruncated);
0296     prof_->Fill(12, tpPtrsMax.size());
0297     hisTracks_->Fill(numTracks);
0298     nEvents_++;
0299   }
0300
0301   void AnalyzerKF::endJob() {
0302     if (nEvents_ == 0)
0303       return;
0304     // effi
0305     effD0_->SetPassedHistogram(*hisEffD0_, "f");
0306     effD0_->SetTotalHistogram(*hisEffD0Total_, "f");
0307     effEta_->SetPassedHistogram(*hisEffEta_, "f");
0308     effEta_->SetTotalHistogram(*hisEffEtaTotal_, "f");
0309     effInv2R_->SetPassedHistogram(*hisEffInv2R_, "f");
0310     effInv2R_->SetTotalHistogram(*hisEffInv2RTotal_, "f");
0311     effPT_->SetPassedHistogram(*hisEffPT_, "f");
0312     effPT_->SetTotalHistogram(*hisEffPTTotal_, "f");
0313     // printout SF summary
0314     const double totalTPs = prof_->GetBinContent(9);
0315     const double numStubs = prof_->GetBinContent(1);
0316     const double numTracks = prof_->GetBinContent(2);
0317     const double totalTracks = prof_->GetBinContent(5);
0318     const double numTracksMatched = prof_->GetBinContent(4);
0319     const double numTPsAll = prof_->GetBinContent(6);
0320     const double numTPsEff = prof_->GetBinContent(7);
0321     const double numTPsEffMax = prof_->GetBinContent(12);
0322     const double errStubs = prof_->GetBinError(1);
0323     const double errTracks = prof_->GetBinError(2);
0324     const double fracFake = (totalTracks - numTracksMatched) / totalTracks;
0325     const double fracDup = (numTracksMatched - numTPsAll) / totalTracks;
0326     const double eff = numTPsEff / totalTPs;
0327     const double errEff = std::sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0328     const double effMax = numTPsEffMax / totalTPs;
0329     const double errEffMax = std::sqrt(effMax * (1. - effMax) / totalTPs / nEvents_);
0330     const int numStates = prof_->GetBinContent(10);
0331     const int numStatesLost = prof_->GetBinContent(11);
0332     const double fracSatest = numStates / (double)(numStates + numStatesLost);
0333     const std::vector<double> nums = {numStubs, numTracks};
0334     const std::vector<double> errs = {errStubs, errTracks};
0335     const int wNums = std::ceil(std::log10(*std::max_element(nums.begin(), nums.end()))) + 5;
0336     const int wErrs = std::ceil(std::log10(*std::max_element(errs.begin(), errs.end()))) + 5;
0337     log_ << "                         KF  SUMMARY                         " << std::endl;
0338     log_ << "number of stubs       per TFP = " << std::setw(wNums) << numStubs << " +- " << std::setw(wErrs) << errStubs
0339          << std::endl;
0340     log_ << "number of tracks      per TFP = " << std::setw(wNums) << numTracks << " +- " << std::setw(wErrs)
0341          << errTracks << std::endl;
0342     log_ << "          tracking efficiency = " << std::setw(wNums) << eff << " +- " << std::setw(wErrs) << errEff
0343          << std::endl;
0344     log_ << "      max tracking efficiency = " << std::setw(wNums) << effMax << " +- " << std::setw(wErrs) << errEffMax
0345          << std::endl;
0346     log_ << "                    fake rate = " << std::setw(wNums) << fracFake << std::endl;
0347     log_ << "               duplicate rate = " << std::setw(wNums) << fracDup << std::endl;
0348     log_ << "    state assessment fraction = " << std::setw(wNums) << fracSatest << std::endl;
0349     log_ << "     number of states per TFP = " << std::setw(wNums) << (numStates + numStatesLost) / setup_->numRegions()
0350          << std::endl;
0351     log_ << "=============================================================";
0352     edm::LogPrint(moduleDescription().moduleName()) << log_.str();
0353   }
0354
0355   //
0356   void AnalyzerKF::associate(const std::vector<TrackKF>& tracks,
0357                              const std::vector<std::vector<StubKF*>>& tracksStubs,
0358                              int region,
0359                              const tt::StubAssociation* ass,
0360                              std::set<TPPtr>& tps,
0361                              int& sum,
0362                              const std::vector<TH1F*>& his,
0363                              const std::vector<TProfile*>& prof,
0364                              bool perfect) const {
0365     for (int frame = 0; frame < static_cast<int>(tracks.size()); frame++) {
0366       const TrackKF& track = tracks[frame];
0367       const std::vector<StubKF*>& stubs = tracksStubs[frame];
0368       std::vector<TTStubRef> ttStubRefs;
0369       ttStubRefs.reserve(stubs.size());
0370       TTBV hitPattern(0, setup_->numLayers());
0371       int layer(-1);
0372       for (StubKF* stub : stubs) {
0373         layer++;
0374         if (!stub)
0375           continue;
0376         hitPattern.set(layer);
0377         ttStubRefs.push_back(stub->frame().first);
0378       }
0379       const std::vector<TPPtr>& tpPtrs = perfect ? ass->associateFinal(ttStubRefs) : ass->associate(ttStubRefs);
0380       if (tpPtrs.empty())
0381         continue;
0382       sum++;
0383       std::copy(tpPtrs.begin(), tpPtrs.end(), std::inserter(tps, tps.begin()));
0384       if (his.empty())
0385         continue;
0386       const double cot = track.cot();
0387       const double z0 = track.zT() - setup_->chosenRofZ() * cot;
0388       const double inv2R = track.inv2R();
0389       const double phi0 = tt::deltaPhi(track.phiT() - setup_->chosenRofPhi() * inv2R +
0390                                        region * dataFormats_->format(Variable::phiT, Process::kf).range());
0391       for (const TPPtr& tpPtr : tpPtrs) {
0392         const double tpPhi0 = tpPtr->phi();
0393         const double tpCot = std::sinh(tpPtr->eta());
0394         const double tpInv2R = -setup_->invPtToDphi() * tpPtr->charge() / tpPtr->pt();
0395         const math::XYZPointD& v = tpPtr->vertex();
0396         const double tpZ0 = v.z() - tpCot * (v.x() * cos(tpPhi0) + v.y() * sin(tpPhi0));
0397         const double dCot = tpCot - cot;
0398         const double dZ0 = tpZ0 - z0;
0399         const double dInv2R = tpInv2R - inv2R;
0400         const double dPhi0 = tt::deltaPhi(tpPhi0 - phi0);
0401         const double dD0 = tpPtr->d0() + track.frame().first->d0();
0402         const std::vector<double> ds = {dPhi0, dInv2R / setup_->invPtToDphi(), dZ0, dCot, dD0};
0403         for (int i = 0; i < (int)ds.size(); i++) {
0404           his[i]->Fill(ds[i]);
0405           prof[i]->Fill(std::abs(tpPtr->eta()), std::abs(ds[i]));
0406         }
0407         if (tpCot < 0) {
0408           his[6]->Fill(dZ0);
0409           his[8]->Fill(dCot);
0410         } else {
0411           his[5]->Fill(dZ0);
0412           his[7]->Fill(dCot);
0413         }
0414       }
0415     }
0416   }
0417
0418 }  // namespace trklet
0419
0420 DEFINE_FWK_MODULE(trklet::AnalyzerKF);