TrackFindingTracklet/test/AnalyzerDRin.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/TrackerTFP/interface/DataFormats.h"
0018 #include "L1Trigger/TrackFindingTracklet/interface/ChannelAssignment.h"
0019
0020 #include <TProfile.h>
0021 #include <TH1F.h>
0022
0023 #include <vector>
0024 #include <deque>
0025 #include <map>
0026 #include <set>
0027 #include <cmath>
0028 #include <numeric>
0029 #include <sstream>
0030
0031 using namespace std;
0032 using namespace edm;
0033 using namespace trackerTFP;
0034 using namespace tt;
0035
0036 namespace trklet {
0037
0038   /*! \class  trklet::AnalyzerDRin
0039    *  \brief  Class to analyze hardware like structured TTStub Collection generated by DRin module
0040    *  \author Thomas Schuh
0041    *  \date   2023, Jan
0042    */
0043   class AnalyzerDRin : public one::EDAnalyzer<one::WatchRuns, one::SharedResources> {
0044   public:
0045     AnalyzerDRin(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 formTracks(const StreamsTrack& streamsTrack,
0055                     const StreamsStub& streamsStubs,
0056                     vector<vector<TTStubRef>>& tracks,
0057                     int channel) const;
0058     //
0059     void associate(const vector<vector<TTStubRef>>& tracks,
0060                    const StubAssociation* ass,
0061                    set<TPPtr>& tps,
0062                    int& sum,
0063                    bool perfect = false) const;
0064
0065     // ED input token of stubs
0066     EDGetTokenT<StreamsStub> edGetTokenAcceptedStubs_;
0067     // ED input token of tracks
0068     EDGetTokenT<StreamsTrack> edGetTokenAcceptedTracks_;
0069     // ED input token of lost stubs
0070     EDGetTokenT<StreamsStub> edGetTokenLostStubs_;
0071     // ED input token of lost tracks
0072     EDGetTokenT<StreamsTrack> edGetTokenLostTracks_;
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     // ChannelAssignment token
0082     ESGetToken<ChannelAssignment, ChannelAssignmentRcd> esGetTokenChannelAssignment_;
0083     // stores, calculates and provides run-time constants
0084     const Setup* setup_ = nullptr;
0085     // helper class to extract structured data from tt::Frames
0086     const DataFormats* dataFormats_ = nullptr;
0087     // helper class to assign tracklet track to channel
0088     const ChannelAssignment* channelAssignment_ = 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
0100     // printout
0101     stringstream log_;
0102   };
0103
0104   AnalyzerDRin::AnalyzerDRin(const ParameterSet& iConfig) : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")) {
0105     usesResource("TFileService");
0106     // book in- and output ED products
0107     const string& label = iConfig.getParameter<string>("LabelDRin");
0108     const string& branchAcceptedStubs = iConfig.getParameter<string>("BranchAcceptedStubs");
0109     const string& branchAcceptedTracks = iConfig.getParameter<string>("BranchAcceptedTracks");
0110     const string& branchLostStubs = iConfig.getParameter<string>("BranchLostStubs");
0111     const string& branchLostTracks = iConfig.getParameter<string>("BranchLostTracks");
0112     edGetTokenAcceptedStubs_ = consumes<StreamsStub>(InputTag(label, branchAcceptedStubs));
0113     edGetTokenAcceptedTracks_ = consumes<StreamsTrack>(InputTag(label, branchAcceptedTracks));
0114     edGetTokenLostStubs_ = consumes<StreamsStub>(InputTag(label, branchLostStubs));
0115     edGetTokenLostTracks_ = consumes<StreamsTrack>(InputTag(label, branchLostTracks));
0116     if (useMCTruth_) {
0117       const auto& inputTagSelecttion = iConfig.getParameter<InputTag>("InputTagSelection");
0118       const auto& inputTagReconstructable = iConfig.getParameter<InputTag>("InputTagReconstructable");
0119       edGetTokenSelection_ = consumes<StubAssociation>(inputTagSelecttion);
0120       edGetTokenReconstructable_ = consumes<StubAssociation>(inputTagReconstructable);
0121     }
0122     // book ES products
0123     esGetTokenSetup_ = esConsumes<Setup, SetupRcd, Transition::BeginRun>();
0124     esGetTokenDataFormats_ = esConsumes<DataFormats, DataFormatsRcd, Transition::BeginRun>();
0125     esGetTokenChannelAssignment_ = esConsumes<ChannelAssignment, ChannelAssignmentRcd, Transition::BeginRun>();
0126     // log config
0127     log_.setf(ios::fixed, ios::floatfield);
0128     log_.precision(4);
0129   }
0130
0131   void AnalyzerDRin::beginRun(const Run& iEvent, const EventSetup& iSetup) {
0132     // helper class to store configurations
0133     setup_ = &iSetup.getData(esGetTokenSetup_);
0134     // helper class to extract structured data from tt::Frames
0135     dataFormats_ = &iSetup.getData(esGetTokenDataFormats_);
0136     // helper class to assign tracklet track to channel
0137     channelAssignment_ = &iSetup.getData(esGetTokenChannelAssignment_);
0138     // book histograms
0139     Service<TFileService> fs;
0140     TFileDirectory dir;
0141     dir = fs->mkdir("DRin");
0142     prof_ = dir.make<TProfile>("Counts", ";", 10, 0.5, 10.5);
0143     prof_->GetXaxis()->SetBinLabel(1, "Stubs");
0144     prof_->GetXaxis()->SetBinLabel(2, "Tracks");
0145     prof_->GetXaxis()->SetBinLabel(3, "Lost Tracks");
0146     prof_->GetXaxis()->SetBinLabel(4, "Matched Tracks");
0147     prof_->GetXaxis()->SetBinLabel(5, "All Tracks");
0148     prof_->GetXaxis()->SetBinLabel(6, "Found TPs");
0149     prof_->GetXaxis()->SetBinLabel(7, "Found selected TPs");
0150     prof_->GetXaxis()->SetBinLabel(8, "Lost TPs");
0151     prof_->GetXaxis()->SetBinLabel(9, "All TPs");
0152     prof_->GetXaxis()->SetBinLabel(10, "Perfect TPs");
0153     // channel occupancy
0154     constexpr int maxOcc = 180;
0155     const int numChannels = channelAssignment_->numNodesDR();
0156     hisChannel_ = dir.make<TH1F>("His Channel Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0157     profChannel_ = dir.make<TProfile>("Prof Channel Occupancy", ";", numChannels, -.5, numChannels - .5);
0158   }
0159
0160   void AnalyzerDRin::analyze(const Event& iEvent, const EventSetup& iSetup) {
0161     // read in ht products
0162     Handle<StreamsStub> handleAcceptedStubs;
0163     iEvent.getByToken<StreamsStub>(edGetTokenAcceptedStubs_, handleAcceptedStubs);
0164     const StreamsStub& acceptedStubs = *handleAcceptedStubs;
0165     Handle<StreamsTrack> handleAcceptedTracks;
0166     iEvent.getByToken<StreamsTrack>(edGetTokenAcceptedTracks_, handleAcceptedTracks);
0167     const StreamsTrack& acceptedTracks = *handleAcceptedTracks;
0168     Handle<StreamsStub> handleLostStubs;
0169     iEvent.getByToken<StreamsStub>(edGetTokenLostStubs_, handleLostStubs);
0170     const StreamsStub& lostStubs = *handleLostStubs;
0171     Handle<StreamsTrack> handleLostTracks;
0172     iEvent.getByToken<StreamsTrack>(edGetTokenLostTracks_, handleLostTracks);
0173     const StreamsTrack& lostTracks = *handleLostTracks;
0174     // read in MCTruth
0175     const StubAssociation* selection = nullptr;
0176     const StubAssociation* reconstructable = nullptr;
0177     if (useMCTruth_) {
0178       Handle<StubAssociation> handleSelection;
0179       iEvent.getByToken<StubAssociation>(edGetTokenSelection_, handleSelection);
0180       selection = handleSelection.product();
0181       prof_->Fill(9, selection->numTPs());
0182       Handle<StubAssociation> handleReconstructable;
0183       iEvent.getByToken<StubAssociation>(edGetTokenReconstructable_, handleReconstructable);
0184       reconstructable = handleReconstructable.product();
0185     }
0186     // analyze ht products and associate found tracks with reconstrucable TrackingParticles
0187     set<TPPtr> tpPtrs;
0188     set<TPPtr> tpPtrsSelection;
0189     set<TPPtr> tpPtrsPerfect;
0190     set<TPPtr> tpPtrsLost;
0191     int allMatched(0);
0192     int allTracks(0);
0193     for (int region = 0; region < setup_->numRegions(); region++) {
0194       const int offset = region * channelAssignment_->numNodesDR();
0195       int nStubs(0);
0196       int nTracks(0);
0197       int nLost(0);
0198       for (int channel = 0; channel < channelAssignment_->numNodesDR(); channel++) {
0199         vector<vector<TTStubRef>> tracks;
0200         formTracks(acceptedTracks, acceptedStubs, tracks, offset + channel);
0201         vector<vector<TTStubRef>> lost;
0202         formTracks(lostTracks, lostStubs, lost, offset + channel);
0203         nTracks += tracks.size();
0204         nStubs += accumulate(tracks.begin(), tracks.end(), 0, [](int sum, const vector<TTStubRef>& track) {
0205           return sum + (int)track.size();
0206         });
0207         nLost += lost.size();
0208         allTracks += tracks.size();
0209         if (!useMCTruth_)
0210           continue;
0211         int tmp(0);
0212         associate(tracks, selection, tpPtrsSelection, tmp);
0213         associate(tracks, selection, tpPtrsPerfect, tmp, true);
0214         associate(lost, selection, tpPtrsLost, tmp);
0215         associate(tracks, reconstructable, tpPtrs, allMatched);
0216         const StreamTrack& stream = acceptedTracks[offset + channel];
0217         const auto end =
0218             find_if(stream.rbegin(), stream.rend(), [](const FrameTrack& frame) { return frame.first.isNonnull(); });
0219         const int size = distance(stream.begin(), end.base()) - 1;
0220         hisChannel_->Fill(size);
0221         profChannel_->Fill(channel, size);
0222       }
0223       prof_->Fill(1, nStubs);
0224       prof_->Fill(2, nTracks);
0225       prof_->Fill(3, nLost);
0226     }
0227     vector<TPPtr> recovered;
0228     recovered.reserve(tpPtrsLost.size());
0229     set_intersection(tpPtrsLost.begin(), tpPtrsLost.end(), tpPtrs.begin(), tpPtrs.end(), back_inserter(recovered));
0230     for (const TPPtr& tpPtr : recovered)
0231       tpPtrsLost.erase(tpPtr);
0232     prof_->Fill(4, allMatched);
0233     prof_->Fill(5, allTracks);
0234     prof_->Fill(6, tpPtrs.size());
0235     prof_->Fill(7, tpPtrsSelection.size());
0236     prof_->Fill(8, tpPtrsLost.size());
0237     prof_->Fill(10, tpPtrsPerfect.size());
0238     nEvents_++;
0239   }
0240
0241   void AnalyzerDRin::endJob() {
0242     if (nEvents_ == 0)
0243       return;
0244     // printout SF summary
0245     const double totalTPs = prof_->GetBinContent(9);
0246     const double numStubs = prof_->GetBinContent(1);
0247     const double numTracks = prof_->GetBinContent(2);
0248     const double numTracksLost = prof_->GetBinContent(3);
0249     const double totalTracks = prof_->GetBinContent(5);
0250     const double numTracksMatched = prof_->GetBinContent(4);
0251     const double numTPsAll = prof_->GetBinContent(6);
0252     const double numTPsEff = prof_->GetBinContent(7);
0253     const double numTPsLost = prof_->GetBinContent(8);
0254     const double numTPsEffPerfect = prof_->GetBinContent(10);
0255     const double errStubs = prof_->GetBinError(1);
0256     const double errTracks = prof_->GetBinError(2);
0257     const double errTracksLost = prof_->GetBinError(3);
0258     const double fracFake = (totalTracks - numTracksMatched) / totalTracks;
0259     const double fracDup = (numTracksMatched - numTPsAll) / totalTracks;
0260     const double eff = numTPsEff / totalTPs;
0261     const double errEff = sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0262     const double effLoss = numTPsLost / totalTPs;
0263     const double errEffLoss = sqrt(effLoss * (1. - effLoss) / totalTPs / nEvents_);
0264     const double effPerfect = numTPsEffPerfect / totalTPs;
0265     const double errEffPerfect = sqrt(effPerfect * (1. - effPerfect) / totalTPs / nEvents_);
0266     const vector<double> nums = {numStubs, numTracks, numTracksLost};
0267     const vector<double> errs = {errStubs, errTracks, errTracksLost};
0268     const int wNums = ceil(log10(*max_element(nums.begin(), nums.end()))) + 5;
0269     const int wErrs = ceil(log10(*max_element(errs.begin(), errs.end()))) + 5;
0270     log_ << "                        DRin  SUMMARY                        " << endl;
0271     log_ << "number of stubs       per TFP = " << setw(wNums) << numStubs << " +- " << setw(wErrs) << errStubs << endl;
0272     log_ << "number of tracks      per TFP = " << setw(wNums) << numTracks << " +- " << setw(wErrs) << errTracks
0273          << endl;
0274     log_ << "number of lost tracks per TFP = " << setw(wNums) << numTracksLost << " +- " << setw(wErrs) << errTracksLost
0275          << endl;
0276     log_ << "  current tracking efficiency = " << setw(wNums) << effPerfect << " +- " << setw(wErrs) << errEffPerfect
0277          << endl;
0278     log_ << "  max     tracking efficiency = " << setw(wNums) << eff << " +- " << setw(wErrs) << errEff << endl;
0279     log_ << "     lost tracking efficiency = " << setw(wNums) << effLoss << " +- " << setw(wErrs) << errEffLoss << endl;
0280     log_ << "                    fake rate = " << setw(wNums) << fracFake << endl;
0281     log_ << "               duplicate rate = " << setw(wNums) << fracDup << endl;
0282     log_ << "=============================================================";
0283     LogPrint("L1Trigger/TrackerTFP") << log_.str();
0284   }
0285
0286   //
0287   void AnalyzerDRin::formTracks(const StreamsTrack& streamsTrack,
0288                                 const StreamsStub& streamsStubs,
0289                                 vector<vector<TTStubRef>>& tracks,
0290                                 int channel) const {
0291     const int offset = channel * setup_->numLayers();
0292     const StreamTrack& streamTrack = streamsTrack[channel];
0293     const int numTracks = accumulate(streamTrack.begin(), streamTrack.end(), 0, [](int sum, const FrameTrack& frame) {
0294       return sum + (frame.first.isNonnull() ? 1 : 0);
0295     });
0296     tracks.reserve(numTracks);
0297     for (int frame = 0; frame < (int)streamTrack.size(); frame++) {
0298       const FrameTrack& frameTrack = streamTrack[frame];
0299       if (frameTrack.first.isNull())
0300         continue;
0301       vector<TTStubRef> ttStubRefs;
0302       ttStubRefs.reserve(setup_->numLayers());
0303       for (int layer = 0; layer < setup_->numLayers(); layer++) {
0304         const FrameStub& stub = streamsStubs[offset + layer][frame];
0305         if (stub.first.isNonnull())
0306           ttStubRefs.push_back(stub.first);
0307       }
0308       tracks.push_back(ttStubRefs);
0309     }
0310   }
0311
0312   //
0313   void AnalyzerDRin::associate(const vector<vector<TTStubRef>>& tracks,
0314                                const StubAssociation* ass,
0315                                set<TPPtr>& tps,
0316                                int& sum,
0317                                bool perfect) const {
0318     for (const vector<TTStubRef>& ttStubRefs : tracks) {
0319       const vector<TPPtr>& tpPtrs = perfect ? ass->associateFinal(ttStubRefs) : ass->associate(ttStubRefs);
0320       if (tpPtrs.empty())
0321         continue;
0322       sum++;
0323       copy(tpPtrs.begin(), tpPtrs.end(), inserter(tps, tps.begin()));
0324     }
0325   }
0326
0327 }  // namespace trklet
0328
0329 DEFINE_FWK_MODULE(trklet::AnalyzerDRin);