RecoPixelVertexing/test/PixelVertexVal.cc

0001 #include "DataFormats/Common/interface/Handle.h"
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003 #include "FWCore/Framework/interface/ESHandle.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/ServiceRegistry/interface/Service.h"
0011
0012 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0013 #include "DataFormats/TrackReco/interface/Track.h"
0014 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0015 #include "DataFormats/VertexReco/interface/Vertex.h"
0016 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0017
0018 #include "HepMC/GenEvent.h"
0019 #include "HepMC/GenVertex.h"
0020 #include "SimDataFormats/Track/interface/SimTrack.h"
0021 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0022 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0023 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0024
0025 #include "TDirectory.h"
0026 #include "TFile.h"
0027 #include "TTree.h"
0028 #include <TH1.h>
0029 #include <cmath>
0030 #include <iostream>
0031 #include <map>
0032 #include <vector>
0033
0034 using namespace std;
0035
0036 class PixelVertexVal : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0037 public:
0038   explicit PixelVertexVal(const edm::ParameterSet &conf);
0039   ~PixelVertexVal() override;
0040   void beginJob() override;
0041   void analyze(const edm::Event &ev, const edm::EventSetup &es) override;
0042
0043 private:
0044   edm::ParameterSet conf_;
0045   int verbose_;
0046   std::string file_;
0047   std::map<std::string, TH1 *> h;
0048   edm::EDGetTokenT<reco::TrackCollection> trackCollectionToken_;
0049   edm::EDGetTokenT<reco::VertexCollection> vertexCollectionToken_;
0050   edm::EDGetTokenT<edm::SimVertexContainer> simVertexContainerToken_;
0051 };
0052
0053 PixelVertexVal::PixelVertexVal(const edm::ParameterSet &conf)
0054     : verbose_(conf.getUntrackedParameter<unsigned int>("Verbosity",
0055                                                         0))  // How noisy?
0056       ,
0057       file_(conf.getUntrackedParameter<std::string>("HistoFile", "pixelVertexHistos.root")),
0058       h(),
0059       trackCollectionToken_(
0060           consumes<reco::TrackCollection>(edm::InputTag(conf.getParameter<std::string>("TrackCollection")))),
0061       vertexCollectionToken_(
0062           consumes<reco::VertexCollection>(edm::InputTag(conf.getParameter<std::string>("VertexCollection")))),
0063       simVertexContainerToken_(consumes<edm::SimVertexContainer>(conf.getParameter<edm::InputTag>("simG4"))) {
0064   usesResource(TFileService::kSharedResource);
0065   edm::LogInfo("PixelVertexVal") << " CTOR";
0066 }
0067
0068 PixelVertexVal::~PixelVertexVal() { edm::LogInfo("PixelVertexVal") << " DTOR"; }
0069
0070 void PixelVertexVal::beginJob() {
0071   // validation histos
0072   edm::Service<TFileService> fs;
0073   h["h_Nbvtx"] = fs->make<TH1F>("h_nbvtx", "nb vertices in event", 16, 0., 16.);
0074   h["h_Nbtrks"] = fs->make<TH1F>("h_Nbtrks", "nb tracks in PV", 100, 0., 100.);
0075   h["h_ResZ"] = fs->make<TH1F>("resz", "residual z", 100, -0.1, 0.1);
0076   h["h_PullZ"] = fs->make<TH1F>("pullz", "pull z", 100, -25., 25.);
0077   h["h_TrkRes"] = fs->make<TH1F>("h_TrkRes", "h_TrkRes", 100, -0.2, 0.2);
0078   h["h_Eff"] = fs->make<TH1F>("h_Etff", "h_Etff", 10, -1., 9.);
0079 }
0080
0081 void PixelVertexVal::analyze(const edm::Event &ev, const edm::EventSetup &es) {
0082   if (verbose_ > 0)
0083     cout << "------------------------------------------------" << endl;
0084   cout << "*** PixelVertexVal, analyze event: " << ev.id() << endl;
0085   edm::Handle<reco::TrackCollection> trackCollection;
0086   ev.getByToken(trackCollectionToken_, trackCollection);
0087   const reco::TrackCollection tracks = *(trackCollection.product());
0088
0089   reco::TrackRefVector trks;
0090
0091   //  if (verbose_ > 0) {
0092   //    std::cout << *(trackCollection.provenance()) << std::endl;
0093   //    cout << "Reconstructed "<< tracks.size() << " tracks" << std::endl;
0094   //  }
0095
0096   edm::Handle<reco::VertexCollection> vertexCollection;
0097   ev.getByToken(vertexCollectionToken_, vertexCollection);
0098   const reco::VertexCollection vertexes = *(vertexCollection.product());
0099   if (verbose_ > 0) {
0100     //    std::cout << *(vertexCollection.provenance()) << std::endl;
0101     cout << "Reconstructed " << vertexes.size() << " vertexes" << std::endl;
0102   }
0103
0104   edm::Handle<edm::SimVertexContainer> simVtcs;
0105   ev.getByToken(simVertexContainerToken_, simVtcs);
0106   if (verbose_ > 0) {
0107     cout << "simulated vertices: " << simVtcs->size() << std::endl;
0108   }
0109
0110   bool hasPV = (!simVtcs->empty());
0111   if (!hasPV)
0112     cout << "Event without PV!, skip" << endl;
0113   float z_PV = hasPV ? (*simVtcs)[0].position().z() : 0.;
0114
0115   int nVertices = vertexes.size();
0116   h["h_Nbvtx"]->Fill(nVertices);
0117   if (nVertices > 0)
0118     h["h_Nbtrks"]->Fill(vertexes[0].tracksSize());
0119
0120   int ivtx_matched = -1;
0121   float min_dist = 0.1;
0122   for (int idx = 0; idx < nVertices; idx++) {
0123     float dz = vertexes[idx].position().z() - z_PV;
0124     if (fabs(dz) < min_dist) {
0125       ivtx_matched = idx;
0126       min_dist = fabs(dz);
0127     }
0128   }
0129   h["h_Eff"]->Fill(ivtx_matched + 0.00001);
0130
0131   if (ivtx_matched == 0) {
0132     const reco::Vertex &pv = vertexes[ivtx_matched];
0133     float dz = pv.position().z() - z_PV;
0134     h["h_ResZ"]->Fill(dz);
0135     h["h_PullZ"]->Fill(dz / pv.zError());
0136
0137     for (reco::Vertex::trackRef_iterator it = pv.tracks_begin(); it != pv.tracks_end(); it++) {
0138       //    for (reco::Vertex::track_iterator it=pv.tracks_begin(); it !=
0139       //    pv.tracks_end(); it++) { for (reco::TrackRefVector::iterator
0140       //    it=pv.tracks_begin(); it != pv.tracks_end(); it++) {
0141       // h["h_TrkRes"]->Fill((*it)->dz());
0142       h["h_TrkRes"]->Fill((*it)->vertex().z() - pv.position().z());
0143     }
0144   }
0145 }
0146
0147 DEFINE_FWK_MODULE(PixelVertexVal);