PatExamples/plugins/PatVertexAnalyzer.cc

0001 #include <iostream>
0002 #include <cmath>
0003 #include <vector>
0004 #include <string>
0005
0006 #include <TH1.h>
0007 #include <TProfile.h>
0008
0009 #include "FWCore/Framework/interface/Frameworkfwd.h"
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0013 #include "FWCore/Utilities/interface/InputTag.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/ServiceRegistry/interface/Service.h"
0016
0017 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0018
0019 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0020 #include "DataFormats/TrackReco/interface/Track.h"
0021 #include "DataFormats/VertexReco/interface/Vertex.h"
0022 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0023 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0024 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0025
0026 class PatVertexAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0027 public:
0028   /// constructor and destructor
0029   PatVertexAnalyzer(const edm::ParameterSet &params);
0030   ~PatVertexAnalyzer() override;
0031
0032   // virtual methods called from base class EDAnalyzer
0033   void beginJob() override;
0034   void analyze(const edm::Event &event, const edm::EventSetup &es) override;
0035
0036 private:
0037   // configuration parameters
0038   edm::EDGetTokenT<reco::VertexCollection> srcToken_;
0039   edm::EDGetTokenT<reco::GenParticleCollection> genParticlesToken_;
0040
0041   TH1 *nVertices_, *nTracks_;
0042   TH1 *x_, *y_, *z_;
0043   TH1 *xErr_, *yErr_, *zErr_;
0044   TH1 *xDelta_, *yDelta_, *zDelta_;
0045   TH1 *xPull_, *yPull_, *zPull_;
0046 };
0047
0048 PatVertexAnalyzer::PatVertexAnalyzer(const edm::ParameterSet &params)
0049     : srcToken_(consumes<reco::VertexCollection>(params.getParameter<edm::InputTag>("src"))),
0050       genParticlesToken_(consumes<reco::GenParticleCollection>(params.getParameter<edm::InputTag>("mc"))) {
0051   usesResource(TFileService::kSharedResource);
0052 }
0053
0054 PatVertexAnalyzer::~PatVertexAnalyzer() {}
0055
0056 void PatVertexAnalyzer::beginJob() {
0057   // retrieve handle to auxiliary service
0058   //  used for storing histograms into ROOT file
0059   edm::Service<TFileService> fs;
0060
0061   nVertices_ = fs->make<TH1F>("nVertices", "number of reconstructed primary vertices", 50, 0, 50);
0062   nTracks_ = fs->make<TH1F>("nTracks", "number of tracks at primary vertex", 100, 0, 300);
0063   x_ = fs->make<TH1F>("pvX", "primary vertex x", 100, -0.1, 0.1);
0064   y_ = fs->make<TH1F>("pvY", "primary vertex y", 100, -0.1, 0.1);
0065   z_ = fs->make<TH1F>("pvZ", "primary vertex z", 100, -30, 30);
0066   xErr_ = fs->make<TH1F>("pvErrorX", "primary vertex x error", 100, 0, 0.005);
0067   yErr_ = fs->make<TH1F>("pvErrorY", "primary vertex y error", 100, 0, 0.005);
0068   zErr_ = fs->make<TH1F>("pvErrorZ", "primary vertex z error", 100, 0, 0.01);
0069   xDelta_ = fs->make<TH1F>("pvDeltaX", "x shift wrt simulated vertex", 100, -0.01, 0.01);
0070   yDelta_ = fs->make<TH1F>("pvDeltaY", "y shift wrt simulated vertex", 100, -0.01, 0.01);
0071   zDelta_ = fs->make<TH1F>("pvDeltaZ", "z shift wrt simulated vertex", 100, -0.02, 0.02);
0072   xPull_ = fs->make<TH1F>("pvPullX", "primary vertex x pull", 100, -5, 5);
0073   yPull_ = fs->make<TH1F>("pvPullY", "primary vertex y pull", 100, -5, 5);
0074   zPull_ = fs->make<TH1F>("pvPullZ", "primary vertex z pull", 100, -5, 5);
0075 }
0076
0077 void PatVertexAnalyzer::analyze(const edm::Event &event, const edm::EventSetup &es) {
0078   // handle to the primary vertex collection
0079   edm::Handle<reco::VertexCollection> pvHandle;
0080   event.getByToken(srcToken_, pvHandle);
0081
0082   // handle to the generator particles (i.e. the MC truth)
0083   edm::Handle<reco::GenParticleCollection> genParticlesHandle;
0084   event.getByToken(genParticlesToken_, genParticlesHandle);
0085
0086   // extract the position of the simulated vertex
0087   math::XYZPoint simPV = (*genParticlesHandle)[2].vertex();
0088
0089   // the number of reconstructed primary vertices
0090   nVertices_->Fill(pvHandle->size());
0091
0092   // if we have at least one, use the first (highest pt^2 sum)
0093   if (!pvHandle->empty()) {
0094     const reco::Vertex &pv = (*pvHandle)[0];
0095
0096     nTracks_->Fill(pv.tracksSize());
0097
0098     x_->Fill(pv.x());
0099     y_->Fill(pv.y());
0100     z_->Fill(pv.z());
0101
0102     xErr_->Fill(pv.xError());
0103     yErr_->Fill(pv.yError());
0104     zErr_->Fill(pv.zError());
0105
0106     xDelta_->Fill(pv.x() - simPV.X());
0107     yDelta_->Fill(pv.y() - simPV.Y());
0108     zDelta_->Fill(pv.z() - simPV.Z());
0109
0110     xPull_->Fill((pv.x() - simPV.X()) / pv.xError());
0111     yPull_->Fill((pv.y() - simPV.Y()) / pv.yError());
0112     zPull_->Fill((pv.z() - simPV.Z()) / pv.zError());
0113
0114     // we could access the tracks using the
0115     // pv.tracks_begin() ... pv.tracks_end() iterators
0116   }
0117 }
0118
0119 #include "FWCore/Framework/interface/MakerMacros.h"
0120
0121 DEFINE_FWK_MODULE(PatVertexAnalyzer);