MCvsRecoVerticesAnalyzer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    MCvsRecoVerticesAnalyzer
// Class:      MCvsRecoVerticesAnalyzer
//
//
// Original Author:  Andrea Venturi
//         Created:  Thu Dec 16 16:32:56 CEST 2010
//
//

// system include files
#include <memory>
#include <numeric>

#include <vector>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/Run.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"

#include "FWCore/Utilities/interface/InputTag.h"

#include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"

#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"

#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"

//
// class decleration
//

class MCvsRecoVerticesAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
public:
  explicit MCvsRecoVerticesAnalyzer(const edm::ParameterSet&);
  ~MCvsRecoVerticesAnalyzer() override;

private:
  void analyze(const edm::Event&, const edm::EventSetup&) override;

  // ----------member data ---------------------------

  const bool m_useweight;
  const bool m_useVisibleVertices;
  const edm::ParameterSet m_histoParameters;
  edm::EDGetTokenT<double> m_doubleToken;
  edm::EDGetTokenT<std::vector<PileupSummaryInfo> > m_vecPileupSummaryInfoToken;
  edm::EDGetTokenT<reco::VertexCollection> m_recoVertexCollectionToken;
  edm::EDGetTokenT<edm::HepMCProduct> m_hepMCProductToken;

  TH2F* m_hrecovsmcnvtx2d;
  TProfile* m_hrecovsmcnvtxprof;
  TProfile* m_hrecovsmcnvtxweightedprof;
  TH2F* m_hrecovsmclumi2d;
  TProfile* m_hrecovsmclumiprof;
  TProfile* m_hrecovsmclumiweightedprof;
  TH1F* m_hdeltazfirst;
  TH1F* m_hdeltazclose;
  TH1F* m_hclosestvtx;

  TH2F* m_hdeltazfirstvsnpu;
  TH2F* m_hdeltazclosevsnpu;
  TH2F* m_hclosestvtxvsnpu;
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
MCvsRecoVerticesAnalyzer::MCvsRecoVerticesAnalyzer(const edm::ParameterSet& iConfig)
    : m_useweight(iConfig.getParameter<bool>("useWeight")),
      m_useVisibleVertices(iConfig.getParameter<bool>("useVisibleVertices")),
      m_histoParameters(iConfig.getUntrackedParameter<edm::ParameterSet>("histoParameters", edm::ParameterSet())),
      m_doubleToken(consumes<double>(iConfig.getParameter<edm::InputTag>("weightProduct"))),
      m_vecPileupSummaryInfoToken(
          consumes<std::vector<PileupSummaryInfo> >(iConfig.getParameter<edm::InputTag>("pileupSummaryCollection"))),
      m_recoVertexCollectionToken(
          consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("pvCollection"))),
      m_hepMCProductToken(consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("mcTruthCollection"))) {
  //now do what ever initialization is needed

  if (m_useVisibleVertices)
    edm::LogInfo("UseVisibleVertices") << "Only visible vertices will be used to compute Npileup";

  usesResource("TFileService");
  edm::Service<TFileService> tfserv;

  m_hrecovsmcnvtx2d = tfserv->make<TH2F>(
      "recovsmcnvtx2d", "Number of reco vertices vs pileup interactions", 60, -0.5, 59.5, 60, -0.5, 59.5);
  m_hrecovsmcnvtx2d->GetXaxis()->SetTitle("Pileup Interactions");
  m_hrecovsmcnvtx2d->GetYaxis()->SetTitle("Reco Vertices");
  m_hrecovsmcnvtxprof =
      tfserv->make<TProfile>("recovsmcnvtxprof", "Mean number of reco vs pileup vertices", 60, -0.5, 59.5);
  m_hrecovsmcnvtxprof->GetXaxis()->SetTitle("Pileup Interactions");
  m_hrecovsmcnvtxprof->GetYaxis()->SetTitle("Reco Vertices");

  m_hrecovsmclumi2d = tfserv->make<TH2F>(
      "recovsmclumi2d", "Number of reco vertices vs ave pileup interactions", 200, 0., 50., 60, -0.5, 59.5);
  m_hrecovsmclumi2d->GetXaxis()->SetTitle("Average Pileup Interactions");
  m_hrecovsmclumi2d->GetYaxis()->SetTitle("Reco Vertices");
  m_hrecovsmclumiprof =
      tfserv->make<TProfile>("recovsmclumiprof", "Mean number of reco vs ave pileup vertices", 200, 0., 50.);
  m_hrecovsmclumiprof->GetXaxis()->SetTitle("Average Pileup Interactions");
  m_hrecovsmclumiprof->GetYaxis()->SetTitle("Reco Vertices");

  if (m_useweight) {
    m_hrecovsmcnvtxweightedprof = tfserv->make<TProfile>(
        "recovsmcnvtxweightedprof", "Mean number of reco vs pileup vertices (1-w) weight", 60, -0.5, 59.5);
    m_hrecovsmcnvtxweightedprof->GetXaxis()->SetTitle("Pileup Interactions");
    m_hrecovsmcnvtxweightedprof->GetYaxis()->SetTitle("Reco Vertices (1-w)");

    m_hrecovsmclumiweightedprof = tfserv->make<TProfile>(
        "recovsmclumiweightedprof", "Mean number of reco vs ave pileup vertices (1-w) weight", 200, 0., 50.);
    m_hrecovsmclumiweightedprof->GetXaxis()->SetTitle("Average Pileup Interactions");
    m_hrecovsmclumiweightedprof->GetYaxis()->SetTitle("Reco Vertices (1-w)");
  }

  m_hdeltazfirst = tfserv->make<TH1F>("deltazfirst",
                                      "Reco-MC vertex z position (first vertex)",
                                      m_histoParameters.getUntrackedParameter<unsigned int>("zBins", 1000),
                                      m_histoParameters.getUntrackedParameter<double>("zMin", -1.),
                                      m_histoParameters.getUntrackedParameter<double>("zMax", 1.));
  m_hdeltazfirst->GetXaxis()->SetTitle("#Delta z (cm)");
  m_hdeltazfirst->GetYaxis()->SetTitle("Events");

  m_hdeltazclose = tfserv->make<TH1F>("deltazclose",
                                      "Reco-MC vertex z position (closest vertex)",
                                      m_histoParameters.getUntrackedParameter<unsigned int>("zBins", 1000),
                                      m_histoParameters.getUntrackedParameter<double>("zMin", -1.),
                                      m_histoParameters.getUntrackedParameter<double>("zMax", 1.));
  m_hdeltazclose->GetXaxis()->SetTitle("#Delta z (cm)");
  m_hdeltazclose->GetYaxis()->SetTitle("Events");

  m_hclosestvtx = tfserv->make<TH1F>("closestvtx", "Closest reco vtx ID", 30, -0.5, 29.5);
  m_hclosestvtx->GetXaxis()->SetTitle("Vtx ID");
  m_hclosestvtx->GetYaxis()->SetTitle("Events");

  m_hdeltazfirstvsnpu = tfserv->make<TH2F>("deltazfirstvsnpu",
                                           "Reco-MC vertex z position (first vertex) vs Npileup",
                                           30,
                                           -0.5,
                                           29.5,
                                           m_histoParameters.getUntrackedParameter<unsigned int>("zBins", 1000),
                                           m_histoParameters.getUntrackedParameter<double>("zMin", -1.),
                                           m_histoParameters.getUntrackedParameter<double>("zMax", 1.));
  m_hdeltazfirstvsnpu->GetXaxis()->SetTitle("pileup Interactions");
  m_hdeltazfirstvsnpu->GetYaxis()->SetTitle("#Delta z (cm)");

  m_hdeltazclosevsnpu = tfserv->make<TH2F>("deltazclosevsnpu",
                                           "Reco-MC vertex z position (closest vertex) v Npileup",
                                           30,
                                           -0.5,
                                           29.5,
                                           m_histoParameters.getUntrackedParameter<unsigned int>("zBins", 1000),
                                           m_histoParameters.getUntrackedParameter<double>("zMin", -1.),
                                           m_histoParameters.getUntrackedParameter<double>("zMax", 1.));
  m_hdeltazclosevsnpu->GetXaxis()->SetTitle("Pileup Interactions");
  m_hdeltazclosevsnpu->GetYaxis()->SetTitle("#Delta z (cm)");

  m_hclosestvtxvsnpu =
      tfserv->make<TH2F>("closestvtxvsnpu", "Closest reco vtx ID vs Npileup", 30, -0.5, 29.5, 30, -0.5, 29.5);
  m_hclosestvtxvsnpu->GetXaxis()->SetTitle("Pileup Interactions");
  m_hclosestvtxvsnpu->GetYaxis()->SetTitle("Vtx ID");
}

