#include <TrackingPFG/PileUp/src/MCvsRecoVerticesAnalyzer.cc>

Inheritance diagram for MCvsRecoVerticesAnalyzer:

Public Member Functions
	MCvsRecoVerticesAnalyzer (const edm::ParameterSet &)
	~MCvsRecoVerticesAnalyzer ()
Private Member Functions
virtual void	analyze (const edm::Event &, const edm::EventSetup &)
virtual void	beginJob ()
virtual void	beginRun (const edm::Run &, const edm::EventSetup &)
virtual void	endJob ()
virtual void	endRun (const edm::Run &, const edm::EventSetup &)
Private Attributes
TH1F *	m_hclosestvtx
TH2F *	m_hclosestvtxvsnpu
TH1F *	m_hdeltazclose
TH2F *	m_hdeltazclosevsnpu
TH1F *	m_hdeltazfirst
TH2F *	m_hdeltazfirstvsnpu
const edm::ParameterSet	m_histoParameters
TH2F *	m_hrecovsmclumi2d
TProfile *	m_hrecovsmclumiprof
TProfile *	m_hrecovsmclumiweightedprof
TH2F *	m_hrecovsmcnvtx2d
TProfile *	m_hrecovsmcnvtxprof
TProfile *	m_hrecovsmcnvtxweightedprof
edm::InputTag	m_mctruthcollection
edm::InputTag	m_pileupcollection
edm::InputTag	m_pvcollection
const bool	m_useVisibleVertices
const bool	m_useweight
edm::InputTag	m_weight

Detailed Description

Description: <one line="" class="" summary>="">

Implementation: <Notes on="" implementation>="">

Definition at line 60 of file MCvsRecoVerticesAnalyzer.cc.

Constructor & Destructor Documentation

MCvsRecoVerticesAnalyzer::MCvsRecoVerticesAnalyzer ( const edm::ParameterSet & iConfig ) [explicit]

Definition at line 112 of file MCvsRecoVerticesAnalyzer.cc.

References edm::ParameterSet::getUntrackedParameter(), m_hclosestvtx, m_hclosestvtxvsnpu, m_hdeltazclose, m_hdeltazclosevsnpu, m_hdeltazfirst, m_hdeltazfirstvsnpu, m_histoParameters, m_hrecovsmclumi2d, m_hrecovsmclumiprof, m_hrecovsmclumiweightedprof, m_hrecovsmcnvtx2d, m_hrecovsmcnvtxprof, m_hrecovsmcnvtxweightedprof, m_useVisibleVertices, and m_useweight.

                                                                                :
  m_pileupcollection(iConfig.getParameter<edm::InputTag>("pileupSummaryCollection")),
  m_mctruthcollection(iConfig.getParameter<edm::InputTag>("mcTruthCollection")),
  m_pvcollection(iConfig.getParameter<edm::InputTag>("pvCollection")),
  m_useweight(iConfig.getParameter<bool>("useWeight")),
  m_weight(iConfig.getParameter<edm::InputTag>("weightProduct")),
  m_useVisibleVertices(iConfig.getParameter<bool>("useVisibleVertices")),
  m_histoParameters(iConfig.getUntrackedParameter<edm::ParameterSet>("histoParameters",edm::ParameterSet()))
{
   //now do what ever initialization is needed

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

  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 ( )

Definition at line 180 of file MCvsRecoVerticesAnalyzer.cc.

{

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

}

Member Function Documentation

void MCvsRecoVerticesAnalyzer::analyze	(	const edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)		`[private, virtual]`

Implements edm::EDAnalyzer.

Definition at line 195 of file MCvsRecoVerticesAnalyzer.cc.

References abs, edm::Event::getByLabel(), m_hclosestvtx, m_hclosestvtxvsnpu, m_hdeltazclose, m_hdeltazclosevsnpu, m_hdeltazfirst, m_hdeltazfirstvsnpu, m_hrecovsmclumi2d, m_hrecovsmclumiprof, m_hrecovsmclumiweightedprof, m_hrecovsmcnvtx2d, m_hrecovsmcnvtxprof, m_hrecovsmcnvtxweightedprof, m_mctruthcollection, m_pileupcollection, m_pvcollection, m_useVisibleVertices, m_useweight, m_weight, histoStyle::weight, and z.

{
  using namespace edm;
  
  double weight = 1.;
  
  if(m_useweight) {
    Handle<double> weightprod;
    iEvent.getByLabel(m_weight,weightprod);
    
    weight = *weightprod;
    
  }
  
  Handle<std::vector<PileupSummaryInfo> > pileupinfos;
  iEvent.getByLabel(m_pileupcollection,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;

  } 
  
  //
  
  Handle<reco::VertexCollection> pvcoll;
  iEvent.getByLabel(m_pvcollection,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);
    }
    //
    
    Handle< HepMCProduct > EvtHandle ;
    iEvent.getByLabel(m_mctruthcollection, 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->size() !=0) {
      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);
    }
    
  }
}