File indexing completed on 2024-04-25 02:14:03
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021
0022
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 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0045 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0046
0047
0048
0049
0050
0051 class VertexTableProducer : public edm::stream::EDProducer<> {
0052 public:
0053 explicit VertexTableProducer(const edm::ParameterSet&);
0054
0055 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0056
0057 private:
0058 void produce(edm::Event&, const edm::EventSetup&) override;
0059
0060
0061
0062 const edm::EDGetTokenT<std::vector<reco::Vertex>> pvs_;
0063 const edm::EDGetTokenT<pat::PackedCandidateCollection> pfc_;
0064 const edm::EDGetTokenT<edm::ValueMap<float>> pvsScore_;
0065 const edm::EDGetTokenT<edm::View<reco::VertexCompositePtrCandidate>> svs_;
0066 const StringCutObjectSelector<reco::Candidate> svCut_;
0067 const StringCutObjectSelector<reco::Vertex> goodPvCut_;
0068 const std::string goodPvCutString_;
0069 const std::string pvName_;
0070 const std::string svName_;
0071 const std::string svDoc_;
0072 const double dlenMin_, dlenSigMin_;
0073 };
0074
0075
0076
0077
0078 VertexTableProducer::VertexTableProducer(const edm::ParameterSet& params)
0079 : pvs_(consumes<std::vector<reco::Vertex>>(params.getParameter<edm::InputTag>("pvSrc"))),
0080 pfc_(consumes<pat::PackedCandidateCollection>(params.getParameter<edm::InputTag>("pfcSrc"))),
0081 pvsScore_(consumes<edm::ValueMap<float>>(params.getParameter<edm::InputTag>("pvSrc"))),
0082 svs_(consumes<edm::View<reco::VertexCompositePtrCandidate>>(params.getParameter<edm::InputTag>("svSrc"))),
0083 svCut_(params.getParameter<std::string>("svCut"), true),
0084 goodPvCut_(params.getParameter<std::string>("goodPvCut"), true),
0085 goodPvCutString_(params.getParameter<std::string>("goodPvCut")),
0086 pvName_(params.getParameter<std::string>("pvName")),
0087 svName_(params.getParameter<std::string>("svName")),
0088 svDoc_(params.getParameter<std::string>("svDoc")),
0089 dlenMin_(params.getParameter<double>("dlenMin")),
0090 dlenSigMin_(params.getParameter<double>("dlenSigMin"))
0091
0092 {
0093 produces<nanoaod::FlatTable>("pv");
0094 produces<nanoaod::FlatTable>("otherPVs");
0095 produces<nanoaod::FlatTable>("svs");
0096 produces<edm::PtrVector<reco::VertexCompositePtrCandidate>>();
0097 }
0098
0099
0100
0101
0102
0103
0104
0105 void VertexTableProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0106 using namespace edm;
0107 const auto& pvsScoreProd = iEvent.get(pvsScore_);
0108 auto pvsIn = iEvent.getHandle(pvs_);
0109
0110
0111 auto pfcIn = iEvent.getHandle(pfc_);
0112
0113 auto pvTable = std::make_unique<nanoaod::FlatTable>(1, pvName_, true);
0114 pvTable->addColumnValue<float>("ndof", (*pvsIn)[0].ndof(), "main primary vertex number of degree of freedom", 8);
0115 pvTable->addColumnValue<float>("x", (*pvsIn)[0].position().x(), "main primary vertex position x coordinate", 10);
0116 pvTable->addColumnValue<float>("y", (*pvsIn)[0].position().y(), "main primary vertex position y coordinate", 10);
0117 pvTable->addColumnValue<float>("z", (*pvsIn)[0].position().z(), "main primary vertex position z coordinate", 16);
0118 pvTable->addColumnValue<float>("chi2", (*pvsIn)[0].normalizedChi2(), "main primary vertex reduced chi2", 8);
0119 int goodPVs = 0;
0120 for (const auto& pv : *pvsIn)
0121 if (goodPvCut_(pv))
0122 goodPVs++;
0123 pvTable->addColumnValue<uint8_t>("npvs", pvsIn->size(), "total number of reconstructed primary vertices");
0124 pvTable->addColumnValue<uint8_t>(
0125 "npvsGood", goodPVs, "number of good reconstructed primary vertices. selection:" + goodPvCutString_);
0126 pvTable->addColumnValue<float>(
0127 "score", pvsScoreProd.get(pvsIn.id(), 0), "main primary vertex score, i.e. sum pt2 of clustered objects", 8);
0128
0129 float pv_sumpt2 = 0.0, pv_sumpx = 0.0, pv_sumpy = 0.0;
0130 for (const auto& obj : *pfcIn) {
0131 if (obj.charge() == 0) {
0132 continue;
0133 }
0134 double dz = fabs(obj.dz((*pvsIn)[0].position()));
0135 bool include_pfc = false;
0136 if (dz < 0.2) {
0137 include_pfc = true;
0138 for (size_t j = 1; j < (*pvsIn).size(); j++) {
0139 double newdz = fabs(obj.dz((*pvsIn)[j].position()));
0140 if (newdz < dz) {
0141 include_pfc = false;
0142 break;
0143 }
0144 }
0145 }
0146 if (include_pfc) {
0147 float pfc_pt = obj.pt();
0148 pv_sumpt2 += pfc_pt * pfc_pt;
0149 pv_sumpx += obj.px();
0150 pv_sumpy += obj.py();
0151 }
0152 }
0153 pvTable->addColumnValue<float>(
0154 "sumpt2", pv_sumpt2, "sum pt2 of pf charged candidates for the main primary vertex", 10);
0155 pvTable->addColumnValue<float>("sumpx", pv_sumpx, "sum px of pf charged candidates for the main primary vertex", 10);
0156 pvTable->addColumnValue<float>("sumpy", pv_sumpy, "sum py of pf charged candidates for the main primary vertex", 10);
0157
0158 auto otherPVsTable =
0159 std::make_unique<nanoaod::FlatTable>((*pvsIn).size() > 4 ? 3 : (*pvsIn).size() - 1, "Other" + pvName_, false);
0160 std::vector<float> pvsz;
0161 std::vector<float> pvscores;
0162 for (size_t i = 1; i < (*pvsIn).size() && i < 4; i++) {
0163 pvsz.push_back((*pvsIn)[i].position().z());
0164 pvscores.push_back(pvsScoreProd.get(pvsIn.id(), i));
0165 }
0166 otherPVsTable->addColumn<float>("z", pvsz, "Z position of other primary vertices, excluding the main PV", 8);
0167 otherPVsTable->addColumn<float>("score", pvscores, "scores of other primary vertices, excluding the main PV", 8);
0168
0169 const auto& svsProd = iEvent.get(svs_);
0170 auto selCandSv = std::make_unique<PtrVector<reco::VertexCompositePtrCandidate>>();
0171 std::vector<float> dlen, dlenSig, pAngle, dxy, dxySig;
0172 std::vector<int16_t> charge;
0173 VertexDistance3D vdist;
0174 VertexDistanceXY vdistXY;
0175
0176 size_t i = 0;
0177 const auto& PV0 = pvsIn->front();
0178 for (const auto& sv : svsProd) {
0179 if (svCut_(sv)) {
0180 Measurement1D dl =
0181 vdist.distance(PV0, VertexState(RecoVertex::convertPos(sv.position()), RecoVertex::convertError(sv.error())));
0182 if (dl.value() > dlenMin_ and dl.significance() > dlenSigMin_) {
0183 dlen.push_back(dl.value());
0184 dlenSig.push_back(dl.significance());
0185 selCandSv->push_back(svsProd.ptrAt(i));
0186 double dx = (PV0.x() - sv.vx()), dy = (PV0.y() - sv.vy()), dz = (PV0.z() - sv.vz());
0187 double pdotv = (dx * sv.px() + dy * sv.py() + dz * sv.pz()) / sv.p() / sqrt(dx * dx + dy * dy + dz * dz);
0188 pAngle.push_back(std::acos(pdotv));
0189 Measurement1D d2d = vdistXY.distance(
0190 PV0, VertexState(RecoVertex::convertPos(sv.position()), RecoVertex::convertError(sv.error())));
0191 dxy.push_back(d2d.value());
0192 dxySig.push_back(d2d.significance());
0193
0194 int sum_charge = 0;
0195 for (unsigned int id = 0; id < sv.numberOfDaughters(); ++id) {
0196 const reco::Candidate* daughter = sv.daughter(id);
0197 sum_charge += daughter->charge();
0198 }
0199 charge.push_back(sum_charge);
0200 }
0201 }
0202 i++;
0203 }
0204
0205 auto svsTable = std::make_unique<nanoaod::FlatTable>(selCandSv->size(), svName_, false);
0206 svsTable->setDoc(svDoc_);
0207
0208 svsTable->addColumn<float>("dlen", dlen, "decay length in cm", 10);
0209 svsTable->addColumn<float>("dlenSig", dlenSig, "decay length significance", 10);
0210 svsTable->addColumn<float>("dxy", dxy, "2D decay length in cm", 10);
0211 svsTable->addColumn<float>("dxySig", dxySig, "2D decay length significance", 10);
0212 svsTable->addColumn<float>("pAngle", pAngle, "pointing angle, i.e. acos(p_SV * (SV - PV)) ", 10);
0213 svsTable->addColumn<int16_t>("charge", charge, "sum of the charge of the SV tracks", 10);
0214
0215 iEvent.put(std::move(pvTable), "pv");
0216 iEvent.put(std::move(otherPVsTable), "otherPVs");
0217 iEvent.put(std::move(svsTable), "svs");
0218 iEvent.put(std::move(selCandSv));
0219 }
0220
0221
0222 void VertexTableProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0223 edm::ParameterSetDescription desc;
0224
0225 desc.add<edm::InputTag>("pvSrc")->setComment(
0226 "std::vector<reco::Vertex> and ValueMap<float> primary vertex input collections");
0227 desc.add<edm::InputTag>("pfcSrc")->setComment("packedPFCandidates input collections");
0228 desc.add<std::string>("goodPvCut")->setComment("selection on the primary vertex");
0229 desc.add<edm::InputTag>("svSrc")->setComment(
0230 "reco::VertexCompositePtrCandidate compatible secondary vertex input collection");
0231 desc.add<std::string>("svCut")->setComment("selection on the secondary vertex");
0232
0233 desc.add<double>("dlenMin")->setComment("minimum value of dl to select secondary vertex");
0234 desc.add<double>("dlenSigMin")->setComment("minimum value of dl significance to select secondary vertex");
0235
0236 desc.add<std::string>("pvName")->setComment("name of the flat table ouput");
0237 desc.add<std::string>("svName")->setComment("name of the flat table ouput");
0238 desc.add<std::string>("svDoc")->setComment("a few words of documentation");
0239
0240 descriptions.addWithDefaultLabel(desc);
0241 }
0242
0243
0244 DEFINE_FWK_MODULE(VertexTableProducer);
