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