CMS 3D CMS Logo

PFJetFilter.cc
Go to the documentation of this file.
4 
11 
12 using namespace reco;
13 using namespace edm;
14 using namespace std;
15 
17  inputTruthLabel_ = consumes<edm::View<reco::Candidate>>(iConfig.getParameter<edm::InputTag>("InputTruthLabel"));
18  inputRecoLabel_ = consumes<edm::View<reco::Candidate>>(iConfig.getParameter<edm::InputTag>("InputRecoLabel"));
19 
20  recPt_cut = iConfig.getParameter<double>("recPt");
21  genPt_cut = iConfig.getParameter<double>("genPt");
22 
23  eta_min = iConfig.getParameter<double>("minEta");
24  eta_max = iConfig.getParameter<double>("maxEta");
25 
26  deltaR_min = iConfig.getParameter<double>("deltaRMin");
27  deltaR_max = iConfig.getParameter<double>("deltaRMax");
28 
29  deltaEt_min = iConfig.getParameter<double>("minDeltaEt");
30  deltaEt_max = iConfig.getParameter<double>("maxDeltaEt");
31 
32  verbose = iConfig.getParameter<bool>("verbose");
33 
34  entry = 0;
35 }
36 
38 
40 
42 
44  // Typedefs to use views
45  typedef edm::View<reco::Candidate> candidateCollection;
46  typedef edm::View<reco::Candidate> candidateCollection;
47 
48  const candidateCollection *truth_candidates;
49  const candidateCollection *reco_candidates;
50 
51  // ==========================================================
52  // Retrieve!
53  // ==========================================================
54 
55  // Get Truth Candidates (GenCandidates, GenJets, etc.)
57  bool isGen = iEvent.getByToken(inputTruthLabel_, truth_hnd);
58  if (!isGen) {
59  std::cout << "Warning : no Gen jets in input !" << std::endl;
60  return false;
61  }
62 
63  truth_candidates = truth_hnd.product();
64 
65  // Get Reco Candidates (PFlow, CaloJet, etc.)
67  bool isReco = iEvent.getByToken(inputRecoLabel_, reco_hnd);
68  if (!isReco) {
69  std::cout << "Warning : no Reco jets in input !" << std::endl;
70  return false;
71  }
72  reco_candidates = reco_hnd.product();
73  if (!truth_candidates || !reco_candidates)
74  return false;
75 
76  bool pass = false;
77 
78  for (unsigned int i = 0; i < reco_candidates->size(); i++) {
79  const reco::Candidate *particle = &(*reco_candidates)[i];
80 
81  // This protection is meant at not being used !
82  assert(particle != nullptr);
83 
84  double rec_pt = particle->pt();
85  double rec_eta = particle->eta();
86  double rec_phi = particle->phi();
87 
88  // skip PFjets with pt < recPt_cut GeV
89  if (rec_pt < recPt_cut)
90  continue;
91 
92  // skip PFjets with eta > maxEta_cut or eta < minEta_cut
93  if (fabs(rec_eta) > eta_max)
94  continue;
95  if (fabs(rec_eta) < eta_min)
96  continue;
97 
98  bool Barrel = false;
99  bool Endcap = false;
100  if (std::abs(rec_eta) < 1.4)
101  Barrel = true;
102  if (std::abs(rec_eta) > 1.4 && std::abs(rec_eta) < 2.6)
103  Endcap = true;
104 
105  // Keep only jets in the barrel or the endcaps, within the tracker
106  // acceptance
107  if (!Barrel && !Endcap)
108  continue;
109 
110  // Find the closets recJet
111  double deltaRmin = 999.;
112  double ptmin = 0.;
113  for (unsigned int j = 0; j < reco_candidates->size(); j++) {
114  if (i == j)
115  continue;
116  const reco::Candidate *other = &(*reco_candidates)[j];
117  double deltaR = algo_->deltaR(particle, other);
118  if (deltaR < deltaRmin && other->pt() > 0.25 * particle->pt() && other->pt() > recPt_cut) {
119  deltaRmin = deltaR;
120  ptmin = other->pt();
121  }
122  if (deltaRmin < deltaR_min)
123  break;
124  }
125  if (deltaRmin < deltaR_min)
126  continue;
127 
128  // Find the closest genJet.
129  const reco::Candidate *gen_particle = algo_->matchByDeltaR(particle, truth_candidates);
130 
131  // Check there is a genJet associated to the recoJet
132  if (gen_particle == nullptr)
133  continue;
134 
135  // check deltaR is small enough
136  double deltaR = algo_->deltaR(particle, gen_particle);
137  if (deltaR > deltaR_max)
138  continue;
139 
140  // double true_E = gen_particle->p();
141  double true_pt = gen_particle->pt();
142  double true_eta = gen_particle->eta();
143  double true_phi = gen_particle->phi();
144 
145  // skip PFjets with pt < genPt_cut GeV
146  if (true_pt < genPt_cut)
147  continue;
148 
149  // Find the pT residual
150  double resPt = (rec_pt - true_pt) / true_pt;
151  double sigma = resolution(true_pt, Barrel);
152  double avera = response(true_pt, Barrel);
153  double nSig = (resPt - avera) / sigma;
154 
155  if (nSig > deltaEt_max || nSig < deltaEt_min) {
156  /* */
157  if (verbose)
158  std::cout << "Entry " << entry << " resPt = " << resPt << " sigma/avera/nSig = " << sigma << "/" << avera << "/"
159  << nSig << " pT (T/R) " << true_pt << "/" << rec_pt << " Eta (T/R) " << true_eta << "/" << rec_eta
160  << " Phi (T/R) " << true_phi << "/" << rec_phi << " deltaRMin/ptmin " << deltaRmin << "/" << ptmin
161  << std::endl;
162  /* */
163  pass = true;
164  }
165 
166  if (pass)
167  break;
168  }
169 
170  entry++;
171  return pass;
172 }
173 
174 double PFJetFilter::resolution(double pt, bool barrel) {
175  double p0 = barrel ? 1.19200e-02 : 8.45341e-03;
176  double p1 = barrel ? 1.06138e+00 : 7.96855e-01;
177  double p2 = barrel ? -2.05929e+00 : -3.12076e-01;
178 
179  double resp = p0 + p1 / sqrt(pt) + p2 / pt;
180  return resp;
181 }
182 
183 double PFJetFilter::response(double pt, bool barrel) {
184  double p0 = barrel ? 1.09906E-1 : 6.91939E+1;
185  double p1 = barrel ? -1.61443E-1 : -6.92733E+1;
186  double p2 = barrel ? 3.45489E+3 : 1.58207E+6;
187 
188  double resp = p0 + p1 * std::exp(-pt / p2);
189  return resp;
190 }
T getParameter(std::string const &) const
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:517
void endJob() override
Definition: PFJetFilter.cc:41
~PFJetFilter() override
Definition: PFJetFilter.cc:37
double response(double, bool)
Definition: PFJetFilter.cc:183
int iEvent
Definition: GenABIO.cc:224
T sqrt(T t)
Definition: SSEVec.h:18
deltaRmin
maximum deltaR of track to jet.
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
PFJetFilter(const edm::ParameterSet &)
Definition: PFJetFilter.cc:16
double p2[4]
Definition: TauolaWrapper.h:90
constexpr auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:28
virtual double eta() const =0
momentum pseudorapidity
T const * product() const
Definition: Handle.h:74
virtual double pt() const =0
transverse momentum
double ptmin
Definition: HydjetWrapper.h:90
fixed size matrix
HLT enums.
double p1[4]
Definition: TauolaWrapper.h:89
bool filter(edm::Event &, const edm::EventSetup &) override
Definition: PFJetFilter.cc:43
void beginJob() override
Definition: PFJetFilter.cc:39
virtual double phi() const =0
momentum azimuthal angle
double resolution(double, bool)
Definition: PFJetFilter.cc:174