File indexing completed on 2024-04-06 12:26:48
0001 #include "RecoMET/METPUSubtraction/plugins/JVFJetIdProducer.h"
0002
0003 #include "FWCore/Utilities/interface/Exception.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include "DataFormats/JetReco/interface/PFJet.h"
0007 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0008 #include "DataFormats/Common/interface/ValueMap.h"
0009 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0010 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0011 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0012 #include "DataFormats/VertexReco/interface/Vertex.h"
0013 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0014 #include "DataFormats/TrackReco/interface/Track.h"
0015 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0016 #include "DataFormats/Common/interface/Handle.h"
0017 #include "DataFormats/Math/interface/deltaR.h"
0018
0019 enum { kNeutralJetPU, kNeutralJetNoPU };
0020
0021 JVFJetIdProducer::JVFJetIdProducer(const edm::ParameterSet& cfg) {
0022 srcJets_ = consumes<reco::PFJetCollection>(cfg.getParameter<edm::InputTag>("srcJets"));
0023
0024 srcPFCandidates_ = consumes<reco::PFCandidateCollection>(cfg.getParameter<edm::InputTag>("srcPFCandidates"));
0025 srcPFCandToVertexAssociations_ =
0026 consumes<PFCandToVertexAssMap>(cfg.getParameter<edm::InputTag>("srcPFCandToVertexAssociations"));
0027 srcHardScatterVertex_ = consumes<reco::VertexCollection>(cfg.getParameter<edm::InputTag>("srcHardScatterVertex"));
0028 minTrackPt_ = cfg.getParameter<double>("minTrackPt");
0029 dZcut_ = cfg.getParameter<double>("dZcut");
0030
0031 JVFcut_ = cfg.getParameter<double>("JVFcut");
0032
0033 std::string neutralJetOption_string = cfg.getParameter<std::string>("neutralJetOption");
0034 if (neutralJetOption_string == "PU")
0035 neutralJetOption_ = kNeutralJetPU;
0036 else if (neutralJetOption_string == "noPU")
0037 neutralJetOption_ = kNeutralJetNoPU;
0038 else
0039 throw cms::Exception("JVFJetIdProducer")
0040 << "Invalid Configuration Parameter 'neutralJetOption' = " << neutralJetOption_string << " !!\n";
0041
0042 verbosity_ = (cfg.exists("verbosity")) ? cfg.getParameter<int>("verbosity") : 0;
0043
0044 produces<edm::ValueMap<double>>("Discriminant");
0045 produces<edm::ValueMap<int>>("Id");
0046 }
0047
0048 JVFJetIdProducer::~JVFJetIdProducer() {
0049
0050 }
0051
0052 namespace {
0053 double computeJVF(const reco::PFJet& jet,
0054 const PFCandToVertexAssMap& pfCandToVertexAssociations,
0055 const reco::VertexCollection& vertices,
0056 double dZ,
0057 double minTrackPt,
0058 int verbosity) {
0059 LogDebug("computeJVF") << "<computeJVF>:" << std::endl
0060 << " jet: Pt = " << jet.pt() << ", eta = " << jet.eta() << ", phi = " << jet.phi()
0061 << std::endl;
0062
0063 double trackSum_isVtxAssociated = 0.;
0064 double trackSum_isNotVtxAssociated = 0.;
0065
0066 std::vector<reco::PFCandidatePtr> pfConsts = jet.getPFConstituents();
0067 for (std::vector<reco::PFCandidatePtr>::const_iterator jetConstituent = pfConsts.begin();
0068 jetConstituent != pfConsts.end();
0069 ++jetConstituent) {
0070 if ((*jetConstituent)->charge() == 0)
0071 continue;
0072
0073 double trackPt = 0.;
0074 if ((*jetConstituent)->gsfTrackRef().isNonnull() && (*jetConstituent)->gsfTrackRef().isAvailable())
0075 trackPt = (*jetConstituent)->gsfTrackRef()->pt();
0076 else if ((*jetConstituent)->trackRef().isNonnull() && (*jetConstituent)->trackRef().isAvailable())
0077 trackPt = (*jetConstituent)->trackRef()->pt();
0078 else
0079 trackPt = (*jetConstituent)->pt();
0080
0081 if (trackPt > minTrackPt) {
0082 int jetConstituent_vtxAssociationType =
0083 noPuUtils::isVertexAssociated((*jetConstituent), pfCandToVertexAssociations, vertices, dZ);
0084
0085 bool jetConstituent_isVtxAssociated = (jetConstituent_vtxAssociationType == noPuUtils::kChHSAssoc);
0086 double jetConstituentPt = (*jetConstituent)->pt();
0087 if (jetConstituent_isVtxAssociated) {
0088 LogDebug("computeJVF") << "associated track: Pt = " << (*jetConstituent)->pt()
0089 << ", eta = " << (*jetConstituent)->eta() << ", phi = " << (*jetConstituent)->phi()
0090 << std::endl
0091 << " (vtxAssociationType = " << jetConstituent_vtxAssociationType << ")" << std::endl;
0092
0093 trackSum_isVtxAssociated += jetConstituentPt;
0094 } else {
0095 LogDebug("computeJVF") << "unassociated track: Pt = " << (*jetConstituent)->pt()
0096 << ", eta = " << (*jetConstituent)->eta() << ", phi = " << (*jetConstituent)->phi()
0097 << std::endl
0098 << " (vtxAssociationType = " << jetConstituent_vtxAssociationType << ")" << std::endl;
0099
0100 trackSum_isNotVtxAssociated += jetConstituentPt;
0101 }
0102 }
0103 }
0104
0105 double trackSum = trackSum_isVtxAssociated + trackSum_isNotVtxAssociated;
0106
0107 double jvf = -1.;
0108 if (std::abs(jet.eta()) < 2.5 && trackSum > 5.) {
0109 jvf = trackSum_isVtxAssociated / trackSum;
0110 }
0111
0112 LogDebug("computeJVF") << "trackSum: associated = " << trackSum_isVtxAssociated
0113 << ", unassociated = " << trackSum_isNotVtxAssociated << std::endl
0114 << " --> JVF = " << jvf << std::endl;
0115
0116 return jvf;
0117 }
0118 }
0119
0120 void JVFJetIdProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0121
0122 edm::Handle<reco::PFJetCollection> jets;
0123 evt.getByToken(srcJets_, jets);
0124
0125
0126 edm::Handle<reco::PFCandidateCollection> pfCandidates;
0127 evt.getByToken(srcPFCandidates_, pfCandidates);
0128
0129
0130 edm::Handle<PFCandToVertexAssMap> pfCandToVertexAssociations;
0131 evt.getByToken(srcPFCandToVertexAssociations_, pfCandToVertexAssociations);
0132
0133 edm::Handle<reco::VertexCollection> hardScatterVertex;
0134 evt.getByToken(srcHardScatterVertex_, hardScatterVertex);
0135
0136 std::vector<double> jetIdDiscriminants;
0137 std::vector<int> jetIdFlags;
0138
0139 size_t numJets = jets->size();
0140 for (size_t iJet = 0; iJet < numJets; ++iJet) {
0141 reco::PFJetRef jet(jets, iJet);
0142
0143 double jetJVF = computeJVF(
0144 *jet, *pfCandToVertexAssociations, *hardScatterVertex, dZcut_, minTrackPt_, verbosity_ && jet->pt() > 20.);
0145 jetIdDiscriminants.push_back(jetJVF);
0146
0147 int jetIdFlag = 0;
0148 if (jetJVF > JVFcut_)
0149 jetIdFlag = 255;
0150 else if (jetJVF < -0.5 && neutralJetOption_ == kNeutralJetNoPU)
0151 jetIdFlag = 255;
0152 jetIdFlags.push_back(jetIdFlag);
0153 }
0154
0155 auto jetIdDiscriminants_ptr = std::make_unique<edm::ValueMap<double>>();
0156 edm::ValueMap<double>::Filler jetIdDiscriminantFiller(*jetIdDiscriminants_ptr);
0157 jetIdDiscriminantFiller.insert(jets, jetIdDiscriminants.begin(), jetIdDiscriminants.end());
0158 jetIdDiscriminantFiller.fill();
0159
0160 auto jetIdFlags_ptr = std::make_unique<edm::ValueMap<int>>();
0161 edm::ValueMap<int>::Filler jetIdFlagFiller(*jetIdFlags_ptr);
0162 jetIdFlagFiller.insert(jets, jetIdFlags.begin(), jetIdFlags.end());
0163 jetIdFlagFiller.fill();
0164
0165 evt.put(std::move(jetIdDiscriminants_ptr), "Discriminant");
0166 evt.put(std::move(jetIdFlags_ptr), "Id");
0167 }
0168
0169 #include "FWCore/Framework/interface/MakerMacros.h"
0170
0171 DEFINE_FWK_MODULE(JVFJetIdProducer);