Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:32:45 

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