TrackerTFP/test/AnalyzerGP.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
0018 #include <TProfile.h>
0019 #include <TH1F.h>
0020
0021 #include <vector>
0022 #include <set>
0023 #include <numeric>
0024 #include <sstream>
0025
0026 using namespace std;
0027 using namespace edm;
0028 using namespace tt;
0029
0030 namespace trackerTFP {
0031
0032   /*! \class  trackerTFP::AnalyzerGP
0033    *  \brief  Class to analyze hardware like structured TTStub Collection generated by Geometric Processor
0034    *  \author Thomas Schuh
0035    *  \date   2020, Apr
0036    */
0037   class AnalyzerGP : public one::EDAnalyzer<one::WatchRuns, one::SharedResources> {
0038   public:
0039     AnalyzerGP(const ParameterSet& iConfig);
0040     void beginJob() override {}
0041     void beginRun(const Run& iEvent, const EventSetup& iSetup) override;
0042     void analyze(const Event& iEvent, const EventSetup& iSetup) override;
0043     void endRun(const Run& iEvent, const EventSetup& iSetup) override {}
0044     void endJob() override;
0045
0046   private:
0047     // ED input token of stubs
0048     EDGetTokenT<StreamsStub> edGetTokenAccepted_;
0049     // ED input token of lost stubs
0050     EDGetTokenT<StreamsStub> edGetTokenLost_;
0051     // ED input token of TTStubRef to TPPtr association for tracking efficiency
0052     EDGetTokenT<StubAssociation> edGetTokenAss_;
0053     // Setup token
0054     ESGetToken<Setup, SetupRcd> esGetToken_;
0055     // stores, calculates and provides run-time constants
0056     const Setup* setup_ = nullptr;
0057     // enables analyze of TPs
0058     bool useMCTruth_;
0059     //
0060     int nEvents_ = 0;
0061
0062     // Histograms
0063
0064     TProfile* prof_;
0065     TProfile* profChannel_;
0066     TH1F* hisChannel_;
0067
0068     // printout
0069     stringstream log_;
0070   };
0071
0072   AnalyzerGP::AnalyzerGP(const ParameterSet& iConfig) : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")) {
0073     usesResource("TFileService");
0074     // book in- and output ED products
0075     const string& label = iConfig.getParameter<string>("LabelGP");
0076     const string& branchAccepted = iConfig.getParameter<string>("BranchAcceptedStubs");
0077     const string& branchLost = iConfig.getParameter<string>("BranchLostStubs");
0078     edGetTokenAccepted_ = consumes<StreamsStub>(InputTag(label, branchAccepted));
0079     edGetTokenLost_ = consumes<StreamsStub>(InputTag(label, branchLost));
0080     if (useMCTruth_) {
0081       const auto& inputTagAss = iConfig.getParameter<InputTag>("InputTagSelection");
0082       edGetTokenAss_ = consumes<StubAssociation>(inputTagAss);
0083     }
0084     // book ES product
0085     esGetToken_ = esConsumes<Setup, SetupRcd, Transition::BeginRun>();
0086     // log config
0087     log_.setf(ios::fixed, ios::floatfield);
0088     log_.precision(4);
0089   }
0090
0091   void AnalyzerGP::beginRun(const Run& iEvent, const EventSetup& iSetup) {
0092     // helper class to store configurations
0093     setup_ = &iSetup.getData(esGetToken_);
0094     // book histograms
0095     Service<TFileService> fs;
0096     TFileDirectory dir;
0097     dir = fs->mkdir("GP");
0098     prof_ = dir.make<TProfile>("Counts", ";", 4, 0.5, 4.5);
0099     prof_->GetXaxis()->SetBinLabel(1, "Stubs");
0100     prof_->GetXaxis()->SetBinLabel(2, "Lost Stubs");
0101     prof_->GetXaxis()->SetBinLabel(3, "Found TPs");
0102     prof_->GetXaxis()->SetBinLabel(4, "Selected TPs");
0103     // channel occupancy
0104     constexpr int maxOcc = 180;
0105     const int numChannels = setup_->numSectors();
0106     hisChannel_ = dir.make<TH1F>("His Channel Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0107     profChannel_ = dir.make<TProfile>("Prof Channel Occupancy", ";", numChannels, -.5, numChannels - .5);
0108   }
0109
0110   void AnalyzerGP::analyze(const Event& iEvent, const EventSetup& iSetup) {
0111     // read in gp products
0112     Handle<StreamsStub> handleAccepted;
0113     iEvent.getByToken<StreamsStub>(edGetTokenAccepted_, handleAccepted);
0114     Handle<StreamsStub> handleLost;
0115     iEvent.getByToken<StreamsStub>(edGetTokenLost_, handleLost);
0116     // read in MCTruth
0117     const StubAssociation* stubAssociation = nullptr;
0118     if (useMCTruth_) {
0119       Handle<StubAssociation> handleAss;
0120       iEvent.getByToken<StubAssociation>(edGetTokenAss_, handleAss);
0121       stubAssociation = handleAss.product();
0122       prof_->Fill(4, stubAssociation->numTPs());
0123     }
0124     // analyze gp products and find still reconstrucable TrackingParticles
0125     set<TPPtr> setTPPtr;
0126     for (int region = 0; region < setup_->numRegions(); region++) {
0127       int nStubs(0);
0128       int nLost(0);
0129       map<TPPtr, vector<TTStubRef>> mapTPsTTStubs;
0130       for (int channel = 0; channel < setup_->numSectors(); channel++) {
0131         const int index = region * setup_->numSectors() + channel;
0132         const StreamStub& accepted = handleAccepted->at(index);
0133         hisChannel_->Fill(accepted.size());
0134         profChannel_->Fill(channel, accepted.size());
0135         for (const FrameStub& frame : accepted) {
0136           if (frame.first.isNull())
0137             continue;
0138           nStubs++;
0139           if (!useMCTruth_)
0140             continue;
0141           const vector<TPPtr>& tpPtrs = stubAssociation->findTrackingParticlePtrs(frame.first);
0142           for (const TPPtr& tpPtr : tpPtrs) {
0143             auto it = mapTPsTTStubs.find(tpPtr);
0144             if (it == mapTPsTTStubs.end()) {
0145               it = mapTPsTTStubs.emplace(tpPtr, vector<TTStubRef>()).first;
0146               it->second.reserve(stubAssociation->findTTStubRefs(tpPtr).size());
0147             }
0148             it->second.push_back(frame.first);
0149           }
0150         }
0151         nLost += handleLost->at(index).size();
0152       }
0153       for (const auto& p : mapTPsTTStubs)
0154         if (setup_->reconstructable(p.second))
0155           setTPPtr.insert(p.first);
0156       prof_->Fill(1, nStubs);
0157       prof_->Fill(2, nLost);
0158     }
0159     prof_->Fill(3, setTPPtr.size());
0160     nEvents_++;
0161   }
0162
0163   void AnalyzerGP::endJob() {
0164     if (nEvents_ == 0)
0165       return;
0166     // printout GP summary
0167     const double numStubs = prof_->GetBinContent(1);
0168     const double numStubsLost = prof_->GetBinContent(2);
0169     const double errStubs = prof_->GetBinError(1);
0170     const double errStubsLost = prof_->GetBinError(2);
0171     const double numTPs = prof_->GetBinContent(3);
0172     const double totalTPs = prof_->GetBinContent(4);
0173     const double eff = numTPs / totalTPs;
0174     const double errEff = sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0175     const vector<double> nums = {numStubs, numStubsLost};
0176     const vector<double> errs = {errStubs, errStubsLost};
0177     const int wNums = ceil(log10(*max_element(nums.begin(), nums.end()))) + 5;
0178     const int wErrs = ceil(log10(*max_element(errs.begin(), errs.end()))) + 5;
0179     log_ << "                         GP  SUMMARY                         " << endl;
0180     log_ << "number of stubs      per TFP = " << setw(wNums) << numStubs << " +- " << setw(wErrs) << errStubs << endl;
0181     log_ << "number of lost stubs per TFP = " << setw(wNums) << numStubsLost << " +- " << setw(wErrs) << errStubsLost
0182          << endl;
0183     log_ << "     max tracking efficiency = " << setw(wNums) << eff << " +- " << setw(wErrs) << errEff << endl;
0184     log_ << "=============================================================";
0185     LogPrint("L1Trigger/TrackerTFP") << log_.str();
0186   }
0187
0188 }  // namespace trackerTFP
0189
0190 DEFINE_FWK_MODULE(trackerTFP::AnalyzerGP);