Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-10-14 01:43:56

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 
0019 #include <TProfile.h>
0020 #include <TH1F.h>
0021 
0022 #include <vector>
0023 #include <deque>
0024 #include <set>
0025 #include <cmath>
0026 #include <numeric>
0027 #include <sstream>
0028 
0029 using namespace std;
0030 using namespace edm;
0031 using namespace tt;
0032 
0033 namespace trackerTFP {
0034 
0035   /*! \class  trackerTFP::AnalyzerHT
0036    *  \brief  Class to analyze hardware like structured TTStub Collection generated by Geometric Processor
0037    *  \author Thomas Schuh
0038    *  \date   2020, Apr
0039    */
0040   class AnalyzerHT : public one::EDAnalyzer<one::WatchRuns, one::SharedResources> {
0041   public:
0042     AnalyzerHT(const ParameterSet& iConfig);
0043     void beginJob() override {}
0044     void beginRun(const Run& iEvent, const EventSetup& iSetup) override;
0045     void analyze(const Event& iEvent, const EventSetup& iSetup) override;
0046     void endRun(const Run& iEvent, const EventSetup& iSetup) override {}
0047     void endJob() override;
0048 
0049   private:
0050     //
0051     void formTracks(const StreamStub& stream, vector<vector<TTStubRef>>& tracks, int qOverPt) const;
0052     //
0053     void associate(const vector<vector<TTStubRef>>& tracks, const StubAssociation* ass, set<TPPtr>& tps, int& sum) const;
0054 
0055     // ED input token of stubs
0056     EDGetTokenT<StreamsStub> edGetTokenAccepted_;
0057     // ED input token of lost stubs
0058     EDGetTokenT<StreamsStub> edGetTokenLost_;
0059     // ED input token of TTStubRef to TPPtr association for tracking efficiency
0060     EDGetTokenT<StubAssociation> edGetTokenSelection_;
0061     // ED input token of TTStubRef to recontructable TPPtr association
0062     EDGetTokenT<StubAssociation> edGetTokenReconstructable_;
0063     // Setup token
0064     ESGetToken<Setup, SetupRcd> esGetTokenSetup_;
0065     // DataFormats token
0066     ESGetToken<DataFormats, DataFormatsRcd> esGetTokenDataFormats_;
0067     // stores, calculates and provides run-time constants
0068     const Setup* setup_ = nullptr;
0069     // helper class to extract structured data from tt::Frames
0070     const DataFormats* dataFormats_ = nullptr;
0071     // enables analyze of TPs
0072     bool useMCTruth_;
0073     //
0074     int nEvents_ = 0;
0075 
0076     // Histograms
0077 
0078     TProfile* prof_;
0079     TProfile* profChannel_;
0080     TH1F* hisChannel_;
0081 
0082     // printout
0083     stringstream log_;
0084   };
0085 
0086   AnalyzerHT::AnalyzerHT(const ParameterSet& iConfig) : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")) {
0087     usesResource("TFileService");
0088     // book in- and output ED products
0089     const string& label = iConfig.getParameter<string>("LabelHT");
0090     const string& branchAccepted = iConfig.getParameter<string>("BranchAcceptedStubs");
0091     const string& branchLost = iConfig.getParameter<string>("BranchLostStubs");
0092     edGetTokenAccepted_ = consumes<StreamsStub>(InputTag(label, branchAccepted));
0093     edGetTokenLost_ = consumes<StreamsStub>(InputTag(label, branchLost));
0094     if (useMCTruth_) {
0095       const auto& inputTagSelecttion = iConfig.getParameter<InputTag>("InputTagSelection");
0096       const auto& inputTagReconstructable = iConfig.getParameter<InputTag>("InputTagReconstructable");
0097       edGetTokenSelection_ = consumes<StubAssociation>(inputTagSelecttion);
0098       edGetTokenReconstructable_ = consumes<StubAssociation>(inputTagReconstructable);
0099     }
0100     // book ES products
0101     esGetTokenSetup_ = esConsumes<Setup, SetupRcd, Transition::BeginRun>();
0102     esGetTokenDataFormats_ = esConsumes<DataFormats, DataFormatsRcd, Transition::BeginRun>();
0103     // log config
0104     log_.setf(ios::fixed, ios::floatfield);
0105     log_.precision(4);
0106   }
0107 
0108   void AnalyzerHT::beginRun(const Run& iEvent, const EventSetup& iSetup) {
0109     // helper class to store configurations
0110     setup_ = &iSetup.getData(esGetTokenSetup_);
0111     // helper class to extract structured data from tt::Frames
0112     dataFormats_ = &iSetup.getData(esGetTokenDataFormats_);
0113     // book histograms
0114     Service<TFileService> fs;
0115     TFileDirectory dir;
0116     dir = fs->mkdir("HT");
0117     prof_ = dir.make<TProfile>("Counts", ";", 9, 0.5, 9.5);
0118     prof_->GetXaxis()->SetBinLabel(1, "Stubs");
0119     prof_->GetXaxis()->SetBinLabel(2, "Tracks");
0120     prof_->GetXaxis()->SetBinLabel(3, "Lost Tracks");
0121     prof_->GetXaxis()->SetBinLabel(4, "Matched Tracks");
0122     prof_->GetXaxis()->SetBinLabel(5, "All Tracks");
0123     prof_->GetXaxis()->SetBinLabel(6, "Found TPs");
0124     prof_->GetXaxis()->SetBinLabel(7, "Found selected TPs");
0125     prof_->GetXaxis()->SetBinLabel(8, "Lost TPs");
0126     prof_->GetXaxis()->SetBinLabel(9, "All TPs");
0127     // binInv2R occupancy
0128     constexpr int maxOcc = 180;
0129     const int numChannel = dataFormats_->numChannel(Process::ht);
0130     hisChannel_ = dir.make<TH1F>("His binInv2R Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0131     profChannel_ = dir.make<TProfile>("Prof binInv2R Occupancy", ";", numChannel, -.5, numChannel - .5);
0132   }
0133 
0134   void AnalyzerHT::analyze(const Event& iEvent, const EventSetup& iSetup) {
0135     // read in ht products
0136     Handle<StreamsStub> handleAccepted;
0137     iEvent.getByToken<StreamsStub>(edGetTokenAccepted_, handleAccepted);
0138     Handle<StreamsStub> handleLost;
0139     iEvent.getByToken<StreamsStub>(edGetTokenLost_, handleLost);
0140     // read in MCTruth
0141     const StubAssociation* selection = nullptr;
0142     const StubAssociation* reconstructable = nullptr;
0143     if (useMCTruth_) {
0144       Handle<StubAssociation> handleSelection;
0145       iEvent.getByToken<StubAssociation>(edGetTokenSelection_, handleSelection);
0146       selection = handleSelection.product();
0147       prof_->Fill(9, selection->numTPs());
0148       Handle<StubAssociation> handleReconstructable;
0149       iEvent.getByToken<StubAssociation>(edGetTokenReconstructable_, handleReconstructable);
0150       reconstructable = handleReconstructable.product();
0151     }
0152     // analyze ht products and associate found tracks with reconstrucable TrackingParticles
0153     set<TPPtr> tpPtrs;
0154     set<TPPtr> tpPtrsSelection;
0155     set<TPPtr> tpPtrsLost;
0156     int allMatched(0);
0157     int allTracks(0);
0158     for (int region = 0; region < setup_->numRegions(); region++) {
0159       int nStubs(0);
0160       int nTracks(0);
0161       int nLost(0);
0162       for (int channel = 0; channel < dataFormats_->numChannel(Process::ht); channel++) {
0163         const int inv2R = dataFormats_->format(Variable::inv2R, Process::ht).toSigned(channel);
0164         const int index = region * dataFormats_->numChannel(Process::ht) + channel;
0165         const StreamStub& accepted = handleAccepted->at(index);
0166         hisChannel_->Fill(accepted.size());
0167         profChannel_->Fill(channel, accepted.size());
0168         nStubs += accepted.size();
0169         vector<vector<TTStubRef>> tracks;
0170         vector<vector<TTStubRef>> lost;
0171         formTracks(accepted, tracks, inv2R);
0172         formTracks(handleLost->at(index), lost, inv2R);
0173         nTracks += tracks.size();
0174         allTracks += tracks.size();
0175         nLost += lost.size();
0176         if (!useMCTruth_)
0177           continue;
0178         int tmp(0);
0179         associate(tracks, selection, tpPtrsSelection, tmp);
0180         associate(lost, selection, tpPtrsLost, tmp);
0181         associate(tracks, reconstructable, tpPtrs, allMatched);
0182       }
0183       prof_->Fill(1, nStubs);
0184       prof_->Fill(2, nTracks);
0185       prof_->Fill(3, nLost);
0186     }
0187     vector<TPPtr> recovered;
0188     recovered.reserve(tpPtrsLost.size());
0189     set_intersection(tpPtrsLost.begin(), tpPtrsLost.end(), tpPtrs.begin(), tpPtrs.end(), back_inserter(recovered));
0190     for (const TPPtr& tpPtr : recovered)
0191       tpPtrsLost.erase(tpPtr);
0192     prof_->Fill(4, allMatched);
0193     prof_->Fill(5, allTracks);
0194     prof_->Fill(6, tpPtrs.size());
0195     prof_->Fill(7, tpPtrsSelection.size());
0196     prof_->Fill(8, tpPtrsLost.size());
0197     nEvents_++;
0198   }
0199 
0200   void AnalyzerHT::endJob() {
0201     if (nEvents_ == 0)
0202       return;
0203     // printout HT summary
0204     const double totalTPs = prof_->GetBinContent(9);
0205     const double numStubs = prof_->GetBinContent(1);
0206     const double numTracks = prof_->GetBinContent(2);
0207     const double numTracksLost = prof_->GetBinContent(3);
0208     const double totalTracks = prof_->GetBinContent(5);
0209     const double numTracksMatched = prof_->GetBinContent(4);
0210     const double numTPsAll = prof_->GetBinContent(6);
0211     const double numTPsEff = prof_->GetBinContent(7);
0212     const double numTPsLost = prof_->GetBinContent(8);
0213     const double errStubs = prof_->GetBinError(1);
0214     const double errTracks = prof_->GetBinError(2);
0215     const double errTracksLost = prof_->GetBinError(3);
0216     const double fracFake = (totalTracks - numTracksMatched) / totalTracks;
0217     const double fracDup = (numTracksMatched - numTPsAll) / totalTracks;
0218     const double eff = numTPsEff / totalTPs;
0219     const double errEff = sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0220     const double effLoss = numTPsLost / totalTPs;
0221     const double errEffLoss = sqrt(effLoss * (1. - effLoss) / totalTPs / nEvents_);
0222     const vector<double> nums = {numStubs, numTracks, numTracksLost};
0223     const vector<double> errs = {errStubs, errTracks, errTracksLost};
0224     const int wNums = ceil(log10(*max_element(nums.begin(), nums.end()))) + 5;
0225     const int wErrs = ceil(log10(*max_element(errs.begin(), errs.end()))) + 5;
0226     log_ << "                         HT  SUMMARY                         " << endl;
0227     log_ << "number of stubs       per TFP = " << setw(wNums) << numStubs << " +- " << setw(wErrs) << errStubs << endl;
0228     log_ << "number of tracks      per TFP = " << setw(wNums) << numTracks << " +- " << setw(wErrs) << errTracks
0229          << endl;
0230     log_ << "number of lost tracks per TFP = " << setw(wNums) << numTracksLost << " +- " << setw(wErrs) << errTracksLost
0231          << endl;
0232     log_ << "     max  tracking efficiency = " << setw(wNums) << eff << " +- " << setw(wErrs) << errEff << endl;
0233     log_ << "     lost tracking efficiency = " << setw(wNums) << effLoss << " +- " << setw(wErrs) << errEffLoss << endl;
0234     log_ << "                    fake rate = " << setw(wNums) << fracFake << endl;
0235     log_ << "               duplicate rate = " << setw(wNums) << fracDup << endl;
0236     log_ << "=============================================================";
0237     LogPrint("L1Trigger/TrackerTFP") << log_.str();
0238   }
0239 
0240   //
0241   void AnalyzerHT::formTracks(const StreamStub& stream, vector<vector<TTStubRef>>& tracks, int inv2R) const {
0242     vector<StubHT> stubs;
0243     stubs.reserve(stream.size());
0244     for (const FrameStub& frame : stream)
0245       stubs.emplace_back(frame, dataFormats_, inv2R);
0246     for (auto it = stubs.begin(); it != stubs.end();) {
0247       const auto start = it;
0248       const int id = it->trackId();
0249       auto different = [id](const StubHT& stub) { return id != stub.trackId(); };
0250       it = find_if(it, stubs.end(), different);
0251       vector<TTStubRef> ttStubRefs;
0252       ttStubRefs.reserve(distance(start, it));
0253       transform(start, it, back_inserter(ttStubRefs), [](const StubHT& stub) { return stub.ttStubRef(); });
0254       tracks.push_back(ttStubRefs);
0255     }
0256   }
0257 
0258   //
0259   void AnalyzerHT::associate(const vector<vector<TTStubRef>>& tracks,
0260                              const StubAssociation* ass,
0261                              set<TPPtr>& tps,
0262                              int& sum) const {
0263     for (const vector<TTStubRef>& ttStubRefs : tracks) {
0264       const vector<TPPtr>& tpPtrs = ass->associate(ttStubRefs);
0265       if (tpPtrs.empty())
0266         continue;
0267       sum++;
0268       copy(tpPtrs.begin(), tpPtrs.end(), inserter(tps, tps.begin()));
0269     }
0270   }
0271 
0272 }  // namespace trackerTFP
0273 
0274 DEFINE_FWK_MODULE(trackerTFP::AnalyzerHT);