NoPileUpPFMEtProducer.cc

Go to the documentation of this file.

#include "RecoMET/METPUSubtraction/plugins/NoPileUpPFMEtProducer.h"

#include "FWCore/Utilities/interface/Exception.h"

//#include "DataFormats/METReco/interface/PFMEtSignCovMatrix.h" //never used so far
#include "RecoMET/METAlgorithms/interface/significanceAlgo.h"
#include "DataFormats/METReco/interface/SigInputObj.h"

#include <cmath>

const double defaultPFMEtResolutionX = 10.;
const double defaultPFMEtResolutionY = 10.;

const double epsilon = 1.e-9;

NoPileUpPFMEtProducer::NoPileUpPFMEtProducer(const edm::ParameterSet& cfg)
    : moduleLabel_(cfg.getParameter<std::string>("@module_label")) {
  srcMEt_ = consumes<reco::PFMETCollection>(cfg.getParameter<edm::InputTag>("srcMEt"));
  srcMEtCov_ = edm::InputTag();  //MM, disabled for the moment until we really need it
  //( cfg.exists("srcMEtCov") ) ?
  //  consumes<edm::Handle<> >(cfg.getParameter<edm::InputTag>("srcMEtCov")) : edm::InputTag();
  srcJetInfo_ = consumes<reco::PUSubMETCandInfoCollection>(cfg.getParameter<edm::InputTag>("srcPUSubMETDataJet"));
  srcJetInfoLeptonMatch_ =
      consumes<reco::PUSubMETCandInfoCollection>(cfg.getParameter<edm::InputTag>("srcPUSubMETDataJetLeptonMatch"));
  srcPFCandInfo_ =
      consumes<reco::PUSubMETCandInfoCollection>(cfg.getParameter<edm::InputTag>("srcPUSubMETDataPFCands"));
  srcPFCandInfoLeptonMatch_ =
      consumes<reco::PUSubMETCandInfoCollection>(cfg.getParameter<edm::InputTag>("srcPUSubMETDataPFCandsLeptonMatch"));
  vInputTag srcLeptonsTags = cfg.getParameter<vInputTag>("srcLeptons");
  for (vInputTag::const_iterator it = srcLeptonsTags.begin(); it != srcLeptonsTags.end(); it++) {
    srcLeptons_.push_back(consumes<edm::View<reco::Candidate> >(*it));
  }

  srcType0Correction_ = consumes<CorrMETData>(cfg.getParameter<edm::InputTag>("srcType0Correction"));

  sfNoPUjets_ = cfg.getParameter<double>("sfNoPUjets");
  sfNoPUjetOffsetEnCorr_ = cfg.getParameter<double>("sfNoPUjetOffsetEnCorr");
  sfPUjets_ = cfg.getParameter<double>("sfPUjets");
  sfNoPUunclChargedCands_ = cfg.getParameter<double>("sfNoPUunclChargedCands");
  sfPUunclChargedCands_ = cfg.getParameter<double>("sfPUunclChargedCands");
  sfUnclNeutralCands_ = cfg.getParameter<double>("sfUnclNeutralCands");
  sfType0Correction_ = cfg.getParameter<double>("sfType0Correction");
  sfLeptonIsoCones_ = cfg.getParameter<double>("sfLeptonIsoCones");

  pfMEtSignInterface_ = new PFMEtSignInterfaceBase(cfg.getParameter<edm::ParameterSet>("resolution"));
  sfMEtCovMin_ = cfg.getParameter<double>("sfMEtCovMin");
  sfMEtCovMax_ = cfg.getParameter<double>("sfMEtCovMax");

  saveInputs_ = (cfg.exists("saveInputs")) ? cfg.getParameter<bool>("saveInputs") : false;

  verbosity_ = (cfg.exists("verbosity")) ? cfg.getParameter<int>("verbosity") : 0;

  produces<reco::PFMETCollection>();

  sfLeptonsName_ = "sumLeptons";
  sfNoPUjetsName_ = "sumNoPUjets";
  sfNoPUjetOffsetEnCorrName_ = "sumNoPUjetOffsetEnCorr";
  sfPUjetsName_ = "sumPUjets";
  sfNoPUunclChargedCandsName_ = "sumNoPUunclChargedCands";
  sfPUunclChargedCandsName_ = "sumPUunclChargedCands";
  sfUnclNeutralCandsName_ = "sumUnclNeutralCands";
  sfType0CorrectionName_ = "type0Correction";
  sfLeptonIsoConesName_ = "sumLeptonIsoCones";

  if (saveInputs_) {
    produces<CommonMETData>(sfLeptonsName_);
    produces<CommonMETData>(sfNoPUjetsName_);
    produces<CommonMETData>(sfNoPUjetOffsetEnCorrName_);
    produces<CommonMETData>(sfPUjetsName_);
    produces<CommonMETData>(sfNoPUunclChargedCandsName_);
    produces<CommonMETData>(sfPUunclChargedCandsName_);
    produces<CommonMETData>(sfUnclNeutralCandsName_);
    produces<CommonMETData>(sfType0CorrectionName_);
    produces<CommonMETData>(sfLeptonIsoConesName_);
  }
  produces<double>("sfNoPU");
}

