Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-06-03 00:12:18 

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 #include <TEfficiency.h>
0022 
0023 #include <vector>
0024 #include <deque>
0025 #include <set>
0026 #include <cmath>
0027 #include <numeric>
0028 #include <sstream>
0029 
0030 namespace trackerTFP {
0031 
0032   /*! \class  trackerTFP::AnalyzerHT
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 AnalyzerHT : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0038   public:
0039     AnalyzerHT(const edm::ParameterSet& iConfig);
0040     void beginJob() override {}
0041     void beginRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) override;
0042     void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0043     void endRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) override {}
0044     void endJob() override;
0045 
0046   private:
0047     //
0048     void formTracks(const tt::StreamStub& stream, std::vector<std::vector<TTStubRef>>& tracks) const;
0049     //
0050     void associate(const std::vector<std::vector<TTStubRef>>& tracks,
0051                    const tt::StubAssociation* ass,
0052                    std::set<TPPtr>& tps,
0053                    int& sum) const;
0054 
0055     // ED input token of stubs
0056     edm::EDGetTokenT<tt::StreamsStub> edGetToken_;
0057     // ED input token of TTStubRef to TPPtr association for tracking efficiency
0058     edm::EDGetTokenT<tt::StubAssociation> edGetTokenSelection_;
0059     // ED input token of TTStubRef to recontructable TPPtr association
0060     edm::EDGetTokenT<tt::StubAssociation> edGetTokenReconstructable_;
0061     // Setup token
0062     edm::ESGetToken<tt::Setup, tt::SetupRcd> esGetTokenSetup_;
0063     // DataFormats token
0064     edm::ESGetToken<DataFormats, DataFormatsRcd> esGetTokenDataFormats_;
0065     // stores, calculates and provides run-time constants
0066     const tt::Setup* setup_ = nullptr;
0067     // helper class to extract structured data from tt::Frames
0068     const DataFormats* dataFormats_ = nullptr;
0069     // enables analyze of TPs
0070     bool useMCTruth_;
0071     //
0072     int nEvents_ = 0;
0073 
0074     // Histograms
0075 
0076     TProfile* prof_;
0077     TProfile* profChannel_;
0078     TH1F* hisChannel_;
0079     TH1F* hisLayers_;
0080     TH1F* hisNumLayers_;
0081     TProfile* profNumLayers_;
0082     TH1F* hisEffEta_;
0083     TH1F* hisEffEtaTotal_;
0084     TEfficiency* effEta_;
0085     TH1F* hisEffZT_;
0086     TH1F* hisEffZTTotal_;
0087     TEfficiency* effZT_;
0088     TH1F* hisEffInv2R_;
0089     TH1F* hisEffInv2RTotal_;
0090     TEfficiency* effInv2R_;
0091     TH1F* hisEffPT_;
0092     TH1F* hisEffPTTotal_;
0093     TEfficiency* effPT_;
0094 
0095     // printout
0096     std::stringstream log_;
0097   };
0098 
0099   AnalyzerHT::AnalyzerHT(const edm::ParameterSet& iConfig) : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")) {
0100     usesResource("TFileService");
0101     // book in- and output ED products
0102     const std::string& label = iConfig.getParameter<std::string>("OutputLabelHT");
0103     const std::string& branch = iConfig.getParameter<std::string>("BranchStubs");
0104     edGetToken_ = consumes<tt::StreamsStub>(edm::InputTag(label, branch));
0105     if (useMCTruth_) {
0106       const auto& inputTagSelecttion = iConfig.getParameter<edm::InputTag>("InputTagSelection");
0107       const auto& inputTagReconstructable = iConfig.getParameter<edm::InputTag>("InputTagReconstructable");
0108       edGetTokenSelection_ = consumes<tt::StubAssociation>(inputTagSelecttion);
0109       edGetTokenReconstructable_ = consumes<tt::StubAssociation>(inputTagReconstructable);
0110     }
0111     // book ES products
0112     esGetTokenSetup_ = esConsumes<edm::Transition::BeginRun>();
0113     esGetTokenDataFormats_ = esConsumes<edm::Transition::BeginRun>();
0114     // log config
0115     log_.setf(std::ios::fixed, std::ios::floatfield);
0116     log_.precision(4);
0117   }
0118 
0119   void AnalyzerHT::beginRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) {
0120     // helper class to store configurations
0121     setup_ = &iSetup.getData(esGetTokenSetup_);
0122     // helper class to extract structured data from tt::Frames
0123     dataFormats_ = &iSetup.getData(esGetTokenDataFormats_);
0124     // book histograms
0125     edm::Service<TFileService> fs;
0126     TFileDirectory dir;
0127     dir = fs->mkdir("HT");
0128     prof_ = dir.make<TProfile>("Counts", ";", 9, 0.5, 9.5);
0129     prof_->GetXaxis()->SetBinLabel(1, "Stubs");
0130     prof_->GetXaxis()->SetBinLabel(2, "Tracks");
0131     prof_->GetXaxis()->SetBinLabel(3, "Truncated Tracks");
0132     prof_->GetXaxis()->SetBinLabel(4, "Matched Tracks");
0133     prof_->GetXaxis()->SetBinLabel(5, "All Tracks");
0134     prof_->GetXaxis()->SetBinLabel(6, "Found TPs");
0135     prof_->GetXaxis()->SetBinLabel(7, "Found selected TPs");
0136     prof_->GetXaxis()->SetBinLabel(8, "Truncated TPs");
0137     prof_->GetXaxis()->SetBinLabel(9, "All TPs");
0138     // Efficiencies
0139     hisEffEtaTotal_ = dir.make<TH1F>("HisTPEtaTotal", ";", 48, -2.4, 2.4);
0140     hisEffEta_ = dir.make<TH1F>("HisTPEta", ";", 48, -2.4, 2.4);
0141     effEta_ = dir.make<TEfficiency>("EffEta", ";", 48, -2.4, 2.4);
0142     const double rangeZT = dataFormats_->format(Variable::zT, Process::dr).range();
0143     const int zTBins = setup_->gpNumBinsZT();
0144     hisEffZTTotal_ = dir.make<TH1F>("HisTPZTTotal", ";", zTBins, -rangeZT / 2, rangeZT / 2);
0145     hisEffZT_ = dir.make<TH1F>("HisTPZT", ";", zTBins, -rangeZT / 2, rangeZT / 2);
0146     effZT_ = dir.make<TEfficiency>("EffZT", ";", zTBins, -rangeZT / 2, rangeZT / 2);
0147     const double rangeInv2R = dataFormats_->format(Variable::inv2R, Process::dr).range();
0148     const int inv2RBins = (setup_->htNumBinsInv2R() + 2) * 2;
0149     hisEffInv2R_ = dir.make<TH1F>("HisTPInv2R", ";", inv2RBins, -rangeInv2R / 2., rangeInv2R / 2.);
0150     hisEffInv2RTotal_ = dir.make<TH1F>("HisTPInv2RTotal", ";", inv2RBins, -rangeInv2R / 2., rangeInv2R / 2.);
0151     effInv2R_ = dir.make<TEfficiency>("EffInv2R", ";", inv2RBins, -rangeInv2R / 2., rangeInv2R / 2.);
0152     hisEffPT_ = dir.make<TH1F>("HisTPPT", ";", 100, 0, 100);
0153     hisEffPTTotal_ = dir.make<TH1F>("HisTPPTTotal", ";", 100, 0, 100);
0154     effPT_ = dir.make<TEfficiency>("EffPT", ";", 100, 0, 100);
0155     // binInv2R occupancy
0156     constexpr int maxOcc = 180;
0157     const int numChannel = dataFormats_->numChannel(Process::ht);
0158     hisChannel_ = dir.make<TH1F>("His binInv2R Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0159     profChannel_ = dir.make<TProfile>("Prof binInv2R Occupancy", ";", numChannel, -.5, numChannel - .5);
0160     // layers
0161     hisLayers_ = dir.make<TH1F>("HisLayers", ";", 8, 0, 8);
0162     hisNumLayers_ = dir.make<TH1F>("HisNumLayers", ";", 9, 0, 9);
0163     profNumLayers_ = dir.make<TProfile>("Prof NumLayers", ";", 32, 0, 2.4);
0164   }
0165 
0166   void AnalyzerHT::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0167     auto fill = [this](const TPPtr& tpPtr, TH1F* hisEta, TH1F* hisZT, TH1F* hisInv2R, TH1F* hisPT) {
0168       const double tpPhi0 = tpPtr->phi();
0169       const double tpCot = sinh(tpPtr->eta());
0170       const math::XYZPointD& v = tpPtr->vertex();
0171       const double tpZ0 = v.z() - tpCot * (v.x() * cos(tpPhi0) + v.y() * sin(tpPhi0));
0172       const double tpZT = tpZ0 + tpCot * setup_->chosenRofZ();
0173       hisEta->Fill(tpPtr->eta());
0174       hisZT->Fill(tpZT);
0175       hisInv2R->Fill(tpPtr->charge() / tpPtr->pt() * setup_->invPtToDphi());
0176       hisPT->Fill(tpPtr->pt());
0177     };
0178     // read in ht products
0179     edm::Handle<tt::StreamsStub> handleStubs;
0180     iEvent.getByToken<tt::StreamsStub>(edGetToken_, handleStubs);
0181     // read in MCTruth
0182     const tt::StubAssociation* selection = nullptr;
0183     const tt::StubAssociation* reconstructable = nullptr;
0184     if (useMCTruth_) {
0185       edm::Handle<tt::StubAssociation> handleSelection;
0186       iEvent.getByToken<tt::StubAssociation>(edGetTokenSelection_, handleSelection);
0187       selection = handleSelection.product();
0188       prof_->Fill(9, selection->numTPs());
0189       edm::Handle<tt::StubAssociation> handleReconstructable;
0190       iEvent.getByToken<tt::StubAssociation>(edGetTokenReconstructable_, handleReconstructable);
0191       reconstructable = handleReconstructable.product();
0192       for (const auto& p : selection->getTrackingParticleToTTStubsMap())
0193         fill(p.first, hisEffEtaTotal_, hisEffZTTotal_, hisEffInv2RTotal_, hisEffPTTotal_);
0194     }
0195     // analyze ht products and associate found tracks with reconstrucable TrackingParticles
0196     std::set<TPPtr> tpPtrs;
0197     std::set<TPPtr> tpPtrsSelection;
0198     int allMatched(0);
0199     int allTracks(0);
0200     for (int region = 0; region < setup_->numRegions(); region++) {
0201       int nStubs(0);
0202       int nTracks(0);
0203       for (int channel = 0; channel < dataFormats_->numChannel(Process::ht); channel++) {
0204         const int index = region * dataFormats_->numChannel(Process::ht) + channel;
0205         const tt::StreamStub& accepted = handleStubs->at(index);
0206         hisChannel_->Fill(accepted.size());
0207         profChannel_->Fill(channel, accepted.size());
0208         nStubs += accepted.size();
0209         std::vector<std::vector<TTStubRef>> tracks;
0210         formTracks(accepted, tracks);
0211         nTracks += tracks.size();
0212         allTracks += tracks.size();
0213         if (!useMCTruth_)
0214           continue;
0215         int tmp(0);
0216         associate(tracks, selection, tpPtrsSelection, tmp);
0217         associate(tracks, reconstructable, tpPtrs, allMatched);
0218       }
0219       prof_->Fill(1, nStubs);
0220       prof_->Fill(2, nTracks);
0221     }
0222     for (const TPPtr& tpPtr : tpPtrsSelection)
0223       fill(tpPtr, hisEffEta_, hisEffZT_, hisEffInv2R_, hisEffPT_);
0224     prof_->Fill(4, allMatched);
0225     prof_->Fill(5, allTracks);
0226     prof_->Fill(6, tpPtrs.size());
0227     prof_->Fill(7, tpPtrsSelection.size());
0228     nEvents_++;
0229   }
0230 
0231   void AnalyzerHT::endJob() {
0232     if (nEvents_ == 0)
0233       return;
0234     // effi
0235     effEta_->SetPassedHistogram(*hisEffEta_, "f");
0236     effEta_->SetTotalHistogram(*hisEffEtaTotal_, "f");
0237     effZT_->SetPassedHistogram(*hisEffZT_, "f");
0238     effZT_->SetTotalHistogram(*hisEffZTTotal_, "f");
0239     effInv2R_->SetPassedHistogram(*hisEffInv2R_, "f");
0240     effInv2R_->SetTotalHistogram(*hisEffInv2RTotal_, "f");
0241     effPT_->SetPassedHistogram(*hisEffPT_, "f");
0242     effPT_->SetTotalHistogram(*hisEffPTTotal_, "f");
0243     // printout HT summary
0244     const double totalTPs = prof_->GetBinContent(9);
0245     const double numStubs = prof_->GetBinContent(1);
0246     const double numTracks = prof_->GetBinContent(2);
0247     const double totalTracks = prof_->GetBinContent(5);
0248     const double numTracksMatched = prof_->GetBinContent(4);
0249     const double numTPsAll = prof_->GetBinContent(6);
0250     const double numTPsEff = prof_->GetBinContent(7);
0251     const double errStubs = prof_->GetBinError(1);
0252     const double errTracks = prof_->GetBinError(2);
0253     const double fracFake = (totalTracks - numTracksMatched) / totalTracks;
0254     const double fracDup = (numTracksMatched - numTPsAll) / totalTracks;
0255     const double eff = numTPsEff / totalTPs;
0256     const double errEff = sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0257     const std::vector<double> nums = {numStubs, numTracks};
0258     const std::vector<double> errs = {errStubs, errTracks};
0259     const int wNums = std::ceil(std::log10(*std::max_element(nums.begin(), nums.end()))) + 5;
0260     const int wErrs = std::ceil(std::log10(*std::max_element(errs.begin(), errs.end()))) + 5;
0261     log_ << "                         HT  SUMMARY                         " << std::endl;
0262     log_ << "number of stubs       per TFP = " << std::setw(wNums) << numStubs << " +- " << std::setw(wErrs) << errStubs
0263          << std::endl;
0264     log_ << "number of tracks      per TFP = " << std::setw(wNums) << numTracks << " +- " << std::setw(wErrs)
0265          << errTracks << std::endl;
0266     log_ << "     max  tracking efficiency = " << std::setw(wNums) << eff << " +- " << std::setw(wErrs) << errEff
0267          << std::endl;
0268     log_ << "                    fake rate = " << std::setw(wNums) << fracFake << std::endl;
0269     log_ << "               duplicate rate = " << std::setw(wNums) << fracDup << std::endl;
0270     log_ << "=============================================================";
0271     edm::LogPrint(moduleDescription().moduleName()) << log_.str();
0272   }
0273 
0274   //
0275   void AnalyzerHT::formTracks(const tt::StreamStub& stream, std::vector<std::vector<TTStubRef>>& tracks) const {
0276     static const DataFormat& layer = dataFormats_->format(Variable::layer, Process::ctb);
0277     auto toTrkId = [this](const StubHT& stub) {
0278       static const DataFormat& phiT = dataFormats_->format(Variable::phiT, Process::ht);
0279       static const DataFormat& zT = dataFormats_->format(Variable::zT, Process::ht);
0280       return (phiT.ttBV(stub.phiT()) + zT.ttBV(stub.zT())).val();
0281     };
0282     std::vector<StubHT> stubs;
0283     stubs.reserve(stream.size());
0284     for (const tt::FrameStub& frame : stream)
0285       stubs.emplace_back(frame, dataFormats_);
0286     for (auto it = stubs.begin(); it != stubs.end();) {
0287       const auto start = it;
0288       const int id = toTrkId(*it);
0289       auto different = [id, toTrkId](const StubHT& stub) { return id != toTrkId(stub); };
0290       it = std::find_if(it, stubs.end(), different);
0291       std::vector<TTStubRef> ttStubRefs;
0292       ttStubRefs.reserve(std::distance(start, it));
0293       std::transform(start, it, std::back_inserter(ttStubRefs), [](const StubHT& stub) { return stub.frame().first; });
0294       tracks.push_back(ttStubRefs);
0295       TTBV hitPattern(0, setup_->numLayers());
0296       for (auto iter = start; iter != it; iter++)
0297         hitPattern.set(iter->layer().val(layer.width()));
0298       const double cot = dataFormats_->format(Variable::zT, Process::ht).floating(start->zT()) / setup_->chosenRofZ();
0299       hisNumLayers_->Fill(hitPattern.count());
0300       profNumLayers_->Fill(abs(sinh(cot)), hitPattern.count());
0301       for (int layer : hitPattern.ids())
0302         hisLayers_->Fill(layer);
0303     }
0304   }
0305 
0306   //
0307   void AnalyzerHT::associate(const std::vector<std::vector<TTStubRef>>& tracks,
0308                              const tt::StubAssociation* ass,
0309                              std::set<TPPtr>& tps,
0310                              int& sum) const {
0311     for (const std::vector<TTStubRef>& ttStubRefs : tracks) {
0312       const std::vector<TPPtr>& tpPtrs = ass->associate(ttStubRefs);
0313       if (tpPtrs.empty())
0314         continue;
0315       sum++;
0316       std::copy(tpPtrs.begin(), tpPtrs.end(), std::inserter(tps, tps.begin()));
0317     }
0318   }
0319 
0320 }  // namespace trackerTFP
0321 
0322 DEFINE_FWK_MODULE(trackerTFP::AnalyzerHT);
[ 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!