TrackerDTC/test/Analyzer.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 "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0014 #include "SimTracker/TrackTriggerAssociation/interface/TTClusterAssociationMap.h"
0015 #include "SimTracker/Common/interface/TrackingParticleSelector.h"
0016 #include "DataFormats/DetId/interface/DetId.h"
0017 #include "DataFormats/Common/interface/Ptr.h"
0018 #include "DataFormats/Common/interface/Handle.h"
0019 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
0020 #include "DataFormats/L1TrackTrigger/interface/TTDTC.h"
0021 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0022 #include "DataFormats/GeometrySurface/interface/Plane.h"
0023 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0024
0025 #include "L1Trigger/TrackTrigger/interface/Setup.h"
0026 #include "L1Trigger/TrackerDTC/interface/LayerEncoding.h"
0027
0028 #include <TProfile.h>
0029 #include <TProfile2D.h>
0030 #include <TH1F.h>
0031 #include <TH2F.h>
0032 #include <TEfficiency.h>
0033
0034 #include <vector>
0035 #include <map>
0036 #include <utility>
0037 #include <set>
0038 #include <algorithm>
0039 #include <cmath>
0040 #include <numeric>
0041 #include <array>
0042 #include <initializer_list>
0043 #include <sstream>
0044
0045 using namespace std;
0046 using namespace edm;
0047 using namespace tt;
0048
0049 namespace trackerDTC {
0050
0051   // mc truth types
0052   typedef TTClusterAssociationMap<Ref_Phase2TrackerDigi_> TTClusterAssMap;
0053   typedef edm::Ptr<TrackingParticle> TPPtr;
0054   // stub resolution plots helper
0055   enum Resolution { R, Phi, Z, NumResolution };
0056   constexpr initializer_list<Resolution> AllResolution = {R, Phi, Z};
0057   constexpr auto NameResolution = {"R", "Phi", "Z"};
0058   inline string name(Resolution r) { return string(*(NameResolution.begin() + r)); }
0059   // max tracking efficiency plots helper
0060   enum Efficiency { Phi0, Pt, InvPt, D0, Z0, Eta, NumEfficiency };
0061   constexpr initializer_list<Efficiency> AllEfficiency = {Phi0, Pt, InvPt, D0, Z0, Eta};
0062   constexpr auto NameEfficiency = {"Phi0", "Pt", "InvPt", "D0", "Z0", "Eta"};
0063   inline string name(Efficiency e) { return string(*(NameEfficiency.begin() + e)); }
0064
0065   /*! \class  trackerDTC::Analyzer
0066    *  \brief  Class to analyze hardware like structured TTStub Collection used by Track Trigger emulators, runs DTC stub emulation, plots performance & stub occupancy
0067    *  \author Thomas Schuh
0068    *  \date   2020, Apr
0069    */
0070   class Analyzer : public one::EDAnalyzer<one::WatchRuns, one::SharedResources> {
0071   public:
0072     Analyzer(const ParameterSet& iConfig);
0073     void beginJob() override {}
0074     void beginRun(const Run& iEvent, const EventSetup& iSetup) override;
0075     void analyze(const Event& iEvent, const EventSetup& iSetup) override;
0076     void endRun(const Run& iEvent, const EventSetup& iSetup) override {}
0077     void endJob() override;
0078
0079   private:
0080     // configuring track particle selector
0081     void configTPSelector();
0082     // book histograms
0083     void bookHistograms();
0084     // associate TPPtr with TTStubRef
0085     void assoc(const Handle<TTStubDetSetVec>&, const Handle<TTClusterAssMap>&, map<TPPtr, set<TTStubRef>>&);
0086     // organize reconstrucable TrackingParticles used for efficiency measurements
0087     void convert(const map<TPPtr, set<TTStubRef>>&, map<TTStubRef, set<TPPtr>>&);
0088     // checks if a stub selection is considered reconstructable
0089     bool reconstructable(const set<TTStubRef>& ttStubRefs) const;
0090     // checks if TrackingParticle is selected for efficiency measurements
0091     bool select(const TrackingParticle& tp) const;
0092     // fills kinematic tp histograms
0093     void fill(const TPPtr& tpPtr, const vector<TH1F*> th1fs) const;
0094     // analyze DTC products and find still reconstrucable TrackingParticles
0095     void analyzeStubs(const TTDTC*, const TTDTC*, const map<TTStubRef, set<TPPtr>>&, map<TPPtr, set<TTStubRef>>&);
0096     // fill stub related histograms
0097     void analyzeStream(const StreamStub& stream, int region, int channel, int& sum, TH2F* th2f);
0098     // returns layerId [1-6, 11-15] of stub
0099     int layerId(const TTStubRef& ttStubRef) const;
0100     // analyze survived TPs
0101     void analyzeTPs(const map<TPPtr, set<TTStubRef>>& mapTPsStubs);
0102     // prints out MC summary
0103     void endJobMC();
0104     // prints out DTC summary
0105     void endJobDTC();
0106
0107     // ED input token of DTC stubs
0108     EDGetTokenT<TTDTC> getTokenTTDTCAccepted_;
0109     // ED input token of lost DTC stubs
0110     EDGetTokenT<TTDTC> getTokenTTDTCLost_;
0111     // ED input token of TT stubs
0112     EDGetTokenT<TTStubDetSetVec> getTokenTTStubDetSetVec_;
0113     // ED input token of TTClsuter
0114     EDGetTokenT<TTClusterDetSetVec> getTokenTTClusterDetSetVec_;
0115     // ED input token of TTCluster to TPPtr association
0116     EDGetTokenT<TTClusterAssMap> getTokenTTClusterAssMap_;
0117     // Setup token
0118     ESGetToken<Setup, SetupRcd> esGetToken_;
0119     // stores, calculates and provides run-time constants
0120     const Setup* setup_ = nullptr;
0121     // selector to partly select TPs for efficiency measurements
0122     TrackingParticleSelector tpSelector_;
0123     //
0124     TrackingParticleSelector tpSelectorLoose_;
0125     // enables analyze of TPs
0126     bool useMCTruth_;
0127     // specifies used TT algorithm
0128     bool hybrid_;
0129     //
0130     int nEvents_ = 0;
0131
0132     // Histograms
0133
0134     TProfile* profMC_;
0135     TProfile* profDTC_;
0136     TProfile* profChannel_;
0137     TH1F* hisChannel_;
0138     TH2F* hisRZStubs_;
0139     TH2F* hisRZStubsLost_;
0140     TH2F* hisRZStubsEff_;
0141     vector<TH1F*> hisResolution_;
0142     vector<TProfile2D*> profResolution_;
0143     vector<TH1F*> hisEff_;
0144     vector<TH1F*> hisEffMC_;
0145     vector<TEfficiency*> eff_;
0146
0147     // printout
0148     stringstream log_;
0149   };
0150
0151   Analyzer::Analyzer(const ParameterSet& iConfig)
0152       : useMCTruth_(iConfig.getParameter<bool>("UseMCTruth")), hybrid_(iConfig.getParameter<bool>("UseHybrid")) {
0153     usesResource("TFileService");
0154     // book in- and output ED products
0155     const auto& inputTagAccepted = iConfig.getParameter<InputTag>("InputTagAccepted");
0156     const auto& inputTagLost = iConfig.getParameter<InputTag>("InputTagLost");
0157     getTokenTTDTCAccepted_ = consumes<TTDTC>(inputTagAccepted);
0158     getTokenTTDTCLost_ = consumes<TTDTC>(inputTagLost);
0159     if (useMCTruth_) {
0160       const auto& inputTagTTStubDetSetVec = iConfig.getParameter<InputTag>("InputTagTTStubDetSetVec");
0161       const auto& inputTagTTClusterDetSetVec = iConfig.getParameter<InputTag>("InputTagTTClusterDetSetVec");
0162       const auto& inputTagTTClusterAssMap = iConfig.getParameter<InputTag>("InputTagTTClusterAssMap");
0163       getTokenTTStubDetSetVec_ = consumes<TTStubDetSetVec>(inputTagTTStubDetSetVec);
0164       getTokenTTClusterDetSetVec_ = consumes<TTClusterDetSetVec>(inputTagTTClusterDetSetVec);
0165       getTokenTTClusterAssMap_ = consumes<TTClusterAssMap>(inputTagTTClusterAssMap);
0166     }
0167     // book ES product
0168     esGetToken_ = esConsumes<Setup, SetupRcd, Transition::BeginRun>();
0169     // log config
0170     log_.setf(ios::fixed, ios::floatfield);
0171     log_.precision(4);
0172   }
0173
0174   void Analyzer::beginRun(const Run& iEvent, const EventSetup& iSetup) {
0175     // helper class to store configurations
0176     setup_ = &iSetup.getData(esGetToken_);
0177     // configuring track particle selector
0178     configTPSelector();
0179     // book histograms
0180     bookHistograms();
0181   }
0182
0183   void Analyzer::analyze(const Event& iEvent, const EventSetup& iSetup) {
0184     // read in TrackingParticle
0185     map<TTStubRef, set<TPPtr>> mapAllStubsTPs;
0186     if (useMCTruth_) {
0187       Handle<TTStubDetSetVec> handleTTStubDetSetVec;
0188       iEvent.getByToken<TTStubDetSetVec>(getTokenTTStubDetSetVec_, handleTTStubDetSetVec);
0189       Handle<TTClusterAssMap> handleTTClusterAssMap;
0190       iEvent.getByToken<TTClusterAssMap>(getTokenTTClusterAssMap_, handleTTClusterAssMap);
0191       // associate TPPtr with TTStubRef
0192       map<TPPtr, set<TTStubRef>> mapAllTPsAllStubs;
0193       assoc(handleTTStubDetSetVec, handleTTClusterAssMap, mapAllTPsAllStubs);
0194       // organize reconstrucable TrackingParticles used for efficiency measurements
0195       convert(mapAllTPsAllStubs, mapAllStubsTPs);
0196       Handle<TTClusterDetSetVec> handleTTClusterDetSetVec;
0197       iEvent.getByToken<TTClusterDetSetVec>(getTokenTTClusterDetSetVec_, handleTTClusterDetSetVec);
0198       int nCluster(0);
0199       for (const auto& detSet : *handleTTClusterDetSetVec)
0200         nCluster += detSet.size();
0201       profMC_->Fill(6, nCluster / (double)setup_->numRegions());
0202     }
0203     // read in dtc products
0204     Handle<TTDTC> handleTTDTCAccepted;
0205     iEvent.getByToken<TTDTC>(getTokenTTDTCAccepted_, handleTTDTCAccepted);
0206     Handle<TTDTC> handleTTDTCLost;
0207     iEvent.getByToken<TTDTC>(getTokenTTDTCLost_, handleTTDTCLost);
0208     map<TPPtr, set<TTStubRef>> mapTPsTTStubs;
0209     // analyze DTC products and find still reconstrucable TrackingParticles
0210     analyzeStubs(handleTTDTCAccepted.product(), handleTTDTCLost.product(), mapAllStubsTPs, mapTPsTTStubs);
0211     // analyze survived TPs
0212     analyzeTPs(mapTPsTTStubs);
0213     nEvents_++;
0214   }
0215
0216   void Analyzer::endJob() {
0217     if (nEvents_ == 0)
0218       return;
0219     // create r-z stub fraction plot
0220     TH2F th2f("", ";;", 400, -300, 300., 400, 0., 120.);
0221     th2f.Add(hisRZStubsLost_);
0222     th2f.Add(hisRZStubs_);
0223     hisRZStubsEff_->Add(hisRZStubsLost_);
0224     hisRZStubsEff_->Divide(&th2f);
0225     // create efficieny plots
0226     if (useMCTruth_) {
0227       for (Efficiency e : AllEfficiency) {
0228         eff_[e]->SetPassedHistogram(*hisEff_[e], "f");
0229         eff_[e]->SetTotalHistogram(*hisEffMC_[e], "f");
0230       }
0231     }
0232     log_ << "'Lost' below refers to truncation losses" << endl;
0233     // printout MC summary
0234     endJobMC();
0235     // printout DTC summary
0236     endJobDTC();
0237     log_ << "=============================================================";
0238     LogPrint("L1Trigger/TrackerDTC") << log_.str();
0239   }
0240
0241   // associate TPPtr with TTStubRef
0242   void Analyzer::assoc(const Handle<TTStubDetSetVec>& handleTTStubDetSetVec,
0243                        const Handle<TTClusterAssMap>& handleTTClusterAssMap,
0244                        map<TPPtr, set<TTStubRef>>& mapTPsStubs) {
0245     int nStubs(0);
0246     int nStubsMatched(0);
0247     for (TTStubDetSetVec::const_iterator ttModule = handleTTStubDetSetVec->begin();
0248          ttModule != handleTTStubDetSetVec->end();
0249          ttModule++) {
0250       nStubs += ttModule->size();
0251       for (TTStubDetSet::const_iterator ttStub = ttModule->begin(); ttStub != ttModule->end(); ttStub++) {
0252         set<TPPtr> tpPtrs;
0253         for (unsigned int iClus = 0; iClus < 2; iClus++) {
0254           const vector<TPPtr>& assocPtrs = handleTTClusterAssMap->findTrackingParticlePtrs(ttStub->clusterRef(iClus));
0255           copy_if(assocPtrs.begin(), assocPtrs.end(), inserter(tpPtrs, tpPtrs.begin()), [](const TPPtr& tpPtr) {
0256             return tpPtr.isNonnull();
0257           });
0258         }
0259         for (const TPPtr& tpPtr : tpPtrs)
0260           mapTPsStubs[tpPtr].emplace(makeRefTo(handleTTStubDetSetVec, ttStub));
0261         if (!tpPtrs.empty())
0262           nStubsMatched++;
0263       }
0264     }
0265     profMC_->Fill(1, nStubs / (double)setup_->numRegions());
0266     profMC_->Fill(2, nStubsMatched / (double)setup_->numRegions());
0267   }
0268
0269   // organize reconstrucable TrackingParticles used for efficiency measurements
0270   void Analyzer::convert(const map<TPPtr, set<TTStubRef>>& mapTPsStubs, map<TTStubRef, set<TPPtr>>& mapStubsTPs) {
0271     int nTPsReco(0);
0272     int nTPsEff(0);
0273     for (const auto& mapTPStubs : mapTPsStubs) {
0274       if (!tpSelectorLoose_(*mapTPStubs.first) || !reconstructable(mapTPStubs.second))
0275         continue;
0276       nTPsReco++;
0277       const bool useForAlgEff = select(*mapTPStubs.first);
0278       if (useForAlgEff) {
0279         nTPsEff++;
0280         fill(mapTPStubs.first, hisEffMC_);
0281         for (const TTStubRef& ttStubRef : mapTPStubs.second)
0282           mapStubsTPs[ttStubRef].insert(mapTPStubs.first);
0283       }
0284     }
0285     profMC_->Fill(3, nTPsReco / (double)setup_->numRegions());
0286     profMC_->Fill(4, nTPsEff / (double)setup_->numRegions());
0287     profMC_->Fill(5, nTPsEff);
0288   }
0289
0290   // checks if a stub selection is considered reconstructable
0291   bool Analyzer::reconstructable(const set<TTStubRef>& ttStubRefs) const {
0292     const TrackerGeometry* trackerGeometry = setup_->trackerGeometry();
0293     const TrackerTopology* trackerTopology = setup_->trackerTopology();
0294     set<int> hitPattern;
0295     set<int> hitPatternPS;
0296     for (const TTStubRef& ttStubRef : ttStubRefs) {
0297       const DetId detId = ttStubRef->getDetId();
0298       const bool barrel = detId.subdetId() == StripSubdetector::TOB;
0299       const bool psModule = trackerGeometry->getDetectorType(detId) == TrackerGeometry::ModuleType::Ph2PSP;
0300       const int layerId = barrel ? trackerTopology->layer(detId) : trackerTopology->tidWheel(detId) + 10;
0301       hitPattern.insert(layerId);
0302       if (psModule)
0303         hitPatternPS.insert(layerId);
0304     }
0305     return (int)hitPattern.size() >= setup_->tpMinLayers() && (int)hitPatternPS.size() >= setup_->tpMinLayersPS();
0306   }
0307
0308   // checks if TrackingParticle is selected for efficiency measurements
0309   bool Analyzer::select(const TrackingParticle& tp) const {
0310     const bool selected = tpSelector_(tp);
0311     const double cot = sinh(tp.eta());
0312     const double s = sin(tp.phi());
0313     const double c = cos(tp.phi());
0314     const TrackingParticle::Point& v = tp.vertex();
0315     const double z0 = v.z() - (v.x() * c + v.y() * s) * cot;
0316     const double d0 = v.x() * s - v.y() * c;
0317     return selected && (fabs(d0) < setup_->tpMaxD0()) && (fabs(z0) < setup_->tpMaxVertZ());
0318   }
0319
0320   // fills kinematic tp histograms
0321   void Analyzer::fill(const TPPtr& tpPtr, const vector<TH1F*> th1fs) const {
0322     const double s = sin(tpPtr->phi());
0323     const double c = cos(tpPtr->phi());
0324     const TrackingParticle::Point& v = tpPtr->vertex();
0325     const vector<double> x = {tpPtr->phi(),
0326                               tpPtr->pt(),
0327                               tpPtr->charge() / tpPtr->pt(),
0328                               v.x() * s - v.y() * c,
0329                               v.z() - (v.x() * c + v.y() * s) * sinh(tpPtr->eta()),
0330                               tpPtr->eta()};
0331     for (Efficiency e : AllEfficiency)
0332       th1fs[e]->Fill(x[e]);
0333   }
0334
0335   // analyze DTC products and find still reconstrucable TrackingParticles
0336   void Analyzer::analyzeStubs(const TTDTC* accepted,
0337                               const TTDTC* lost,
0338                               const map<TTStubRef, set<TPPtr>>& mapStubsTPs,
0339                               map<TPPtr, set<TTStubRef>>& mapTPsStubs) {
0340     for (int region = 0; region < setup_->numRegions(); region++) {
0341       int nStubs(0);
0342       int nLost(0);
0343       for (int channel = 0; channel < setup_->numDTCsPerTFP(); channel++) {
0344         const StreamStub& stream = accepted->stream(region, channel);
0345         hisChannel_->Fill(stream.size());
0346         profChannel_->Fill(region * setup_->numDTCsPerTFP() + channel, stream.size());
0347         for (const FrameStub& frame : stream) {
0348           if (frame.first.isNull())
0349             continue;
0350           const auto it = mapStubsTPs.find(frame.first);
0351           if (it == mapStubsTPs.end())
0352             continue;
0353           for (const TPPtr& tp : it->second)
0354             mapTPsStubs[tp].insert(frame.first);
0355         }
0356         analyzeStream(stream, region, channel, nStubs, hisRZStubs_);
0357         analyzeStream(lost->stream(region, channel), region, channel, nLost, hisRZStubsLost_);
0358       }
0359       profDTC_->Fill(1, nStubs);
0360       profDTC_->Fill(2, nLost);
0361     }
0362   }
0363
0364   // fill stub related histograms
0365   void Analyzer::analyzeStream(const StreamStub& stream, int region, int channel, int& sum, TH2F* th2f) {
0366     for (const FrameStub& frame : stream) {
0367       if (frame.first.isNull())
0368         continue;
0369       sum++;
0370       const GlobalPoint& pos = setup_->stubPos(hybrid_, frame, region);
0371       const GlobalPoint& ttPos = setup_->stubPos(frame.first);
0372       const vector<double> resolutions = {
0373           ttPos.perp() - pos.perp(), deltaPhi(ttPos.phi() - pos.phi()), ttPos.z() - pos.z()};
0374       for (Resolution r : AllResolution) {
0375         hisResolution_[r]->Fill(resolutions[r]);
0376         profResolution_[r]->Fill(ttPos.z(), ttPos.perp(), abs(resolutions[r]));
0377       }
0378       th2f->Fill(ttPos.z(), ttPos.perp());
0379     }
0380   }
0381
0382   // returns layerId [1-6, 11-15] of stub
0383   int Analyzer::layerId(const TTStubRef& ttStubRef) const {
0384     const TrackerTopology* trackerTopology = setup_->trackerTopology();
0385     const DetId detId = ttStubRef->getDetId() + setup_->offsetDetIdDSV();
0386     const bool barrel = detId.subdetId() == StripSubdetector::TOB;
0387     return barrel ? trackerTopology->layer(detId) : trackerTopology->tidWheel(detId) + setup_->offsetLayerDisks();
0388   }
0389
0390   // analyze survived TPs
0391   void Analyzer::analyzeTPs(const map<TPPtr, set<TTStubRef>>& mapTPsStubs) {
0392     int nTPs(0);
0393     for (const auto& mapTPStubs : mapTPsStubs) {
0394       if (!reconstructable(mapTPStubs.second))
0395         continue;
0396       nTPs++;
0397       fill(mapTPStubs.first, hisEff_);
0398     }
0399     profDTC_->Fill(3, nTPs);
0400   }
0401
0402   // prints out MC summary
0403   void Analyzer::endJobMC() {
0404     const double numStubs = profMC_->GetBinContent(1);
0405     const double numStubsMatched = profMC_->GetBinContent(2);
0406     const double numTPsReco = profMC_->GetBinContent(3);
0407     const double numTPsEff = profMC_->GetBinContent(4);
0408     const double errStubs = profMC_->GetBinError(1);
0409     const double errStubsMatched = profMC_->GetBinError(2);
0410     const double errTPsReco = profMC_->GetBinError(3);
0411     const double errTPsEff = profMC_->GetBinError(4);
0412     const double numCluster = profMC_->GetBinContent(6);
0413     const double errCluster = profMC_->GetBinError(6);
0414     const vector<double> nums = {numStubs, numStubsMatched, numTPsReco, numTPsEff, numCluster};
0415     const vector<double> errs = {errStubs, errStubsMatched, errTPsReco, errTPsEff, errCluster};
0416     const int wNums = ceil(log10(*max_element(nums.begin(), nums.end()))) + 5;
0417     const int wErrs = ceil(log10(*max_element(errs.begin(), errs.end()))) + 5;
0418     log_ << "=============================================================" << endl;
0419     log_ << "                         MC  SUMMARY                         " << endl;
0420     log_ << "number of cluster       per TFP = " << setw(wNums) << numCluster << " +- " << setw(wErrs) << errCluster
0421          << endl;
0422     log_ << "number of stubs         per TFP = " << setw(wNums) << numStubs << " +- " << setw(wErrs) << errStubs
0423          << endl;
0424     log_ << "number of matched stubs per TFP = " << setw(wNums) << numStubsMatched << " +- " << setw(wErrs)
0425          << errStubsMatched << endl;
0426     log_ << "number of TPs           per TFP = " << setw(wNums) << numTPsReco << " +- " << setw(wErrs) << errTPsReco
0427          << endl;
0428     log_ << "number of TPs for eff   per TFP = " << setw(wNums) << numTPsEff << " +- " << setw(wErrs) << errTPsEff
0429          << endl;
0430   }
0431
0432   // prints out DTC summary
0433   void Analyzer::endJobDTC() {
0434     const double numStubs = profDTC_->GetBinContent(1);
0435     const double numStubsLost = profDTC_->GetBinContent(2);
0436     const double numTPs = profDTC_->GetBinContent(3);
0437     const double errStubs = profDTC_->GetBinError(1);
0438     const double errStubsLost = profDTC_->GetBinError(2);
0439     const double totalTPs = profMC_->GetBinContent(5);
0440     const double eff = numTPs / totalTPs;
0441     const double errEff = sqrt(eff * (1. - eff) / totalTPs / nEvents_);
0442     const vector<double> nums = {numStubs, numStubsLost};
0443     const vector<double> errs = {errStubs, errStubsLost};
0444     const int wNums = ceil(log10(*max_element(nums.begin(), nums.end()))) + 5;
0445     const int wErrs = ceil(log10(*max_element(errs.begin(), errs.end()))) + 5;
0446     log_ << "=============================================================" << endl;
0447     log_ << "                         DTC SUMMARY                         " << endl;
0448     log_ << "number of stubs      per TFP = " << setw(wNums) << numStubs << " +- " << setw(wErrs) << errStubs << endl;
0449     log_ << "number of lost stubs per TFP = " << setw(wNums) << numStubsLost << " +- " << setw(wErrs) << errStubsLost
0450          << endl;
0451     log_ << "     max tracking efficiency = " << setw(wNums) << eff << " +- " << setw(wErrs) << errEff << endl;
0452   }
0453
0454   // configuring track particle selector
0455   void Analyzer::configTPSelector() {
0456     const double ptMin = hybrid_ ? setup_->hybridMinPtStub() : setup_->minPt();
0457     constexpr double ptMax = 9999999999.;
0458     const double etaMax = setup_->tpMaxEta();
0459     const double tip = setup_->tpMaxVertR();
0460     const double lip = setup_->tpMaxVertZ();
0461     constexpr int minHit = 0;
0462     constexpr bool signalOnly = true;
0463     constexpr bool intimeOnly = true;
0464     constexpr bool chargedOnly = true;
0465     constexpr bool stableOnly = false;
0466     tpSelector_ = TrackingParticleSelector(
0467         ptMin, ptMax, -etaMax, etaMax, tip, lip, minHit, signalOnly, intimeOnly, chargedOnly, stableOnly);
0468     tpSelectorLoose_ =
0469         TrackingParticleSelector(ptMin, ptMax, -etaMax, etaMax, tip, lip, minHit, false, false, false, stableOnly);
0470   }
0471
0472   // book histograms
0473   void Analyzer::bookHistograms() {
0474     Service<TFileService> fs;
0475     TFileDirectory dir;
0476     // mc
0477     dir = fs->mkdir("MC");
0478     profMC_ = dir.make<TProfile>("Counts", ";", 6, 0.5, 6.5);
0479     profMC_->GetXaxis()->SetBinLabel(1, "Stubs");
0480     profMC_->GetXaxis()->SetBinLabel(2, "Matched Stubs");
0481     profMC_->GetXaxis()->SetBinLabel(3, "reco TPs");
0482     profMC_->GetXaxis()->SetBinLabel(4, "eff TPs");
0483     profMC_->GetXaxis()->SetBinLabel(5, "total eff TPs");
0484     profMC_->GetXaxis()->SetBinLabel(6, "Cluster");
0485     constexpr array<int, NumEfficiency> binsEff{{9 * 8, 10, 16, 10, 30, 24}};
0486     constexpr array<pair<double, double>, NumEfficiency> rangesEff{
0487         {{-M_PI, M_PI}, {0., 100.}, {-1. / 3., 1. / 3.}, {-5., 5.}, {-15., 15.}, {-2.4, 2.4}}};
0488     if (useMCTruth_) {
0489       hisEffMC_.reserve(NumEfficiency);
0490       for (Efficiency e : AllEfficiency)
0491         hisEffMC_.emplace_back(
0492             dir.make<TH1F>(("HisTP" + name(e)).c_str(), ";", binsEff[e], rangesEff[e].first, rangesEff[e].second));
0493     }
0494     // dtc
0495     dir = fs->mkdir("DTC");
0496     profDTC_ = dir.make<TProfile>("Counts", ";", 3, 0.5, 3.5);
0497     profDTC_->GetXaxis()->SetBinLabel(1, "Stubs");
0498     profDTC_->GetXaxis()->SetBinLabel(2, "Lost Stubs");
0499     profDTC_->GetXaxis()->SetBinLabel(3, "TPs");
0500     // channel occupancy
0501     constexpr int maxOcc = 180;
0502     const int numChannels = setup_->numDTCs() * setup_->numOverlappingRegions();
0503     hisChannel_ = dir.make<TH1F>("His Channel Occupancy", ";", maxOcc, -.5, maxOcc - .5);
0504     profChannel_ = dir.make<TProfile>("Prof Channel Occupancy", ";", numChannels, -.5, numChannels - .5);
0505     // max tracking efficiencies
0506     if (useMCTruth_) {
0507       dir = fs->mkdir("DTC/Effi");
0508       hisEff_.reserve(NumEfficiency);
0509       for (Efficiency e : AllEfficiency)
0510         hisEff_.emplace_back(
0511             dir.make<TH1F>(("HisTP" + name(e)).c_str(), ";", binsEff[e], rangesEff[e].first, rangesEff[e].second));
0512       eff_.reserve(NumEfficiency);
0513       for (Efficiency e : AllEfficiency)
0514         eff_.emplace_back(
0515             dir.make<TEfficiency>(("Eff" + name(e)).c_str(), ";", binsEff[e], rangesEff[e].first, rangesEff[e].second));
0516     }
0517     // lost stub fraction in r-z
0518     dir = fs->mkdir("DTC/Loss");
0519     constexpr int bins = 400;
0520     constexpr double maxZ = 300.;
0521     constexpr double maxR = 120.;
0522     hisRZStubs_ = dir.make<TH2F>("RZ Stubs", ";;", bins, -maxZ, maxZ, bins, 0., maxR);
0523     hisRZStubsLost_ = dir.make<TH2F>("RZ Stubs Lost", ";;", bins, -maxZ, maxZ, bins, 0., maxR);
0524     hisRZStubsEff_ = dir.make<TH2F>("RZ Stubs Eff", ";;", bins, -maxZ, maxZ, bins, 0., maxR);
0525     // stub parameter resolutions
0526     dir = fs->mkdir("DTC/Res");
0527     constexpr array<double, NumResolution> ranges{{.2, .0001, .5}};
0528     constexpr int binsHis = 100;
0529     hisResolution_.reserve(NumResolution);
0530     profResolution_.reserve(NumResolution);
0531     for (Resolution r : AllResolution) {
0532       hisResolution_.emplace_back(dir.make<TH1F>(("HisRes" + name(r)).c_str(), ";", binsHis, -ranges[r], ranges[r]));
0533       profResolution_.emplace_back(
0534           dir.make<TProfile2D>(("ProfRes" + name(r)).c_str(), ";;", bins, -maxZ, maxZ, bins, 0., maxR));
0535     }
0536   }
0537
0538 }  // namespace trackerDTC
0539
0540 DEFINE_FWK_MODULE(trackerDTC::Analyzer);