Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-12-04 23:57:27 

0001 // -*- C++ -*-
0002 //
0003 // Package:    TkAlTools/JetHTAnalyzer
0004 // Class:      JetHTAnalyzer
0005 //
0006 /**\class JetHTAnalyzer JetHTAnalyzer.cc Alignment/OfflineValidation/plugins/JetHTAnalyzer.cc
0007 
0008  Description: [one line class summary]
0009 
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Marco Musich
0015 //         Created:  Fri, 29 Mar 2019 14:54:59 GMT
0016 //
0017 //      Updated by:  Jussi Viinikainen
0018 //
0019 
0020 // system include files
0021 #include <algorithm>  // std::sort
0022 #include <memory>
0023 #include <string>
0024 #include <vector>  // std::vector
0025 
0026 // CMSSW includes
0027 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0028 #include "DataFormats/TrackReco/interface/Track.h"
0029 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0030 #include "DataFormats/VertexReco/interface/Vertex.h"
0031 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0032 #include "FWCore/Framework/interface/Event.h"
0033 #include "FWCore/Framework/interface/EventSetup.h"
0034 #include "FWCore/Framework/interface/Frameworkfwd.h"
0035 #include "FWCore/Framework/interface/MakerMacros.h"
0036 #include "FWCore/Framework/interface/MakerMacros.h"
0037 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0038 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0039 #include "FWCore/ServiceRegistry/interface/Service.h"
0040 #include "FWCore/Utilities/interface/InputTag.h"
0041 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0042 #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h"
0043 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0044 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0045 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0046 #include "Alignment/OfflineValidation/interface/SmartSelectionMonitor.h"
0047 #include "DataFormats/Common/interface/TriggerResults.h"  // Classes needed to print trigger results
0048 #include "FWCore/Common/interface/TriggerNames.h"
0049 
0050 // ROOT includes
0051 #include "TRandom.h"
0052 #include "TTree.h"
0053 #include "TString.h"
0054 #include "TMath.h"
0055 
0056 //
0057 // class declaration
0058 //
0059 using reco::TrackCollection;
0060 
0061 class JetHTAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0062 public:
0063   explicit JetHTAnalyzer(const edm::ParameterSet&);
0064   ~JetHTAnalyzer() override = default;
0065 
0066   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0067   static bool mysorter(reco::Track i, reco::Track j) { return (i.pt() > j.pt()); }
0068 
0069 private:
0070   void beginJob() override;
0071   void analyze(const edm::Event&, const edm::EventSetup&) override;
0072   void endJob() override;
0073 
0074   // ----------member data ---------------------------
0075   const edm::InputTag pvsTag_;
0076   const edm::EDGetTokenT<reco::VertexCollection> pvsToken_;
0077 
0078   const edm::InputTag tracksTag_;
0079   const edm::EDGetTokenT<reco::TrackCollection> tracksToken_;
0080 
0081   const edm::InputTag triggerTag_;
0082   const edm::EDGetTokenT<edm::TriggerResults> triggerToken_;
0083 
0084   const int printTriggerTable_;
0085   const double minVtxNdf_;
0086   const double minVtxWgt_;
0087 
0088   const std::vector<double> profilePtBorders_;
0089   const std::vector<int> iovList_;
0090 
0091   // output histograms
0092   edm::Service<TFileService> outfile_;
0093   TH1F* h_ntrks;
0094   TH1F* h_probePt;
0095   TH1F* h_probeEta;
0096   TH1F* h_probePhi;
0097   TH1F* h_probeDxy;
0098   TH1F* h_probeDz;
0099   TH1F* h_probeDxyErr;
0100   TH1F* h_probeDzErr;
0101 
0102   SmartSelectionMonitor mon;
0103 
0104   // for the conversions
0105   static constexpr double cmToum = 10000;
0106 };
0107 
0108 //
0109 // Constructor
0110 //
0111 JetHTAnalyzer::JetHTAnalyzer(const edm::ParameterSet& iConfig)
0112     : pvsTag_(iConfig.getParameter<edm::InputTag>("vtxCollection")),
0113       pvsToken_(consumes<reco::VertexCollection>(pvsTag_)),
0114       tracksTag_(iConfig.getParameter<edm::InputTag>("trackCollection")),
0115       tracksToken_(consumes<reco::TrackCollection>(tracksTag_)),
0116       triggerTag_(iConfig.getParameter<edm::InputTag>("triggerResults")),
0117       triggerToken_(consumes<edm::TriggerResults>(triggerTag_)),
0118       printTriggerTable_(iConfig.getUntrackedParameter<int>("printTriggerTable")),
0119       minVtxNdf_(iConfig.getUntrackedParameter<double>("minVertexNdf")),
0120       minVtxWgt_(iConfig.getUntrackedParameter<double>("minVertexMeanWeight")),
0121       profilePtBorders_(iConfig.getUntrackedParameter<std::vector<double>>("profilePtBorders")),
0122       iovList_(iConfig.getUntrackedParameter<std::vector<int>>("iovList")) {
0123   // Specify that TFileService is used by the class
0124   usesResource(TFileService::kSharedResource);
0125 }
0126 
0127 //
0128 // member functions
0129 //
0130 
0131 // ------------ method called for each event  ------------
0132 void JetHTAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0133   using namespace edm;
0134 
0135   const auto& vertices = iEvent.get(pvsToken_);
0136   const reco::VertexCollection& pvtx = vertices;
0137 
0138   edm::Handle<reco::TrackCollection> tracks = iEvent.getHandle(tracksToken_);
0139 
0140   // Find the IOV of the current event so that we can tag a histogram with the IOV
0141   const int runNumber = iEvent.id().run();
0142   std::string iovString = "iovNotFound";
0143 
0144   // Find the current IOV from the IOV list
0145   if (runNumber >= iovList_.at(0)) {  // If run number is smaller than the first item in the list, it is not in any IOV
0146     for (std::vector<int>::size_type i = 1; i < iovList_.size(); i++) {
0147       if (iovList_.at(i) > runNumber) {
0148         iovString = Form("iov%d-%d", iovList_.at(i - 1), iovList_.at(i) - 1);
0149         break;
0150       }
0151     }
0152   }
0153 
0154   // Print the triggers to console
0155   if (printTriggerTable_) {
0156     const auto& triggerResults = iEvent.get(triggerToken_);
0157 
0158     const edm::TriggerNames& triggerNames = iEvent.triggerNames(triggerResults);
0159     for (unsigned i = 0; i < triggerNames.size(); i++) {
0160       const std::string& hltName = triggerNames.triggerName(i);
0161       bool decision = triggerResults.accept(triggerNames.triggerIndex(hltName));
0162       std::cout << hltName << " " << decision << std::endl;
0163     }
0164   }
0165 
0166   int counter = 0;
0167   for (reco::VertexCollection::const_iterator pvIt = pvtx.begin(); pvIt != pvtx.end(); pvIt++) {
0168     const reco::Vertex& iPV = *pvIt;
0169     counter++;
0170 
0171     if (iPV.isFake())
0172       continue;
0173     reco::Vertex::trackRef_iterator trki;
0174 
0175     const math::XYZPoint pos_(iPV.x(), iPV.y(), iPV.z());
0176 
0177     // vertex selection as in bs code
0178     if (iPV.ndof() < minVtxNdf_ || (iPV.ndof() + 3.) / iPV.tracksSize() < 2 * minVtxWgt_)
0179       continue;
0180 
0181     reco::TrackCollection allTracks;
0182     for (trki = iPV.tracks_begin(); trki != iPV.tracks_end(); ++trki) {
0183       if (trki->isNonnull()) {
0184         reco::TrackRef trk_now(tracks, (*trki).key());
0185         allTracks.push_back(*trk_now);
0186       }
0187     }
0188 
0189     // order with decreasing pt
0190     std::sort(allTracks.begin(), allTracks.end(), mysorter);
0191     uint ntrks = allTracks.size();
0192     h_ntrks->Fill(ntrks);
0193 
0194     for (uint tracksIt = 0; tracksIt < ntrks; tracksIt++) {
0195       auto tk = allTracks.at(tracksIt);
0196 
0197       double dxyRes = tk.dxy(pos_) * cmToum;
0198       double dzRes = tk.dz(pos_) * cmToum;
0199 
0200       double dxy_err = tk.dxyError() * cmToum;
0201       double dz_err = tk.dzError() * cmToum;
0202 
0203       float trackphi = tk.phi();
0204       float tracketa = tk.eta();
0205       float trackpt = tk.pt();
0206 
0207       h_probePt->Fill(trackpt);
0208       h_probeEta->Fill(tracketa);
0209       h_probePhi->Fill(trackphi);
0210 
0211       h_probeDxy->Fill(dxyRes);
0212       h_probeDz->Fill(dzRes);
0213       h_probeDxyErr->Fill(dxy_err);
0214       h_probeDzErr->Fill(dz_err);
0215 
0216       mon.fillHisto("dxy", "all", dxyRes, 1.);
0217       mon.fillHisto("dz", "all", dzRes, 1.);
0218       mon.fillHisto("dxyerr", "all", dxy_err, 1.);
0219       mon.fillHisto("dzerr", "all", dz_err, 1.);
0220 
0221       mon.fillProfile("dxyErrVsPt", "all", trackpt, dxy_err, 1.);
0222       mon.fillProfile("dzErrVsPt", "all", trackpt, dz_err, 1.);
0223 
0224       mon.fillProfile("dxyErrVsPhi", "all", trackphi, dxy_err, 1.);
0225       mon.fillProfile("dzErrVsPhi", "all", trackphi, dz_err, 1.);
0226 
0227       mon.fillProfile("dxyErrVsEta", "all", tracketa, dxy_err, 1.);
0228       mon.fillProfile("dzErrVsEta", "all", tracketa, dz_err, 1.);
0229 
0230       // Integrated pT bins
0231       for (std::vector<double>::size_type i = 0; i < profilePtBorders_.size(); i++) {
0232         if (trackpt < profilePtBorders_.at(i))
0233           break;
0234         mon.fillProfile("dxyErrVsPtWide", "all", i, dxy_err, 1.);
0235         mon.fillProfile("dzErrVsPtWide", "all", i, dz_err, 1.);
0236       }
0237 
0238       // Fill IOV specific histograms
0239       mon.fillHisto("dxy", iovString, dxyRes, 1.);
0240       mon.fillHisto("dz", iovString, dzRes, 1.);
0241       mon.fillHisto("dxyerr", iovString, dxy_err, 1.);
0242       mon.fillHisto("dzerr", iovString, dz_err, 1.);
0243 
0244       mon.fillProfile("dxyErrVsPt", iovString, trackpt, dxy_err, 1.);
0245       mon.fillProfile("dzErrVsPt", iovString, trackpt, dz_err, 1.);
0246 
0247       mon.fillProfile("dxyErrVsPhi", iovString, trackphi, dxy_err, 1.);
0248       mon.fillProfile("dzErrVsPhi", iovString, trackphi, dz_err, 1.);
0249 
0250       mon.fillProfile("dxyErrVsEta", iovString, tracketa, dxy_err, 1.);
0251       mon.fillProfile("dzErrVsEta", iovString, tracketa, dz_err, 1.);
0252 
0253       // Integrated pT bins
0254       for (std::vector<double>::size_type i = 0; i < profilePtBorders_.size(); i++) {
0255         if (trackpt < profilePtBorders_.at(i))
0256           break;
0257         mon.fillProfile("dxyErrVsPtWide", iovString, i, dxy_err, 1.);
0258         mon.fillProfile("dzErrVsPtWide", iovString, i, dz_err, 1.);
0259       }
0260 
0261       if (std::abs(tracketa) < 1.) {
0262         mon.fillHisto("dxy", "central", dxyRes, 1.);
0263         mon.fillHisto("dz", "central", dzRes, 1.);
0264         mon.fillHisto("dxyerr", "central", dxy_err, 1.);
0265         mon.fillHisto("dzerr", "central", dz_err, 1.);
0266 
0267         mon.fillProfile("dxyErrVsPt", "central", trackpt, dxy_err, 1.);
0268         mon.fillProfile("dzErrVsPt", "central", trackpt, dz_err, 1.);
0269 
0270         mon.fillProfile("dxyErrVsPhi", "central", trackphi, dxy_err, 1.);
0271         mon.fillProfile("dzErrVsPhi", "central", trackphi, dz_err, 1.);
0272 
0273         // Integrated pT bins
0274         for (std::vector<double>::size_type i = 0; i < profilePtBorders_.size(); i++) {
0275           if (trackpt < profilePtBorders_.at(i))
0276             break;
0277           mon.fillProfile("dxyErrVsPtWide", "central", i, dxy_err, 1.);
0278           mon.fillProfile("dzErrVsPtWide", "central", i, dz_err, 1.);
0279         }
0280       }
0281 
0282     }  // loop on tracks in vertex
0283   }    // loop on vertices
0284 
0285   mon.fillHisto("nvtx", "all", counter, 1.);
0286 }
0287 
0288 // ------------ method called once each job just before starting event loop  ------------
0289 void JetHTAnalyzer::beginJob() {
0290   h_ntrks = outfile_->make<TH1F>("h_ntrks", "n. trks;n. of tracks/vertex;n. vertices", 100, 0, 100);
0291   h_probePt = outfile_->make<TH1F>("h_probePt", "p_{T} of probe track;track p_{T} (GeV); tracks", 100, 0., 500.);
0292   h_probeEta = outfile_->make<TH1F>("h_probeEta", "#eta of the probe track;track #eta;tracks", 54, -2.8, 2.8);
0293   h_probePhi = outfile_->make<TH1F>("h_probePhi", "#phi of probe track;track #phi (rad);tracks", 100, -3.15, 3.15);
0294 
0295   h_probeDxy =
0296       outfile_->make<TH1F>("h_probeDxy", "d_{xy}(PV) of the probe track;track d_{xy}(PV);tracks", 200, -100, 100);
0297   h_probeDz = outfile_->make<TH1F>("h_probeDz", "d_{z}(PV) of the probe track;track d_{z}(PV);tracks", 200, -100, 100);
0298   h_probeDxyErr = outfile_->make<TH1F>(
0299       "h_probeDxyErr", "error on d_{xy}(PV) of the probe track;track error on d_{xy}(PV);tracks", 100, 0., 100);
0300   h_probeDzErr = outfile_->make<TH1F>(
0301       "h_probeDzErr", "error on d_{z}(PV)  of the probe track;track error on d_{z}(PV);tracks", 100, 0., 100);
0302 
0303   mon.addHistogram(new TH1F("nvtx", ";Vertices;Events", 50, 0, 50));
0304   mon.addHistogram(new TH1F("dxy", ";d_{xy};tracks", 100, -100, 100));
0305   mon.addHistogram(new TH1F("dz", ";d_{z};tracks", 100, -100, 100));
0306   mon.addHistogram(new TH1F("dxyerr", ";d_{xy} error;tracks", 100, 0., 200));
0307   mon.addHistogram(new TH1F("dzerr", ";d_{z} error;tracks", 100, 0., 200));
0308   mon.addHistogram(new TProfile("dxyErrVsPt", ";track p_{T};d_{xy} error", 100, 0., 200, 0., 100.));
0309   mon.addHistogram(new TProfile("dzErrVsPt", ";track p_{T};d_{z} error", 100, 0., 200, 0., 100.));
0310   mon.addHistogram(
0311       new TProfile("dxyErrVsPhi", ";track #varphi;d_{xy} error", 100, -TMath::Pi(), TMath::Pi(), 0., 100.));
0312   mon.addHistogram(new TProfile("dzErrVsPhi", ";track #varphi;d_{z} error", 100, -TMath::Pi(), TMath::Pi(), 0., 100.));
0313   mon.addHistogram(new TProfile("dxyErrVsEta", ";track #eta;d_{xy} error", 100, -2.5, 2.5, 0., 100.));
0314   mon.addHistogram(new TProfile("dzErrVsEta", ";track #eta;d_{z} error", 100, -2.5, 2.5, 0., 100.));
0315 
0316   // Variable size histogram depending on the given pT bin borders
0317   int nBins = profilePtBorders_.size();
0318   mon.addHistogram(
0319       new TProfile("dxyErrVsPtWide", ";track p_{T} wide bin;d_{xy} error", nBins, -0.5, nBins - 0.5, 0.0, 100.0));
0320   mon.addHistogram(
0321       new TProfile("dzErrVsPtWide", ";track p_{T} wide bin;d_{z} error", nBins, -0.5, nBins - 0.5, 0.0, 100.0));
0322 }
0323 
0324 // ------------ method called once each job just after ending the event loop  ------------
0325 void JetHTAnalyzer::endJob() { mon.Write(); }
0326 
0327 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0328 void JetHTAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0329   edm::ParameterSetDescription desc;
0330   desc.setComment("JetHT validation analyzer plugin.");
0331   desc.add<edm::InputTag>("vtxCollection", edm::InputTag("offlinePrimaryVerticesFromRefittedTrks"));
0332   desc.add<edm::InputTag>("triggerResults", edm::InputTag("TriggerResults", "", "HLT"));
0333   desc.add<edm::InputTag>("trackCollection", edm::InputTag("TrackRefitter"));
0334   desc.add<int>("printTriggerTable", false);
0335   desc.add<double>("minVertexNdf", 10.);
0336   desc.add<double>("minVertexMeanWeight", 0.5);
0337   desc.add<std::vector<double>>("profilePtBorders", {3, 5, 10, 20, 50, 100});
0338   desc.add<std::vector<double>>("iovList", {0, 500000});
0339   descriptions.addWithDefaultLabel(desc);
0340 }
0341 
0342 //define this as a plug-in
0343 DEFINE_FWK_MODULE(JetHTAnalyzer);
[ 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!