Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-11-23 23:27:32

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   edm::Handle<edm::ValueMap<float>> pvsScoreIn;
0116   edm::Handle<std::vector<reco::Vertex>> pvsIn;
0117   iEvent.getByToken(pvs_, pvsIn);
0118   iEvent.getByToken(pvsScore_, pvsScoreIn);
0119   auto pvTable = std::make_unique<nanoaod::FlatTable>(1, pvName_, true);
0120   pvTable->addColumnValue<float>("ndof", (*pvsIn)[0].ndof(), "main primary vertex number of degree of freedom", 8);
0121   pvTable->addColumnValue<float>("x", (*pvsIn)[0].position().x(), "main primary vertex position x coordinate", 10);
0122   pvTable->addColumnValue<float>("y", (*pvsIn)[0].position().y(), "main primary vertex position y coordinate", 10);
0123   pvTable->addColumnValue<float>("z", (*pvsIn)[0].position().z(), "main primary vertex position z coordinate", 16);
0124   pvTable->addColumnValue<float>("chi2", (*pvsIn)[0].normalizedChi2(), "main primary vertex reduced chi2", 8);
0125   int goodPVs = 0;
0126   for (const auto& pv : *pvsIn)
0127     if (goodPvCut_(pv))
0128       goodPVs++;
0129   pvTable->addColumnValue<int>("npvs", pvsIn->size(), "total number of reconstructed primary vertices");
0130   pvTable->addColumnValue<int>(
0131       "npvsGood", goodPVs, "number of good reconstructed primary vertices. selection:" + goodPvCutString_);
0132   pvTable->addColumnValue<float>(
0133       "score", pvsScoreIn->get(pvsIn.id(), 0), "main primary vertex score, i.e. sum pt2 of clustered objects", 8);
0134 
0135   auto otherPVsTable =
0136       std::make_unique<nanoaod::FlatTable>((*pvsIn).size() > 4 ? 3 : (*pvsIn).size() - 1, "Other" + pvName_, false);
0137   std::vector<float> pvsz;
0138   std::vector<float> pvscores;
0139   for (size_t i = 1; i < (*pvsIn).size() && i < 4; i++) {
0140     pvsz.push_back((*pvsIn)[i].position().z());
0141     pvscores.push_back((*pvsScoreIn).get(pvsIn.id(), i));
0142   }
0143   otherPVsTable->addColumn<float>("z", pvsz, "Z position of other primary vertices, excluding the main PV", 8);
0144   otherPVsTable->addColumn<float>("score", pvscores, "scores of other primary vertices, excluding the main PV", 8);
0145 
0146   edm::Handle<edm::View<reco::VertexCompositePtrCandidate>> svsIn;
0147   iEvent.getByToken(svs_, svsIn);
0148   auto selCandSv = std::make_unique<PtrVector<reco::Candidate>>();
0149   std::vector<float> dlen, dlenSig, pAngle, dxy, dxySig;
0150   std::vector<int> charge;
0151   VertexDistance3D vdist;
0152   VertexDistanceXY vdistXY;
0153 
0154   size_t i = 0;
0155   const auto& PV0 = pvsIn->front();
0156   for (const auto& sv : *svsIn) {
0157     if (svCut_(sv)) {
0158       Measurement1D dl =
0159           vdist.distance(PV0, VertexState(RecoVertex::convertPos(sv.position()), RecoVertex::convertError(sv.error())));
0160       if (dl.value() > dlenMin_ and dl.significance() > dlenSigMin_) {
0161         dlen.push_back(dl.value());
0162         dlenSig.push_back(dl.significance());
0163         edm::Ptr<reco::Candidate> c = svsIn->ptrAt(i);
0164         selCandSv->push_back(c);
0165         double dx = (PV0.x() - sv.vx()), dy = (PV0.y() - sv.vy()), dz = (PV0.z() - sv.vz());
0166         double pdotv = (dx * sv.px() + dy * sv.py() + dz * sv.pz()) / sv.p() / sqrt(dx * dx + dy * dy + dz * dz);
0167         pAngle.push_back(std::acos(pdotv));
0168         Measurement1D d2d = vdistXY.distance(
0169             PV0, VertexState(RecoVertex::convertPos(sv.position()), RecoVertex::convertError(sv.error())));
0170         dxy.push_back(d2d.value());
0171         dxySig.push_back(d2d.significance());
0172 
0173         int sum_charge = 0;
0174         for (unsigned int id = 0; id < sv.numberOfDaughters(); ++id) {
0175           const reco::Candidate* daughter = sv.daughter(id);
0176           sum_charge += daughter->charge();
0177         }
0178         charge.push_back(sum_charge);
0179       }
0180     }
0181     i++;
0182   }
0183 
0184   auto svsTable = std::make_unique<nanoaod::FlatTable>(selCandSv->size(), svName_, false);
0185   // For SV we fill from here only stuff that cannot be created with the SimpleFlatTableProducer
0186   svsTable->addColumn<float>("dlen", dlen, "decay length in cm", 10);
0187   svsTable->addColumn<float>("dlenSig", dlenSig, "decay length significance", 10);
0188   svsTable->addColumn<float>("dxy", dxy, "2D decay length in cm", 10);
0189   svsTable->addColumn<float>("dxySig", dxySig, "2D decay length significance", 10);
0190   svsTable->addColumn<float>("pAngle", pAngle, "pointing angle, i.e. acos(p_SV * (SV - PV)) ", 10);
0191   svsTable->addColumn<int>("charge", charge, "sum of the charge of the SV tracks", 10);
0192 
0193   iEvent.put(std::move(pvTable), "pv");
0194   iEvent.put(std::move(otherPVsTable), "otherPVs");
0195   iEvent.put(std::move(svsTable), "svs");
0196   iEvent.put(std::move(selCandSv));
0197 }
0198 
0199 // ------------ method called once each stream before processing any runs, lumis or events  ------------
0200 void VertexTableProducer::beginStream(edm::StreamID) {}
0201 
0202 // ------------ method called once each stream after processing all runs, lumis and events  ------------
0203 void VertexTableProducer::endStream() {}
0204 
0205 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0206 void VertexTableProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0207   //The following says we do not know what parameters are allowed so do no validation
0208   // Please change this to state exactly what you do use, even if it is no parameters
0209   edm::ParameterSetDescription desc;
0210   desc.setUnknown();
0211   descriptions.addDefault(desc);
0212 }
0213 
0214 //define this as a plug-in
0215 DEFINE_FWK_MODULE(VertexTableProducer);