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();
229 std::vector<reco::Candidate::LorentzVector>
leptons;
230 std::vector<metsig::SigInputObj> metSignObjectsLeptons;
238 leptons.push_back(lepton->p4());
240 sumLeptonP4s += lepton->p4();
243 LogDebug(
"produce") <<
" sum(leptons): Pt = " << sumLeptonP4s.pt() <<
", eta = " << sumLeptonP4s.eta()
244 <<
", phi = " << sumLeptonP4s.phi() <<
"," 245 <<
" mass = " << sumLeptonP4s.mass() << std::endl;
260 std::vector<CommonMETData> sumJetsPlusPFCandidates_leptons(
leptons.size());
261 for (std::vector<CommonMETData>::iterator sumJetsPlusPFCandidates = sumJetsPlusPFCandidates_leptons.begin();
262 sumJetsPlusPFCandidates != sumJetsPlusPFCandidates_leptons.end();
263 ++sumJetsPlusPFCandidates) {
266 for (reco::PUSubMETCandInfoCollection::const_iterator
jet = jets_leptons.begin();
jet != jets_leptons.end(); ++
jet) {
268 assert(leptonIdx_dRmin >= 0 && leptonIdx_dRmin < (
int)sumJetsPlusPFCandidates_leptons.size());
270 LogDebug(
"produce") <<
"jet-to-lepton match:" 271 <<
" jetPt = " <<
jet->p4().pt() <<
", jetEta = " <<
jet->p4().eta()
272 <<
", jetPhi = " <<
jet->p4().phi() <<
" leptonPt = " <<
leptons[leptonIdx_dRmin].pt()
273 <<
", leptonEta = " <<
leptons[leptonIdx_dRmin].eta()
274 <<
", leptonPhi = " <<
leptons[leptonIdx_dRmin].phi() << std::endl;
276 sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mex +=
jet->p4().px();
277 sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mey +=
jet->p4().py();
278 sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].sumet +=
jet->p4().pt();
280 for (reco::PUSubMETCandInfoCollection::const_iterator pfCandidate = pfCandidates_leptons.begin();
281 pfCandidate != pfCandidates_leptons.end();
283 bool isWithinJet_lepton =
false;
284 if (pfCandidate->isWithinJet()) {
285 for (reco::PUSubMETCandInfoCollection::const_iterator
jet = jets_leptons.begin();
jet != jets_leptons.end();
287 double dR2 =
deltaR2(pfCandidate->p4(),
jet->p4());
289 isWithinJet_lepton =
true;
292 if (!isWithinJet_lepton) {
294 assert(leptonIdx_dRmin >= 0 && leptonIdx_dRmin < (
int)sumJetsPlusPFCandidates_leptons.size());
295 LogDebug(
"produce") <<
"pfCandidate-to-lepton match:" 296 <<
" pfCandidatePt = " << pfCandidate->p4().pt()
297 <<
", pfCandidateEta = " << pfCandidate->p4().eta()
298 <<
", pfCandidatePhi = " << pfCandidate->p4().phi()
299 <<
" leptonPt = " <<
leptons[leptonIdx_dRmin].pt()
300 <<
", leptonEta = " <<
leptons[leptonIdx_dRmin].eta()
301 <<
", leptonPhi = " <<
leptons[leptonIdx_dRmin].phi() << std::endl;
303 sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mex += pfCandidate->p4().px();
304 sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mey += pfCandidate->p4().py();
305 sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].sumet += pfCandidate->p4().pt();
307 LogDebug(
"produce") <<
" pfCandidate is within jet --> skipping." << std::endl;
310 auto sumLeptons = std::make_unique<CommonMETData>();
312 auto sumLeptonIsoCones = std::make_unique<CommonMETData>();
315 for (std::vector<CommonMETData>::iterator sumJetsPlusPFCandidates = sumJetsPlusPFCandidates_leptons.begin();
316 sumJetsPlusPFCandidates != sumJetsPlusPFCandidates_leptons.end();
317 ++sumJetsPlusPFCandidates) {
318 if (sumJetsPlusPFCandidates->sumet >
leptons[leptonIdx].pt()) {
319 double leptonEnFrac =
leptons[leptonIdx].pt() / sumJetsPlusPFCandidates->sumet;
320 assert(leptonEnFrac >= 0.0 && leptonEnFrac <= 1.0);
321 sumLeptons->mex += (leptonEnFrac * sumJetsPlusPFCandidates->mex);
322 sumLeptons->mey += (leptonEnFrac * sumJetsPlusPFCandidates->mey);
323 sumLeptons->sumet += (leptonEnFrac * sumJetsPlusPFCandidates->sumet);
324 double leptonIsoConeEnFrac = 1.0 - leptonEnFrac;
325 assert(leptonIsoConeEnFrac >= 0.0 && leptonIsoConeEnFrac <= 1.0);
326 sumLeptonIsoCones->mex += (leptonIsoConeEnFrac * sumJetsPlusPFCandidates->mex);
327 sumLeptonIsoCones->mey += (leptonIsoConeEnFrac * sumJetsPlusPFCandidates->mey);
328 sumLeptonIsoCones->sumet += (leptonIsoConeEnFrac * sumJetsPlusPFCandidates->sumet);
330 sumLeptons->mex += sumJetsPlusPFCandidates->mex;
331 sumLeptons->mey += sumJetsPlusPFCandidates->mey;
332 sumLeptons->sumet += sumJetsPlusPFCandidates->sumet;
340 auto sumNoPUjets = std::make_unique<CommonMETData>();
342 std::vector<metsig::SigInputObj> metSignObjectsNoPUjets;
343 auto sumNoPUjetOffsetEnCorr = std::make_unique<CommonMETData>();
345 std::vector<metsig::SigInputObj> metSignObjectsNoPUjetOffsetEnCorr;
346 auto sumPUjets = std::make_unique<CommonMETData>();
348 std::vector<metsig::SigInputObj> metSignObjectsPUjets;
349 for (reco::PUSubMETCandInfoCollection::const_iterator
jet = jets_cleaned.begin();
jet != jets_cleaned.end(); ++
jet) {
350 if (
jet->passesLooseJetId()) {
353 metSignObjectsNoPUjets.push_back(
jet->metSignObj());
354 float jetp =
jet->p4().P();
355 float jetcorr =
jet->offsetEnCorr();
356 sumNoPUjetOffsetEnCorr->mex += jetcorr *
jet->p4().px() / jetp;
357 sumNoPUjetOffsetEnCorr->mey += jetcorr *
jet->p4().py() / jetp;
358 sumNoPUjetOffsetEnCorr->mez += jetcorr *
jet->p4().pz() / jetp;
359 sumNoPUjetOffsetEnCorr->sumet += jetcorr *
jet->p4().pt() / jetp;
361 jet->metSignObj().get_type(),
363 jet->metSignObj().get_phi(),
364 (
jet->offsetEnCorr() /
jet->p4().E()) *
jet->metSignObj().get_sigma_e(),
365 jet->metSignObj().get_sigma_tan());
366 metSignObjectsNoPUjetOffsetEnCorr.push_back(pfMEtSignObjectOffsetEnCorr);
369 metSignObjectsPUjets.push_back(
jet->metSignObj());
374 auto sumNoPUunclChargedCands = std::make_unique<CommonMETData>();
376 std::vector<metsig::SigInputObj> metSignObjectsNoPUunclChargedCands;
377 auto sumPUunclChargedCands = std::make_unique<CommonMETData>();
379 std::vector<metsig::SigInputObj> metSignObjectsPUunclChargedCands;
380 auto sumUnclNeutralCands = std::make_unique<CommonMETData>();
382 std::vector<metsig::SigInputObj> metSignObjectsUnclNeutralCands;
383 for (reco::PUSubMETCandInfoCollection::const_iterator pfCandidate = pfCandidates_cleaned.begin();
384 pfCandidate != pfCandidates_cleaned.end();
386 if (pfCandidate->passesLooseJetId()) {
387 if (!pfCandidate->isWithinJet()) {
390 metSignObjectsNoPUunclChargedCands.push_back(pfCandidate->metSignObj());
393 metSignObjectsPUunclChargedCands.push_back(pfCandidate->metSignObj());
396 metSignObjectsUnclNeutralCands.push_back(pfCandidate->metSignObj());
404 auto type0Correction_output = std::make_unique<CommonMETData>();
406 type0Correction_output->mex = type0Correction_input->
mex;
407 type0Correction_output->mey = type0Correction_input->
mey;
419 double noPileUpScaleFactor =
420 (sumPUunclChargedCands->sumet > 0.)
421 ? (sumPUunclChargedCands->sumet / (sumNoPUunclChargedCands->sumet + sumPUunclChargedCands->sumet))
423 LogDebug(
"produce") <<
"noPileUpScaleFactor = " << noPileUpScaleFactor << std::endl;
425 double noPileUpMEtPx =
426 -(sumLeptons->mex + sumNoPUjets->mex + sumNoPUunclChargedCands->mex +
427 noPileUpScaleFactor *
432 noPileUpMEtPx -= (noPileUpScaleFactor *
sfLeptonIsoCones_ * sumLeptonIsoCones->mex);
435 double noPileUpMEtPy =
436 -(sumLeptons->mey + sumNoPUjets->mey + sumNoPUunclChargedCands->mey +
437 noPileUpScaleFactor *
442 noPileUpMEtPy -= (noPileUpScaleFactor *
sfLeptonIsoCones_ * sumLeptonIsoCones->mey);
445 double noPileUpMEtPt =
sqrt(noPileUpMEtPx * noPileUpMEtPx + noPileUpMEtPy * noPileUpMEtPy);
449 noPileUpMEt.setP4(noPileUpMEtP4);
452 std::vector<metsig::SigInputObj> metSignObjects_scaled;
462 metSignObjectsPUunclChargedCands,
467 metSignObjectsUnclNeutralCands,
472 noPileUpMEt.setSignificanceMatrix(pfMEtCov_recomputed);
474 LogDebug(
"produce") <<
"<NoPileUpPFMEtProducer::produce>:" << std::endl
476 <<
" PFMET: Pt = " << pfMEt_original.
pt() <<
", phi = " << pfMEt_original.
phi() <<
" " 477 <<
"(Px = " << pfMEt_original.
px() <<
", Py = " << pfMEt_original.
py() <<
")" << std::endl
478 <<
" Cov:" << std::endl
479 <<
" " << pfMEtCov(0, 0) <<
" " << pfMEtCov(0, 1) <<
"\n " << pfMEtCov(1, 0) <<
" " 480 << pfMEtCov(1, 1) << std::endl
481 <<
" no-PU MET: Pt = " << noPileUpMEt.pt() <<
", phi = " << noPileUpMEt.phi() <<
" " 482 <<
"(Px = " << noPileUpMEt.px() <<
", Py = " << noPileUpMEt.py() <<
")" << std::endl
483 <<
" Cov:" << std::endl
484 <<
" " << (noPileUpMEt.getSignificanceMatrix())(0, 0) <<
" " 485 << (noPileUpMEt.getSignificanceMatrix())(0, 1) << std::endl
486 << (noPileUpMEt.getSignificanceMatrix())(1, 0) <<
" " 487 << (noPileUpMEt.getSignificanceMatrix())(1, 1) << std::endl;
490 auto noPileUpMEtCollection = std::make_unique<reco::PFMETCollection>();
491 noPileUpMEtCollection->push_back(noPileUpMEt);
506 evt.
put(std::make_unique<double>(noPileUpScaleFactor),
"sfNoPU");
double sfNoPUunclChargedCands_
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
reco::METCovMatrix getSignificanceMatrix(void) const
double pt() const final
transverse momentum
edm::EDGetTokenT< reco::PFMETCollection > srcMEt_
std::string sfLeptonIsoConesName_
edm::EDGetTokenT< reco::PUSubMETCandInfoCollection > srcPFCandInfoLeptonMatch_
ROOT::Math::SMatrix< double, 2 > METCovMatrix
bool getByToken(EDGetToken token, Handle< PROD > &result) const
double sfType0Correction_
std::string sfNoPUjetsName_
reco::PUSubMETCandInfoCollection cleanJets(const reco::PUSubMETCandInfoCollection &, const std::vector< reco::Candidate::LorentzVector > &, double, bool)
edm::EDGetTokenT< reco::PUSubMETCandInfoCollection > srcJetInfo_
std::vector< edm::EDGetTokenT< reco::CandidateView > > srcLeptons_
double px() const final
x coordinate of momentum vector
int findBestMatchingLepton(const std::vector< reco::Candidate::LorentzVector > &leptons, const reco::Candidate::LorentzVector &p4_ref)
void addToCommonMETData(CommonMETData &metData, const reco::Candidate::LorentzVector &p4)
std::string sfPUunclChargedCandsName_
std::string sfPUjetsName_
std::string sfUnclNeutralCandsName_
Abs< T >::type abs(const T &t)
std::vector< reco::PUSubMETCandInfo > PUSubMETCandInfoCollection
double py() const final
y coordinate of momentum vector
double sfNoPUjetOffsetEnCorr_
edm::EDGetTokenT< CorrMETData > srcType0Correction_
void finalizeCommonMETData(CommonMETData &metData)
std::string sfNoPUunclChargedCandsName_
std::string sfType0CorrectionName_
edm::EDGetTokenT< reco::PUSubMETCandInfoCollection > srcJetInfoLeptonMatch_
std::string sfLeptonsName_
math::XYZTLorentzVector LorentzVector
Lorentz vector.
edm::EDGetTokenT< reco::PUSubMETCandInfoCollection > srcPFCandInfo_
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
reco::PUSubMETCandInfoCollection cleanPFCandidates(const reco::PUSubMETCandInfoCollection &, const std::vector< reco::Candidate::LorentzVector > &, double, bool)
metsig::SigInputObj compResolution(const T *particle) const
const_iterator begin() const
NoPileUpMEtUtilities utils_
void scaleAndAddPFMEtSignObjects(std::vector< metsig::SigInputObj > &metSignObjects_scaled, const std::vector< metsig::SigInputObj > &metSignObjects, double sf, double sfMin, double sfMax)
double phi() const final
momentum azimuthal angle
const_iterator end() const
reco::METCovMatrix computePFMEtSignificance(const std::vector< metsig::SigInputObj > &metSignObjects)
double sfPUunclChargedCands_
void initializeCommonMETData(CommonMETData &metData)
double sfUnclNeutralCands_
PFMEtSignInterfaceBase * pfMEtSignInterface_
std::string sfNoPUjetOffsetEnCorrName_