Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:26:48 

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     for (reco::CandidateView::const_iterator lepton = leptons_i->begin(); lepton != leptons_i->end(); ++lepton) {
0238       leptons.push_back(lepton->p4());
0239       metSignObjectsLeptons.push_back(pfMEtSignInterface_->compResolution(&(*lepton)));
0240       sumLeptonP4s += lepton->p4();
0241     }
0242   }
0243   LogDebug("produce") << " sum(leptons): Pt = " << sumLeptonP4s.pt() << ", eta = " << sumLeptonP4s.eta()
0244                       << ", phi = " << sumLeptonP4s.phi() << ","
0245                       << " mass = " << sumLeptonP4s.mass() << std::endl;
0246 
0247   // get jet and PFCandidate information
0248   edm::Handle<reco::PUSubMETCandInfoCollection> jets;
0249   evt.getByToken(srcJetInfo_, jets);
0250   edm::Handle<reco::PUSubMETCandInfoCollection> jetsLeptonMatch;
0251   evt.getByToken(srcJetInfoLeptonMatch_, jetsLeptonMatch);
0252   edm::Handle<reco::PUSubMETCandInfoCollection> pfCandidates;
0253   evt.getByToken(srcPFCandInfo_, pfCandidates);
0254   edm::Handle<reco::PUSubMETCandInfoCollection> pfCandidatesLeptonMatch;
0255   evt.getByToken(srcPFCandInfoLeptonMatch_, pfCandidatesLeptonMatch);
0256 
0257   reco::PUSubMETCandInfoCollection jets_leptons = utils_.cleanJets(*jetsLeptonMatch, leptons, 0.5, true);
0258   reco::PUSubMETCandInfoCollection pfCandidates_leptons =
0259       utils_.cleanPFCandidates(*pfCandidatesLeptonMatch, leptons, 0.3, true);
0260   std::vector<CommonMETData> sumJetsPlusPFCandidates_leptons(leptons.size());
0261   for (std::vector<CommonMETData>::iterator sumJetsPlusPFCandidates = sumJetsPlusPFCandidates_leptons.begin();
0262        sumJetsPlusPFCandidates != sumJetsPlusPFCandidates_leptons.end();
0263        ++sumJetsPlusPFCandidates) {
0264     initializeCommonMETData(*sumJetsPlusPFCandidates);
0265   }
0266   for (reco::PUSubMETCandInfoCollection::const_iterator jet = jets_leptons.begin(); jet != jets_leptons.end(); ++jet) {
0267     int leptonIdx_dRmin = findBestMatchingLepton(leptons, jet->p4());
0268     assert(leptonIdx_dRmin >= 0 && leptonIdx_dRmin < (int)sumJetsPlusPFCandidates_leptons.size());
0269 
0270     LogDebug("produce") << "jet-to-lepton match:"
0271                         << " jetPt = " << jet->p4().pt() << ", jetEta = " << jet->p4().eta()
0272                         << ", jetPhi = " << jet->p4().phi() << " leptonPt = " << leptons[leptonIdx_dRmin].pt()
0273                         << ", leptonEta = " << leptons[leptonIdx_dRmin].eta()
0274                         << ", leptonPhi = " << leptons[leptonIdx_dRmin].phi() << std::endl;
0275 
0276     sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mex += jet->p4().px();
0277     sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mey += jet->p4().py();
0278     sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].sumet += jet->p4().pt();
0279   }
0280   for (reco::PUSubMETCandInfoCollection::const_iterator pfCandidate = pfCandidates_leptons.begin();
0281        pfCandidate != pfCandidates_leptons.end();
0282        ++pfCandidate) {
0283     bool isWithinJet_lepton = false;
0284     if (pfCandidate->isWithinJet()) {
0285       for (reco::PUSubMETCandInfoCollection::const_iterator jet = jets_leptons.begin(); jet != jets_leptons.end();
0286            ++jet) {
0287         double dR2 = deltaR2(pfCandidate->p4(), jet->p4());
0288         if (dR2 < 0.5 * 0.5)
0289           isWithinJet_lepton = true;
0290       }
0291     }
0292     if (!isWithinJet_lepton) {
0293       int leptonIdx_dRmin = findBestMatchingLepton(leptons, pfCandidate->p4());
0294       assert(leptonIdx_dRmin >= 0 && leptonIdx_dRmin < (int)sumJetsPlusPFCandidates_leptons.size());
0295       LogDebug("produce") << "pfCandidate-to-lepton match:"
0296                           << " pfCandidatePt = " << pfCandidate->p4().pt()
0297                           << ", pfCandidateEta = " << pfCandidate->p4().eta()
0298                           << ", pfCandidatePhi = " << pfCandidate->p4().phi()
0299                           << " leptonPt = " << leptons[leptonIdx_dRmin].pt()
0300                           << ", leptonEta = " << leptons[leptonIdx_dRmin].eta()
0301                           << ", leptonPhi = " << leptons[leptonIdx_dRmin].phi() << std::endl;
0302 
0303       sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mex += pfCandidate->p4().px();
0304       sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].mey += pfCandidate->p4().py();
0305       sumJetsPlusPFCandidates_leptons[leptonIdx_dRmin].sumet += pfCandidate->p4().pt();
0306     } else {
0307       LogDebug("produce") << " pfCandidate is within jet --> skipping." << std::endl;
0308     }
0309   }
0310   auto sumLeptons = std::make_unique<CommonMETData>();
0311   initializeCommonMETData(*sumLeptons);
0312   auto sumLeptonIsoCones = std::make_unique<CommonMETData>();
0313   initializeCommonMETData(*sumLeptonIsoCones);
0314   int leptonIdx = 0;
0315   for (std::vector<CommonMETData>::iterator sumJetsPlusPFCandidates = sumJetsPlusPFCandidates_leptons.begin();
0316        sumJetsPlusPFCandidates != sumJetsPlusPFCandidates_leptons.end();
0317        ++sumJetsPlusPFCandidates) {
0318     if (sumJetsPlusPFCandidates->sumet > leptons[leptonIdx].pt()) {
0319       double leptonEnFrac = leptons[leptonIdx].pt() / sumJetsPlusPFCandidates->sumet;
0320       assert(leptonEnFrac >= 0.0 && leptonEnFrac <= 1.0);
0321       sumLeptons->mex += (leptonEnFrac * sumJetsPlusPFCandidates->mex);
0322       sumLeptons->mey += (leptonEnFrac * sumJetsPlusPFCandidates->mey);
0323       sumLeptons->sumet += (leptonEnFrac * sumJetsPlusPFCandidates->sumet);
0324       double leptonIsoConeEnFrac = 1.0 - leptonEnFrac;
0325       assert(leptonIsoConeEnFrac >= 0.0 && leptonIsoConeEnFrac <= 1.0);
0326       sumLeptonIsoCones->mex += (leptonIsoConeEnFrac * sumJetsPlusPFCandidates->mex);
0327       sumLeptonIsoCones->mey += (leptonIsoConeEnFrac * sumJetsPlusPFCandidates->mey);
0328       sumLeptonIsoCones->sumet += (leptonIsoConeEnFrac * sumJetsPlusPFCandidates->sumet);
0329     } else {
0330       sumLeptons->mex += sumJetsPlusPFCandidates->mex;
0331       sumLeptons->mey += sumJetsPlusPFCandidates->mey;
0332       sumLeptons->sumet += sumJetsPlusPFCandidates->sumet;
0333     }
0334     ++leptonIdx;
0335   }
0336 
0337   reco::PUSubMETCandInfoCollection jets_cleaned = utils_.cleanJets(*jets, leptons, 0.5, false);
0338   reco::PUSubMETCandInfoCollection pfCandidates_cleaned = utils_.cleanPFCandidates(*pfCandidates, leptons, 0.3, false);
0339 
0340   auto sumNoPUjets = std::make_unique<CommonMETData>();
0341   initializeCommonMETData(*sumNoPUjets);
0342   std::vector<metsig::SigInputObj> metSignObjectsNoPUjets;
0343   auto sumNoPUjetOffsetEnCorr = std::make_unique<CommonMETData>();
0344   initializeCommonMETData(*sumNoPUjetOffsetEnCorr);
0345   std::vector<metsig::SigInputObj> metSignObjectsNoPUjetOffsetEnCorr;
0346   auto sumPUjets = std::make_unique<CommonMETData>();
0347   initializeCommonMETData(*sumPUjets);
0348   std::vector<metsig::SigInputObj> metSignObjectsPUjets;
0349   for (reco::PUSubMETCandInfoCollection::const_iterator jet = jets_cleaned.begin(); jet != jets_cleaned.end(); ++jet) {
0350     if (jet->passesLooseJetId()) {
0351       if (jet->type() == reco::PUSubMETCandInfo::kHS) {
0352         addToCommonMETData(*sumNoPUjets, jet->p4());
0353         metSignObjectsNoPUjets.push_back(jet->metSignObj());
0354         float jetp = jet->p4().P();
0355         float jetcorr = jet->offsetEnCorr();
0356         sumNoPUjetOffsetEnCorr->mex += jetcorr * jet->p4().px() / jetp;
0357         sumNoPUjetOffsetEnCorr->mey += jetcorr * jet->p4().py() / jetp;
0358         sumNoPUjetOffsetEnCorr->mez += jetcorr * jet->p4().pz() / jetp;
0359         sumNoPUjetOffsetEnCorr->sumet += jetcorr * jet->p4().pt() / jetp;
0360         metsig::SigInputObj pfMEtSignObjectOffsetEnCorr(
0361             jet->metSignObj().get_type(),
0362             jet->offsetEnCorr(),
0363             jet->metSignObj().get_phi(),
0364             (jet->offsetEnCorr() / jet->p4().E()) * jet->metSignObj().get_sigma_e(),
0365             jet->metSignObj().get_sigma_tan());
0366         metSignObjectsNoPUjetOffsetEnCorr.push_back(pfMEtSignObjectOffsetEnCorr);
0367       } else {
0368         addToCommonMETData(*sumPUjets, jet->p4());
0369         metSignObjectsPUjets.push_back(jet->metSignObj());
0370       }
0371     }
0372   }
0373 
0374   auto sumNoPUunclChargedCands = std::make_unique<CommonMETData>();
0375   initializeCommonMETData(*sumNoPUunclChargedCands);
0376   std::vector<metsig::SigInputObj> metSignObjectsNoPUunclChargedCands;
0377   auto sumPUunclChargedCands = std::make_unique<CommonMETData>();
0378   initializeCommonMETData(*sumPUunclChargedCands);
0379   std::vector<metsig::SigInputObj> metSignObjectsPUunclChargedCands;
0380   auto sumUnclNeutralCands = std::make_unique<CommonMETData>();
0381   initializeCommonMETData(*sumUnclNeutralCands);
0382   std::vector<metsig::SigInputObj> metSignObjectsUnclNeutralCands;
0383   for (reco::PUSubMETCandInfoCollection::const_iterator pfCandidate = pfCandidates_cleaned.begin();
0384        pfCandidate != pfCandidates_cleaned.end();
0385        ++pfCandidate) {
0386     if (pfCandidate->passesLooseJetId()) {
0387       if (!pfCandidate->isWithinJet()) {
0388         if (pfCandidate->type() == reco::PUSubMETCandInfo::kChHS) {
0389           addToCommonMETData(*sumNoPUunclChargedCands, pfCandidate->p4());
0390           metSignObjectsNoPUunclChargedCands.push_back(pfCandidate->metSignObj());
0391         } else if (pfCandidate->type() == reco::PUSubMETCandInfo::kChPU) {
0392           addToCommonMETData(*sumPUunclChargedCands, pfCandidate->p4());
0393           metSignObjectsPUunclChargedCands.push_back(pfCandidate->metSignObj());
0394         } else if (pfCandidate->type() == reco::PUSubMETCandInfo::kNeutral) {
0395           addToCommonMETData(*sumUnclNeutralCands, pfCandidate->p4());
0396           metSignObjectsUnclNeutralCands.push_back(pfCandidate->metSignObj());
0397         }
0398       }
0399     }
0400   }
0401 
0402   edm::Handle<CorrMETData> type0Correction_input;
0403   evt.getByToken(srcType0Correction_, type0Correction_input);
0404   auto type0Correction_output = std::make_unique<CommonMETData>();
0405   initializeCommonMETData(*type0Correction_output);
0406   type0Correction_output->mex = type0Correction_input->mex;
0407   type0Correction_output->mey = type0Correction_input->mey;
0408 
0409   finalizeCommonMETData(*sumLeptons);
0410   finalizeCommonMETData(*sumNoPUjetOffsetEnCorr);
0411   finalizeCommonMETData(*sumNoPUjets);
0412   finalizeCommonMETData(*sumPUjets);
0413   finalizeCommonMETData(*sumNoPUunclChargedCands);
0414   finalizeCommonMETData(*sumPUunclChargedCands);
0415   finalizeCommonMETData(*sumUnclNeutralCands);
0416   finalizeCommonMETData(*type0Correction_output);
0417   finalizeCommonMETData(*sumLeptonIsoCones);
0418 
0419   double noPileUpScaleFactor =
0420       (sumPUunclChargedCands->sumet > 0.)
0421           ? (sumPUunclChargedCands->sumet / (sumNoPUunclChargedCands->sumet + sumPUunclChargedCands->sumet))
0422           : 1.;
0423   LogDebug("produce") << "noPileUpScaleFactor = " << noPileUpScaleFactor << std::endl;
0424 
0425   double noPileUpMEtPx =
0426       -(sumLeptons->mex + sumNoPUjets->mex + sumNoPUunclChargedCands->mex +
0427         noPileUpScaleFactor *
0428             (sfNoPUjetOffsetEnCorr_ * sumNoPUjetOffsetEnCorr->mex + sfUnclNeutralCands_ * sumUnclNeutralCands->mex +
0429              sfPUunclChargedCands_ * sumPUunclChargedCands->mex + sfPUjets_ * sumPUjets->mex)) +
0430       noPileUpScaleFactor * sfType0Correction_ * type0Correction_output->mex;
0431   if (sfLeptonIsoCones_ >= 0.)
0432     noPileUpMEtPx -= (noPileUpScaleFactor * sfLeptonIsoCones_ * sumLeptonIsoCones->mex);
0433   else
0434     noPileUpMEtPx -= (std::abs(sfLeptonIsoCones_) * sumLeptonIsoCones->mex);
0435   double noPileUpMEtPy =
0436       -(sumLeptons->mey + sumNoPUjets->mey + sumNoPUunclChargedCands->mey +
0437         noPileUpScaleFactor *
0438             (sfNoPUjetOffsetEnCorr_ * sumNoPUjetOffsetEnCorr->mey + sfUnclNeutralCands_ * sumUnclNeutralCands->mey +
0439              sfPUunclChargedCands_ * sumPUunclChargedCands->mey + sfPUjets_ * sumPUjets->mey)) +
0440       noPileUpScaleFactor * sfType0Correction_ * type0Correction_output->mey;
0441   if (sfLeptonIsoCones_ >= 0.)
0442     noPileUpMEtPy -= (noPileUpScaleFactor * sfLeptonIsoCones_ * sumLeptonIsoCones->mey);
0443   else
0444     noPileUpMEtPy -= (std::abs(sfLeptonIsoCones_) * sumLeptonIsoCones->mey);
0445   double noPileUpMEtPt = sqrt(noPileUpMEtPx * noPileUpMEtPx + noPileUpMEtPy * noPileUpMEtPy);
0446   reco::Candidate::LorentzVector noPileUpMEtP4(noPileUpMEtPx, noPileUpMEtPy, 0., noPileUpMEtPt);
0447 
0448   reco::PFMET noPileUpMEt(pfMEt_original);
0449   noPileUpMEt.setP4(noPileUpMEtP4);
0450   //noPileUpMEt.setSignificanceMatrix(pfMEtCov);
0451 
0452   std::vector<metsig::SigInputObj> metSignObjects_scaled;
0453   scaleAndAddPFMEtSignObjects(metSignObjects_scaled, metSignObjectsLeptons, 1.0, sfMEtCovMin_, sfMEtCovMax_);
0454   scaleAndAddPFMEtSignObjects(
0455       metSignObjects_scaled, metSignObjectsNoPUjetOffsetEnCorr, sfNoPUjetOffsetEnCorr_, sfMEtCovMin_, sfMEtCovMax_);
0456   scaleAndAddPFMEtSignObjects(metSignObjects_scaled, metSignObjectsNoPUjets, sfNoPUjets_, sfMEtCovMin_, sfMEtCovMax_);
0457   scaleAndAddPFMEtSignObjects(
0458       metSignObjects_scaled, metSignObjectsPUjets, noPileUpScaleFactor * sfPUjets_, sfMEtCovMin_, sfMEtCovMax_);
0459   scaleAndAddPFMEtSignObjects(
0460       metSignObjects_scaled, metSignObjectsNoPUunclChargedCands, sfNoPUunclChargedCands_, sfMEtCovMin_, sfMEtCovMax_);
0461   scaleAndAddPFMEtSignObjects(metSignObjects_scaled,
0462                               metSignObjectsPUunclChargedCands,
0463                               noPileUpScaleFactor * sfPUunclChargedCands_,
0464                               sfMEtCovMin_,
0465                               sfMEtCovMax_);
0466   scaleAndAddPFMEtSignObjects(metSignObjects_scaled,
0467                               metSignObjectsUnclNeutralCands,
0468                               noPileUpScaleFactor * sfUnclNeutralCands_,
0469                               sfMEtCovMin_,
0470                               sfMEtCovMax_);
0471   reco::METCovMatrix pfMEtCov_recomputed = computePFMEtSignificance(metSignObjects_scaled);
0472   noPileUpMEt.setSignificanceMatrix(pfMEtCov_recomputed);
0473 
0474   LogDebug("produce") << "<NoPileUpPFMEtProducer::produce>:" << std::endl
0475                       << " moduleLabel = " << moduleLabel_ << std::endl
0476                       << " PFMET: Pt = " << pfMEt_original.pt() << ", phi = " << pfMEt_original.phi() << " "
0477                       << "(Px = " << pfMEt_original.px() << ", Py = " << pfMEt_original.py() << ")" << std::endl
0478                       << " Cov:" << std::endl
0479                       << " " << pfMEtCov(0, 0) << "  " << pfMEtCov(0, 1) << "\n " << pfMEtCov(1, 0) << "  "
0480                       << pfMEtCov(1, 1) << std::endl
0481                       << " no-PU MET: Pt = " << noPileUpMEt.pt() << ", phi = " << noPileUpMEt.phi() << " "
0482                       << "(Px = " << noPileUpMEt.px() << ", Py = " << noPileUpMEt.py() << ")" << std::endl
0483                       << " Cov:" << std::endl
0484                       << " " << (noPileUpMEt.getSignificanceMatrix())(0, 0) << "  "
0485                       << (noPileUpMEt.getSignificanceMatrix())(0, 1) << std::endl
0486                       << (noPileUpMEt.getSignificanceMatrix())(1, 0) << "  "
0487                       << (noPileUpMEt.getSignificanceMatrix())(1, 1) << std::endl;
0488 
0489   // add no-PU MET object to the event
0490   auto noPileUpMEtCollection = std::make_unique<reco::PFMETCollection>();
0491   noPileUpMEtCollection->push_back(noPileUpMEt);
0492 
0493   evt.put(std::move(noPileUpMEtCollection));
0494   if (saveInputs_) {
0495     evt.put(std::move(sumLeptons), sfLeptonsName_);
0496     evt.put(std::move(sumNoPUjetOffsetEnCorr), sfNoPUjetOffsetEnCorrName_);
0497     evt.put(std::move(sumNoPUjets), sfNoPUjetsName_);
0498     evt.put(std::move(sumPUjets), sfPUjetsName_);
0499     evt.put(std::move(sumNoPUunclChargedCands), sfNoPUunclChargedCandsName_);
0500     evt.put(std::move(sumPUunclChargedCands), sfPUunclChargedCandsName_);
0501     evt.put(std::move(sumUnclNeutralCands), sfUnclNeutralCandsName_);
0502     evt.put(std::move(type0Correction_output), sfType0CorrectionName_);
0503     evt.put(std::move(sumLeptonIsoCones), sfLeptonIsoConesName_);
0504   }
0505 
0506   evt.put(std::make_unique<double>(noPileUpScaleFactor), "sfNoPU");
0507 }
0508 
0509 #include "FWCore/Framework/interface/MakerMacros.h"
0510 
0511 DEFINE_FWK_MODULE(NoPileUpPFMEtProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!