METPUSubtraction/plugins/JVFJetIdProducer.cc

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   // nothing to be done yet...
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         //isVertexAssociated_fast(pfCandidateRef, pfCandToVertexAssociations_reversed, *hardScatterVertex, dZcut_, numWarnings_, maxWarnings_);
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 }  // namespace
0119
0120 void JVFJetIdProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0121   // get jets
0122   edm::Handle<reco::PFJetCollection> jets;
0123   evt.getByToken(srcJets_, jets);
0124
0125   // get PFCandidates
0126   edm::Handle<reco::PFCandidateCollection> pfCandidates;
0127   evt.getByToken(srcPFCandidates_, pfCandidates);
0128
0129   // get PFCandidate-to-vertex associations and "the" hard-scatter vertex
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);