File indexing completed on 2022-12-04 23:57:27
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021 #include <algorithm> // std::sort
0022 #include <memory>
0023 #include <string>
0024 #include <vector> // std::vector
0025
0026
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 #include "TRandom.h"
0052 #include "TTree.h"
0053 #include "TString.h"
0054 #include "TMath.h"
0055
0056
0057
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
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
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
0105 static constexpr double cmToum = 10000;
0106 };
0107
0108
0109
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
0124 usesResource(TFileService::kSharedResource);
0125 }
0126
0127
0128
0129
0130
0131
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
0141 const int runNumber = iEvent.id().run();
0142 std::string iovString = "iovNotFound";
0143
0144
0145 if (runNumber >= iovList_.at(0)) {
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
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
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
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
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
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
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
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 }
0283 }
0284
0285 mon.fillHisto("nvtx", "all", counter, 1.);
0286 }
0287
0288
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
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
0325 void JetHTAnalyzer::endJob() { mon.Write(); }
0326
0327
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
0343 DEFINE_FWK_MODULE(JetHTAnalyzer);