NanoAOD/plugins/VertexTableProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    PhysicsTools/NanoAOD
0004 // Class:      VertexTableProducer
0005 //
0006 /**\class VertexTableProducer VertexTableProducer.cc PhysicsTools/VertexTableProducer/plugins/VertexTableProducer.cc
0007
0008  Description: [one line class summary]
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Andrea Rizzi
0015 //         Created:  Mon, 28 Aug 2017 09:26:39 GMT
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/stream/EDProducer.h"
0025
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/Utilities/interface/StreamID.h"
0031
0032 #include "DataFormats/VertexReco/interface/Vertex.h"
0033 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0034
0035 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0036
0037 #include "DataFormats/NanoAOD/interface/FlatTable.h"
0038 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0039 #include "RecoVertex/VertexTools/interface/VertexDistanceXY.h"
0040 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0041 #include "RecoVertex/VertexPrimitives/interface/VertexState.h"
0042 #include "DataFormats/Common/interface/ValueMap.h"
0043
0044 //
0045 // class declaration
0046 //
0047
0048 class VertexTableProducer : public edm::stream::EDProducer<> {
0049 public:
0050   explicit VertexTableProducer(const edm::ParameterSet&);
0051   ~VertexTableProducer() override;
0052
0053   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0054
0055 private:
0056   void beginStream(edm::StreamID) override;
0057   void produce(edm::Event&, const edm::EventSetup&) override;
0058   void endStream() override;
0059
0060   //virtual void beginRun(edm::Run const&, edm::EventSetup const&) override;
0061   //virtual void endRun(edm::Run const&, edm::EventSetup const&) override;
0062   //virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override;
0063   //virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override;
0064
0065   // ----------member data ---------------------------
0066
0067   const edm::EDGetTokenT<std::vector<reco::Vertex>> pvs_;
0068   const edm::EDGetTokenT<edm::ValueMap<float>> pvsScore_;
0069   const edm::EDGetTokenT<edm::View<reco::VertexCompositePtrCandidate>> svs_;
0070   const StringCutObjectSelector<reco::Candidate> svCut_;
0071   const StringCutObjectSelector<reco::Vertex> goodPvCut_;
0072   const std::string goodPvCutString_;
0073   const std::string pvName_;
0074   const std::string svName_;
0075   const std::string svDoc_;
0076   const double dlenMin_, dlenSigMin_;
0077 };
0078
0079 //
0080 // constructors and destructor
0081 //
0082 VertexTableProducer::VertexTableProducer(const edm::ParameterSet& params)
0083     : pvs_(consumes<std::vector<reco::Vertex>>(params.getParameter<edm::InputTag>("pvSrc"))),
0084       pvsScore_(consumes<edm::ValueMap<float>>(params.getParameter<edm::InputTag>("pvSrc"))),
0085       svs_(consumes<edm::View<reco::VertexCompositePtrCandidate>>(params.getParameter<edm::InputTag>("svSrc"))),
0086       svCut_(params.getParameter<std::string>("svCut"), true),
0087       goodPvCut_(params.getParameter<std::string>("goodPvCut"), true),
0088       goodPvCutString_(params.getParameter<std::string>("goodPvCut")),
0089       pvName_(params.getParameter<std::string>("pvName")),
0090       svName_(params.getParameter<std::string>("svName")),
0091       svDoc_(params.getParameter<std::string>("svDoc")),
0092       dlenMin_(params.getParameter<double>("dlenMin")),
0093       dlenSigMin_(params.getParameter<double>("dlenSigMin"))
0094
0095 {
0096   produces<nanoaod::FlatTable>("pv");
0097   produces<nanoaod::FlatTable>("otherPVs");
0098   produces<nanoaod::FlatTable>("svs");
0099   produces<edm::PtrVector<reco::Candidate>>();
0100 }
0101
0102 VertexTableProducer::~VertexTableProducer() {
0103   // do anything here that needs to be done at destruction time
0104   // (e.g. close files, deallocate resources etc.)
0105 }
0106
0107 //
0108 // member functions
0109 //
0110
0111 // ------------ method called to produce the data  ------------
0112
0113 void VertexTableProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0114   using namespace edm;
0115   const auto& pvsScoreProd = iEvent.get(pvsScore_);
0116   auto pvsIn = iEvent.getHandle(pvs_);
0117
0118   auto pvTable = std::make_unique<nanoaod::FlatTable>(1, pvName_, true);
0119   pvTable->addColumnValue<float>("ndof", (*pvsIn)[0].ndof(), "main primary vertex number of degree of freedom", 8);
0120   pvTable->addColumnValue<float>("x", (*pvsIn)[0].position().x(), "main primary vertex position x coordinate", 10);
0121   pvTable->addColumnValue<float>("y", (*pvsIn)[0].position().y(), "main primary vertex position y coordinate", 10);
0122   pvTable->addColumnValue<float>("z", (*pvsIn)[0].position().z(), "main primary vertex position z coordinate", 16);
0123   pvTable->addColumnValue<float>("chi2", (*pvsIn)[0].normalizedChi2(), "main primary vertex reduced chi2", 8);
0124   int goodPVs = 0;
0125   for (const auto& pv : *pvsIn)
0126     if (goodPvCut_(pv))
0127       goodPVs++;
0128   pvTable->addColumnValue<uint8_t>("npvs", pvsIn->size(), "total number of reconstructed primary vertices");
0129   pvTable->addColumnValue<uint8_t>(
0130       "npvsGood", goodPVs, "number of good reconstructed primary vertices. selection:" + goodPvCutString_);
0131   pvTable->addColumnValue<float>(
0132       "score", pvsScoreProd.get(pvsIn.id(), 0), "main primary vertex score, i.e. sum pt2 of clustered objects", 8);
0133
0134   auto otherPVsTable =
0135       std::make_unique<nanoaod::FlatTable>((*pvsIn).size() > 4 ? 3 : (*pvsIn).size() - 1, "Other" + pvName_, false);
0136   std::vector<float> pvsz;
0137   std::vector<float> pvscores;
0138   for (size_t i = 1; i < (*pvsIn).size() && i < 4; i++) {
0139     pvsz.push_back((*pvsIn)[i].position().z());
0140     pvscores.push_back(pvsScoreProd.get(pvsIn.id(), i));
0141   }
0142   otherPVsTable->addColumn<float>("z", pvsz, "Z position of other primary vertices, excluding the main PV", 8);
0143   otherPVsTable->addColumn<float>("score", pvscores, "scores of other primary vertices, excluding the main PV", 8);
0144
0145   const auto& svsProd = iEvent.get(svs_);
0146   auto selCandSv = std::make_unique<PtrVector<reco::Candidate>>();
0147   std::vector<float> dlen, dlenSig, pAngle, dxy, dxySig;
0148   std::vector<int16_t> charge;
0149   VertexDistance3D vdist;
0150   VertexDistanceXY vdistXY;
0151
0152   size_t i = 0;
0153   const auto& PV0 = pvsIn->front();
0154   for (const auto& sv : svsProd) {
0155     if (svCut_(sv)) {
0156       Measurement1D dl =
0157           vdist.distance(PV0, VertexState(RecoVertex::convertPos(sv.position()), RecoVertex::convertError(sv.error())));
0158       if (dl.value() > dlenMin_ and dl.significance() > dlenSigMin_) {
0159         dlen.push_back(dl.value());
0160         dlenSig.push_back(dl.significance());
0161         edm::Ptr<reco::Candidate> c = svsProd.ptrAt(i);
0162         selCandSv->push_back(c);
0163         double dx = (PV0.x() - sv.vx()), dy = (PV0.y() - sv.vy()), dz = (PV0.z() - sv.vz());
0164         double pdotv = (dx * sv.px() + dy * sv.py() + dz * sv.pz()) / sv.p() / sqrt(dx * dx + dy * dy + dz * dz);
0165         pAngle.push_back(std::acos(pdotv));
0166         Measurement1D d2d = vdistXY.distance(
0167             PV0, VertexState(RecoVertex::convertPos(sv.position()), RecoVertex::convertError(sv.error())));
0168         dxy.push_back(d2d.value());
0169         dxySig.push_back(d2d.significance());
0170
0171         int sum_charge = 0;
0172         for (unsigned int id = 0; id < sv.numberOfDaughters(); ++id) {
0173           const reco::Candidate* daughter = sv.daughter(id);
0174           sum_charge += daughter->charge();
0175         }
0176         charge.push_back(sum_charge);
0177       }
0178     }
0179     i++;
0180   }
0181
0182   auto svsTable = std::make_unique<nanoaod::FlatTable>(selCandSv->size(), svName_, false);
0183   svsTable->setDoc(svDoc_);
0184   // For SV we fill from here only stuff that cannot be created with the SimpleFlatTableProducer
0185   svsTable->addColumn<float>("dlen", dlen, "decay length in cm", 10);
0186   svsTable->addColumn<float>("dlenSig", dlenSig, "decay length significance", 10);
0187   svsTable->addColumn<float>("dxy", dxy, "2D decay length in cm", 10);
0188   svsTable->addColumn<float>("dxySig", dxySig, "2D decay length significance", 10);
0189   svsTable->addColumn<float>("pAngle", pAngle, "pointing angle, i.e. acos(p_SV * (SV - PV)) ", 10);
0190   svsTable->addColumn<int16_t>("charge", charge, "sum of the charge of the SV tracks", 10);
0191
0192   iEvent.put(std::move(pvTable), "pv");
0193   iEvent.put(std::move(otherPVsTable), "otherPVs");
0194   iEvent.put(std::move(svsTable), "svs");
0195   iEvent.put(std::move(selCandSv));
0196 }
0197
0198 // ------------ method called once each stream before processing any runs, lumis or events  ------------
0199 void VertexTableProducer::beginStream(edm::StreamID) {}
0200
0201 // ------------ method called once each stream after processing all runs, lumis and events  ------------
0202 void VertexTableProducer::endStream() {}
0203
0204 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0205 void VertexTableProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0206   edm::ParameterSetDescription desc;
0207
0208   desc.add<edm::InputTag>("pvSrc")->setComment(
0209       "std::vector<reco::Vertex> and ValueMap<float> primary vertex input collections");
0210   desc.add<std::string>("goodPvCut")->setComment("selection on the primary vertex");
0211   desc.add<edm::InputTag>("svSrc")->setComment(
0212       "reco::VertexCompositePtrCandidate compatible secondary vertex input collection");
0213   desc.add<std::string>("svCut")->setComment("selection on the secondary vertex");
0214
0215   desc.add<double>("dlenMin")->setComment("minimum value of dl to select secondary vertex");
0216   desc.add<double>("dlenSigMin")->setComment("minimum value of dl significance to select secondary vertex");
0217
0218   desc.add<std::string>("pvName")->setComment("name of the flat table ouput");
0219   desc.add<std::string>("svName")->setComment("name of the flat table ouput");
0220   desc.add<std::string>("svDoc")->setComment("a few words of documentation");
0221
0222   descriptions.addWithDefaultLabel(desc);
0223 }
0224
0225 //define this as a plug-in
0226 DEFINE_FWK_MODULE(VertexTableProducer);