RecoParticleFlow/plugins/PFJetFilter.cc

0001 #include "DataFormats/Candidate/interface/Candidate.h"
0002 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0003 #include "DataFormats/Common/interface/Handle.h"
0004 #include "DataFormats/JetReco/interface/PFJet.h"
0005 #include "FWCore/Framework/interface/global/EDFilter.h"
0006 #include "FWCore/Framework/interface/ESHandle.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/ServiceRegistry/interface/Service.h"
0012 #include "FWCore/Utilities/interface/InputTag.h"
0013 #include "RecoParticleFlow/Benchmark/interface/PFBenchmarkAlgo.h"
0014 #include <atomic>
0015
0016 class PFJetFilter : public edm::global::EDFilter<> {
0017 public:
0018   explicit PFJetFilter(const edm::ParameterSet &);
0019   ~PFJetFilter() override;
0020
0021 private:
0022   void beginJob() override;
0023   bool filter(edm::StreamID, edm::Event &, const edm::EventSetup &) const override;
0024   void endJob() override;
0025
0026   double resolution(double, bool) const;
0027   double response(double, bool) const;
0028
0029   double recPt_cut;
0030   double genPt_cut;
0031   double deltaEt_min;
0032   double deltaEt_max;
0033   double deltaR_min;
0034   double deltaR_max;
0035   double eta_min;
0036   double eta_max;
0037   edm::EDGetTokenT<edm::View<reco::Candidate>> inputTruthLabel_;
0038   edm::EDGetTokenT<edm::View<reco::Candidate>> inputRecoLabel_;
0039
0040   mutable std::atomic<unsigned int> entry;
0041   bool verbose;
0042 };
0043
0044 #include "FWCore/Framework/interface/MakerMacros.h"
0045 DEFINE_FWK_MODULE(PFJetFilter);
0046
0047 using namespace reco;
0048 using namespace edm;
0049 using namespace std;
0050
0051 PFJetFilter::PFJetFilter(const edm::ParameterSet &iConfig) {
0052   inputTruthLabel_ = consumes<edm::View<reco::Candidate>>(iConfig.getParameter<edm::InputTag>("InputTruthLabel"));
0053   inputRecoLabel_ = consumes<edm::View<reco::Candidate>>(iConfig.getParameter<edm::InputTag>("InputRecoLabel"));
0054
0055   recPt_cut = iConfig.getParameter<double>("recPt");
0056   genPt_cut = iConfig.getParameter<double>("genPt");
0057
0058   eta_min = iConfig.getParameter<double>("minEta");
0059   eta_max = iConfig.getParameter<double>("maxEta");
0060
0061   deltaR_min = iConfig.getParameter<double>("deltaRMin");
0062   deltaR_max = iConfig.getParameter<double>("deltaRMax");
0063
0064   deltaEt_min = iConfig.getParameter<double>("minDeltaEt");
0065   deltaEt_max = iConfig.getParameter<double>("maxDeltaEt");
0066
0067   verbose = iConfig.getParameter<bool>("verbose");
0068
0069   entry = 0;
0070 }
0071
0072 PFJetFilter::~PFJetFilter() {}
0073
0074 void PFJetFilter::beginJob() {}
0075
0076 void PFJetFilter::endJob() {}
0077
0078 bool PFJetFilter::filter(edm::StreamID, edm::Event &iEvent, const edm::EventSetup &iESetup) const {
0079   // Typedefs to use views
0080   typedef edm::View<reco::Candidate> candidateCollection;
0081   typedef edm::View<reco::Candidate> candidateCollection;
0082
0083   const candidateCollection *truth_candidates;
0084   const candidateCollection *reco_candidates;
0085
0086   // ==========================================================
0087   // Retrieve!
0088   // ==========================================================
0089
0090   // Get Truth Candidates (GenCandidates, GenJets, etc.)
0091   Handle<candidateCollection> truth_hnd;
0092   bool isGen = iEvent.getByToken(inputTruthLabel_, truth_hnd);
0093   if (!isGen) {
0094     std::cout << "Warning : no Gen jets in input !" << std::endl;
0095     return false;
0096   }
0097
0098   truth_candidates = truth_hnd.product();
0099
0100   // Get Reco Candidates (PFlow, CaloJet, etc.)
0101   Handle<candidateCollection> reco_hnd;
0102   bool isReco = iEvent.getByToken(inputRecoLabel_, reco_hnd);
0103   if (!isReco) {
0104     std::cout << "Warning : no Reco jets in input !" << std::endl;
0105     return false;
0106   }
0107   reco_candidates = reco_hnd.product();
0108   if (!truth_candidates || !reco_candidates)
0109     return false;
0110
0111   bool pass = false;
0112
0113   for (unsigned int i = 0; i < reco_candidates->size(); i++) {
0114     const reco::Candidate *particle = &(*reco_candidates)[i];
0115
0116     // This protection is meant at not being used !
0117     assert(particle != nullptr);
0118
0119     double rec_pt = particle->pt();
0120     double rec_eta = particle->eta();
0121     double rec_phi = particle->phi();
0122
0123     // skip PFjets with pt < recPt_cut GeV
0124     if (rec_pt < recPt_cut)
0125       continue;
0126
0127     // skip PFjets with eta > maxEta_cut or eta < minEta_cut
0128     if (fabs(rec_eta) > eta_max)
0129       continue;
0130     if (fabs(rec_eta) < eta_min)
0131       continue;
0132
0133     bool Barrel = false;
0134     bool Endcap = false;
0135     if (std::abs(rec_eta) < 1.4)
0136       Barrel = true;
0137     if (std::abs(rec_eta) > 1.4 && std::abs(rec_eta) < 2.6)
0138       Endcap = true;
0139
0140     // Keep only jets in the barrel or the endcaps, within the tracker
0141     // acceptance
0142     if (!Barrel && !Endcap)
0143       continue;
0144
0145     // Find the closets recJet
0146     double deltaRmin = 999.;
0147     double ptmin = 0.;
0148     for (unsigned int j = 0; j < reco_candidates->size(); j++) {
0149       if (i == j)
0150         continue;
0151       const reco::Candidate *other = &(*reco_candidates)[j];
0152       double deltaR = PFBenchmarkAlgo::deltaR(particle, other);
0153       if (deltaR < deltaRmin && other->pt() > 0.25 * particle->pt() && other->pt() > recPt_cut) {
0154         deltaRmin = deltaR;
0155         ptmin = other->pt();
0156       }
0157       if (deltaRmin < deltaR_min)
0158         break;
0159     }
0160     if (deltaRmin < deltaR_min)
0161       continue;
0162
0163     // Find the closest genJet.
0164     const reco::Candidate *gen_particle = PFBenchmarkAlgo::matchByDeltaR(particle, truth_candidates);
0165
0166     // Check there is a genJet associated to the recoJet
0167     if (gen_particle == nullptr)
0168       continue;
0169
0170     // check deltaR is small enough
0171     double deltaR = PFBenchmarkAlgo::deltaR(particle, gen_particle);
0172     if (deltaR > deltaR_max)
0173       continue;
0174
0175     // double true_E = gen_particle->p();
0176     double true_pt = gen_particle->pt();
0177     double true_eta = gen_particle->eta();
0178     double true_phi = gen_particle->phi();
0179
0180     // skip PFjets with pt < genPt_cut GeV
0181     if (true_pt < genPt_cut)
0182       continue;
0183
0184     // Find the pT residual
0185     double resPt = (rec_pt - true_pt) / true_pt;
0186     double sigma = resolution(true_pt, Barrel);
0187     double avera = response(true_pt, Barrel);
0188     double nSig = (resPt - avera) / sigma;
0189
0190     if (nSig > deltaEt_max || nSig < deltaEt_min) {
0191       /* */
0192       if (verbose)
0193         std::cout << "Entry " << entry++ << " resPt = " << resPt << " sigma/avera/nSig = " << sigma << "/" << avera
0194                   << "/" << nSig << " pT (T/R) " << true_pt << "/" << rec_pt << " Eta (T/R) " << true_eta << "/"
0195                   << rec_eta << " Phi (T/R) " << true_phi << "/" << rec_phi << " deltaRMin/ptmin " << deltaRmin << "/"
0196                   << ptmin << std::endl;
0197       /* */
0198       pass = true;
0199     }
0200
0201     if (pass)
0202       break;
0203   }
0204
0205   return pass;
0206 }
0207
0208 double PFJetFilter::resolution(double pt, bool barrel) const {
0209   double p0 = barrel ? 1.19200e-02 : 8.45341e-03;
0210   double p1 = barrel ? 1.06138e+00 : 7.96855e-01;
0211   double p2 = barrel ? -2.05929e+00 : -3.12076e-01;
0212
0213   double resp = p0 + p1 / sqrt(pt) + p2 / pt;
0214   return resp;
0215 }
0216
0217 double PFJetFilter::response(double pt, bool barrel) const {
0218   double p0 = barrel ? 1.09906E-1 : 6.91939E+1;
0219   double p1 = barrel ? -1.61443E-1 : -6.92733E+1;
0220   double p2 = barrel ? 3.45489E+3 : 1.58207E+6;
0221
0222   double resp = p0 + p1 * std::exp(-pt / p2);
0223   return resp;
0224 }