NoPileUpPFMEtProducer::~NoPileUpPFMEtProducer() { delete pfMEtSignInterface_; }

void initializeCommonMETData(CommonMETData& metData) {
  metData.met = 0.;
  metData.mex = 0.;
  metData.mey = 0.;
  metData.mez = 0.;
  metData.sumet = 0.;
  metData.phi = 0.;
}

void addToCommonMETData(CommonMETData& metData, const reco::Candidate::LorentzVector& p4) {
  metData.mex += p4.px();
  metData.mey += p4.py();
  metData.mez += p4.pz();
  metData.sumet += p4.pt();
}

void finalizeCommonMETData(CommonMETData& metData) {
  metData.met = sqrt(metData.mex * metData.mex + metData.mey * metData.mey);
  metData.phi = atan2(metData.mey, metData.mex);
}

int findBestMatchingLepton(const std::vector<reco::Candidate::LorentzVector>& leptons,
                           const reco::Candidate::LorentzVector& p4_ref) {
  int leptonIdx_dR2min = -1;
  double dR2min = 1.e+3;
  int leptonIdx = 0;
  for (std::vector<reco::Candidate::LorentzVector>::const_iterator lepton = leptons.begin(); lepton != leptons.end();
       ++lepton) {
    double dR2 = deltaR2(*lepton, p4_ref);
    if (leptonIdx_dR2min == -1 || dR2 < dR2min) {
      leptonIdx_dR2min = leptonIdx;
      dR2min = dR2;
    }
    ++leptonIdx;
  }
  assert(leptonIdx_dR2min >= 0 && leptonIdx_dR2min < (int)leptons.size());
  return leptonIdx_dR2min;
}

void scaleAndAddPFMEtSignObjects(std::vector<metsig::SigInputObj>& metSignObjects_scaled,
                                 const std::vector<metsig::SigInputObj>& metSignObjects,
                                 double sf,
                                 double sfMin,
                                 double sfMax) {
  double sf_value = sf;
  if (sf_value > sfMax)
    sf_value = sfMax;
  if (sf_value < sfMin)
    sf_value = sfMin;
  for (std::vector<metsig::SigInputObj>::const_iterator metSignObject = metSignObjects.begin();
       metSignObject != metSignObjects.end();
       ++metSignObject) {
    metsig::SigInputObj metSignObject_scaled;
    metSignObject_scaled.set(metSignObject->get_type(),
                             sf_value * metSignObject->get_energy(),
                             metSignObject->get_phi(),
                             sf_value * metSignObject->get_sigma_e(),
                             metSignObject->get_sigma_tan());
    metSignObjects_scaled.push_back(metSignObject_scaled);
  }
}

