Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:25:37

0001 #include "RecoMET/METPUSubtraction/plugins/NoPileUpPFMEtProducer.h"
0002 
0003 #include "FWCore/Utilities/interface/Exception.h"
0004 
0005 //#include "DataFormats/METReco/interface/PFMEtSignCovMatrix.h" //never used so far
0006 #include "RecoMET/METAlgorithms/interface/significanceAlgo.h"
0007 #include "DataFormats/METReco/interface/SigInputObj.h"
0008 
0009 #include <cmath>
0010 
0011 const double defaultPFMEtResolutionX = 10.;
0012 const double defaultPFMEtResolutionY = 10.;
0013 
0014 const double epsilon = 1.e-9;
0015 
0016 NoPileUpPFMEtProducer::NoPileUpPFMEtProducer(const edm::ParameterSet& cfg)
0017     : moduleLabel_(cfg.getParameter<std::string>("@module_label")) {
0018   srcMEt_ = consumes<reco::PFMETCollection>(cfg.getParameter<edm::InputTag>("srcMEt"));
0019   srcMEtCov_ = edm::InputTag();  //MM, disabled for the moment until we really need it
0020   //( cfg.exists("srcMEtCov") ) ?
0021   //  consumes<edm::Handle<> >(cfg.getParameter<edm::InputTag>("srcMEtCov")) : edm::InputTag();
0022   srcJetInfo_ = consumes<reco::PUSubMETCandInfoCollection>(cfg.getParameter<edm::InputTag>("srcPUSubMETDataJet"));
0023   srcJetInfoLeptonMatch_ =
0024       consumes<reco::PUSubMETCandInfoCollection>(cfg.getParameter<edm::InputTag>("srcPUSubMETDataJetLeptonMatch"));
0025   srcPFCandInfo_ =
0026       consumes<reco::PUSubMETCandInfoCollection>(cfg.getParameter<edm::InputTag>("srcPUSubMETDataPFCands"));
0027   srcPFCandInfoLeptonMatch_ =
0028       consumes<reco::PUSubMETCandInfoCollection>(cfg.getParameter<edm::InputTag>("srcPUSubMETDataPFCandsLeptonMatch"));
0029   vInputTag srcLeptonsTags = cfg.getParameter<vInputTag>("srcLeptons");
0030   for (vInputTag::const_iterator it = srcLeptonsTags.begin(); it != srcLeptonsTags.end(); it++) {
0031     srcLeptons_.push_back(consumes<edm::View<reco::Candidate> >(*it));
0032   }
0033 
0034   srcType0Correction_ = consumes<CorrMETData>(cfg.getParameter<edm::InputTag>("srcType0Correction"));
0035 
0036   sfNoPUjets_ = cfg.getParameter<double>("sfNoPUjets");
0037   sfNoPUjetOffsetEnCorr_ = cfg.getParameter<double>("sfNoPUjetOffsetEnCorr");
0038   sfPUjets_ = cfg.getParameter<double>("sfPUjets");
0039   sfNoPUunclChargedCands_ = cfg.getParameter<double>("sfNoPUunclChargedCands");
0040   sfPUunclChargedCands_ = cfg.getParameter<double>("sfPUunclChargedCands");
0041   sfUnclNeutralCands_ = cfg.getParameter<double>("sfUnclNeutralCands");
0042   sfType0Correction_ = cfg.getParameter<double>("sfType0Correction");
0043   sfLeptonIsoCones_ = cfg.getParameter<double>("sfLeptonIsoCones");
0044 
0045   pfMEtSignInterface_ = new PFMEtSignInterfaceBase(cfg.getParameter<edm::ParameterSet>("resolution"));
0046   sfMEtCovMin_ = cfg.getParameter<double>("sfMEtCovMin");
0047   sfMEtCovMax_ = cfg.getParameter<double>("sfMEtCovMax");
0048 
0049   saveInputs_ = (cfg.exists("saveInputs")) ? cfg.getParameter<bool>("saveInputs") : false;
0050 
0051   verbosity_ = (cfg.exists("verbosity")) ? cfg.getParameter<int>("verbosity") : 0;
0052 
0053   produces<reco::PFMETCollection>();
0054 
0055   sfLeptonsName_ = "sumLeptons";
0056   sfNoPUjetsName_ = "sumNoPUjets";
0057   sfNoPUjetOffsetEnCorrName_ = "sumNoPUjetOffsetEnCorr";
0058   sfPUjetsName_ = "sumPUjets";
0059   sfNoPUunclChargedCandsName_ = "sumNoPUunclChargedCands";
0060   sfPUunclChargedCandsName_ = "sumPUunclChargedCands";
0061   sfUnclNeutralCandsName_ = "sumUnclNeutralCands";
0062   sfType0CorrectionName_ = "type0Correction";
0063   sfLeptonIsoConesName_ = "sumLeptonIsoCones";
0064 
0065   if (saveInputs_) {
0066     produces<CommonMETData>(sfLeptonsName_);
0067     produces<CommonMETData>(sfNoPUjetsName_);
0068     produces<CommonMETData>(sfNoPUjetOffsetEnCorrName_);
0069     produces<CommonMETData>(sfPUjetsName_);
0070     produces<CommonMETData>(sfNoPUunclChargedCandsName_);
0071     produces<CommonMETData>(sfPUunclChargedCandsName_);
0072     produces<CommonMETData>(sfUnclNeutralCandsName_);
0073     produces<CommonMETData>(sfType0CorrectionName_);
0074     produces<CommonMETData>(sfLeptonIsoConesName_);
0075   }
0076   produces<double>("sfNoPU");
0077 }
0078 
0079 NoPileUpPFMEtProducer::~NoPileUpPFMEtProducer() { delete pfMEtSignInterface_; }
0080 
0081 void initializeCommonMETData(CommonMETData& metData) {
0082   metData.met = 0.;
0083   metData.mex = 0.;
0084   metData.mey = 0.;
0085   metData.mez = 0.;
0086   metData.sumet = 0.;
0087   metData.phi = 0.;
0088 }
0089 
0090 void addToCommonMETData(CommonMETData& metData, const reco::Candidate::LorentzVector& p4) {
0091   metData.mex += p4.px();
0092   metData.mey += p4.py();
0093   metData.mez += p4.pz();
0094   metData.sumet += p4.pt();
0095 }
0096 
0097 void finalizeCommonMETData(CommonMETData& metData) {
0098   metData.met = sqrt(metData.mex * metData.mex + metData.mey * metData.mey);
0099   metData.phi = atan2(metData.mey, metData.mex);
0100 }
0101 
0102 int findBestMatchingLepton(const std::vector<reco::Candidate::LorentzVector>& leptons,
0103                            const reco::Candidate::LorentzVector& p4_ref) {
0104   int leptonIdx_dR2min = -1;
0105   double dR2min = 1.e+3;
0106   int leptonIdx = 0;
0107   for (std::vector<reco::Candidate::LorentzVector>::const_iterator lepton = leptons.begin(); lepton != leptons.end();
0108        ++lepton) {
0109     double dR2 = deltaR2(*lepton, p4_ref);
0110     if (leptonIdx_dR2min == -1 || dR2 < dR2min) {
0111       leptonIdx_dR2min = leptonIdx;
0112       dR2min = dR2;
0113     }
0114     ++leptonIdx;
0115   }
0116   assert(leptonIdx_dR2min >= 0 && leptonIdx_dR2min < (int)leptons.size());
0117   return leptonIdx_dR2min;
0118 }
0119 
0120 void scaleAndAddPFMEtSignObjects(std::vector<metsig::SigInputObj>& metSignObjects_scaled,
0121                                  const std::vector<metsig::SigInputObj>& metSignObjects,
0122                                  double sf,
0123                                  double sfMin,
0124                                  double sfMax) {
0125   double sf_value = sf;
0126   if (sf_value > sfMax)
0127     sf_value = sfMax;
0128   if (sf_value < sfMin)
0129     sf_value = sfMin;
0130   for (std::vector<metsig::SigInputObj>::const_iterator metSignObject = metSignObjects.begin();
0131        metSignObject != metSignObjects.end();
0132        ++metSignObject) {
0133     metsig::SigInputObj metSignObject_scaled;
0134     metSignObject_scaled.set(metSignObject->get_type(),
0135                              sf_value * metSignObject->get_energy(),
0136                              metSignObject->get_phi(),
0137                              sf_value * metSignObject->get_sigma_e(),
0138                              metSignObject->get_sigma_tan());
0139     metSignObjects_scaled.push_back(metSignObject_scaled);
0140   }
0141 }
0142 
0143 reco::METCovMatrix computePFMEtSignificance(const std::vector<metsig::SigInputObj>& metSignObjects) {
0144   reco::METCovMatrix pfMEtCov;
0145   if (metSignObjects.size() >= 2) {
0146     metsig::significanceAlgo pfMEtSignAlgorithm;
0147     pfMEtSignAlgorithm.addObjects(metSignObjects);
0148     pfMEtCov = pfMEtSignAlgorithm.getSignifMatrix();
0149   }
0150 
0151   double det = 0;
0152   pfMEtCov.Det(det);
0153   if (std::abs(det) < epsilon) {
0154     edm::LogWarning("computePFMEtSignificance") << "Inversion of PFMEt covariance matrix failed, det = " << det
0155                                                 << " --> replacing covariance matrix by resolution defaults !!";
0156     pfMEtCov(0, 0) = defaultPFMEtResolutionX * defaultPFMEtResolutionX;
0157     pfMEtCov(0, 1) = 0.;
0158     pfMEtCov(1, 0) = 0.;
0159     pfMEtCov(1, 1) = defaultPFMEtResolutionY * defaultPFMEtResolutionY;
0160   }
0161 
0162   return pfMEtCov;
0163 }
0164 
0165 void printP4(const std::string& label_part1, int idx, const std::string& label_part2, const reco::Candidate& candidate) {
0166   std::cout << label_part1 << " #" << idx << label_part2 << ": Pt = " << candidate.pt() << ", eta = " << candidate.eta()
0167             << ", phi = " << candidate.phi() << " (charge = " << candidate.charge() << ")" << std::endl;
0168 }
0169 
0170 void printCommonMETData(const std::string& label, const CommonMETData& metData) {
0171   std::cout << label << ": Px = " << metData.mex << ", Py = " << metData.mey << ", sumEt = " << metData.sumet
0172             << std::endl;
0173 }
0174 
0175 void printMVAMEtJetInfo(const std::string& label, int idx, const reco::PUSubMETCandInfo& jet) {
0176   std::cout << label << " #" << idx << " (";
0177   if (jet.type() == reco::PUSubMETCandInfo::kHS)
0178     std::cout << "no-PU";
0179   else if (jet.type() == reco::PUSubMETCandInfo::kPU)
0180     std::cout << "PU";
0181   std::cout << "): Pt = " << jet.p4().pt() << ", eta = " << jet.p4().eta() << ", phi = " << jet.p4().phi();
0182   std::cout << " id. flags: anti-noise = " << jet.passesLooseJetId() << std::endl;
0183   std::cout << std::endl;
0184 }
0185 
0186 void printMVAMEtPFCandInfo(const std::string& label, int idx, const reco::PUSubMETCandInfo& pfCand) {
0187   std::cout << label << " #" << idx << " (";
0188   if (pfCand.type() == reco::PUSubMETCandInfo::kChHS)
0189     std::cout << "no-PU charged";
0190   else if (pfCand.type() == reco::PUSubMETCandInfo::kChPU)
0191     std::cout << "PU charged";
0192   else if (pfCand.type() == reco::PUSubMETCandInfo::kNeutral)
0193     std::cout << "neutral";
0194   std::cout << "): Pt = " << pfCand.p4().pt() << ", eta = " << pfCand.p4().eta() << ", phi = " << pfCand.p4().phi();
0195   std::string isWithinJet_string;
0196   if (pfCand.isWithinJet())
0197     isWithinJet_string = "true";
0198   else
0199     isWithinJet_string = "false";
0200   std::cout << " (isWithinJet = " << isWithinJet_string << ")";
0201   if (pfCand.isWithinJet())
0202     std::cout << " Jet id. flags: anti-noise = " << pfCand.passesLooseJetId() << std::endl;
0203   std::cout << std::endl;
0204 }
0205 
0206 void NoPileUpPFMEtProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0207   LogDebug("produce") << " moduleLabel = " << moduleLabel_ << std::endl;
0208 
0209   // get original MET
0210   edm::Handle<reco::PFMETCollection> pfMETs;
0211   evt.getByToken(srcMEt_, pfMETs);
0212   if (!(pfMETs->size() == 1))
0213     throw cms::Exception("NoPileUpPFMEtProducer::produce") << "Failed to find unique MET object !!\n";
0214   const reco::PFMET& pfMEt_original = pfMETs->front();
0215 
0216   // get MET covariance matrix
0217   reco::METCovMatrix pfMEtCov;
0218   if (!srcMEtCov_.label().empty()) {
0219     //MM manual bypass to pfMET as this case has neer been presented
0220     // edm::Handle<PFMEtSignCovMatrix> pfMEtCovHandle;
0221     // evt.getByToken(srcMEtCov_, pfMEtCovHandle);
0222     // pfMEtCov = (*pfMEtCovHandle);
0223     pfMEtCov = pfMEt_original.getSignificanceMatrix();
0224   } else {
0225     pfMEtCov = pfMEt_original.getSignificanceMatrix();
0226   }
0227 
0228   // get lepton momenta
0229   std::vector<reco::Candidate::LorentzVector> leptons;
0230   std::vector<metsig::SigInputObj> metSignObjectsLeptons;
0231   reco::Candidate::LorentzVector sumLeptonP4s;
0232   for (std::vector<edm::EDGetTokenT<edm::View<reco::Candidate> > >::const_iterator srcLeptons_i = srcLeptons_.begin();
0233        srcLeptons_i != srcLeptons_.end();
0234        ++srcLeptons_i) {
0235     edm::Handle<reco::CandidateView> leptons_i;
0236     evt.getByToken(*srcLeptons_i, leptons_i);
0237     int leptonIdx = 0;
0238     for (reco::CandidateView::const_iterator lepton = leptons_i->begin(); lepton != leptons_i->end(); ++lepton) {
0239       leptons.push_back(lepton->p4());
0240       metSignObjectsLeptons.push_back(pfMEtSignInterface_->compResolution(&(*lepton)));
0241       sumLeptonP4s += lepton->p4();
0242       ++leptonIdx;
0243     }
0244   }
0245   LogDebug("produce") << " sum(leptons): Pt = " << sumLeptonP4s.pt() << ", eta = " << sumLeptonP4s.eta()
0246                       << ", phi = " << sumLeptonP4s.phi() << ","
0247                       << " mass = " << sumLeptonP4s.mass() << std::endl;
0248 
0249   // get jet and PFCandidate information
0250   edm::Handle<reco::PUSubMETCandInfoCollection> jets;
0251   evt.getByToken(srcJetInfo_, jets);
0252   edm::Handle<reco::PUSubMETCandInfoCollection> jetsLeptonMatch;
0253   evt.getByToken(srcJetInfoLeptonMatch_, jetsLeptonMatch);
0254   edm::Handle<reco::PUSubMETCandInfoCollection> pfCandidates;
0255   evt.getByToken(srcPFCandInfo_, pfCandidates);
0256   edm::Handle<reco::PUSubMETCandInfoCollection> pfCandidatesLeptonMatch;
0257   evt.getByToken(srcPFCandInfoLeptonMatch_, pfCandidatesLeptonMatch);
0258 
0259   reco::PUSubMETCandInfoCollection jets_leptons = utils_.cleanJets(*jetsLeptonMatch, leptons, 0.5, true);
0260   reco::PUSubMETCandInfoCollection pfCandidates_leptons =
0261       utils_.cleanPFCandidates(*pfCandidatesLeptonMatch, leptons, 0.3, true);
0262   std::vector<CommonMETData> sumJetsPlusPFCandidates_leptons(leptons.size());
0263   for (std::vector<CommonMETData>::iterator sumJetsPlusPFCandidates = sumJetsPlusPFCandidates_leptons.begin();
0264        sumJetsPlusPFCandidates != sumJetsPlusPFCandidates_leptons.end();
0265        ++sumJetsPlusPFCandidates) {
0266     initializeCommonMETData(*sumJetsPlusPFCandidates);
0267   }
0268   for (reco::PUSubMETCandInfoCollection::const_iterator jet = jets_leptons.begin(); jet != jets_leptons.end(); ++jet) {
0269     int leptonIdx_dRmin = findBestMatchingLepton(leptons, jet->p4());
0270     assert(leptonIdx_dRmin >= 0 && leptonIdx_dRmin < (int)sumJetsPlusPFCandidates_leptons.size());
0271 
0272     LogDebug("produce") << "jet-to-lepton match:"
0273                         << " jetPt = " << jet->p4().pt() << ", jetEta = " << jet->p4().eta()
0274                         << ", jetPhi = " << jet->p4().phi() << " leptonPt = " << leptons[leptonIdx_dRmin].pt()
0275                         << ", leptonEta = " << leptons[leptonIdx_dRmin].eta()
0276                         << ", leptonPhi = " << leptons[leptonIdx_dRmin].phi() << std::endl;
0277 
0278     sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mex += jet->p4().px();
0279     sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mey += jet->p4().py();
0280     sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].sumet += jet->p4().pt();
0281   }
0282   for (reco::PUSubMETCandInfoCollection::const_iterator pfCandidate = pfCandidates_leptons.begin();
0283        pfCandidate != pfCandidates_leptons.end();
0284        ++pfCandidate) {
0285     bool isWithinJet_lepton = false;
0286     if (pfCandidate->isWithinJet()) {
0287       for (reco::PUSubMETCandInfoCollection::const_iterator jet = jets_leptons.begin(); jet != jets_leptons.end();
0288            ++jet) {
0289         double dR2 = deltaR2(pfCandidate->p4(), jet->p4());
0290         if (dR2 < 0.5 * 0.5)
0291           isWithinJet_lepton = true;
0292       }
0293     }
0294     if (!isWithinJet_lepton) {
0295       int leptonIdx_dRmin = findBestMatchingLepton(leptons, pfCandidate->p4());
0296       assert(leptonIdx_dRmin >= 0 && leptonIdx_dRmin < (int)sumJetsPlusPFCandidates_leptons.size());
0297       LogDebug("produce") << "pfCandidate-to-lepton match:"
0298                           << " pfCandidatePt = " << pfCandidate->p4().pt()
0299                           << ", pfCandidateEta = " << pfCandidate->p4().eta()
0300                           << ", pfCandidatePhi = " << pfCandidate->p4().phi()
0301                           << " leptonPt = " << leptons[leptonIdx_dRmin].pt()
0302                           << ", leptonEta = " << leptons[leptonIdx_dRmin].eta()
0303                           << ", leptonPhi = " << leptons[leptonIdx_dRmin].phi() << std::endl;
0304 
0305       sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mex += pfCandidate->p4().px();
0306       sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mey += pfCandidate->p4().py();
0307       sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].sumet += pfCandidate->p4().pt();
0308     } else {
0309       LogDebug("produce") << " pfCandidate is within jet --> skipping." << std::endl;
0310     }
0311   }
0312   auto sumLeptons = std::make_unique<CommonMETData>();
0313   initializeCommonMETData(*sumLeptons);
0314   auto sumLeptonIsoCones = std::make_unique<CommonMETData>();
0315   initializeCommonMETData(*sumLeptonIsoCones);
0316   int leptonIdx = 0;
0317   for (std::vector<CommonMETData>::iterator sumJetsPlusPFCandidates = sumJetsPlusPFCandidates_leptons.begin();
0318        sumJetsPlusPFCandidates != sumJetsPlusPFCandidates_leptons.end();
0319        ++sumJetsPlusPFCandidates) {
0320     if (sumJetsPlusPFCandidates->sumet > leptons[leptonIdx].pt()) {
0321       double leptonEnFrac = leptons[leptonIdx].pt() / sumJetsPlusPFCandidates->sumet;
0322       assert(leptonEnFrac >= 0.0 && leptonEnFrac <= 1.0);
0323       sumLeptons->mex += (leptonEnFrac * sumJetsPlusPFCandidates->mex);
0324       sumLeptons->mey += (leptonEnFrac * sumJetsPlusPFCandidates->mey);
0325       sumLeptons->sumet += (leptonEnFrac * sumJetsPlusPFCandidates->sumet);
0326       double leptonIsoConeEnFrac = 1.0 - leptonEnFrac;
0327       assert(leptonIsoConeEnFrac >= 0.0 && leptonIsoConeEnFrac <= 1.0);
0328       sumLeptonIsoCones->mex += (leptonIsoConeEnFrac * sumJetsPlusPFCandidates->mex);
0329       sumLeptonIsoCones->mey += (leptonIsoConeEnFrac * sumJetsPlusPFCandidates->mey);
0330       sumLeptonIsoCones->sumet += (leptonIsoConeEnFrac * sumJetsPlusPFCandidates->sumet);
0331     } else {
0332       sumLeptons->mex += sumJetsPlusPFCandidates->mex;
0333       sumLeptons->mey += sumJetsPlusPFCandidates->mey;
0334       sumLeptons->sumet += sumJetsPlusPFCandidates->sumet;
0335     }
0336     ++leptonIdx;
0337   }
0338 
0339   reco::PUSubMETCandInfoCollection jets_cleaned = utils_.cleanJets(*jets, leptons, 0.5, false);
0340   reco::PUSubMETCandInfoCollection pfCandidates_cleaned = utils_.cleanPFCandidates(*pfCandidates, leptons, 0.3, false);
0341 
0342   auto sumNoPUjets = std::make_unique<CommonMETData>();
0343   initializeCommonMETData(*sumNoPUjets);
0344   std::vector<metsig::SigInputObj> metSignObjectsNoPUjets;
0345   auto sumNoPUjetOffsetEnCorr = std::make_unique<CommonMETData>();
0346   initializeCommonMETData(*sumNoPUjetOffsetEnCorr);
0347   std::vector<metsig::SigInputObj> metSignObjectsNoPUjetOffsetEnCorr;
0348   auto sumPUjets = std::make_unique<CommonMETData>();
0349   initializeCommonMETData(*sumPUjets);
0350   std::vector<metsig::SigInputObj> metSignObjectsPUjets;
0351   int jetIdx = 0;
0352   for (reco::PUSubMETCandInfoCollection::const_iterator jet = jets_cleaned.begin(); jet != jets_cleaned.end(); ++jet) {
0353     if (jet->passesLooseJetId()) {
0354       if (jet->type() == reco::PUSubMETCandInfo::kHS) {
0355         addToCommonMETData(*sumNoPUjets, jet->p4());
0356         metSignObjectsNoPUjets.push_back(jet->metSignObj());
0357         float jetp = jet->p4().P();
0358         float jetcorr = jet->offsetEnCorr();
0359         sumNoPUjetOffsetEnCorr->mex += jetcorr * jet->p4().px() / jetp;
0360         sumNoPUjetOffsetEnCorr->mey += jetcorr * jet->p4().py() / jetp;
0361         sumNoPUjetOffsetEnCorr->mez += jetcorr * jet->p4().pz() / jetp;
0362         sumNoPUjetOffsetEnCorr->sumet += jetcorr * jet->p4().pt() / jetp;
0363         metsig::SigInputObj pfMEtSignObjectOffsetEnCorr(
0364             jet->metSignObj().get_type(),
0365             jet->offsetEnCorr(),
0366             jet->metSignObj().get_phi(),
0367             (jet->offsetEnCorr() / jet->p4().E()) * jet->metSignObj().get_sigma_e(),
0368             jet->metSignObj().get_sigma_tan());
0369         metSignObjectsNoPUjetOffsetEnCorr.push_back(pfMEtSignObjectOffsetEnCorr);
0370       } else {
0371         addToCommonMETData(*sumPUjets, jet->p4());
0372         metSignObjectsPUjets.push_back(jet->metSignObj());
0373       }
0374     }
0375     ++jetIdx;
0376   }
0377 
0378   auto sumNoPUunclChargedCands = std::make_unique<CommonMETData>();
0379   initializeCommonMETData(*sumNoPUunclChargedCands);
0380   std::vector<metsig::SigInputObj> metSignObjectsNoPUunclChargedCands;
0381   auto sumPUunclChargedCands = std::make_unique<CommonMETData>();
0382   initializeCommonMETData(*sumPUunclChargedCands);
0383   std::vector<metsig::SigInputObj> metSignObjectsPUunclChargedCands;
0384   auto sumUnclNeutralCands = std::make_unique<CommonMETData>();
0385   initializeCommonMETData(*sumUnclNeutralCands);
0386   std::vector<metsig::SigInputObj> metSignObjectsUnclNeutralCands;
0387   int pfCandIdx = 0;
0388   for (reco::PUSubMETCandInfoCollection::const_iterator pfCandidate = pfCandidates_cleaned.begin();
0389        pfCandidate != pfCandidates_cleaned.end();
0390        ++pfCandidate) {
0391     if (pfCandidate->passesLooseJetId()) {
0392       if (!pfCandidate->isWithinJet()) {
0393         if (pfCandidate->type() == reco::PUSubMETCandInfo::kChHS) {
0394           addToCommonMETData(*sumNoPUunclChargedCands, pfCandidate->p4());
0395           metSignObjectsNoPUunclChargedCands.push_back(pfCandidate->metSignObj());
0396         } else if (pfCandidate->type() == reco::PUSubMETCandInfo::kChPU) {
0397           addToCommonMETData(*sumPUunclChargedCands, pfCandidate->p4());
0398           metSignObjectsPUunclChargedCands.push_back(pfCandidate->metSignObj());
0399         } else if (pfCandidate->type() == reco::PUSubMETCandInfo::kNeutral) {
0400           addToCommonMETData(*sumUnclNeutralCands, pfCandidate->p4());
0401           metSignObjectsUnclNeutralCands.push_back(pfCandidate->metSignObj());
0402         }
0403       }
0404     }
0405     ++pfCandIdx;
0406   }
0407 
0408   edm::Handle<CorrMETData> type0Correction_input;
0409   evt.getByToken(srcType0Correction_, type0Correction_input);
0410   auto type0Correction_output = std::make_unique<CommonMETData>();
0411   initializeCommonMETData(*type0Correction_output);
0412   type0Correction_output->mex = type0Correction_input->mex;
0413   type0Correction_output->mey = type0Correction_input->mey;
0414 
0415   finalizeCommonMETData(*sumLeptons);
0416   finalizeCommonMETData(*sumNoPUjetOffsetEnCorr);
0417   finalizeCommonMETData(*sumNoPUjets);
0418   finalizeCommonMETData(*sumPUjets);
0419   finalizeCommonMETData(*sumNoPUunclChargedCands);
0420   finalizeCommonMETData(*sumPUunclChargedCands);
0421   finalizeCommonMETData(*sumUnclNeutralCands);
0422   finalizeCommonMETData(*type0Correction_output);
0423   finalizeCommonMETData(*sumLeptonIsoCones);
0424 
0425   double noPileUpScaleFactor =
0426       (sumPUunclChargedCands->sumet > 0.)
0427           ? (sumPUunclChargedCands->sumet / (sumNoPUunclChargedCands->sumet + sumPUunclChargedCands->sumet))
0428           : 1.;
0429   LogDebug("produce") << "noPileUpScaleFactor = " << noPileUpScaleFactor << std::endl;
0430 
0431   double noPileUpMEtPx =
0432       -(sumLeptons->mex + sumNoPUjets->mex + sumNoPUunclChargedCands->mex +
0433         noPileUpScaleFactor *
0434             (sfNoPUjetOffsetEnCorr_ * sumNoPUjetOffsetEnCorr->mex + sfUnclNeutralCands_ * sumUnclNeutralCands->mex +
0435              sfPUunclChargedCands_ * sumPUunclChargedCands->mex + sfPUjets_ * sumPUjets->mex)) +
0436       noPileUpScaleFactor * sfType0Correction_ * type0Correction_output->mex;
0437   if (sfLeptonIsoCones_ >= 0.)
0438     noPileUpMEtPx -= (noPileUpScaleFactor * sfLeptonIsoCones_ * sumLeptonIsoCones->mex);
0439   else
0440     noPileUpMEtPx -= (std::abs(sfLeptonIsoCones_) * sumLeptonIsoCones->mex);
0441   double noPileUpMEtPy =
0442       -(sumLeptons->mey + sumNoPUjets->mey + sumNoPUunclChargedCands->mey +
0443         noPileUpScaleFactor *
0444             (sfNoPUjetOffsetEnCorr_ * sumNoPUjetOffsetEnCorr->mey + sfUnclNeutralCands_ * sumUnclNeutralCands->mey +
0445              sfPUunclChargedCands_ * sumPUunclChargedCands->mey + sfPUjets_ * sumPUjets->mey)) +
0446       noPileUpScaleFactor * sfType0Correction_ * type0Correction_output->mey;
0447   if (sfLeptonIsoCones_ >= 0.)
0448     noPileUpMEtPy -= (noPileUpScaleFactor * sfLeptonIsoCones_ * sumLeptonIsoCones->mey);
0449   else
0450     noPileUpMEtPy -= (std::abs(sfLeptonIsoCones_) * sumLeptonIsoCones->mey);
0451   double noPileUpMEtPt = sqrt(noPileUpMEtPx * noPileUpMEtPx + noPileUpMEtPy * noPileUpMEtPy);
0452   reco::Candidate::LorentzVector noPileUpMEtP4(noPileUpMEtPx, noPileUpMEtPy, 0., noPileUpMEtPt);
0453 
0454   reco::PFMET noPileUpMEt(pfMEt_original);
0455   noPileUpMEt.setP4(noPileUpMEtP4);
0456   //noPileUpMEt.setSignificanceMatrix(pfMEtCov);
0457 
0458   std::vector<metsig::SigInputObj> metSignObjects_scaled;
0459   scaleAndAddPFMEtSignObjects(metSignObjects_scaled, metSignObjectsLeptons, 1.0, sfMEtCovMin_, sfMEtCovMax_);
0460   scaleAndAddPFMEtSignObjects(
0461       metSignObjects_scaled, metSignObjectsNoPUjetOffsetEnCorr, sfNoPUjetOffsetEnCorr_, sfMEtCovMin_, sfMEtCovMax_);
0462   scaleAndAddPFMEtSignObjects(metSignObjects_scaled, metSignObjectsNoPUjets, sfNoPUjets_, sfMEtCovMin_, sfMEtCovMax_);
0463   scaleAndAddPFMEtSignObjects(
0464       metSignObjects_scaled, metSignObjectsPUjets, noPileUpScaleFactor * sfPUjets_, sfMEtCovMin_, sfMEtCovMax_);
0465   scaleAndAddPFMEtSignObjects(
0466       metSignObjects_scaled, metSignObjectsNoPUunclChargedCands, sfNoPUunclChargedCands_, sfMEtCovMin_, sfMEtCovMax_);
0467   scaleAndAddPFMEtSignObjects(metSignObjects_scaled,
0468                               metSignObjectsPUunclChargedCands,
0469                               noPileUpScaleFactor * sfPUunclChargedCands_,
0470                               sfMEtCovMin_,
0471                               sfMEtCovMax_);
0472   scaleAndAddPFMEtSignObjects(metSignObjects_scaled,
0473                               metSignObjectsUnclNeutralCands,
0474                               noPileUpScaleFactor * sfUnclNeutralCands_,
0475                               sfMEtCovMin_,
0476                               sfMEtCovMax_);
0477   reco::METCovMatrix pfMEtCov_recomputed = computePFMEtSignificance(metSignObjects_scaled);
0478   noPileUpMEt.setSignificanceMatrix(pfMEtCov_recomputed);
0479 
0480   LogDebug("produce") << "<NoPileUpPFMEtProducer::produce>:" << std::endl
0481                       << " moduleLabel = " << moduleLabel_ << std::endl
0482                       << " PFMET: Pt = " << pfMEt_original.pt() << ", phi = " << pfMEt_original.phi() << " "
0483                       << "(Px = " << pfMEt_original.px() << ", Py = " << pfMEt_original.py() << ")" << std::endl
0484                       << " Cov:" << std::endl
0485                       << " " << pfMEtCov(0, 0) << "  " << pfMEtCov(0, 1) << "\n " << pfMEtCov(1, 0) << "  "
0486                       << pfMEtCov(1, 1) << std::endl
0487                       << " no-PU MET: Pt = " << noPileUpMEt.pt() << ", phi = " << noPileUpMEt.phi() << " "
0488                       << "(Px = " << noPileUpMEt.px() << ", Py = " << noPileUpMEt.py() << ")" << std::endl
0489                       << " Cov:" << std::endl
0490                       << " " << (noPileUpMEt.getSignificanceMatrix())(0, 0) << "  "
0491                       << (noPileUpMEt.getSignificanceMatrix())(0, 1) << std::endl
0492                       << (noPileUpMEt.getSignificanceMatrix())(1, 0) << "  "
0493                       << (noPileUpMEt.getSignificanceMatrix())(1, 1) << std::endl;
0494 
0495   // add no-PU MET object to the event
0496   auto noPileUpMEtCollection = std::make_unique<reco::PFMETCollection>();
0497   noPileUpMEtCollection->push_back(noPileUpMEt);
0498 
0499   evt.put(std::move(noPileUpMEtCollection));
0500   if (saveInputs_) {
0501     evt.put(std::move(sumLeptons), sfLeptonsName_);
0502     evt.put(std::move(sumNoPUjetOffsetEnCorr), sfNoPUjetOffsetEnCorrName_);
0503     evt.put(std::move(sumNoPUjets), sfNoPUjetsName_);
0504     evt.put(std::move(sumPUjets), sfPUjetsName_);
0505     evt.put(std::move(sumNoPUunclChargedCands), sfNoPUunclChargedCandsName_);
0506     evt.put(std::move(sumPUunclChargedCands), sfPUunclChargedCandsName_);
0507     evt.put(std::move(sumUnclNeutralCands), sfUnclNeutralCandsName_);
0508     evt.put(std::move(type0Correction_output), sfType0CorrectionName_);
0509     evt.put(std::move(sumLeptonIsoCones), sfLeptonIsoConesName_);
0510   }
0511 
0512   evt.put(std::make_unique<double>(noPileUpScaleFactor), "sfNoPU");
0513 }
0514 
0515 #include "FWCore/Framework/interface/MakerMacros.h"
0516 
0517 DEFINE_FWK_MODULE(NoPileUpPFMEtProducer);