CMS 3D CMS Logo
NoPileUpPFMEtProducer.cc
Go to the documentation of this file.
1 #include "RecoMET/METPUSubtraction/plugins/NoPileUpPFMEtProducer.h"
2 
3 #include "FWCore/Utilities/interface/Exception.h"
4 
5 //#include "DataFormats/METReco/interface/PFMEtSignCovMatrix.h" //never used so far
6 #include "RecoMET/METAlgorithms/interface/significanceAlgo.h"
7 #include "DataFormats/METReco/interface/SigInputObj.h"
8 
9 #include <cmath>
10 
11 const double defaultPFMEtResolutionX = 10.;
12 const double defaultPFMEtResolutionY = 10.;
13 
14 const double epsilon = 1.e-9;
15 
16 NoPileUpPFMEtProducer::NoPileUpPFMEtProducer(const edm::ParameterSet& cfg)
17  : moduleLabel_(cfg.getParameter<std::string>("@module_label")) {
18  srcMEt_ = consumes<reco::PFMETCollection>(cfg.getParameter<edm::InputTag>("srcMEt"));
19  srcMEtCov_ = edm::InputTag(); //MM, disabled for the moment until we really need it
20  //( cfg.exists("srcMEtCov") ) ?
21  // consumes<edm::Handle<> >(cfg.getParameter<edm::InputTag>("srcMEtCov")) : edm::InputTag();
22  srcJetInfo_ = consumes<reco::PUSubMETCandInfoCollection>(cfg.getParameter<edm::InputTag>("srcPUSubMETDataJet"));
23  srcJetInfoLeptonMatch_ =
24  consumes<reco::PUSubMETCandInfoCollection>(cfg.getParameter<edm::InputTag>("srcPUSubMETDataJetLeptonMatch"));
25  srcPFCandInfo_ =
26  consumes<reco::PUSubMETCandInfoCollection>(cfg.getParameter<edm::InputTag>("srcPUSubMETDataPFCands"));
27  srcPFCandInfoLeptonMatch_ =
28  consumes<reco::PUSubMETCandInfoCollection>(cfg.getParameter<edm::InputTag>("srcPUSubMETDataPFCandsLeptonMatch"));
29  vInputTag srcLeptonsTags = cfg.getParameter<vInputTag>("srcLeptons");
30  for (vInputTag::const_iterator it = srcLeptonsTags.begin(); it != srcLeptonsTags.end(); it++) {
31  srcLeptons_.push_back(consumes<edm::View<reco::Candidate> >(*it));
32  }
33 
34  srcType0Correction_ = consumes<CorrMETData>(cfg.getParameter<edm::InputTag>("srcType0Correction"));
35 
36  sfNoPUjets_ = cfg.getParameter<double>("sfNoPUjets");
37  sfNoPUjetOffsetEnCorr_ = cfg.getParameter<double>("sfNoPUjetOffsetEnCorr");
38  sfPUjets_ = cfg.getParameter<double>("sfPUjets");
39  sfNoPUunclChargedCands_ = cfg.getParameter<double>("sfNoPUunclChargedCands");
40  sfPUunclChargedCands_ = cfg.getParameter<double>("sfPUunclChargedCands");
41  sfUnclNeutralCands_ = cfg.getParameter<double>("sfUnclNeutralCands");
42  sfType0Correction_ = cfg.getParameter<double>("sfType0Correction");
43  sfLeptonIsoCones_ = cfg.getParameter<double>("sfLeptonIsoCones");
44 
45  pfMEtSignInterface_ = new PFMEtSignInterfaceBase(cfg.getParameter<edm::ParameterSet>("resolution"));
46  sfMEtCovMin_ = cfg.getParameter<double>("sfMEtCovMin");
47  sfMEtCovMax_ = cfg.getParameter<double>("sfMEtCovMax");
48 
49  saveInputs_ = (cfg.exists("saveInputs")) ? cfg.getParameter<bool>("saveInputs") : false;
50 
51  verbosity_ = (cfg.exists("verbosity")) ? cfg.getParameter<int>("verbosity") : 0;
52 
53  produces<reco::PFMETCollection>();
54 
55  sfLeptonsName_ = "sumLeptons";
56  sfNoPUjetsName_ = "sumNoPUjets";
57  sfNoPUjetOffsetEnCorrName_ = "sumNoPUjetOffsetEnCorr";
58  sfPUjetsName_ = "sumPUjets";
59  sfNoPUunclChargedCandsName_ = "sumNoPUunclChargedCands";
60  sfPUunclChargedCandsName_ = "sumPUunclChargedCands";
61  sfUnclNeutralCandsName_ = "sumUnclNeutralCands";
62  sfType0CorrectionName_ = "type0Correction";
63  sfLeptonIsoConesName_ = "sumLeptonIsoCones";
64 
65  if (saveInputs_) {
66  produces<CommonMETData>(sfLeptonsName_);
67  produces<CommonMETData>(sfNoPUjetsName_);
68  produces<CommonMETData>(sfNoPUjetOffsetEnCorrName_);
69  produces<CommonMETData>(sfPUjetsName_);
70  produces<CommonMETData>(sfNoPUunclChargedCandsName_);
71  produces<CommonMETData>(sfPUunclChargedCandsName_);
72  produces<CommonMETData>(sfUnclNeutralCandsName_);
73  produces<CommonMETData>(sfType0CorrectionName_);
74  produces<CommonMETData>(sfLeptonIsoConesName_);
75  }
76  produces<double>("sfNoPU");
77 }
78 
79 NoPileUpPFMEtProducer::~NoPileUpPFMEtProducer() { delete pfMEtSignInterface_; }
80 
81 void initializeCommonMETData(CommonMETData& metData) {
82  metData.met = 0.;
83  metData.mex = 0.;
84  metData.mey = 0.;
85  metData.mez = 0.;
86  metData.sumet = 0.;
87  metData.phi = 0.;
88 }
89 
90 void addToCommonMETData(CommonMETData& metData, const reco::Candidate::LorentzVector& p4) {
91  metData.mex += p4.px();
92  metData.mey += p4.py();
93  metData.mez += p4.pz();
94  metData.sumet += p4.pt();
95 }
96 
97 void finalizeCommonMETData(CommonMETData& metData) {
98  metData.met = sqrt(metData.mex * metData.mex + metData.mey * metData.mey);
99  metData.phi = atan2(metData.mey, metData.mex);
100 }
101 
102 int findBestMatchingLepton(const std::vector<reco::Candidate::LorentzVector>& leptons,
103  const reco::Candidate::LorentzVector& p4_ref) {
104  int leptonIdx_dR2min = -1;
105  double dR2min = 1.e+3;
106  int leptonIdx = 0;
107  for (std::vector<reco::Candidate::LorentzVector>::const_iterator lepton = leptons.begin(); lepton != leptons.end();
108  ++lepton) {
109  double dR2 = deltaR2(*lepton, p4_ref);
110  if (leptonIdx_dR2min == -1 || dR2 < dR2min) {
111  leptonIdx_dR2min = leptonIdx;
112  dR2min = dR2;
113  }
114  ++leptonIdx;
115  }
116  assert(leptonIdx_dR2min >= 0 && leptonIdx_dR2min < (int)leptons.size());
117  return leptonIdx_dR2min;
118 }
119 
120 void scaleAndAddPFMEtSignObjects(std::vector<metsig::SigInputObj>& metSignObjects_scaled,
121  const std::vector<metsig::SigInputObj>& metSignObjects,
122  double sf,
123  double sfMin,
124  double sfMax) {
125  double sf_value = sf;
126  if (sf_value > sfMax)
127  sf_value = sfMax;
128  if (sf_value < sfMin)
129  sf_value = sfMin;
130  for (std::vector<metsig::SigInputObj>::const_iterator metSignObject = metSignObjects.begin();
131  metSignObject != metSignObjects.end();
132  ++metSignObject) {
133  metsig::SigInputObj metSignObject_scaled;
134  metSignObject_scaled.set(metSignObject->get_type(),
135  sf_value * metSignObject->get_energy(),
136  metSignObject->get_phi(),
137  sf_value * metSignObject->get_sigma_e(),
138  metSignObject->get_sigma_tan());
139  metSignObjects_scaled.push_back(metSignObject_scaled);
140  }
141 }
142 
143 reco::METCovMatrix computePFMEtSignificance(const std::vector<metsig::SigInputObj>& metSignObjects) {
144  reco::METCovMatrix pfMEtCov;
145  if (metSignObjects.size() >= 2) {
146  metsig::significanceAlgo pfMEtSignAlgorithm;
147  pfMEtSignAlgorithm.addObjects(metSignObjects);
148  pfMEtCov = pfMEtSignAlgorithm.getSignifMatrix();
149  }
150 
151  double det = 0;
152  pfMEtCov.Det(det);
153  if (std::abs(det) < epsilon) {
154  edm::LogWarning("computePFMEtSignificance") << "Inversion of PFMEt covariance matrix failed, det = " << det
155  << " --> replacing covariance matrix by resolution defaults !!";
156  pfMEtCov(0, 0) = defaultPFMEtResolutionX * defaultPFMEtResolutionX;
157  pfMEtCov(0, 1) = 0.;
158  pfMEtCov(1, 0) = 0.;
159  pfMEtCov(1, 1) = defaultPFMEtResolutionY * defaultPFMEtResolutionY;
160  }
161 
162  return pfMEtCov;
163 }
164 
165 void printP4(const std::string& label_part1, int idx, const std::string& label_part2, const reco::Candidate& candidate) {
166  std::cout << label_part1 << " #" << idx << label_part2 << ": Pt = " << candidate.pt() << ", eta = " << candidate.eta()
167  << ", phi = " << candidate.phi() << " (charge = " << candidate.charge() << ")" << std::endl;
168 }
169 
170 void printCommonMETData(const std::string& label, const CommonMETData& metData) {
171  std::cout << label << ": Px = " << metData.mex << ", Py = " << metData.mey << ", sumEt = " << metData.sumet
172  << std::endl;
173 }
174 
175 void printMVAMEtJetInfo(const std::string& label, int idx, const reco::PUSubMETCandInfo& jet) {
176  std::cout << label << " #" << idx << " (";
177  if (jet.type() == reco::PUSubMETCandInfo::kHS)
178  std::cout << "no-PU";
179  else if (jet.type() == reco::PUSubMETCandInfo::kPU)
180  std::cout << "PU";
181  std::cout << "): Pt = " << jet.p4().pt() << ", eta = " << jet.p4().eta() << ", phi = " << jet.p4().phi();
182  std::cout << " id. flags: anti-noise = " << jet.passesLooseJetId() << std::endl;
183  std::cout << std::endl;
184 }
185 
186 void printMVAMEtPFCandInfo(const std::string& label, int idx, const reco::PUSubMETCandInfo& pfCand) {
187  std::cout << label << " #" << idx << " (";
188  if (pfCand.type() == reco::PUSubMETCandInfo::kChHS)
189  std::cout << "no-PU charged";
190  else if (pfCand.type() == reco::PUSubMETCandInfo::kChPU)
191  std::cout << "PU charged";
192  else if (pfCand.type() == reco::PUSubMETCandInfo::kNeutral)
193  std::cout << "neutral";
194  std::cout << "): Pt = " << pfCand.p4().pt() << ", eta = " << pfCand.p4().eta() << ", phi = " << pfCand.p4().phi();
195  std::string isWithinJet_string;
196  if (pfCand.isWithinJet())
197  isWithinJet_string = "true";
198  else
199  isWithinJet_string = "false";
200  std::cout << " (isWithinJet = " << isWithinJet_string << ")";
201  if (pfCand.isWithinJet())
202  std::cout << " Jet id. flags: anti-noise = " << pfCand.passesLooseJetId() << std::endl;
203  std::cout << std::endl;
204 }
205 
206 void NoPileUpPFMEtProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
207  LogDebug("produce") << " moduleLabel = " << moduleLabel_ << std::endl;
208 
209  // get original MET
210  edm::Handle<reco::PFMETCollection> pfMETs;
211  evt.getByToken(srcMEt_, pfMETs);
212  if (!(pfMETs->size() == 1))
213  throw cms::Exception("NoPileUpPFMEtProducer::produce") << "Failed to find unique MET object !!\n";
214  const reco::PFMET& pfMEt_original = pfMETs->front();
215 
216  // get MET covariance matrix
217  reco::METCovMatrix pfMEtCov;
218  if (!srcMEtCov_.label().empty()) {
219  //MM manual bypass to pfMET as this case has neer been presented
220  // edm::Handle<PFMEtSignCovMatrix> pfMEtCovHandle;
221  // evt.getByToken(srcMEtCov_, pfMEtCovHandle);
222  // pfMEtCov = (*pfMEtCovHandle);
223  pfMEtCov = pfMEt_original.getSignificanceMatrix();
224  } else {
225  pfMEtCov = pfMEt_original.getSignificanceMatrix();
226  }
227 
228  // get lepton momenta
229  std::vector<reco::Candidate::LorentzVector> leptons;
230  std::vector<metsig::SigInputObj> metSignObjectsLeptons;
231  reco::Candidate::LorentzVector sumLeptonP4s;
232  for (std::vector<edm::EDGetTokenT<edm::View<reco::Candidate> > >::const_iterator srcLeptons_i = srcLeptons_.begin();
233  srcLeptons_i != srcLeptons_.end();
234  ++srcLeptons_i) {
235  edm::Handle<reco::CandidateView> leptons_i;
236  evt.getByToken(*srcLeptons_i, leptons_i);
237  int leptonIdx = 0;
238  for (reco::CandidateView::const_iterator lepton = leptons_i->begin(); lepton != leptons_i->end(); ++lepton) {
239  leptons.push_back(lepton->p4());
240  metSignObjectsLeptons.push_back(pfMEtSignInterface_->compResolution(&(*lepton)));
241  sumLeptonP4s += lepton->p4();
242  ++leptonIdx;
243  }
244  }
245  LogDebug("produce") << " sum(leptons): Pt = " << sumLeptonP4s.pt() << ", eta = " << sumLeptonP4s.eta()
246  << ", phi = " << sumLeptonP4s.phi() << ","
247  << " mass = " << sumLeptonP4s.mass() << std::endl;
248 
249  // get jet and PFCandidate information
250  edm::Handle<reco::PUSubMETCandInfoCollection> jets;
251  evt.getByToken(srcJetInfo_, jets);
252  edm::Handle<reco::PUSubMETCandInfoCollection> jetsLeptonMatch;
253  evt.getByToken(srcJetInfoLeptonMatch_, jetsLeptonMatch);
254  edm::Handle<reco::PUSubMETCandInfoCollection> pfCandidates;
255  evt.getByToken(srcPFCandInfo_, pfCandidates);
256  edm::Handle<reco::PUSubMETCandInfoCollection> pfCandidatesLeptonMatch;
257  evt.getByToken(srcPFCandInfoLeptonMatch_, pfCandidatesLeptonMatch);
258 
259  reco::PUSubMETCandInfoCollection jets_leptons = utils_.cleanJets(*jetsLeptonMatch, leptons, 0.5, true);
260  reco::PUSubMETCandInfoCollection pfCandidates_leptons =
261  utils_.cleanPFCandidates(*pfCandidatesLeptonMatch, leptons, 0.3, true);
262  std::vector<CommonMETData> sumJetsPlusPFCandidates_leptons(leptons.size());
263  for (std::vector<CommonMETData>::iterator sumJetsPlusPFCandidates = sumJetsPlusPFCandidates_leptons.begin();
264  sumJetsPlusPFCandidates != sumJetsPlusPFCandidates_leptons.end();
265  ++sumJetsPlusPFCandidates) {
266  initializeCommonMETData(*sumJetsPlusPFCandidates);
267  }
268  for (reco::PUSubMETCandInfoCollection::const_iterator jet = jets_leptons.begin(); jet != jets_leptons.end(); ++jet) {
269  int leptonIdx_dRmin = findBestMatchingLepton(leptons, jet->p4());
270  assert(leptonIdx_dRmin >= 0 && leptonIdx_dRmin < (int)sumJetsPlusPFCandidates_leptons.size());
271 
272  LogDebug("produce") << "jet-to-lepton match:"
273  << " jetPt = " << jet->p4().pt() << ", jetEta = " << jet->p4().eta()
274  << ", jetPhi = " << jet->p4().phi() << " leptonPt = " << leptons[leptonIdx_dRmin].pt()
275  << ", leptonEta = " << leptons[leptonIdx_dRmin].eta()
276  << ", leptonPhi = " << leptons[leptonIdx_dRmin].phi() << std::endl;
277 
278  sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mex += jet->p4().px();
279  sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mey += jet->p4().py();
280  sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].sumet += jet->p4().pt();
281  }
282  for (reco::PUSubMETCandInfoCollection::const_iterator pfCandidate = pfCandidates_leptons.begin();
283  pfCandidate != pfCandidates_leptons.end();
284  ++pfCandidate) {
285  bool isWithinJet_lepton = false;
286  if (pfCandidate->isWithinJet()) {
287  for (reco::PUSubMETCandInfoCollection::const_iterator jet = jets_leptons.begin(); jet != jets_leptons.end();
288  ++jet) {
289  double dR2 = deltaR2(pfCandidate->p4(), jet->p4());
290  if (dR2 < 0.5 * 0.5)
291  isWithinJet_lepton = true;
292  }
293  }
294  if (!isWithinJet_lepton) {
295  int leptonIdx_dRmin = findBestMatchingLepton(leptons, pfCandidate->p4());
296  assert(leptonIdx_dRmin >= 0 && leptonIdx_dRmin < (int)sumJetsPlusPFCandidates_leptons.size());
297  LogDebug("produce") << "pfCandidate-to-lepton match:"
298  << " pfCandidatePt = " << pfCandidate->p4().pt()
299  << ", pfCandidateEta = " << pfCandidate->p4().eta()
300  << ", pfCandidatePhi = " << pfCandidate->p4().phi()
301  << " leptonPt = " << leptons[leptonIdx_dRmin].pt()
302  << ", leptonEta = " << leptons[leptonIdx_dRmin].eta()
303  << ", leptonPhi = " << leptons[leptonIdx_dRmin].phi() << std::endl;
304 
305  sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mex += pfCandidate->p4().px();
306  sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mey += pfCandidate->p4().py();
307  sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].sumet += pfCandidate->p4().pt();
308  } else {
309  LogDebug("produce") << " pfCandidate is within jet --> skipping." << std::endl;
310  }
311  }
312  auto sumLeptons = std::make_unique<CommonMETData>();
313  initializeCommonMETData(*sumLeptons);
314  auto sumLeptonIsoCones = std::make_unique<CommonMETData>();
315  initializeCommonMETData(*sumLeptonIsoCones);
316  int leptonIdx = 0;
317  for (std::vector<CommonMETData>::iterator sumJetsPlusPFCandidates = sumJetsPlusPFCandidates_leptons.begin();
318  sumJetsPlusPFCandidates != sumJetsPlusPFCandidates_leptons.end();
319  ++sumJetsPlusPFCandidates) {
320  if (sumJetsPlusPFCandidates->sumet > leptons[leptonIdx].pt()) {
321  double leptonEnFrac = leptons[leptonIdx].pt() / sumJetsPlusPFCandidates->sumet;
322  assert(leptonEnFrac >= 0.0 && leptonEnFrac <= 1.0);
323  sumLeptons->mex += (leptonEnFrac * sumJetsPlusPFCandidates->mex);
324  sumLeptons->mey += (leptonEnFrac * sumJetsPlusPFCandidates->mey);
325  sumLeptons->sumet += (leptonEnFrac * sumJetsPlusPFCandidates->sumet);
326  double leptonIsoConeEnFrac = 1.0 - leptonEnFrac;
327  assert(leptonIsoConeEnFrac >= 0.0 && leptonIsoConeEnFrac <= 1.0);
328  sumLeptonIsoCones->mex += (leptonIsoConeEnFrac * sumJetsPlusPFCandidates->mex);
329  sumLeptonIsoCones->mey += (leptonIsoConeEnFrac * sumJetsPlusPFCandidates->mey);
330  sumLeptonIsoCones->sumet += (leptonIsoConeEnFrac * sumJetsPlusPFCandidates->sumet);
331  } else {
332  sumLeptons->mex += sumJetsPlusPFCandidates->mex;
333  sumLeptons->mey += sumJetsPlusPFCandidates->mey;
334  sumLeptons->sumet += sumJetsPlusPFCandidates->sumet;
335  }
336  ++leptonIdx;
337  }
338 
339  reco::PUSubMETCandInfoCollection jets_cleaned = utils_.cleanJets(*jets, leptons, 0.5, false);
340  reco::PUSubMETCandInfoCollection pfCandidates_cleaned = utils_.cleanPFCandidates(*pfCandidates, leptons, 0.3, false);
341 
342  auto sumNoPUjets = std::make_unique<CommonMETData>();
343  initializeCommonMETData(*sumNoPUjets);
344  std::vector<metsig::SigInputObj> metSignObjectsNoPUjets;
345  auto sumNoPUjetOffsetEnCorr = std::make_unique<CommonMETData>();
346  initializeCommonMETData(*sumNoPUjetOffsetEnCorr);
347  std::vector<metsig::SigInputObj> metSignObjectsNoPUjetOffsetEnCorr;
348  auto sumPUjets = std::make_unique<CommonMETData>();
349  initializeCommonMETData(*sumPUjets);
350  std::vector<metsig::SigInputObj> metSignObjectsPUjets;
351  int jetIdx = 0;
352  for (reco::PUSubMETCandInfoCollection::const_iterator jet = jets_cleaned.begin(); jet != jets_cleaned.end(); ++jet) {
353  if (jet->passesLooseJetId()) {
354  if (jet->type() == reco::PUSubMETCandInfo::kHS) {
355  addToCommonMETData(*sumNoPUjets, jet->p4());
356  metSignObjectsNoPUjets.push_back(jet->metSignObj());
357  float jetp = jet->p4().P();
358  float jetcorr = jet->offsetEnCorr();
359  sumNoPUjetOffsetEnCorr->mex += jetcorr * jet->p4().px() / jetp;
360  sumNoPUjetOffsetEnCorr->mey += jetcorr * jet->p4().py() / jetp;
361  sumNoPUjetOffsetEnCorr->mez += jetcorr * jet->p4().pz() / jetp;
362  sumNoPUjetOffsetEnCorr->sumet += jetcorr * jet->p4().pt() / jetp;
363  metsig::SigInputObj pfMEtSignObjectOffsetEnCorr(
364  jet->metSignObj().get_type(),
365  jet->offsetEnCorr(),
366  jet->metSignObj().get_phi(),
367  (jet->offsetEnCorr() / jet->p4().E()) * jet->metSignObj().get_sigma_e(),
368  jet->metSignObj().get_sigma_tan());
369  metSignObjectsNoPUjetOffsetEnCorr.push_back(pfMEtSignObjectOffsetEnCorr);
370  } else {
371  addToCommonMETData(*sumPUjets, jet->p4());
372  metSignObjectsPUjets.push_back(jet->metSignObj());
373  }
374  }
375  ++jetIdx;
376  }
377 
378  auto sumNoPUunclChargedCands = std::make_unique<CommonMETData>();
379  initializeCommonMETData(*sumNoPUunclChargedCands);
380  std::vector<metsig::SigInputObj> metSignObjectsNoPUunclChargedCands;
381  auto sumPUunclChargedCands = std::make_unique<CommonMETData>();
382  initializeCommonMETData(*sumPUunclChargedCands);
383  std::vector<metsig::SigInputObj> metSignObjectsPUunclChargedCands;
384  auto sumUnclNeutralCands = std::make_unique<CommonMETData>();
385  initializeCommonMETData(*sumUnclNeutralCands);
386  std::vector<metsig::SigInputObj> metSignObjectsUnclNeutralCands;
387  int pfCandIdx = 0;
388  for (reco::PUSubMETCandInfoCollection::const_iterator pfCandidate = pfCandidates_cleaned.begin();
389  pfCandidate != pfCandidates_cleaned.end();
390  ++pfCandidate) {
391  if (pfCandidate->passesLooseJetId()) {
392  if (!pfCandidate->isWithinJet()) {
393  if (pfCandidate->type() == reco::PUSubMETCandInfo::kChHS) {
394  addToCommonMETData(*sumNoPUunclChargedCands, pfCandidate->p4());
395  metSignObjectsNoPUunclChargedCands.push_back(pfCandidate->metSignObj());
396  } else if (pfCandidate->type() == reco::PUSubMETCandInfo::kChPU) {
397  addToCommonMETData(*sumPUunclChargedCands, pfCandidate->p4());
398  metSignObjectsPUunclChargedCands.push_back(pfCandidate->metSignObj());
399  } else if (pfCandidate->type() == reco::PUSubMETCandInfo::kNeutral) {
400  addToCommonMETData(*sumUnclNeutralCands, pfCandidate->p4());
401  metSignObjectsUnclNeutralCands.push_back(pfCandidate->metSignObj());
402  }
403  }
404  }
405  ++pfCandIdx;
406  }
407 
408  edm::Handle<CorrMETData> type0Correction_input;
409  evt.getByToken(srcType0Correction_, type0Correction_input);
410  auto type0Correction_output = std::make_unique<CommonMETData>();
411  initializeCommonMETData(*type0Correction_output);
412  type0Correction_output->mex = type0Correction_input->mex;
413  type0Correction_output->mey = type0Correction_input->mey;
414 
415  finalizeCommonMETData(*sumLeptons);
416  finalizeCommonMETData(*sumNoPUjetOffsetEnCorr);
417  finalizeCommonMETData(*sumNoPUjets);
418  finalizeCommonMETData(*sumPUjets);
419  finalizeCommonMETData(*sumNoPUunclChargedCands);
420  finalizeCommonMETData(*sumPUunclChargedCands);
421  finalizeCommonMETData(*sumUnclNeutralCands);
422  finalizeCommonMETData(*type0Correction_output);
423  finalizeCommonMETData(*sumLeptonIsoCones);
424 
425  double noPileUpScaleFactor =
426  (sumPUunclChargedCands->sumet > 0.)
427  ? (sumPUunclChargedCands->sumet / (sumNoPUunclChargedCands->sumet + sumPUunclChargedCands->sumet))
428  : 1.;
429  LogDebug("produce") << "noPileUpScaleFactor = " << noPileUpScaleFactor << std::endl;
430 
431  double noPileUpMEtPx =
432  -(sumLeptons->mex + sumNoPUjets->mex + sumNoPUunclChargedCands->mex +
433  noPileUpScaleFactor *
434  (sfNoPUjetOffsetEnCorr_ * sumNoPUjetOffsetEnCorr->mex + sfUnclNeutralCands_ * sumUnclNeutralCands->mex +
435  sfPUunclChargedCands_ * sumPUunclChargedCands->mex + sfPUjets_ * sumPUjets->mex)) +
436  noPileUpScaleFactor * sfType0Correction_ * type0Correction_output->mex;
437  if (sfLeptonIsoCones_ >= 0.)
438  noPileUpMEtPx -= (noPileUpScaleFactor * sfLeptonIsoCones_ * sumLeptonIsoCones->mex);
439  else
440  noPileUpMEtPx -= (std::abs(sfLeptonIsoCones_) * sumLeptonIsoCones->mex);
441  double noPileUpMEtPy =
442  -(sumLeptons->mey + sumNoPUjets->mey + sumNoPUunclChargedCands->mey +
443  noPileUpScaleFactor *
444  (sfNoPUjetOffsetEnCorr_ * sumNoPUjetOffsetEnCorr->mey + sfUnclNeutralCands_ * sumUnclNeutralCands->mey +
445  sfPUunclChargedCands_ * sumPUunclChargedCands->mey + sfPUjets_ * sumPUjets->mey)) +
446  noPileUpScaleFactor * sfType0Correction_ * type0Correction_output->mey;
447  if (sfLeptonIsoCones_ >= 0.)
448  noPileUpMEtPy -= (noPileUpScaleFactor * sfLeptonIsoCones_ * sumLeptonIsoCones->mey);
449  else
450  noPileUpMEtPy -= (std::abs(sfLeptonIsoCones_) * sumLeptonIsoCones->mey);
451  double noPileUpMEtPt = sqrt(noPileUpMEtPx * noPileUpMEtPx + noPileUpMEtPy * noPileUpMEtPy);
452  reco::Candidate::LorentzVector noPileUpMEtP4(noPileUpMEtPx, noPileUpMEtPy, 0., noPileUpMEtPt);
453 
454  reco::PFMET noPileUpMEt(pfMEt_original);
455  noPileUpMEt.setP4(noPileUpMEtP4);
456  //noPileUpMEt.setSignificanceMatrix(pfMEtCov);
457 
458  std::vector<metsig::SigInputObj> metSignObjects_scaled;
459  scaleAndAddPFMEtSignObjects(metSignObjects_scaled, metSignObjectsLeptons, 1.0, sfMEtCovMin_, sfMEtCovMax_);
460  scaleAndAddPFMEtSignObjects(
461  metSignObjects_scaled, metSignObjectsNoPUjetOffsetEnCorr, sfNoPUjetOffsetEnCorr_, sfMEtCovMin_, sfMEtCovMax_);
462  scaleAndAddPFMEtSignObjects(metSignObjects_scaled, metSignObjectsNoPUjets, sfNoPUjets_, sfMEtCovMin_, sfMEtCovMax_);
463  scaleAndAddPFMEtSignObjects(
464  metSignObjects_scaled, metSignObjectsPUjets, noPileUpScaleFactor * sfPUjets_, sfMEtCovMin_, sfMEtCovMax_);
465  scaleAndAddPFMEtSignObjects(
466  metSignObjects_scaled, metSignObjectsNoPUunclChargedCands, sfNoPUunclChargedCands_, sfMEtCovMin_, sfMEtCovMax_);
467  scaleAndAddPFMEtSignObjects(metSignObjects_scaled,
468  metSignObjectsPUunclChargedCands,
469  noPileUpScaleFactor * sfPUunclChargedCands_,
470  sfMEtCovMin_,
471  sfMEtCovMax_);
472  scaleAndAddPFMEtSignObjects(metSignObjects_scaled,
473  metSignObjectsUnclNeutralCands,
474  noPileUpScaleFactor * sfUnclNeutralCands_,
475  sfMEtCovMin_,
476  sfMEtCovMax_);
477  reco::METCovMatrix pfMEtCov_recomputed = computePFMEtSignificance(metSignObjects_scaled);
478  noPileUpMEt.setSignificanceMatrix(pfMEtCov_recomputed);
479 
480  LogDebug("produce") << "<NoPileUpPFMEtProducer::produce>:" << std::endl
481  << " moduleLabel = " << moduleLabel_ << std::endl
482  << " PFMET: Pt = " << pfMEt_original.pt() << ", phi = " << pfMEt_original.phi() << " "
483  << "(Px = " << pfMEt_original.px() << ", Py = " << pfMEt_original.py() << ")" << std::endl
484  << " Cov:" << std::endl
485  << " " << pfMEtCov(0, 0) << " " << pfMEtCov(0, 1) << "\n " << pfMEtCov(1, 0) << " "
486  << pfMEtCov(1, 1) << std::endl
487  << " no-PU MET: Pt = " << noPileUpMEt.pt() << ", phi = " << noPileUpMEt.phi() << " "
488  << "(Px = " << noPileUpMEt.px() << ", Py = " << noPileUpMEt.py() << ")" << std::endl
489  << " Cov:" << std::endl
490  << " " << (noPileUpMEt.getSignificanceMatrix())(0, 0) << " "
491  << (noPileUpMEt.getSignificanceMatrix())(0, 1) << std::endl
492  << (noPileUpMEt.getSignificanceMatrix())(1, 0) << " "
493  << (noPileUpMEt.getSignificanceMatrix())(1, 1) << std::endl;
494 
495  // add no-PU MET object to the event
496  auto noPileUpMEtCollection = std::make_unique<reco::PFMETCollection>();
497  noPileUpMEtCollection->push_back(noPileUpMEt);
498 
499  evt.put(std::move(noPileUpMEtCollection));
500  if (saveInputs_) {
501  evt.put(std::move(sumLeptons), sfLeptonsName_);
502  evt.put(std::move(sumNoPUjetOffsetEnCorr), sfNoPUjetOffsetEnCorrName_);
503  evt.put(std::move(sumNoPUjets), sfNoPUjetsName_);
504  evt.put(std::move(sumPUjets), sfPUjetsName_);
505  evt.put(std::move(sumNoPUunclChargedCands), sfNoPUunclChargedCandsName_);
506  evt.put(std::move(sumPUunclChargedCands), sfPUunclChargedCandsName_);
507  evt.put(std::move(sumUnclNeutralCands), sfUnclNeutralCandsName_);
508  evt.put(std::move(type0Correction_output), sfType0CorrectionName_);
509  evt.put(std::move(sumLeptonIsoCones), sfLeptonIsoConesName_);
510  }
511 
512  evt.put(std::make_unique<double>(noPileUpScaleFactor), "sfNoPU");
513 }
514 
515 #include "FWCore/Framework/interface/MakerMacros.h"
516 
517 DEFINE_FWK_MODULE(NoPileUpPFMEtProducer);
NoPileUpPFMEtProducer::sfUnclNeutralCandsName_
std::string sfUnclNeutralCandsName_
Definition: NoPileUpPFMEtProducer.h:66
CommonMETData::met
double met
Definition: CommonMETData.h:14
significanceAlgo.h
reco::PUSubMETCandInfo::kHS
Definition: PUSubMETData.h:63
NoPileUpPFMEtProducer::sfPUjets_
double sfPUjets_
Definition: NoPileUpPFMEtProducer.h:53
edm::View::begin
const_iterator begin() const
NoPileUpPFMEtProducer::vInputTag
std::vector< edm::InputTag > vInputTag
Definition: NoPileUpPFMEtProducer.h:46
zmumugammaAnalyzer_cfi.pfCandidates
pfCandidates
Definition: zmumugammaAnalyzer_cfi.py:11
NoPileUpPFMEtProducer::sfLeptonsName_
std::string sfLeptonsName_
Definition: NoPileUpPFMEtProducer.h:60
NoPileUpPFMEtProducer::srcJetInfo_
edm::EDGetTokenT< reco::PUSubMETCandInfoCollection > srcJetInfo_
Definition: NoPileUpPFMEtProducer.h:42
NoPileUpPFMEtProducer::pfMEtSignInterface_
PFMEtSignInterfaceBase * pfMEtSignInterface_
Definition: NoPileUpPFMEtProducer.h:70
NoPileUpPFMEtProducer::sfNoPUjetOffsetEnCorr_
double sfNoPUjetOffsetEnCorr_
Definition: NoPileUpPFMEtProducer.h:52
NoPileUpPFMEtProducer::sfNoPUunclChargedCands_
double sfNoPUunclChargedCands_
Definition: NoPileUpPFMEtProducer.h:54
CorrMETData::mey
double mey
Definition: CorrMETData.h:16
printMVAMEtPFCandInfo
void printMVAMEtPFCandInfo(const std::string &label, int idx, const reco::PUSubMETCandInfo &pfCand)
Definition: NoPileUpPFMEtProducer.cc:186
NoPileUpPFMEtProducer::srcMEtCov_
edm::InputTag srcMEtCov_
Definition: NoPileUpPFMEtProducer.h:41
reco::Candidate::eta
virtual double eta() const =0
momentum pseudorapidity
edm::EDGetTokenT
Definition: EDGetToken.h:33
reco::MET::getSignificanceMatrix
reco::METCovMatrix getSignificanceMatrix(void) const
Definition: MET.cc:128
gather_cfg.cout
cout
Definition: gather_cfg.py:144
metsig::SigInputObj::set
void set(const std::string &m_type, const double &m_energy, const double &m_phi, const double &m_sigma_e, const double &m_sigma_tan)
Definition: SigInputObj.h:42
NoPileUpMEtUtilities::cleanPFCandidates
reco::PUSubMETCandInfoCollection cleanPFCandidates(const reco::PUSubMETCandInfoCollection &, const std::vector< reco::Candidate::LorentzVector > &, double, bool)
Definition: NoPileUpMEtUtilities.cc:81
NoPileUpPFMEtProducer::moduleLabel_
std::string moduleLabel_
Definition: NoPileUpPFMEtProducer.h:38
reco::Candidate::pt
virtual double pt() const =0
transverse momentum
cms::cuda::assert
assert(be >=bs)
defaultPFMEtResolutionX
const double defaultPFMEtResolutionX
Definition: NoPileUpPFMEtProducer.cc:11
defaultPFMEtResolutionY
const double defaultPFMEtResolutionY
Definition: NoPileUpPFMEtProducer.cc:12
CommonMETData::mez
double mez
Definition: CommonMETData.h:17
NoPileUpPFMEtProducer::saveInputs_
bool saveInputs_
Definition: NoPileUpPFMEtProducer.h:74
NoPileUpPFMEtProducer::~NoPileUpPFMEtProducer
~NoPileUpPFMEtProducer() override
Definition: NoPileUpPFMEtProducer.cc:79
CorrMETData::mex
double mex
Definition: CorrMETData.h:15
singleTopDQM_cfi.jets
jets
Definition: singleTopDQM_cfi.py:42
findBestMatchingLepton
int findBestMatchingLepton(const std::vector< reco::Candidate::LorentzVector > &leptons, const reco::Candidate::LorentzVector &p4_ref)
Definition: NoPileUpPFMEtProducer.cc:102
CommonMETData
Structure containing data common to all types of MET.
Definition: CommonMETData.h:12
reco::LeafCandidate::pt
double pt() const final
transverse momentum
Definition: LeafCandidate.h:146
reco::PUSubMETCandInfo::kPU
Definition: PUSubMETData.h:63
edm::Handle< reco::PFMETCollection >
training_settings.idx
idx
Definition: training_settings.py:16
initializeCommonMETData
void initializeCommonMETData(CommonMETData &metData)
Definition: NoPileUpPFMEtProducer.cc:81
PFMEtSignInterfaceBase::compResolution
metsig::SigInputObj compResolution(const T *particle) const
Definition: PFMEtSignInterfaceBase.h:45
NoPileUpPFMEtProducer.h
NoPileUpPFMEtProducer::sfNoPUjets_
double sfNoPUjets_
Definition: NoPileUpPFMEtProducer.h:51
edm::InputTag::label
std::string const & label() const
Definition: InputTag.h:36
NoPileUpPFMEtProducer::sfLeptonIsoCones_
double sfLeptonIsoCones_
Definition: NoPileUpPFMEtProducer.h:58
metsig::significanceAlgo::addObjects
const void addObjects(const std::vector< metsig::SigInputObj > &EventVec)
Definition: significanceAlgo.cc:109
reco::MET::setSignificanceMatrix
void setSignificanceMatrix(const reco::METCovMatrix &matrix)
Definition: MET.cc:142
MakerMacros.h
reco::PUSubMETCandInfo
Definition: PUSubMETData.h:21
NoPileUpMEtUtilities::cleanJets
reco::PUSubMETCandInfoCollection cleanJets(const reco::PUSubMETCandInfoCollection &, const std::vector< reco::Candidate::LorentzVector > &, double, bool)
Definition: NoPileUpMEtUtilities.cc:28
reco::PUSubMETCandInfo::type
int type() const
Definition: PUSubMETData.h:33
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
CommonMETData::mex
double mex
Definition: CommonMETData.h:15
reco::PUSubMETCandInfo::p4
const reco::Candidate::LorentzVector & p4() const
Definition: PUSubMETData.h:30
reco::PUSubMETCandInfo::kChHS
Definition: PUSubMETData.h:63
NoPileUpPFMEtProducer::sfType0CorrectionName_
std::string sfType0CorrectionName_
Definition: NoPileUpPFMEtProducer.h:67
reco::LeafCandidate::py
double py() const final
y coordinate of momentum vector
Definition: LeafCandidate.h:142
reco::Candidate::charge
virtual int charge() const =0
electric charge
SigInputObj.h
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
NoPileUpPFMEtProducer::sfMEtCovMin_
double sfMEtCovMin_
Definition: NoPileUpPFMEtProducer.h:71
reco::PUSubMETCandInfo::kChPU
Definition: PUSubMETData.h:63
reco::PUSubMETCandInfo::passesLooseJetId
float passesLooseJetId() const
Definition: PUSubMETData.h:38
printCommonMETData
void printCommonMETData(const std::string &label, const CommonMETData &metData)
Definition: NoPileUpPFMEtProducer.cc:170
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:528
NoPileUpPFMEtProducer::sfLeptonIsoConesName_
std::string sfLeptonIsoConesName_
Definition: NoPileUpPFMEtProducer.h:68
NoPileUpPFMEtProducer::sfNoPUjetOffsetEnCorrName_
std::string sfNoPUjetOffsetEnCorrName_
Definition: NoPileUpPFMEtProducer.h:62
addToCommonMETData
void addToCommonMETData(CommonMETData &metData, const reco::Candidate::LorentzVector &p4)
Definition: NoPileUpPFMEtProducer.cc:90
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
edm::LogWarning
Definition: MessageLogger.h:141
CommonMETData::sumet
double sumet
Definition: CommonMETData.h:18
CommonMETData::phi
double phi
Definition: CommonMETData.h:19
edm::View
Definition: CaloClusterFwd.h:14
reco::PFMET
Definition: PFMET.h:18
NoPileUpPFMEtProducer::produce
void produce(edm::Event &, const edm::EventSetup &) override
Definition: NoPileUpPFMEtProducer.cc:206
HLT_2018_cff.InputTag
InputTag
Definition: HLT_2018_cff.py:79016
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:670
edm::ParameterSet
Definition: ParameterSet.h:36
electrons_cff.jetIdx
jetIdx
Definition: electrons_cff.py:356
NoPileUpPFMEtProducer::sfPUjetsName_
std::string sfPUjetsName_
Definition: NoPileUpPFMEtProducer.h:63
NoPileUpPFMEtProducer::srcType0Correction_
edm::EDGetTokenT< CorrMETData > srcType0Correction_
Definition: NoPileUpPFMEtProducer.h:49
HLT_2018_cff.leptons
leptons
Definition: HLT_2018_cff.py:24820
epsilon
const double epsilon
Definition: NoPileUpPFMEtProducer.cc:14
NoPileUpPFMEtProducer::NoPileUpPFMEtProducer
NoPileUpPFMEtProducer(const edm::ParameterSet &)
Definition: NoPileUpPFMEtProducer.cc:16
p4
double p4[4]
Definition: TauolaWrapper.h:92
edm::Event::put
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:132
NoPileUpPFMEtProducer::sfPUunclChargedCands_
double sfPUunclChargedCands_
Definition: NoPileUpPFMEtProducer.h:55
NoPileUpPFMEtProducer::srcJetInfoLeptonMatch_
edm::EDGetTokenT< reco::PUSubMETCandInfoCollection > srcJetInfoLeptonMatch_
Definition: NoPileUpPFMEtProducer.h:43
edm::EventSetup
Definition: EventSetup.h:57
looper.cfg
cfg
Definition: looper.py:297
reco::Candidate
Definition: Candidate.h:27
computePFMEtSignificance
reco::METCovMatrix computePFMEtSignificance(const std::vector< metsig::SigInputObj > &metSignObjects)
Definition: NoPileUpPFMEtProducer.cc:143
reco::LeafCandidate::setP4
void setP4(const LorentzVector &p4) final
set 4-momentum
Definition: LeafCandidate.h:158
NoPileUpPFMEtProducer::sfNoPUunclChargedCandsName_
std::string sfNoPUunclChargedCandsName_
Definition: NoPileUpPFMEtProducer.h:64
eostools.move
def move(src, dest)
Definition: eostools.py:511
NoPileUpPFMEtProducer::sfType0Correction_
double sfType0Correction_
Definition: NoPileUpPFMEtProducer.h:57
std
Definition: JetResolutionObject.h:76
reco::PUSubMETCandInfo::kNeutral
Definition: PUSubMETData.h:63
reco::LeafCandidate::phi
double phi() const final
momentum azimuthal angle
Definition: LeafCandidate.h:148
NoPileUpPFMEtProducer::verbosity_
int verbosity_
Definition: NoPileUpPFMEtProducer.h:76
metsig::significanceAlgo
Definition: significanceAlgo.h:86
HLTMuonOfflineAnalyzer_cfi.deltaR2
deltaR2
Definition: HLTMuonOfflineAnalyzer_cfi.py:105
metsig::jet
Definition: SignAlgoResolutions.h:47
NoPileUpPFMEtProducer::srcMEt_
edm::EDGetTokenT< reco::PFMETCollection > srcMEt_
Definition: NoPileUpPFMEtProducer.h:40
Exception
Definition: hltDiff.cc:246
printP4
void printP4(const std::string &label_part1, int idx, const std::string &label_part2, const reco::Candidate &candidate)
Definition: NoPileUpPFMEtProducer.cc:165
NoPileUpPFMEtProducer::sfNoPUjetsName_
std::string sfNoPUjetsName_
Definition: NoPileUpPFMEtProducer.h:61
metsig::significanceAlgo::getSignifMatrix
reco::METCovMatrix getSignifMatrix() const
Definition: significanceAlgo.h:99
Exception.h
edm::View::const_iterator
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:86
PFMEtSignInterfaceBase
Definition: PFMEtSignInterfaceBase.h:39
reco::PUSubMETCandInfo::isWithinJet
bool isWithinJet() const
Definition: PUSubMETData.h:37
NoPileUpPFMEtProducer::sfMEtCovMax_
double sfMEtCovMax_
Definition: NoPileUpPFMEtProducer.h:72
NoPileUpPFMEtProducer::srcPFCandInfoLeptonMatch_
edm::EDGetTokenT< reco::PUSubMETCandInfoCollection > srcPFCandInfoLeptonMatch_
Definition: NoPileUpPFMEtProducer.h:45
NoPileUpPFMEtProducer::srcLeptons_
std::vector< edm::EDGetTokenT< reco::CandidateView > > srcLeptons_
Definition: NoPileUpPFMEtProducer.h:47
reco::PUSubMETCandInfoCollection
std::vector< reco::PUSubMETCandInfo > PUSubMETCandInfoCollection
Definition: PUSubMETDataFwd.h:9
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
edm::View::end
const_iterator end() const
finalizeCommonMETData
void finalizeCommonMETData(CommonMETData &metData)
Definition: NoPileUpPFMEtProducer.cc:97
NoPileUpPFMEtProducer::sfUnclNeutralCands_
double sfUnclNeutralCands_
Definition: NoPileUpPFMEtProducer.h:56
NoPileUpPFMEtProducer::srcPFCandInfo_
edm::EDGetTokenT< reco::PUSubMETCandInfoCollection > srcPFCandInfo_
Definition: NoPileUpPFMEtProducer.h:44
reco::Candidate::LorentzVector
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Candidate.h:36
NoPileUpPFMEtProducer
Definition: NoPileUpPFMEtProducer.h:30
edm::Event
Definition: Event.h:73
printMVAMEtJetInfo
void printMVAMEtJetInfo(const std::string &label, int idx, const reco::PUSubMETCandInfo &jet)
Definition: NoPileUpPFMEtProducer.cc:175
reco::Candidate::phi
virtual double phi() const =0
momentum azimuthal angle
reco::LeafCandidate::px
double px() const final
x coordinate of momentum vector
Definition: LeafCandidate.h:140
CommonMETData::mey
double mey
Definition: CommonMETData.h:16
NoPileUpPFMEtProducer::sfPUunclChargedCandsName_
std::string sfPUunclChargedCandsName_
Definition: NoPileUpPFMEtProducer.h:65
edm::InputTag
Definition: InputTag.h:15
metsig::SigInputObj
Definition: SigInputObj.h:29
label
const char * label
Definition: PFTauDecayModeTools.cc:11
reco::METCovMatrix
ROOT::Math::SMatrix< double, 2 > METCovMatrix
Definition: MET.h:39
NoPileUpPFMEtProducer::utils_
NoPileUpMEtUtilities utils_
Definition: NoPileUpPFMEtProducer.h:78
scaleAndAddPFMEtSignObjects
void scaleAndAddPFMEtSignObjects(std::vector< metsig::SigInputObj > &metSignObjects_scaled, const std::vector< metsig::SigInputObj > &metSignObjects, double sf, double sfMin, double sfMax)
Definition: NoPileUpPFMEtProducer.cc:120
Generated for CMSSW Reference Manual by doxygen 1.8.16