MCvsRecoVerticesAnalyzer::~MCvsRecoVerticesAnalyzer() {
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
}

//
// member functions
//

// ------------ method called to for each event  ------------
void MCvsRecoVerticesAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  double weight = 1.;

  if (m_useweight) {
    edm::Handle<double> weightprod;
    iEvent.getByToken(m_doubleToken, weightprod);

    weight = *weightprod;
  }

  edm::Handle<std::vector<PileupSummaryInfo> > pileupinfos;
  iEvent.getByToken(m_vecPileupSummaryInfoToken, pileupinfos);

  // look for the intime PileupSummaryInfo

  std::vector<PileupSummaryInfo>::const_iterator pileupinfo;

  for (pileupinfo = pileupinfos->begin(); pileupinfo != pileupinfos->end(); ++pileupinfo) {
    if (pileupinfo->getBunchCrossing() == 0)
      break;
  }

  //

  edm::Handle<reco::VertexCollection> pvcoll;
  iEvent.getByToken(m_recoVertexCollectionToken, pvcoll);

  //

  if (pileupinfo->getBunchCrossing() != 0) {
    edm::LogError("NoInTimePileUpInfo") << "Cannot find the in-time pileup info " << pileupinfo->getBunchCrossing();

  } else {
    int npileup = pileupinfo->getPU_NumInteractions();

    if (m_useVisibleVertices)
      npileup = pileupinfo->getPU_zpositions().size();

    m_hrecovsmcnvtx2d->Fill(npileup, pvcoll->size(), weight);
    m_hrecovsmcnvtxprof->Fill(npileup, pvcoll->size(), weight);

    m_hrecovsmclumi2d->Fill(pileupinfo->getTrueNumInteractions(), pvcoll->size(), weight);
    m_hrecovsmclumiprof->Fill(pileupinfo->getTrueNumInteractions(), pvcoll->size(), weight);

    if (m_useweight) {
      m_hrecovsmcnvtxweightedprof->Fill(npileup, pvcoll->size(), 1. - weight);
      m_hrecovsmclumiweightedprof->Fill(pileupinfo->getTrueNumInteractions(), pvcoll->size(), 1. - weight);
    }
    //

    edm::Handle<edm::HepMCProduct> EvtHandle;
    iEvent.getByToken(m_hepMCProductToken, EvtHandle);

    const HepMC::GenEvent* Evt = EvtHandle->GetEvent();

    // compute the difference between the main interaction vertex z position and the first vertex of the collection

    if (!pvcoll->empty()) {
      if (!(*pvcoll)[0].isFake()) {
        // get the first vertex
        if (Evt->vertices_begin() != Evt->vertices_end()) {
          m_hdeltazfirst->Fill((*pvcoll)[0].z() - (*Evt->vertices_begin())->point3d().z() / 10., weight);
          m_hdeltazfirstvsnpu->Fill(npileup, (*pvcoll)[0].z() - (*Evt->vertices_begin())->point3d().z() / 10., weight);
        }
      }
    }

    // compute the difference between the main interaction vertex z position and the closest reco vertex

    double minabsdist = -1.;
    double mindist = -999.;
    int closestvtx = -1;

    for (unsigned int ivtx = 0; ivtx < pvcoll->size(); ++ivtx) {
      if (closestvtx < 0 ||
          minabsdist > std::abs((*pvcoll)[ivtx].z() - (*Evt->vertices_begin())->point3d().z() / 10.)) {
        mindist = (*pvcoll)[ivtx].z() - (*Evt->vertices_begin())->point3d().z() / 10.;
        closestvtx = ivtx;
        minabsdist = std::abs(mindist);
      }
    }
    if (closestvtx >= 0) {
      m_hdeltazclose->Fill(mindist, weight);
      m_hdeltazclosevsnpu->Fill(npileup, mindist, weight);
      m_hclosestvtx->Fill(closestvtx, weight);
      m_hclosestvtxvsnpu->Fill(npileup, closestvtx, weight);
    }
  }
}

//define this as a plug-in
DEFINE_FWK_MODULE(MCvsRecoVerticesAnalyzer);