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::AnalyzerCTB
0033    *  \brief  Class to analyze hardware like structured TTStub Collection generated by Clean Track Builder
0034    *  \author Thomas Schuh
0035    *  \date   2020, April
0036    */
0037   class AnalyzerCTB : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0038   public:
0039     AnalyzerCTB(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::StreamsTrack& streamsTrack,
0049                     const tt::StreamsStub& streamsStubs,
0050                     std::vector<std::vector<TTStubRef>>& tracks,
0051                     int channel) const;
0052     //
0053     void associate(const std::vector<std::vector<TTStubRef>>& tracks,
0054                    const tt::StubAssociation* ass,
0055                    std::set<TPPtr>& tps,
0056                    int& sum) const;
0057     // ED input token of stubs
0058     edm::EDGetTokenT<tt::StreamsStub> edGetTokenStubs_;
0059     // ED input token of tracks
0060     edm::EDGetTokenT<tt::StreamsTrack> edGetTokenTracks_;
0061     // ED input token of TTStubRef to TPPtr association for tracking efficiency
0062     edm::EDGetTokenT<tt::StubAssociation> edGetTokenSelection_;
0063     // ED input token of TTStubRef to recontructable TPPtr association
0064     edm::EDGetTokenT<tt::StubAssociation> edGetTokenReconstructable_;
0065     // Setup token
0066     edm::ESGetToken<tt::Setup, tt::SetupRcd> esGetTokenSetup_;
0067     // DataFormats token
0068     edm::ESGetToken<DataFormats, DataFormatsRcd> esGetTokenDataFormats_;
0069     // stores, calculates and provides run-time constants
0070     const tt::Setup* setup_ = nullptr;
0071     // helper class to extract structured data from tt::Frames
0072     const DataFormats* dataFormats_ = nullptr;
0073     // enables analyze of TPs
0074     bool useMCTruth_;
0075     //
0076     int nEvents_ = 0;
0077 
0078     // Histograms
0079 
0080     TProfile* prof_;
0081     TProfile* profChan_;
0082     TProfile* profStubs_;
0083     TProfile* profTracks_;
0084     TH1F* hisChan_;
0085     TH1F* hisStubs_;
0086     TH1F* hisTracks_;
0087     TH1F* hisLayers_;
0088     TH1F* hisNumLayers_;
0089     TProfile* profNumLayers_;
0090     TH1F* hisEffZT_;
0091     TH1F* hisEffZTTotal_;
0092     TEfficiency* effZT_;
0093     TH1F* hisEffInv2R_;
0094     TH1F* hisEffInv2RTotal_;
0095     TEfficiency* effInv2R_;
0096 
0097     // printout
0098     std::stringstream log_;
0099   };
0100 
0101   AnalyzerCTB::AnalyzerCTB(const edm::ParameterSet& iConfig) : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")) {
0102     usesResource("TFileService");
0103     // book in- and output ED products
0104     const std::string& label = iConfig.getParameter<std::string>("OutputLabelCTB");
0105     const std::string& branchStubs = iConfig.getParameter<std::string>("BranchStubs");
0106     const std::string& branchTracks = iConfig.getParameter<std::string>("BranchTracks");
0107     edGetTokenStubs_ = consumes<tt::StreamsStub>(edm::InputTag(label, branchStubs));
0108     edGetTokenTracks_ = consumes<tt::StreamsTrack>(edm::InputTag(label, branchTracks));
0109     if (useMCTruth_) {
0110       const auto& inputTagSelecttion = iConfig.getParameter<edm::InputTag>("InputTagSelection");
0111       const auto& inputTagReconstructable = iConfig.getParameter<edm::InputTag>("InputTagReconstructable");
0112       edGetTokenSelection_ = consumes<tt::StubAssociation>(inputTagSelecttion);
0113       edGetTokenReconstructable_ = consumes<tt::StubAssociation>(inputTagReconstructable);
0114     }
0115     // book ES products
0116     esGetTokenSetup_ = esConsumes<edm::Transition::BeginRun>();
0117     esGetTokenDataFormats_ = esConsumes<edm::Transition::BeginRun>();
0118     // log config
0119     log_.setf(std::ios::fixed, std::ios::floatfield);
0120     log_.precision(4);
0121   }
0122 
0123   void AnalyzerCTB::beginRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) {
0124     // helper class to store configurations
0125     setup_ = &iSetup.getData(esGetTokenSetup_);
0126     // helper class to extract structured data from tt::Frames
0127     dataFormats_ = &iSetup.getData(esGetTokenDataFormats_);
0128     // book histograms
0129     edm::Service<TFileService> fs;
0130     TFileDirectory dir;
0131     dir = fs->mkdir("CTB");
0132     prof_ = dir.make<TProfile>("Counts", ";", 9, 0.5, 9.5);
0133     prof_->GetXaxis()->SetBinLabel(1, "Stubs");
0134     prof_->GetXaxis()->SetBinLabel(2, "Tracks");
0135     prof_->GetXaxis()->SetBinLabel(4, "Matched Tracks");
0136     prof_->GetXaxis()->SetBinLabel(5, "All Tracks");
0137     prof_->GetXaxis()->SetBinLabel(6, "Found TPs");
0138     prof_->GetXaxis()->SetBinLabel(7, "Found selected TPs");
0139     prof_->GetXaxis()->SetBinLabel(9, "All TPs");
0140     // channel occupancy
0141     constexpr int maxOcc = 180;
0142     const int numChannelsTracks = dataFormats_->numChannel(Process::ctb);
0143     const int numChannelsStubs = numChannelsTracks * setup_->numLayers();
0144     hisChan_ = dir.make<TH1F>("His Channel Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0145     profChan_ = dir.make<TProfile>("Prof Channel Occupancy", ";", numChannelsTracks, -.5, numChannelsTracks - .5);
0146     // stub occupancy
0147     hisStubs_ = dir.make<TH1F>("His Channel Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0148     profStubs_ = dir.make<TProfile>("Prof Channel Occupancy", ";", numChannelsStubs, -.5, numChannelsStubs - .5);
0149     // track occupancy
0150     hisTracks_ = dir.make<TH1F>("His Track Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0151     profTracks_ = dir.make<TProfile>("Prof Track Occupancy", ";", numChannelsTracks, -.5, numChannelsTracks - .5);
0152     // layers
0153     hisLayers_ = dir.make<TH1F>("HisLayers", ";", 8, 0, 8);
0154     hisNumLayers_ = dir.make<TH1F>("HisNumLayers", ";", 9, 0, 9);
0155     profNumLayers_ = dir.make<TProfile>("Prof NumLayers", ";", 32, 0, 2.4);
0156     // Efficiencies
0157     const double rangeZT = dataFormats_->format(Variable::zT, Process::dr).range();
0158     const int zTBins = setup_->gpNumBinsZT();
0159     hisEffZTTotal_ = dir.make<TH1F>("HisTPZTTotal", ";", zTBins, -rangeZT / 2, rangeZT / 2);
0160     hisEffZT_ = dir.make<TH1F>("HisTPZT", ";", zTBins, -rangeZT / 2, rangeZT / 2);
0161     effZT_ = dir.make<TEfficiency>("EffZT", ";", zTBins, -rangeZT / 2, rangeZT / 2);
0162     const double rangeInv2R = dataFormats_->format(Variable::inv2R, Process::dr).range();
0163     const int inv2RBins = (setup_->htNumBinsInv2R() + 2) * 2;
0164     hisEffInv2R_ = dir.make<TH1F>("HisTPInv2R", ";", inv2RBins, -rangeInv2R / 2., rangeInv2R / 2.);
0165     hisEffInv2RTotal_ = dir.make<TH1F>("HisTPInv2RTotal", ";", inv2RBins, -rangeInv2R / 2., rangeInv2R / 2.);
0166     effInv2R_ = dir.make<TEfficiency>("EffInv2R", ";", inv2RBins, -rangeInv2R / 2., rangeInv2R / 2.);
0167   }
0168 
0169   void AnalyzerCTB::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0170     auto fill = [this](const TPPtr& tpPtr, TH1F* hisZT, TH1F* hisInv2R) {
0171       const double tpPhi0 = tpPtr->phi();
0172       const double tpCot = sinh(tpPtr->eta());
0173       const math::XYZPointD& v = tpPtr->vertex();
0174       const double tpZ0 = v.z() - tpCot * (v.x() * cos(tpPhi0) + v.y() * sin(tpPhi0));
0175       const double tpZT = tpZ0 + tpCot * setup_->chosenRofZ();
0176       const double tpInv2R = tpPtr->charge() / tpPtr->pt() * setup_->invPtToDphi();
0177       hisZT->Fill(tpZT);
0178       hisInv2R->Fill(tpInv2R);
0179     };
0180     // read in ht products
0181     edm::Handle<tt::StreamsStub> handleStubs;
0182     iEvent.getByToken<tt::StreamsStub>(edGetTokenStubs_, handleStubs);
0183     const tt::StreamsStub& acceptedStubs = *handleStubs;
0184     edm::Handle<tt::StreamsTrack> handleTracks;
0185     iEvent.getByToken<tt::StreamsTrack>(edGetTokenTracks_, handleTracks);
0186     const tt::StreamsTrack& acceptedTracks = *handleTracks;
0187     // read in MCTruth
0188     const tt::StubAssociation* selection = nullptr;
0189     const tt::StubAssociation* reconstructable = nullptr;
0190     if (useMCTruth_) {
0191       edm::Handle<tt::StubAssociation> handleSelection;
0192       iEvent.getByToken<tt::StubAssociation>(edGetTokenSelection_, handleSelection);
0193       selection = handleSelection.product();
0194       prof_->Fill(9, selection->numTPs());
0195       edm::Handle<tt::StubAssociation> handleReconstructable;
0196       iEvent.getByToken<tt::StubAssociation>(edGetTokenReconstructable_, handleReconstructable);
0197       reconstructable = handleReconstructable.product();
0198       for (const auto& p : selection->getTrackingParticleToTTStubsMap())
0199         fill(p.first, hisEffZTTotal_, hisEffInv2RTotal_);
0200     }
0201     // analyze ht products and associate found tracks with reconstrucable TrackingParticles
0202     std::set<TPPtr> tpPtrs;
0203     std::set<TPPtr> tpPtrsSelection;
0204     int allMatched(0);
0205     int allTracks(0);
0206     for (int region = 0; region < setup_->numRegions(); region++) {
0207       const int offsetTrack = region * dataFormats_->numChannel(Process::ctb);
0208       int nStubs(0);
0209       int nTracks(0);
0210       for (int channel = 0; channel < dataFormats_->numChannel(Process::ctb); channel++) {
0211         const int indexTrack = offsetTrack + channel;
0212         const int size = acceptedTracks[indexTrack].size();
0213         hisChan_->Fill(size);
0214         profChan_->Fill(channel, size);
0215         const int offsetStub = indexTrack * setup_->numLayers();
0216         for (int layer = 0; layer < setup_->numLayers(); layer++) {
0217           const tt::StreamStub& stream = acceptedStubs[offsetStub + layer];
0218           const int nStubs = std::accumulate(stream.begin(), stream.end(), 0, [](int sum, const tt::FrameStub& frame) {
0219             return sum + (frame.first.isNonnull() ? 1 : 0);
0220           });
0221           hisStubs_->Fill(nStubs);
0222           profStubs_->Fill(channel * setup_->numLayers() + layer, nStubs);
0223         }
0224         std::vector<std::vector<TTStubRef>> tracks;
0225         formTracks(acceptedTracks, acceptedStubs, tracks, indexTrack);
0226         hisTracks_->Fill(tracks.size());
0227         profTracks_->Fill(channel, tracks.size());
0228         nTracks += tracks.size();
0229         nStubs += std::accumulate(tracks.begin(), tracks.end(), 0, [](int sum, const std::vector<TTStubRef>& track) {
0230           return sum + track.size();
0231         });
0232         allTracks += tracks.size();
0233         if (!useMCTruth_)
0234           continue;
0235         int tmp(0);
0236         associate(tracks, selection, tpPtrsSelection, tmp);
0237         associate(tracks, reconstructable, tpPtrs, allMatched);
0238       }
0239       prof_->Fill(1, nStubs);
0240       prof_->Fill(2, nTracks);
0241     }
0242     for (const TPPtr& tpPtr : tpPtrsSelection)
0243       fill(tpPtr, hisEffZT_, hisEffInv2R_);
0244     prof_->Fill(4, allMatched);
0245     prof_->Fill(5, allTracks);
0246     prof_->Fill(6, tpPtrs.size());
0247     prof_->Fill(7, tpPtrsSelection.size());
0248     nEvents_++;
0249   }
0250 
0251   void AnalyzerCTB::endJob() {
0252     if (nEvents_ == 0)
0253       return;
0254     // effi
0255     effZT_->SetPassedHistogram(*hisEffZT_, "f");
0256     effZT_->SetTotalHistogram(*hisEffZTTotal_, "f");
0257     effInv2R_->SetPassedHistogram(*hisEffInv2R_, "f");
0258     effInv2R_->SetTotalHistogram(*hisEffInv2RTotal_, "f");
0259     // printout TB summary
0260     const double totalTPs = prof_->GetBinContent(9);
0261     const double numStubs = prof_->GetBinContent(1);
0262     const double numTracks = prof_->GetBinContent(2);
0263     const double totalTracks = prof_->GetBinContent(5);
0264     const double numTracksMatched = prof_->GetBinContent(4);
0265     const double numTPsAll = prof_->GetBinContent(6);
0266     const double numTPsEff = prof_->GetBinContent(7);
0267     const double errStubs = prof_->GetBinError(1);
0268     const double errTracks = prof_->GetBinError(2);
0269     const double fracFake = (totalTracks - numTracksMatched) / totalTracks;
0270     const double fracDup = (numTracksMatched - numTPsAll) / totalTracks;
0271     const double eff = numTPsEff / totalTPs;
0272     const double errEff = sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0273     const std::vector<double> nums = {numStubs, numTracks};
0274     const std::vector<double> errs = {errStubs, errTracks};
0275     const int wNums = std::ceil(std::log10(*std::max_element(nums.begin(), nums.end()))) + 5;
0276     const int wErrs = std::ceil(std::log10(*std::max_element(errs.begin(), errs.end()))) + 5;
0277     log_ << "                         CTB SUMMARY                         " << std::endl;
0278     log_ << "number of stubs       per TFP = " << std::setw(wNums) << numStubs << " +- " << std::setw(wErrs) << errStubs
0279          << std::endl;
0280     log_ << "number of tracks      per TFP = " << std::setw(wNums) << numTracks << " +- " << std::setw(wErrs)
0281          << errTracks << std::endl;
0282     log_ << "     max  tracking efficiency = " << std::setw(wNums) << eff << " +- " << std::setw(wErrs) << errEff
0283          << std::endl;
0284     log_ << "                    fake rate = " << std::setw(wNums) << fracFake << std::endl;
0285     log_ << "               duplicate rate = " << std::setw(wNums) << fracDup << std::endl;
0286     log_ << "=============================================================";
0287     edm::LogPrint(moduleDescription().moduleName()) << log_.str();
0288   }
0289 
0290   //
0291   void AnalyzerCTB::formTracks(const tt::StreamsTrack& streamsTrack,
0292                                const tt::StreamsStub& streamsStubs,
0293                                std::vector<std::vector<TTStubRef>>& tracks,
0294                                int channel) const {
0295     const int offset = channel * setup_->numLayers();
0296     const tt::StreamTrack& streamTrack = streamsTrack[channel];
0297     const int numTracks =
0298         std::accumulate(streamTrack.begin(), streamTrack.end(), 0, [](int sum, const tt::FrameTrack& frame) {
0299           return sum + (frame.first.isNonnull() ? 1 : 0);
0300         });
0301     tracks.reserve(numTracks);
0302     for (int frame = 0; frame < static_cast<int>(streamTrack.size()); frame++) {
0303       const tt::FrameTrack& frameTrack = streamTrack[frame];
0304       if (frameTrack.first.isNull())
0305         continue;
0306       const auto end = std::find_if(std::next(streamTrack.begin(), frame + 1),
0307                                     streamTrack.end(),
0308                                     [](const tt::FrameTrack& frame) { return frame.first.isNonnull(); });
0309       const int size = std::distance(std::next(streamTrack.begin(), frame), end);
0310       int numStubs(0);
0311       for (int layer = 0; layer < setup_->numLayers(); layer++) {
0312         const tt::StreamStub& stream = streamsStubs[offset + layer];
0313         numStubs += std::accumulate(
0314             stream.begin() + frame, stream.begin() + frame + size, 0, [](int sum, const tt::FrameStub& frame) {
0315               return sum + (frame.first.isNonnull() ? 1 : 0);
0316             });
0317       }
0318       std::vector<TTStubRef> stubs;
0319       stubs.reserve(numStubs);
0320       int numLayers(0);
0321       for (int layer = 0; layer < setup_->numLayers(); layer++) {
0322         bool any(false);
0323         for (int f = frame; f < frame + size; f++) {
0324           const tt::FrameStub& stub = streamsStubs[offset + layer][f];
0325           if (stub.first.isNonnull()) {
0326             any = true;
0327             stubs.push_back(stub.first);
0328           }
0329         }
0330         if (any) {
0331           hisLayers_->Fill(layer);
0332           numLayers++;
0333         }
0334       }
0335       const double cot = TrackCTB(frameTrack, dataFormats_).zT() / setup_->chosenRofZ();
0336       hisNumLayers_->Fill(numLayers);
0337       profNumLayers_->Fill(abs(sinh(cot)), numLayers);
0338       tracks.push_back(stubs);
0339     }
0340   }
0341 
0342   //
0343   void AnalyzerCTB::associate(const std::vector<std::vector<TTStubRef>>& tracks,
0344                               const tt::StubAssociation* ass,
0345                               std::set<TPPtr>& tps,
0346                               int& sum) const {
0347     for (const std::vector<TTStubRef>& ttStubRefs : tracks) {
0348       const std::vector<TPPtr>& tpPtrs = ass->associate(ttStubRefs);
0349       if (tpPtrs.empty())
0350         continue;
0351       sum++;
0352       std::copy(tpPtrs.begin(), tpPtrs.end(), std::inserter(tps, tps.begin()));
0353     }
0354   }
0355 
0356 }  // namespace trackerTFP
0357 
0358 DEFINE_FWK_MODULE(trackerTFP::AnalyzerCTB);
[ 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!