reco::METCovMatrix computePFMEtSignificance(const std::vector<metsig::SigInputObj>& metSignObjects) {
  reco::METCovMatrix pfMEtCov;
  if (metSignObjects.size() >= 2) {
    metsig::significanceAlgo pfMEtSignAlgorithm;
    pfMEtSignAlgorithm.addObjects(metSignObjects);
    pfMEtCov = pfMEtSignAlgorithm.getSignifMatrix();
  }

  double det = 0;
  pfMEtCov.Det(det);
  if (std::abs(det) < epsilon) {
    edm::LogWarning("computePFMEtSignificance") << "Inversion of PFMEt covariance matrix failed, det = " << det
                                                << " --> replacing covariance matrix by resolution defaults !!";
    pfMEtCov(0, 0) = defaultPFMEtResolutionX * defaultPFMEtResolutionX;
    pfMEtCov(0, 1) = 0.;
    pfMEtCov(1, 0) = 0.;
    pfMEtCov(1, 1) = defaultPFMEtResolutionY * defaultPFMEtResolutionY;
  }

  return pfMEtCov;
}

void printP4(const std::string& label_part1, int idx, const std::string& label_part2, const reco::Candidate& candidate) {
  std::cout << label_part1 << " #" << idx << label_part2 << ": Pt = " << candidate.pt() << ", eta = " << candidate.eta()
            << ", phi = " << candidate.phi() << " (charge = " << candidate.charge() << ")" << std::endl;
}

void printCommonMETData(const std::string& label, const CommonMETData& metData) {
  std::cout << label << ": Px = " << metData.mex << ", Py = " << metData.mey << ", sumEt = " << metData.sumet
            << std::endl;
}

void printMVAMEtJetInfo(const std::string& label, int idx, const reco::PUSubMETCandInfo& jet) {
  std::cout << label << " #" << idx << " (";
  if (jet.type() == reco::PUSubMETCandInfo::kHS)
    std::cout << "no-PU";
  else if (jet.type() == reco::PUSubMETCandInfo::kPU)
    std::cout << "PU";
  std::cout << "): Pt = " << jet.p4().pt() << ", eta = " << jet.p4().eta() << ", phi = " << jet.p4().phi();
  std::cout << " id. flags: anti-noise = " << jet.passesLooseJetId() << std::endl;
  std::cout << std::endl;
}

void printMVAMEtPFCandInfo(const std::string& label, int idx, const reco::PUSubMETCandInfo& pfCand) {
  std::cout << label << " #" << idx << " (";
  if (pfCand.type() == reco::PUSubMETCandInfo::kChHS)
    std::cout << "no-PU charged";
  else if (pfCand.type() == reco::PUSubMETCandInfo::kChPU)
    std::cout << "PU charged";
  else if (pfCand.type() == reco::PUSubMETCandInfo::kNeutral)
    std::cout << "neutral";
  std::cout << "): Pt = " << pfCand.p4().pt() << ", eta = " << pfCand.p4().eta() << ", phi = " << pfCand.p4().phi();
  std::string isWithinJet_string;
  if (pfCand.isWithinJet())
    isWithinJet_string = "true";
  else
    isWithinJet_string = "false";
  std::cout << " (isWithinJet = " << isWithinJet_string << ")";
  if (pfCand.isWithinJet())
    std::cout << " Jet id. flags: anti-noise = " << pfCand.passesLooseJetId() << std::endl;
  std::cout << std::endl;
}

void NoPileUpPFMEtProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
  LogDebug("produce") << " moduleLabel = " << moduleLabel_ << std::endl;

  // get original MET
  edm::Handle<reco::PFMETCollection> pfMETs;
  evt.getByToken(srcMEt_, pfMETs);
  if (!(pfMETs->size() == 1))
    throw cms::Exception("NoPileUpPFMEtProducer::produce") << "Failed to find unique MET object !!\n";
  const reco::PFMET& pfMEt_original = pfMETs->front();

  // get MET covariance matrix
  reco::METCovMatrix pfMEtCov;
  if (!srcMEtCov_.label().empty()) {
    //MM manual bypass to pfMET as this case has neer been presented
    // edm::Handle<PFMEtSignCovMatrix> pfMEtCovHandle;
    // evt.getByToken(srcMEtCov_, pfMEtCovHandle);
    // pfMEtCov = (*pfMEtCovHandle);
    pfMEtCov = pfMEt_original.getSignificanceMatrix();
  } else {
    pfMEtCov = pfMEt_original.getSignificanceMatrix();
  }

  // get lepton momenta
  std::vector<reco::Candidate::LorentzVector> leptons;
  std::vector<metsig::SigInputObj> metSignObjectsLeptons;
  reco::Candidate::LorentzVector sumLeptonP4s;
  for (std::vector<edm::EDGetTokenT<edm::View<reco::Candidate> > >::const_iterator srcLeptons_i = srcLeptons_.begin();
       srcLeptons_i != srcLeptons_.end();
       ++srcLeptons_i) {
    edm::Handle<reco::CandidateView> leptons_i;
    evt.getByToken(*srcLeptons_i, leptons_i);
    int leptonIdx = 0;
    for (reco::CandidateView::const_iterator lepton = leptons_i->begin(); lepton != leptons_i->end(); ++lepton) {
      leptons.push_back(lepton->p4());
      metSignObjectsLeptons.push_back(pfMEtSignInterface_->compResolution(&(*lepton)));
      sumLeptonP4s += lepton->p4();
      ++leptonIdx;
    }
  }
  LogDebug("produce") << " sum(leptons): Pt = " << sumLeptonP4s.pt() << ", eta = " << sumLeptonP4s.eta()
                      << ", phi = " << sumLeptonP4s.phi() << ","
                      << " mass = " << sumLeptonP4s.mass() << std::endl;

  // get jet and PFCandidate information
  edm::Handle<reco::PUSubMETCandInfoCollection> jets;
  evt.getByToken(srcJetInfo_, jets);
  edm::Handle<reco::PUSubMETCandInfoCollection> jetsLeptonMatch;
  evt.getByToken(srcJetInfoLeptonMatch_, jetsLeptonMatch);
  edm::Handle<reco::PUSubMETCandInfoCollection> pfCandidates;
  evt.getByToken(srcPFCandInfo_, pfCandidates);
  edm::Handle<reco::PUSubMETCandInfoCollection> pfCandidatesLeptonMatch;
  evt.getByToken(srcPFCandInfoLeptonMatch_, pfCandidatesLeptonMatch);

  reco::PUSubMETCandInfoCollection jets_leptons = utils_.cleanJets(*jetsLeptonMatch, leptons, 0.5, true);
  reco::PUSubMETCandInfoCollection pfCandidates_leptons =
      utils_.cleanPFCandidates(*pfCandidatesLeptonMatch, leptons, 0.3, true);
  std::vector<CommonMETData> sumJetsPlusPFCandidates_leptons(leptons.size());
  for (std::vector<CommonMETData>::iterator sumJetsPlusPFCandidates = sumJetsPlusPFCandidates_leptons.begin();
       sumJetsPlusPFCandidates != sumJetsPlusPFCandidates_leptons.end();
       ++sumJetsPlusPFCandidates) {
    initializeCommonMETData(*sumJetsPlusPFCandidates);
  }
  for (reco::PUSubMETCandInfoCollection::const_iterator jet = jets_leptons.begin(); jet != jets_leptons.end(); ++jet) {
    int leptonIdx_dRmin = findBestMatchingLepton(leptons, jet->p4());
    assert(leptonIdx_dRmin >= 0 && leptonIdx_dRmin < (int)sumJetsPlusPFCandidates_leptons.size());

    LogDebug("produce") << "jet-to-lepton match:"
                        << " jetPt = " << jet->p4().pt() << ", jetEta = " << jet->p4().eta()
                        << ", jetPhi = " << jet->p4().phi() << " leptonPt = " << leptons[leptonIdx_dRmin].pt()
                        << ", leptonEta = " << leptons[leptonIdx_dRmin].eta()
                        << ", leptonPhi = " << leptons[leptonIdx_dRmin].phi() << std::endl;

    sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mex += jet->p4().px();
    sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mey += jet->p4().py();
    sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].sumet += jet->p4().pt();
  }
  for (reco::PUSubMETCandInfoCollection::const_iterator pfCandidate = pfCandidates_leptons.begin();
       pfCandidate != pfCandidates_leptons.end();
       ++pfCandidate) {
    bool isWithinJet_lepton = false;
    if (pfCandidate->isWithinJet()) {
      for (reco::PUSubMETCandInfoCollection::const_iterator jet = jets_leptons.begin(); jet != jets_leptons.end();
           ++jet) {
        double dR2 = deltaR2(pfCandidate->p4(), jet->p4());
        if (dR2 < 0.5 * 0.5)
          isWithinJet_lepton = true;
      }
    }
    if (!isWithinJet_lepton) {
      int leptonIdx_dRmin = findBestMatchingLepton(leptons, pfCandidate->p4());
      assert(leptonIdx_dRmin >= 0 && leptonIdx_dRmin < (int)sumJetsPlusPFCandidates_leptons.size());
      LogDebug("produce") << "pfCandidate-to-lepton match:"
                          << " pfCandidatePt = " << pfCandidate->p4().pt()
                          << ", pfCandidateEta = " << pfCandidate->p4().eta()
                          << ", pfCandidatePhi = " << pfCandidate->p4().phi()
                          << " leptonPt = " << leptons[leptonIdx_dRmin].pt()
                          << ", leptonEta = " << leptons[leptonIdx_dRmin].eta()
                          << ", leptonPhi = " << leptons[leptonIdx_dRmin].phi() << std::endl;

      sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mex += pfCandidate->p4().px();
      sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mey += pfCandidate->p4().py();
      sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].sumet += pfCandidate->p4().pt();
    } else {
      LogDebug("produce") << " pfCandidate is within jet --> skipping." << std::endl;
    }
  }
  auto sumLeptons = std::make_unique<CommonMETData>();
  initializeCommonMETData(*sumLeptons);
  auto sumLeptonIsoCones = std::make_unique<CommonMETData>();
  initializeCommonMETData(*sumLeptonIsoCones);
  int leptonIdx = 0;
  for (std::vector<CommonMETData>::iterator sumJetsPlusPFCandidates = sumJetsPlusPFCandidates_leptons.begin();
       sumJetsPlusPFCandidates != sumJetsPlusPFCandidates_leptons.end();
       ++sumJetsPlusPFCandidates) {
    if (sumJetsPlusPFCandidates->sumet > leptons[leptonIdx].pt()) {
      double leptonEnFrac = leptons[leptonIdx].pt() / sumJetsPlusPFCandidates->sumet;
      assert(leptonEnFrac >= 0.0 && leptonEnFrac <= 1.0);
      sumLeptons->mex += (leptonEnFrac * sumJetsPlusPFCandidates->mex);
      sumLeptons->mey += (leptonEnFrac * sumJetsPlusPFCandidates->mey);
      sumLeptons->sumet += (leptonEnFrac * sumJetsPlusPFCandidates->sumet);
      double leptonIsoConeEnFrac = 1.0 - leptonEnFrac;
      assert(leptonIsoConeEnFrac >= 0.0 && leptonIsoConeEnFrac <= 1.0);
      sumLeptonIsoCones->mex += (leptonIsoConeEnFrac * sumJetsPlusPFCandidates->mex);
      sumLeptonIsoCones->mey += (leptonIsoConeEnFrac * sumJetsPlusPFCandidates->mey);
      sumLeptonIsoCones->sumet += (leptonIsoConeEnFrac * sumJetsPlusPFCandidates->sumet);
    } else {
      sumLeptons->mex += sumJetsPlusPFCandidates->mex;
      sumLeptons->mey += sumJetsPlusPFCandidates->mey;
      sumLeptons->sumet += sumJetsPlusPFCandidates->sumet;
    }
    ++leptonIdx;
  }

  reco::PUSubMETCandInfoCollection jets_cleaned = utils_.cleanJets(*jets, leptons, 0.5, false);
  reco::PUSubMETCandInfoCollection pfCandidates_cleaned = utils_.cleanPFCandidates(*pfCandidates, leptons, 0.3, false);

  auto sumNoPUjets = std::make_unique<CommonMETData>();
  initializeCommonMETData(*sumNoPUjets);
  std::vector<metsig::SigInputObj> metSignObjectsNoPUjets;
  auto sumNoPUjetOffsetEnCorr = std::make_unique<CommonMETData>();
  initializeCommonMETData(*sumNoPUjetOffsetEnCorr);
  std::vector<metsig::SigInputObj> metSignObjectsNoPUjetOffsetEnCorr;
  auto sumPUjets = std::make_unique<CommonMETData>();
  initializeCommonMETData(*sumPUjets);
  std::vector<metsig::SigInputObj> metSignObjectsPUjets;
  int jetIdx = 0;
  for (reco::PUSubMETCandInfoCollection::const_iterator jet = jets_cleaned.begin(); jet != jets_cleaned.end(); ++jet) {
    if (jet->passesLooseJetId()) {
      if (jet->type() == reco::PUSubMETCandInfo::kHS) {
        addToCommonMETData(*sumNoPUjets, jet->p4());
        metSignObjectsNoPUjets.push_back(jet->metSignObj());
        float jetp = jet->p4().P();
        float jetcorr = jet->offsetEnCorr();
        sumNoPUjetOffsetEnCorr->mex += jetcorr * jet->p4().px() / jetp;
        sumNoPUjetOffsetEnCorr->mey += jetcorr * jet->p4().py() / jetp;
        sumNoPUjetOffsetEnCorr->mez += jetcorr * jet->p4().pz() / jetp;
        sumNoPUjetOffsetEnCorr->sumet += jetcorr * jet->p4().pt() / jetp;
        metsig::SigInputObj pfMEtSignObjectOffsetEnCorr(
            jet->metSignObj().get_type(),
            jet->offsetEnCorr(),
            jet->metSignObj().get_phi(),
            (jet->offsetEnCorr() / jet->p4().E()) * jet->metSignObj().get_sigma_e(),
            jet->metSignObj().get_sigma_tan());
        metSignObjectsNoPUjetOffsetEnCorr.push_back(pfMEtSignObjectOffsetEnCorr);
      } else {
        addToCommonMETData(*sumPUjets, jet->p4());
        metSignObjectsPUjets.push_back(jet->metSignObj());
      }
    }
    ++jetIdx;
  }

  auto sumNoPUunclChargedCands = std::make_unique<CommonMETData>();
  initializeCommonMETData(*sumNoPUunclChargedCands);
  std::vector<metsig::SigInputObj> metSignObjectsNoPUunclChargedCands;
  auto sumPUunclChargedCands = std::make_unique<CommonMETData>();
  initializeCommonMETData(*sumPUunclChargedCands);
  std::vector<metsig::SigInputObj> metSignObjectsPUunclChargedCands;
  auto sumUnclNeutralCands = std::make_unique<CommonMETData>();
  initializeCommonMETData(*sumUnclNeutralCands);
  std::vector<metsig::SigInputObj> metSignObjectsUnclNeutralCands;
  int pfCandIdx = 0;
  for (reco::PUSubMETCandInfoCollection::const_iterator pfCandidate = pfCandidates_cleaned.begin();
       pfCandidate != pfCandidates_cleaned.end();
       ++pfCandidate) {
    if (pfCandidate->passesLooseJetId()) {
      if (!pfCandidate->isWithinJet()) {
        if (pfCandidate->type() == reco::PUSubMETCandInfo::kChHS) {
          addToCommonMETData(*sumNoPUunclChargedCands, pfCandidate->p4());
          metSignObjectsNoPUunclChargedCands.push_back(pfCandidate->metSignObj());
        } else if (pfCandidate->type() == reco::PUSubMETCandInfo::kChPU) {
          addToCommonMETData(*sumPUunclChargedCands, pfCandidate->p4());
          metSignObjectsPUunclChargedCands.push_back(pfCandidate->metSignObj());
        } else if (pfCandidate->type() == reco::PUSubMETCandInfo::kNeutral) {
          addToCommonMETData(*sumUnclNeutralCands, pfCandidate->p4());
          metSignObjectsUnclNeutralCands.push_back(pfCandidate->metSignObj());
        }
      }
    }
    ++pfCandIdx;
  }

  edm::Handle<CorrMETData> type0Correction_input;
  evt.getByToken(srcType0Correction_, type0Correction_input);
  auto type0Correction_output = std::make_unique<CommonMETData>();
  initializeCommonMETData(*type0Correction_output);
  type0Correction_output->mex = type0Correction_input->mex;
  type0Correction_output->mey = type0Correction_input->mey;

  finalizeCommonMETData(*sumLeptons);
  finalizeCommonMETData(*sumNoPUjetOffsetEnCorr);
  finalizeCommonMETData(*sumNoPUjets);
  finalizeCommonMETData(*sumPUjets);
  finalizeCommonMETData(*sumNoPUunclChargedCands);
  finalizeCommonMETData(*sumPUunclChargedCands);
  finalizeCommonMETData(*sumUnclNeutralCands);
  finalizeCommonMETData(*type0Correction_output);
  finalizeCommonMETData(*sumLeptonIsoCones);

  double noPileUpScaleFactor =
      (sumPUunclChargedCands->sumet > 0.)
          ? (sumPUunclChargedCands->sumet / (sumNoPUunclChargedCands->sumet + sumPUunclChargedCands->sumet))
          : 1.;
  LogDebug("produce") << "noPileUpScaleFactor = " << noPileUpScaleFactor << std::endl;

  double noPileUpMEtPx =
      -(sumLeptons->mex + sumNoPUjets->mex + sumNoPUunclChargedCands->mex +
        noPileUpScaleFactor *
            (sfNoPUjetOffsetEnCorr_ * sumNoPUjetOffsetEnCorr->mex + sfUnclNeutralCands_ * sumUnclNeutralCands->mex +
             sfPUunclChargedCands_ * sumPUunclChargedCands->mex + sfPUjets_ * sumPUjets->mex)) +
      noPileUpScaleFactor * sfType0Correction_ * type0Correction_output->mex;
  if (sfLeptonIsoCones_ >= 0.)
    noPileUpMEtPx -= (noPileUpScaleFactor * sfLeptonIsoCones_ * sumLeptonIsoCones->mex);
  else
    noPileUpMEtPx -= (std::abs(sfLeptonIsoCones_) * sumLeptonIsoCones->mex);
  double noPileUpMEtPy =
      -(sumLeptons->mey + sumNoPUjets->mey + sumNoPUunclChargedCands->mey +
        noPileUpScaleFactor *
            (sfNoPUjetOffsetEnCorr_ * sumNoPUjetOffsetEnCorr->mey + sfUnclNeutralCands_ * sumUnclNeutralCands->mey +
             sfPUunclChargedCands_ * sumPUunclChargedCands->mey + sfPUjets_ * sumPUjets->mey)) +
      noPileUpScaleFactor * sfType0Correction_ * type0Correction_output->mey;
  if (sfLeptonIsoCones_ >= 0.)
    noPileUpMEtPy -= (noPileUpScaleFactor * sfLeptonIsoCones_ * sumLeptonIsoCones->mey);
  else
    noPileUpMEtPy -= (std::abs(sfLeptonIsoCones_) * sumLeptonIsoCones->mey);
  double noPileUpMEtPt = sqrt(noPileUpMEtPx * noPileUpMEtPx + noPileUpMEtPy * noPileUpMEtPy);
  reco::Candidate::LorentzVector noPileUpMEtP4(noPileUpMEtPx, noPileUpMEtPy, 0., noPileUpMEtPt);

  reco::PFMET noPileUpMEt(pfMEt_original);
  noPileUpMEt.setP4(noPileUpMEtP4);
  //noPileUpMEt.setSignificanceMatrix(pfMEtCov);

  std::vector<metsig::SigInputObj> metSignObjects_scaled;
  scaleAndAddPFMEtSignObjects(metSignObjects_scaled, metSignObjectsLeptons, 1.0, sfMEtCovMin_, sfMEtCovMax_);
  scaleAndAddPFMEtSignObjects(
      metSignObjects_scaled, metSignObjectsNoPUjetOffsetEnCorr, sfNoPUjetOffsetEnCorr_, sfMEtCovMin_, sfMEtCovMax_);
  scaleAndAddPFMEtSignObjects(metSignObjects_scaled, metSignObjectsNoPUjets, sfNoPUjets_, sfMEtCovMin_, sfMEtCovMax_);
  scaleAndAddPFMEtSignObjects(
      metSignObjects_scaled, metSignObjectsPUjets, noPileUpScaleFactor * sfPUjets_, sfMEtCovMin_, sfMEtCovMax_);
  scaleAndAddPFMEtSignObjects(
      metSignObjects_scaled, metSignObjectsNoPUunclChargedCands, sfNoPUunclChargedCands_, sfMEtCovMin_, sfMEtCovMax_);
  scaleAndAddPFMEtSignObjects(metSignObjects_scaled,
                              metSignObjectsPUunclChargedCands,
                              noPileUpScaleFactor * sfPUunclChargedCands_,
                              sfMEtCovMin_,
                              sfMEtCovMax_);
  scaleAndAddPFMEtSignObjects(metSignObjects_scaled,
                              metSignObjectsUnclNeutralCands,
                              noPileUpScaleFactor * sfUnclNeutralCands_,
                              sfMEtCovMin_,
                              sfMEtCovMax_);
  reco::METCovMatrix pfMEtCov_recomputed = computePFMEtSignificance(metSignObjects_scaled);
  noPileUpMEt.setSignificanceMatrix(pfMEtCov_recomputed);

  LogDebug("produce") << "<NoPileUpPFMEtProducer::produce>:" << std::endl
                      << " moduleLabel = " << moduleLabel_ << std::endl
                      << " PFMET: Pt = " << pfMEt_original.pt() << ", phi = " << pfMEt_original.phi() << " "
                      << "(Px = " << pfMEt_original.px() << ", Py = " << pfMEt_original.py() << ")" << std::endl
                      << " Cov:" << std::endl
                      << " " << pfMEtCov(0, 0) << "  " << pfMEtCov(0, 1) << "\n " << pfMEtCov(1, 0) << "  "
                      << pfMEtCov(1, 1) << std::endl
                      << " no-PU MET: Pt = " << noPileUpMEt.pt() << ", phi = " << noPileUpMEt.phi() << " "
                      << "(Px = " << noPileUpMEt.px() << ", Py = " << noPileUpMEt.py() << ")" << std::endl
                      << " Cov:" << std::endl
                      << " " << (noPileUpMEt.getSignificanceMatrix())(0, 0) << "  "
                      << (noPileUpMEt.getSignificanceMatrix())(0, 1) << std::endl
                      << (noPileUpMEt.getSignificanceMatrix())(1, 0) << "  "
                      << (noPileUpMEt.getSignificanceMatrix())(1, 1) << std::endl;

  // add no-PU MET object to the event
  auto noPileUpMEtCollection = std::make_unique<reco::PFMETCollection>();
  noPileUpMEtCollection->push_back(noPileUpMEt);

  evt.put(std::move(noPileUpMEtCollection));
  if (saveInputs_) {
    evt.put(std::move(sumLeptons), sfLeptonsName_);
    evt.put(std::move(sumNoPUjetOffsetEnCorr), sfNoPUjetOffsetEnCorrName_);
    evt.put(std::move(sumNoPUjets), sfNoPUjetsName_);
    evt.put(std::move(sumPUjets), sfPUjetsName_);
    evt.put(std::move(sumNoPUunclChargedCands), sfNoPUunclChargedCandsName_);
    evt.put(std::move(sumPUunclChargedCands), sfPUunclChargedCandsName_);
    evt.put(std::move(sumUnclNeutralCands), sfUnclNeutralCandsName_);
    evt.put(std::move(type0Correction_output), sfType0CorrectionName_);
    evt.put(std::move(sumLeptonIsoCones), sfLeptonIsoConesName_);
  }

  evt.put(std::make_unique<double>(noPileUpScaleFactor), "sfNoPU");
}

#include "FWCore/Framework/interface/MakerMacros.h"

DEFINE_FWK_MODULE(NoPileUpPFMEtProducer);