Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:56:13 

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 #include <TEfficiency.h>
0024 
0025 #include <vector>
0026 #include <deque>
0027 #include <set>
0028 #include <cmath>
0029 #include <numeric>
0030 #include <sstream>
0031 
0032 using namespace std;
0033 using namespace edm;
0034 using namespace tt;
0035 
0036 namespace trackerTFP {
0037 
0038   /*! \class  trackerTFP::AnalyzerKF
0039    *  \brief  Class to analyze hardware like structured TTTrack Collection generated by Kalman Filter
0040    *  \author Thomas Schuh
0041    *  \date   2020, Sep
0042    */
0043   class AnalyzerKF : public one::EDAnalyzer<one::WatchRuns, one::SharedResources> {
0044   public:
0045     AnalyzerKF(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 TTTracks& ttTracks,
0055                    const StubAssociation* ass,
0056                    set<TPPtr>& tps,
0057                    int& sum,
0058                    const vector<TH1F*>& his,
0059                    TProfile* prof) const;
0060 
0061     // ED input token of accepted Tracks
0062     EDGetTokenT<StreamsStub> edGetTokenAcceptedStubs_;
0063     // ED input token of accepted Stubs
0064     EDGetTokenT<StreamsTrack> edGetTokenAcceptedTracks_;
0065     // ED input token of lost Stubs
0066     EDGetTokenT<StreamsStub> edGetTokenLostStubs_;
0067     // ED input token of lost Tracks
0068     EDGetTokenT<StreamsTrack> edGetTokenLostTracks_;
0069     // ED input token for number of accepted States
0070     EDGetTokenT<int> edGetTokenNumAcceptedStates_;
0071     // ED input token for number of lost States
0072     EDGetTokenT<int> edGetTokenNumLostStates_;
0073     // ED input token of TTStubRef to TPPtr association for tracking efficiency
0074     EDGetTokenT<StubAssociation> edGetTokenSelection_;
0075     // ED input token of TTStubRef to recontructable TPPtr association
0076     EDGetTokenT<StubAssociation> edGetTokenReconstructable_;
0077     // Setup token
0078     ESGetToken<Setup, SetupRcd> esGetTokenSetup_;
0079     // DataFormats token
0080     ESGetToken<DataFormats, DataFormatsRcd> esGetTokenDataFormats_;
0081     // LayerEncoding token
0082     ESGetToken<LayerEncoding, LayerEncodingRcd> esGetTokenLayerEncoding_;
0083     // stores, calculates and provides run-time constants
0084     const Setup* setup_ = nullptr;
0085     //
0086     const DataFormats* dataFormats_ = nullptr;
0087     //
0088     const LayerEncoding* layerEncoding_ = nullptr;
0089     // enables analyze of TPs
0090     bool useMCTruth_;
0091     //
0092     int nEvents_ = 0;
0093 
0094     // Histograms
0095 
0096     TProfile* prof_;
0097     TProfile* profChannel_;
0098     TH1F* hisChannel_;
0099     vector<TH1F*> hisRes_;
0100     TProfile* profResZ0_;
0101     TH1F* hisEffEta_;
0102     TH1F* hisEffEtaTotal_;
0103     TEfficiency* effEta_;
0104     TH1F* hisEffInv2R_;
0105     TH1F* hisEffInv2RTotal_;
0106     TEfficiency* effInv2R_;
0107     TH1F* hisChi2_;
0108     TH1F* hisPhi_;
0109 
0110     // printout
0111     stringstream log_;
0112   };
0113 
0114   AnalyzerKF::AnalyzerKF(const ParameterSet& iConfig)
0115       : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")), hisRes_(4) {
0116     usesResource("TFileService");
0117     // book in- and output ED products
0118     const string& label = iConfig.getParameter<string>("LabelKF");
0119     const string& branchAcceptedStubs = iConfig.getParameter<string>("BranchAcceptedStubs");
0120     const string& branchAcceptedTracks = iConfig.getParameter<string>("BranchAcceptedTracks");
0121     const string& branchLostStubs = iConfig.getParameter<string>("BranchLostStubs");
0122     const string& branchLostTracks = iConfig.getParameter<string>("BranchLostTracks");
0123     edGetTokenAcceptedStubs_ = consumes<StreamsStub>(InputTag(label, branchAcceptedStubs));
0124     edGetTokenAcceptedTracks_ = consumes<StreamsTrack>(InputTag(label, branchAcceptedTracks));
0125     edGetTokenLostStubs_ = consumes<StreamsStub>(InputTag(label, branchLostStubs));
0126     edGetTokenLostTracks_ = consumes<StreamsTrack>(InputTag(label, branchLostTracks));
0127     edGetTokenNumAcceptedStates_ = consumes<int>(InputTag(label, branchAcceptedTracks));
0128     ;
0129     edGetTokenNumLostStates_ = consumes<int>(InputTag(label, branchLostTracks));
0130     ;
0131     if (useMCTruth_) {
0132       const auto& inputTagSelecttion = iConfig.getParameter<InputTag>("InputTagSelection");
0133       const auto& inputTagReconstructable = iConfig.getParameter<InputTag>("InputTagReconstructable");
0134       edGetTokenSelection_ = consumes<StubAssociation>(inputTagSelecttion);
0135       edGetTokenReconstructable_ = consumes<StubAssociation>(inputTagReconstructable);
0136     }
0137     // book ES products
0138     esGetTokenSetup_ = esConsumes<Setup, SetupRcd, Transition::BeginRun>();
0139     esGetTokenDataFormats_ = esConsumes<DataFormats, DataFormatsRcd, Transition::BeginRun>();
0140     esGetTokenLayerEncoding_ = esConsumes<LayerEncoding, LayerEncodingRcd, Transition::BeginRun>();
0141     // log config
0142     log_.setf(ios::fixed, ios::floatfield);
0143     log_.precision(4);
0144   }
0145 
0146   void AnalyzerKF::beginRun(const Run& iEvent, const EventSetup& iSetup) {
0147     // helper class to store configurations
0148     setup_ = &iSetup.getData(esGetTokenSetup_);
0149     dataFormats_ = &iSetup.getData(esGetTokenDataFormats_);
0150     layerEncoding_ = &iSetup.getData(esGetTokenLayerEncoding_);
0151     // book histograms
0152     Service<TFileService> fs;
0153     TFileDirectory dir;
0154     dir = fs->mkdir("KF");
0155     prof_ = dir.make<TProfile>("Counts", ";", 11, 0.5, 11.5);
0156     prof_->GetXaxis()->SetBinLabel(1, "Stubs");
0157     prof_->GetXaxis()->SetBinLabel(2, "Tracks");
0158     prof_->GetXaxis()->SetBinLabel(3, "Lost Tracks");
0159     prof_->GetXaxis()->SetBinLabel(4, "Matched Tracks");
0160     prof_->GetXaxis()->SetBinLabel(5, "All Tracks");
0161     prof_->GetXaxis()->SetBinLabel(6, "Found TPs");
0162     prof_->GetXaxis()->SetBinLabel(7, "Found selected TPs");
0163     prof_->GetXaxis()->SetBinLabel(8, "Lost TPs");
0164     prof_->GetXaxis()->SetBinLabel(9, "All TPs");
0165     prof_->GetXaxis()->SetBinLabel(10, "states");
0166     prof_->GetXaxis()->SetBinLabel(11, "lost states");
0167     // channel occupancy
0168     constexpr int maxOcc = 180;
0169     const int numChannels = dataFormats_->numChannel(Process::kf);
0170     hisChannel_ = dir.make<TH1F>("His Channel Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0171     profChannel_ = dir.make<TProfile>("Prof Channel Occupancy", ";", numChannels, -.5, numChannels - .5);
0172     // resoultions
0173     static const vector<string> names = {"phiT", "inv2R", "zT", "cot"};
0174     static const vector<double> ranges = {.01, .1, 5, .1};
0175     for (int i = 0; i < 4; i++) {
0176       const double range = ranges[i];
0177       hisRes_[i] = dir.make<TH1F>(("HisRes" + names[i]).c_str(), ";", 100, -range, range);
0178     }
0179     profResZ0_ = dir.make<TProfile>("ProfResZ0", ";", 32, 0, 2.5);
0180     // Efficiencies
0181     hisEffEtaTotal_ = dir.make<TH1F>("HisTPEtaTotal", ";", 128, -2.5, 2.5);
0182     hisEffEta_ = dir.make<TH1F>("HisTPEta", ";", 128, -2.5, 2.5);
0183     effEta_ = dir.make<TEfficiency>("EffEta", ";", 128, -2.5, 2.5);
0184     const double rangeInv2R = dataFormats_->format(Variable::inv2R, Process::dr).range();
0185     hisEffInv2R_ = dir.make<TH1F>("HisTPInv2R", ";", 32, -rangeInv2R / 2., rangeInv2R / 2.);
0186     hisEffInv2RTotal_ = dir.make<TH1F>("HisTPInv2RTotal", ";", 32, -rangeInv2R / 2., rangeInv2R / 2.);
0187     effInv2R_ = dir.make<TEfficiency>("EffInv2R", ";", 32, -rangeInv2R / 2., rangeInv2R / 2.);
0188     // chi2
0189     hisChi2_ = dir.make<TH1F>("HisChi2", ";", 100, -.5, 99.5);
0190     const double rangePhi = dataFormats_->format(Variable::phi0, Process::dr).range();
0191     hisPhi_ = dir.make<TH1F>("HisPhi", ";", 100, -rangePhi, rangePhi);
0192   }
0193 
0194   void AnalyzerKF::analyze(const Event& iEvent, const EventSetup& iSetup) {
0195     auto fill = [this](const TPPtr& tpPtr, TH1F* hisEta, TH1F* hisInv2R) {
0196       hisEta->Fill(tpPtr->eta());
0197       hisInv2R->Fill(tpPtr->charge() / tpPtr->pt() * setup_->invPtToDphi());
0198     };
0199     // read in kf products
0200     Handle<StreamsStub> handleAcceptedStubs;
0201     iEvent.getByToken<StreamsStub>(edGetTokenAcceptedStubs_, handleAcceptedStubs);
0202     const StreamsStub& acceptedStubs = *handleAcceptedStubs;
0203     Handle<StreamsTrack> handleAcceptedTracks;
0204     iEvent.getByToken<StreamsTrack>(edGetTokenAcceptedTracks_, handleAcceptedTracks);
0205     Handle<StreamsStub> handleLostStubs;
0206     iEvent.getByToken<StreamsStub>(edGetTokenLostStubs_, handleLostStubs);
0207     const StreamsStub& lostStubs = *handleLostStubs;
0208     Handle<StreamsTrack> handleLostTracks;
0209     iEvent.getByToken<StreamsTrack>(edGetTokenLostTracks_, handleLostTracks);
0210     Handle<int> handleNumAcceptedStates;
0211     iEvent.getByToken<int>(edGetTokenNumAcceptedStates_, handleNumAcceptedStates);
0212     Handle<int> handleNumLostStates;
0213     iEvent.getByToken<int>(edGetTokenNumLostStates_, handleNumLostStates);
0214     // read in MCTruth
0215     const StubAssociation* selection = nullptr;
0216     const StubAssociation* reconstructable = nullptr;
0217     if (useMCTruth_) {
0218       Handle<StubAssociation> handleSelection;
0219       iEvent.getByToken<StubAssociation>(edGetTokenSelection_, handleSelection);
0220       selection = handleSelection.product();
0221       prof_->Fill(9, selection->numTPs());
0222       Handle<StubAssociation> handleReconstructable;
0223       iEvent.getByToken<StubAssociation>(edGetTokenReconstructable_, handleReconstructable);
0224       reconstructable = handleReconstructable.product();
0225       for (const auto& p : selection->getTrackingParticleToTTStubsMap())
0226         fill(p.first, hisEffEtaTotal_, hisEffInv2RTotal_);
0227     }
0228     // analyze kf products and associate found tracks with reconstrucable TrackingParticles
0229     set<TPPtr> tpPtrs;
0230     set<TPPtr> tpPtrsSelection;
0231     set<TPPtr> tpPtrsLost;
0232     int allMatched(0);
0233     int allTracks(0);
0234     auto consume = [this](const StreamTrack& tracks, const StreamsStub& streams, int channel, TTTracks& ttTracks) {
0235       const int offset = channel * setup_->numLayers();
0236       int pos(0);
0237       for (const FrameTrack& frameTrack : tracks) {
0238         vector<StubKF> stubs;
0239         stubs.reserve(setup_->numLayers());
0240         for (int layer = 0; layer < setup_->numLayers(); layer++) {
0241           const FrameStub& frameStub = streams[offset + layer][pos];
0242           if (frameStub.first.isNonnull())
0243             stubs.emplace_back(frameStub, dataFormats_, layer);
0244         }
0245         TrackKF track(frameTrack, dataFormats_);
0246         ttTracks.emplace_back(track.ttTrack(stubs));
0247         pos++;
0248       }
0249     };
0250     for (int region = 0; region < setup_->numRegions(); region++) {
0251       int nStubsRegion(0);
0252       int nTracksRegion(0);
0253       int nLostRegion(0);
0254       for (int channel = 0; channel < dataFormats_->numChannel(Process::kf); channel++) {
0255         const int index = region * dataFormats_->numChannel(Process::kf) + channel;
0256         const StreamTrack& accepted = handleAcceptedTracks->at(index);
0257         const StreamTrack& lost = handleLostTracks->at(index);
0258         hisChannel_->Fill(accepted.size());
0259         profChannel_->Fill(channel, accepted.size());
0260         TTTracks tracks;
0261         const int nTracks = accumulate(accepted.begin(), accepted.end(), 0, [](int sum, const FrameTrack& frame) {
0262           return sum + (frame.first.isNonnull() ? 1 : 0);
0263         });
0264         nTracksRegion += nTracks;
0265         tracks.reserve(nTracks);
0266         consume(accepted, acceptedStubs, index, tracks);
0267         for (const TTTrack<Ref_Phase2TrackerDigi_>& ttTrack : tracks)
0268           hisPhi_->Fill(ttTrack.momentum().phi());
0269         nStubsRegion += accumulate(tracks.begin(), tracks.end(), 0, [](int sum, const auto& ttTrack) {
0270           return sum + (int)ttTrack.getStubRefs().size();
0271         });
0272         TTTracks tracksLost;
0273         const int nLost = accumulate(lost.begin(), lost.end(), 0, [](int sum, const FrameTrack& frame) {
0274           return sum + (frame.first.isNonnull() ? 1 : 0);
0275         });
0276         nLostRegion += nLost;
0277         tracksLost.reserve(nLost);
0278         consume(lost, lostStubs, index, tracksLost);
0279         allTracks += nTracks;
0280         if (!useMCTruth_)
0281           continue;
0282         int tmp(0);
0283         associate(tracks, selection, tpPtrsSelection, tmp, hisRes_, profResZ0_);
0284         associate(tracksLost, selection, tpPtrsLost, tmp, vector<TH1F*>(), nullptr);
0285         associate(tracks, reconstructable, tpPtrs, allMatched, vector<TH1F*>(), nullptr);
0286       }
0287       prof_->Fill(1, nStubsRegion);
0288       prof_->Fill(2, nTracksRegion);
0289       prof_->Fill(3, nLostRegion);
0290     }
0291     for (const TPPtr& tpPtr : tpPtrsSelection)
0292       fill(tpPtr, hisEffEta_, hisEffInv2R_);
0293     deque<TPPtr> tpPtrsRealLost;
0294     set_difference(tpPtrsLost.begin(), tpPtrsLost.end(), tpPtrs.begin(), tpPtrs.end(), back_inserter(tpPtrsRealLost));
0295     prof_->Fill(4, allMatched);
0296     prof_->Fill(5, allTracks);
0297     prof_->Fill(6, tpPtrs.size());
0298     prof_->Fill(7, tpPtrsSelection.size());
0299     prof_->Fill(8, tpPtrsRealLost.size());
0300     prof_->Fill(10, *handleNumAcceptedStates);
0301     prof_->Fill(11, *handleNumLostStates);
0302     nEvents_++;
0303   }
0304 
0305   void AnalyzerKF::endJob() {
0306     if (nEvents_ == 0)
0307       return;
0308     // effi
0309     effEta_->SetPassedHistogram(*hisEffEta_, "f");
0310     effEta_->SetTotalHistogram(*hisEffEtaTotal_, "f");
0311     effInv2R_->SetPassedHistogram(*hisEffInv2R_, "f");
0312     effInv2R_->SetTotalHistogram(*hisEffInv2RTotal_, "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 numTracksLost = prof_->GetBinContent(3);
0318     const double totalTracks = prof_->GetBinContent(5);
0319     const double numTracksMatched = prof_->GetBinContent(4);
0320     const double numTPsAll = prof_->GetBinContent(6);
0321     const double numTPsEff = prof_->GetBinContent(7);
0322     const double numTPsLost = prof_->GetBinContent(8);
0323     const double errStubs = prof_->GetBinError(1);
0324     const double errTracks = prof_->GetBinError(2);
0325     const double errTracksLost = prof_->GetBinError(3);
0326     const double fracFake = (totalTracks - numTracksMatched) / totalTracks;
0327     const double fracDup = (numTracksMatched - numTPsAll) / totalTracks;
0328     const double eff = numTPsEff / totalTPs;
0329     const double errEff = sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0330     const double effLoss = numTPsLost / totalTPs;
0331     const double errEffLoss = sqrt(effLoss * (1. - effLoss) / totalTPs / nEvents_);
0332     const int numStates = prof_->GetBinContent(10);
0333     const int numStatesLost = prof_->GetBinContent(11);
0334     const double fracSatest = numStates / (double)(numStates + numStatesLost);
0335     const vector<double> nums = {numStubs, numTracks, numTracksLost};
0336     const vector<double> errs = {errStubs, errTracks, errTracksLost};
0337     const int wNums = ceil(log10(*max_element(nums.begin(), nums.end()))) + 5;
0338     const int wErrs = ceil(log10(*max_element(errs.begin(), errs.end()))) + 5;
0339     log_ << "                         KF  SUMMARY                         " << endl;
0340     log_ << "number of stubs       per TFP = " << setw(wNums) << numStubs << " +- " << setw(wErrs) << errStubs << endl;
0341     log_ << "number of tracks      per TFP = " << setw(wNums) << numTracks << " +- " << setw(wErrs) << errTracks
0342          << endl;
0343     log_ << "number of lost tracks per TFP = " << setw(wNums) << numTracksLost << " +- " << setw(wErrs) << errTracksLost
0344          << endl;
0345     log_ << "          tracking efficiency = " << setw(wNums) << eff << " +- " << setw(wErrs) << errEff << endl;
0346     log_ << "     lost tracking efficiency = " << setw(wNums) << effLoss << " +- " << setw(wErrs) << errEffLoss << endl;
0347     log_ << "                    fake rate = " << setw(wNums) << fracFake << endl;
0348     log_ << "               duplicate rate = " << setw(wNums) << fracDup << endl;
0349     log_ << "    state assessment fraction = " << setw(wNums) << fracSatest << endl;
0350     log_ << "=============================================================";
0351     LogPrint("L1Trigger/TrackerTFP") << log_.str();
0352   }
0353 
0354   //
0355   void AnalyzerKF::associate(const TTTracks& ttTracks,
0356                              const StubAssociation* ass,
0357                              set<TPPtr>& tps,
0358                              int& sum,
0359                              const vector<TH1F*>& his,
0360                              TProfile* prof) const {
0361     for (const TTTrack<Ref_Phase2TrackerDigi_>& ttTrack : ttTracks) {
0362       const vector<TTStubRef>& ttStubRefs = ttTrack.getStubRefs();
0363       const vector<TPPtr>& tpPtrs = ass->associateFinal(ttStubRefs);
0364       if (tpPtrs.empty())
0365         continue;
0366       sum++;
0367       copy(tpPtrs.begin(), tpPtrs.end(), inserter(tps, tps.begin()));
0368       if (his.empty())
0369         continue;
0370       for (const TPPtr& tpPtr : tpPtrs) {
0371         const double phi0 = tpPtr->phi();
0372         const double cot = sinh(tpPtr->eta());
0373         const double inv2R = setup_->invPtToDphi() * tpPtr->charge() / tpPtr->pt();
0374         const math::XYZPointD& v = tpPtr->vertex();
0375         const double z0 = v.z() - cot * (v.x() * cos(phi0) + v.y() * sin(phi0));
0376         const double dCot = cot - ttTrack.tanL();
0377         const double dZ0 = z0 - ttTrack.z0();
0378         const double dInv2R = inv2R - ttTrack.rInv();
0379         const double dPhi0 = deltaPhi(phi0 - ttTrack.phi());
0380         const vector<double> ds = {dPhi0, dInv2R, dZ0, dCot};
0381         for (int i = 0; i < (int)ds.size(); i++)
0382           his[i]->Fill(ds[i]);
0383         prof->Fill(abs(tpPtr->eta()), abs(dZ0));
0384       }
0385     }
0386   }
0387 
0388 }  // namespace trackerTFP
0389 
0390 DEFINE_FWK_MODULE(trackerTFP::AnalyzerKF);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!