RecoVertex/src/MCVerticesAnalyzer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    MCVerticesAnalyzer
0004 // Class:      MCVerticesAnalyzer
0005 //
0006 /**\class MCVerticesAnalyzer MCVerticesAnalyzer.cc TrackingPFG/PileUp/src/MCVerticesAnalyzer.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 "TH1F.h"
0047 #include "TH2F.h"
0048 #include "TProfile.h"
0049
0050 //
0051 // class decleration
0052 //
0053
0054 class MCVerticesAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0055 public:
0056   explicit MCVerticesAnalyzer(const edm::ParameterSet&);
0057   ~MCVerticesAnalyzer() override;
0058
0059 private:
0060   void analyze(const edm::Event&, const edm::EventSetup&) override;
0061
0062   // ----------member data ---------------------------
0063
0064   const bool m_useweight;
0065
0066   edm::EDGetTokenT<double> m_doubleToken;
0067   edm::EDGetTokenT<std::vector<PileupSummaryInfo> > m_vecPileupSummaryInfoToken;
0068   edm::EDGetTokenT<edm::HepMCProduct> m_hepMCProductToken;
0069
0070   TH1F* m_hnvtx;
0071   TH2F* m_hnvtxvsbx;
0072   TH1F* m_hlumi;
0073   TH2F* m_hnvtxvslumi;
0074   TH1F* m_hnvtxweight;
0075   TProfile* m_hnvtxweightprof;
0076   TH1F* m_hmainvtxx;
0077   TH1F* m_hmainvtxy;
0078   TH1F* m_hmainvtxz;
0079   TH1F* m_hpileupvtxz;
0080 };
0081
0082 //
0083 // constants, enums and typedefs
0084 //
0085
0086 //
0087 // static data member definitions
0088 //
0089
0090 //
0091 // constructors and destructor
0092 //
0093 MCVerticesAnalyzer::MCVerticesAnalyzer(const edm::ParameterSet& iConfig)
0094     : m_useweight(iConfig.getParameter<bool>("useWeight")),
0095       m_doubleToken(consumes<double>(iConfig.getParameter<edm::InputTag>("weightProduct"))),
0096       m_vecPileupSummaryInfoToken(
0097           consumes<std::vector<PileupSummaryInfo> >(iConfig.getParameter<edm::InputTag>("pileupSummaryCollection"))),
0098       m_hepMCProductToken(consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("mcTruthCollection"))) {
0099   //now do what ever initialization is needed
0100
0101   usesResource("TFileService");
0102   edm::Service<TFileService> tfserv;
0103
0104   m_hnvtx = tfserv->make<TH1F>("nvtx", "Number of pileup vertices", 400, -0.5, 399.5);
0105   m_hnvtx->GetXaxis()->SetTitle("Number of Interactions");
0106
0107   m_hnvtxvsbx = tfserv->make<TH2F>("nvtxvsbx", "Number of pileup vertices vs BX", 9, -4.5, 4.5, 400, -0.5, 399.5);
0108   m_hnvtxvsbx->GetXaxis()->SetTitle("BX number");
0109   m_hnvtxvsbx->GetYaxis()->SetTitle("Number of Interactions");
0110
0111   m_hlumi = tfserv->make<TH1F>("lumi", "BX luminosity*xsect", 200, 0., 50.);
0112   m_hlumi->GetXaxis()->SetTitle("Average Number of Interactions");
0113
0114   m_hnvtxvslumi = tfserv->make<TH2F>("nvtxvslumi", "Npileup vs BX luminosity*xsect", 200, 0., 50., 400, -0.5, 399.5);
0115   m_hnvtxvslumi->GetXaxis()->SetTitle("Average Number of Interactions");
0116   m_hnvtxvslumi->GetYaxis()->SetTitle("Number of Interactions");
0117
0118   if (m_useweight) {
0119     m_hnvtxweight = tfserv->make<TH1F>("nvtxweight", "Number of pileup vertices (1-w)", 400, -0.5, 399.5);
0120     m_hnvtxweight->GetXaxis()->SetTitle("Number of Interactions");
0121     m_hnvtxweightprof =
0122         tfserv->make<TProfile>("nvtxweightprof", "Mean (1-w) vs Number of pileup interactions", 400, -0.5, 399.5);
0123     m_hnvtxweightprof->GetXaxis()->SetTitle("Number of Interactions");
0124   }
0125
0126   m_hmainvtxx = tfserv->make<TH1F>("mainvtxx", "Main vertex x position", 200, -.5, .5);
0127   m_hmainvtxx->GetXaxis()->SetTitle("X (cm)");
0128   m_hmainvtxy = tfserv->make<TH1F>("mainvtxy", "Main vertex y position", 200, -.5, .5);
0129   m_hmainvtxy->GetXaxis()->SetTitle("Y (cm)");
0130   m_hmainvtxz = tfserv->make<TH1F>("mainvtxz", "Main vertex z position", 600, -30., 30.);
0131   m_hmainvtxz->GetXaxis()->SetTitle("Z (cm)");
0132   m_hpileupvtxz = tfserv->make<TH1F>("pileupvtxz", "PileUp vertices z position", 600, -30., 30.);
0133   m_hpileupvtxz->GetXaxis()->SetTitle("Z (cm)");
0134 }
0135
0136 MCVerticesAnalyzer::~MCVerticesAnalyzer() {
0137   // do anything here that needs to be done at desctruction time
0138   // (e.g. close files, deallocate resources etc.)
0139 }
0140
0141 //
0142 // member functions
0143 //
0144
0145 // ------------ method called to for each event  ------------
0146 void MCVerticesAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0147   double weight = 1.;
0148
0149   if (m_useweight) {
0150     edm::Handle<double> weightprod;
0151     iEvent.getByToken(m_doubleToken, weightprod);
0152
0153     weight = *weightprod;
0154   }
0155
0156   edm::Handle<std::vector<PileupSummaryInfo> > pileupinfos;
0157   iEvent.getByToken(m_vecPileupSummaryInfoToken, pileupinfos);
0158
0159   //
0160
0161   if (pileupinfos.isValid()) {
0162     // look for the intime PileupSummaryInfo
0163
0164     std::vector<PileupSummaryInfo>::const_iterator pileupinfo;
0165
0166     for (pileupinfo = pileupinfos->begin(); pileupinfo != pileupinfos->end(); ++pileupinfo) {
0167       m_hnvtxvsbx->Fill(pileupinfo->getBunchCrossing(), pileupinfo->getPU_NumInteractions(), weight);
0168     }
0169
0170     for (pileupinfo = pileupinfos->begin(); pileupinfo != pileupinfos->end(); ++pileupinfo) {
0171       if (pileupinfo->getBunchCrossing() == 0)
0172         break;
0173     }
0174
0175     //
0176
0177     if (pileupinfo->getBunchCrossing() != 0) {
0178       edm::LogError("NoInTimePileUpInfo") << "Cannot find the in-time pileup info " << pileupinfo->getBunchCrossing();
0179     } else {
0180       m_hlumi->Fill(pileupinfo->getTrueNumInteractions(), weight);
0181       m_hnvtx->Fill(pileupinfo->getPU_NumInteractions(), weight);
0182       m_hnvtxvslumi->Fill(pileupinfo->getTrueNumInteractions(), pileupinfo->getPU_NumInteractions(), weight);
0183
0184       if (m_useweight) {
0185         m_hnvtxweight->Fill(pileupinfo->getPU_NumInteractions(), 1. - weight);
0186         m_hnvtxweightprof->Fill(pileupinfo->getPU_NumInteractions(), 1. - weight);
0187       }
0188
0189       const std::vector<float>& zpositions = pileupinfo->getPU_zpositions();
0190
0191       for (std::vector<float>::const_iterator zpos = zpositions.begin(); zpos != zpositions.end(); ++zpos) {
0192         m_hpileupvtxz->Fill(*zpos, weight);
0193       }
0194     }
0195   }
0196   // main interaction part
0197
0198   edm::Handle<edm::HepMCProduct> EvtHandle;
0199   iEvent.getByToken(m_hepMCProductToken, EvtHandle);
0200
0201   if (EvtHandle.isValid()) {
0202     const HepMC::GenEvent* Evt = EvtHandle->GetEvent();
0203
0204     // get the first vertex
0205
0206     if (Evt->vertices_begin() != Evt->vertices_end()) {
0207       m_hmainvtxx->Fill((*Evt->vertices_begin())->point3d().x() / 10., weight);
0208       m_hmainvtxy->Fill((*Evt->vertices_begin())->point3d().y() / 10., weight);
0209       m_hmainvtxz->Fill((*Evt->vertices_begin())->point3d().z() / 10., weight);
0210     }
0211   }
0212 }
0213
0214 //define this as a plug-in
0215 DEFINE_FWK_MODULE(MCVerticesAnalyzer);