RecoVertex/src/MCvsRecoVerticesAnalyzer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    MCvsRecoVerticesAnalyzer
0004 // Class:      MCvsRecoVerticesAnalyzer
0005 //
0006 /**\class MCvsRecoVerticesAnalyzer MCvsRecoVerticesAnalyzer.cc TrackingPFG/PileUp/src/MCvsRecoVerticesAnalyzer.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Andrea Venturi
0015 //         Created:  Thu Dec 16 16:32:56 CEST 2010
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021 #include <numeric>
0022
0023 #include <vector>
0024
0025 // user include files
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0028
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/Run.h"
0031 #include "FWCore/Framework/interface/MakerMacros.h"
0032
0033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0034
0035 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0036
0037 #include "FWCore/ServiceRegistry/interface/Service.h"
0038 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0039
0040 #include "FWCore/Utilities/interface/InputTag.h"
0041
0042 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0043 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0044 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0045
0046 #include "DataFormats/VertexReco/interface/Vertex.h"
0047 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0048
0049 #include "TH1F.h"
0050 #include "TH2F.h"
0051 #include "TProfile.h"
0052
0053 //
0054 // class decleration
0055 //
0056
0057 class MCvsRecoVerticesAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0058 public:
0059   explicit MCvsRecoVerticesAnalyzer(const edm::ParameterSet&);
0060   ~MCvsRecoVerticesAnalyzer() override;
0061
0062 private:
0063   void analyze(const edm::Event&, const edm::EventSetup&) override;
0064
0065   // ----------member data ---------------------------
0066
0067   const bool m_useweight;
0068   const bool m_useVisibleVertices;
0069   const edm::ParameterSet m_histoParameters;
0070   edm::EDGetTokenT<double> m_doubleToken;
0071   edm::EDGetTokenT<std::vector<PileupSummaryInfo> > m_vecPileupSummaryInfoToken;
0072   edm::EDGetTokenT<reco::VertexCollection> m_recoVertexCollectionToken;
0073   edm::EDGetTokenT<edm::HepMCProduct> m_hepMCProductToken;
0074
0075   TH2F* m_hrecovsmcnvtx2d;
0076   TProfile* m_hrecovsmcnvtxprof;
0077   TProfile* m_hrecovsmcnvtxweightedprof;
0078   TH2F* m_hrecovsmclumi2d;
0079   TProfile* m_hrecovsmclumiprof;
0080   TProfile* m_hrecovsmclumiweightedprof;
0081   TH1F* m_hdeltazfirst;
0082   TH1F* m_hdeltazclose;
0083   TH1F* m_hclosestvtx;
0084
0085   TH2F* m_hdeltazfirstvsnpu;
0086   TH2F* m_hdeltazclosevsnpu;
0087   TH2F* m_hclosestvtxvsnpu;
0088 };
0089
0090 //
0091 // constants, enums and typedefs
0092 //
0093
0094 //
0095 // static data member definitions
0096 //
0097
0098 //
0099 // constructors and destructor
0100 //
0101 MCvsRecoVerticesAnalyzer::MCvsRecoVerticesAnalyzer(const edm::ParameterSet& iConfig)
0102     : m_useweight(iConfig.getParameter<bool>("useWeight")),
0103       m_useVisibleVertices(iConfig.getParameter<bool>("useVisibleVertices")),
0104       m_histoParameters(iConfig.getUntrackedParameter<edm::ParameterSet>("histoParameters", edm::ParameterSet())),
0105       m_doubleToken(consumes<double>(iConfig.getParameter<edm::InputTag>("weightProduct"))),
0106       m_vecPileupSummaryInfoToken(
0107           consumes<std::vector<PileupSummaryInfo> >(iConfig.getParameter<edm::InputTag>("pileupSummaryCollection"))),
0108       m_recoVertexCollectionToken(
0109           consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("pvCollection"))),
0110       m_hepMCProductToken(consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("mcTruthCollection"))) {
0111   //now do what ever initialization is needed
0112
0113   if (m_useVisibleVertices)
0114     edm::LogInfo("UseVisibleVertices") << "Only visible vertices will be used to compute Npileup";
0115
0116   usesResource("TFileService");
0117   edm::Service<TFileService> tfserv;
0118
0119   m_hrecovsmcnvtx2d = tfserv->make<TH2F>(
0120       "recovsmcnvtx2d", "Number of reco vertices vs pileup interactions", 60, -0.5, 59.5, 60, -0.5, 59.5);
0121   m_hrecovsmcnvtx2d->GetXaxis()->SetTitle("Pileup Interactions");
0122   m_hrecovsmcnvtx2d->GetYaxis()->SetTitle("Reco Vertices");
0123   m_hrecovsmcnvtxprof =
0124       tfserv->make<TProfile>("recovsmcnvtxprof", "Mean number of reco vs pileup vertices", 60, -0.5, 59.5);
0125   m_hrecovsmcnvtxprof->GetXaxis()->SetTitle("Pileup Interactions");
0126   m_hrecovsmcnvtxprof->GetYaxis()->SetTitle("Reco Vertices");
0127
0128   m_hrecovsmclumi2d = tfserv->make<TH2F>(
0129       "recovsmclumi2d", "Number of reco vertices vs ave pileup interactions", 200, 0., 50., 60, -0.5, 59.5);
0130   m_hrecovsmclumi2d->GetXaxis()->SetTitle("Average Pileup Interactions");
0131   m_hrecovsmclumi2d->GetYaxis()->SetTitle("Reco Vertices");
0132   m_hrecovsmclumiprof =
0133       tfserv->make<TProfile>("recovsmclumiprof", "Mean number of reco vs ave pileup vertices", 200, 0., 50.);
0134   m_hrecovsmclumiprof->GetXaxis()->SetTitle("Average Pileup Interactions");
0135   m_hrecovsmclumiprof->GetYaxis()->SetTitle("Reco Vertices");
0136
0137   if (m_useweight) {
0138     m_hrecovsmcnvtxweightedprof = tfserv->make<TProfile>(
0139         "recovsmcnvtxweightedprof", "Mean number of reco vs pileup vertices (1-w) weight", 60, -0.5, 59.5);
0140     m_hrecovsmcnvtxweightedprof->GetXaxis()->SetTitle("Pileup Interactions");
0141     m_hrecovsmcnvtxweightedprof->GetYaxis()->SetTitle("Reco Vertices (1-w)");
0142
0143     m_hrecovsmclumiweightedprof = tfserv->make<TProfile>(
0144         "recovsmclumiweightedprof", "Mean number of reco vs ave pileup vertices (1-w) weight", 200, 0., 50.);
0145     m_hrecovsmclumiweightedprof->GetXaxis()->SetTitle("Average Pileup Interactions");
0146     m_hrecovsmclumiweightedprof->GetYaxis()->SetTitle("Reco Vertices (1-w)");
0147   }
0148
0149   m_hdeltazfirst = tfserv->make<TH1F>("deltazfirst",
0150                                       "Reco-MC vertex z position (first vertex)",
0151                                       m_histoParameters.getUntrackedParameter<unsigned int>("zBins", 1000),
0152                                       m_histoParameters.getUntrackedParameter<double>("zMin", -1.),
0153                                       m_histoParameters.getUntrackedParameter<double>("zMax", 1.));
0154   m_hdeltazfirst->GetXaxis()->SetTitle("#Delta z (cm)");
0155   m_hdeltazfirst->GetYaxis()->SetTitle("Events");
0156
0157   m_hdeltazclose = tfserv->make<TH1F>("deltazclose",
0158                                       "Reco-MC vertex z position (closest vertex)",
0159                                       m_histoParameters.getUntrackedParameter<unsigned int>("zBins", 1000),
0160                                       m_histoParameters.getUntrackedParameter<double>("zMin", -1.),
0161                                       m_histoParameters.getUntrackedParameter<double>("zMax", 1.));
0162   m_hdeltazclose->GetXaxis()->SetTitle("#Delta z (cm)");
0163   m_hdeltazclose->GetYaxis()->SetTitle("Events");
0164
0165   m_hclosestvtx = tfserv->make<TH1F>("closestvtx", "Closest reco vtx ID", 30, -0.5, 29.5);
0166   m_hclosestvtx->GetXaxis()->SetTitle("Vtx ID");
0167   m_hclosestvtx->GetYaxis()->SetTitle("Events");
0168
0169   m_hdeltazfirstvsnpu = tfserv->make<TH2F>("deltazfirstvsnpu",
0170                                            "Reco-MC vertex z position (first vertex) vs Npileup",
0171                                            30,
0172                                            -0.5,
0173                                            29.5,
0174                                            m_histoParameters.getUntrackedParameter<unsigned int>("zBins", 1000),
0175                                            m_histoParameters.getUntrackedParameter<double>("zMin", -1.),
0176                                            m_histoParameters.getUntrackedParameter<double>("zMax", 1.));
0177   m_hdeltazfirstvsnpu->GetXaxis()->SetTitle("pileup Interactions");
0178   m_hdeltazfirstvsnpu->GetYaxis()->SetTitle("#Delta z (cm)");
0179
0180   m_hdeltazclosevsnpu = tfserv->make<TH2F>("deltazclosevsnpu",
0181                                            "Reco-MC vertex z position (closest vertex) v Npileup",
0182                                            30,
0183                                            -0.5,
0184                                            29.5,
0185                                            m_histoParameters.getUntrackedParameter<unsigned int>("zBins", 1000),
0186                                            m_histoParameters.getUntrackedParameter<double>("zMin", -1.),
0187                                            m_histoParameters.getUntrackedParameter<double>("zMax", 1.));
0188   m_hdeltazclosevsnpu->GetXaxis()->SetTitle("Pileup Interactions");
0189   m_hdeltazclosevsnpu->GetYaxis()->SetTitle("#Delta z (cm)");
0190
0191   m_hclosestvtxvsnpu =
0192       tfserv->make<TH2F>("closestvtxvsnpu", "Closest reco vtx ID vs Npileup", 30, -0.5, 29.5, 30, -0.5, 29.5);
0193   m_hclosestvtxvsnpu->GetXaxis()->SetTitle("Pileup Interactions");
0194   m_hclosestvtxvsnpu->GetYaxis()->SetTitle("Vtx ID");
0195 }
0196
0197 MCvsRecoVerticesAnalyzer::~MCvsRecoVerticesAnalyzer() {
0198   // do anything here that needs to be done at desctruction time
0199   // (e.g. close files, deallocate resources etc.)
0200 }
0201
0202 //
0203 // member functions
0204 //
0205
0206 // ------------ method called to for each event  ------------
0207 void MCvsRecoVerticesAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0208   double weight = 1.;
0209
0210   if (m_useweight) {
0211     edm::Handle<double> weightprod;
0212     iEvent.getByToken(m_doubleToken, weightprod);
0213
0214     weight = *weightprod;
0215   }
0216
0217   edm::Handle<std::vector<PileupSummaryInfo> > pileupinfos;
0218   iEvent.getByToken(m_vecPileupSummaryInfoToken, pileupinfos);
0219
0220   // look for the intime PileupSummaryInfo
0221
0222   std::vector<PileupSummaryInfo>::const_iterator pileupinfo;
0223
0224   for (pileupinfo = pileupinfos->begin(); pileupinfo != pileupinfos->end(); ++pileupinfo) {
0225     if (pileupinfo->getBunchCrossing() == 0)
0226       break;
0227   }
0228
0229   //
0230
0231   edm::Handle<reco::VertexCollection> pvcoll;
0232   iEvent.getByToken(m_recoVertexCollectionToken, pvcoll);
0233
0234   //
0235
0236   if (pileupinfo->getBunchCrossing() != 0) {
0237     edm::LogError("NoInTimePileUpInfo") << "Cannot find the in-time pileup info " << pileupinfo->getBunchCrossing();
0238
0239   } else {
0240     int npileup = pileupinfo->getPU_NumInteractions();
0241
0242     if (m_useVisibleVertices)
0243       npileup = pileupinfo->getPU_zpositions().size();
0244
0245     m_hrecovsmcnvtx2d->Fill(npileup, pvcoll->size(), weight);
0246     m_hrecovsmcnvtxprof->Fill(npileup, pvcoll->size(), weight);
0247
0248     m_hrecovsmclumi2d->Fill(pileupinfo->getTrueNumInteractions(), pvcoll->size(), weight);
0249     m_hrecovsmclumiprof->Fill(pileupinfo->getTrueNumInteractions(), pvcoll->size(), weight);
0250
0251     if (m_useweight) {
0252       m_hrecovsmcnvtxweightedprof->Fill(npileup, pvcoll->size(), 1. - weight);
0253       m_hrecovsmclumiweightedprof->Fill(pileupinfo->getTrueNumInteractions(), pvcoll->size(), 1. - weight);
0254     }
0255     //
0256
0257     edm::Handle<edm::HepMCProduct> EvtHandle;
0258     iEvent.getByToken(m_hepMCProductToken, EvtHandle);
0259
0260     const HepMC::GenEvent* Evt = EvtHandle->GetEvent();
0261
0262     // compute the difference between the main interaction vertex z position and the first vertex of the collection
0263
0264     if (!pvcoll->empty()) {
0265       if (!(*pvcoll)[0].isFake()) {
0266         // get the first vertex
0267         if (Evt->vertices_begin() != Evt->vertices_end()) {
0268           m_hdeltazfirst->Fill((*pvcoll)[0].z() - (*Evt->vertices_begin())->point3d().z() / 10., weight);
0269           m_hdeltazfirstvsnpu->Fill(npileup, (*pvcoll)[0].z() - (*Evt->vertices_begin())->point3d().z() / 10., weight);
0270         }
0271       }
0272     }
0273
0274     // compute the difference between the main interaction vertex z position and the closest reco vertex
0275
0276     double minabsdist = -1.;
0277     double mindist = -999.;
0278     int closestvtx = -1;
0279
0280     for (unsigned int ivtx = 0; ivtx < pvcoll->size(); ++ivtx) {
0281       if (closestvtx < 0 ||
0282           minabsdist > std::abs((*pvcoll)[ivtx].z() - (*Evt->vertices_begin())->point3d().z() / 10.)) {
0283         mindist = (*pvcoll)[ivtx].z() - (*Evt->vertices_begin())->point3d().z() / 10.;
0284         closestvtx = ivtx;
0285         minabsdist = std::abs(mindist);
0286       }
0287     }
0288     if (closestvtx >= 0) {
0289       m_hdeltazclose->Fill(mindist, weight);
0290       m_hdeltazclosevsnpu->Fill(npileup, mindist, weight);
0291       m_hclosestvtx->Fill(closestvtx, weight);
0292       m_hclosestvtxvsnpu->Fill(npileup, closestvtx, weight);
0293     }
0294   }
0295 }
0296
0297 //define this as a plug-in
0298 DEFINE_FWK_MODULE(MCvsRecoVerticesAnalyzer);