METAlgorithms/src/METSignificance.cc

0001
0002 // -*- C++ -*-
0003 //
0004 // Package:    METAlgorithms
0005 // Class:      METSignificance
0006 //
0007 /**\class METSignificance METSignificance.cc RecoMET/METAlgorithms/src/METSignificance.cc
0008 Description: [one line class summary]
0009 Implementation:
0010 [Notes on implementation]
0011 */
0012 //
0013 // Original Author:  Nathan Mirman (Cornell University)
0014 //         Created:  Thu May 30 16:39:52 CDT 2013
0015 //
0016 //
0017
0018 #include "RecoMET/METAlgorithms/interface/METSignificance.h"
0019 #include <unordered_set>
0020
0021 namespace {
0022   struct ptr_hash {
0023     std::size_t operator()(const reco::CandidatePtr& k) const {
0024       if (k.refCore().isTransient())
0025         return (unsigned long)k.refCore().productPtr() ^ k.key();
0026       else
0027         return k.refCore().id().processIndex() ^ k.refCore().id().productIndex() ^ k.key();
0028     }
0029   };
0030 }  // namespace
0031
0032 metsig::METSignificance::METSignificance(const edm::ParameterSet& iConfig) {
0033   edm::ParameterSet cfgParams = iConfig.getParameter<edm::ParameterSet>("parameters");
0034
0035   double dRmatch = cfgParams.getParameter<double>("dRMatch");
0036   dR2match_ = dRmatch * dRmatch;
0037
0038   jetThreshold_ = cfgParams.getParameter<double>("jetThreshold");
0039   jetEtas_ = cfgParams.getParameter<std::vector<double> >("jeta");
0040   jetParams_ = cfgParams.getParameter<std::vector<double> >("jpar");
0041   pjetParams_ = cfgParams.getParameter<std::vector<double> >("pjpar");
0042 }
0043
0044 metsig::METSignificance::~METSignificance() {}
0045
0046 reco::METCovMatrix metsig::METSignificance::getCovariance(const edm::View<reco::Jet>& jets,
0047                                                           const std::vector<edm::Handle<reco::CandidateView> >& leptons,
0048                                                           const edm::Handle<edm::View<reco::Candidate> >& pfCandidatesH,
0049                                                           double rho,
0050                                                           JME::JetResolution& resPtObj,
0051                                                           JME::JetResolution& resPhiObj,
0052                                                           JME::JetResolutionScaleFactor& resSFObj,
0053                                                           bool isRealData,
0054                                                           double& sumPtUnclustered,
0055                                                           edm::ValueMap<float> const* weights) {
0056   //pfcandidates
0057   const edm::View<reco::Candidate>& pfCandidates = *pfCandidatesH;
0058
0059   // metsig covariance
0060   double cov_xx = 0;
0061   double cov_xy = 0;
0062   double cov_yy = 0;
0063
0064   // for lepton and jet subtraction
0065   std::unordered_set<reco::CandidatePtr, ptr_hash> footprint;
0066
0067   // subtract leptons out of sumPtUnclustered
0068   for (const auto& lep_i : leptons) {
0069     for (const auto& lep : lep_i->ptrs()) {
0070       if (lep->pt() > 10) {
0071         for (unsigned int n = 0; n < lep->numberOfSourceCandidatePtrs(); n++)
0072           footprint.insert(lep->sourceCandidatePtr(n));
0073       }
0074     }
0075   }
0076
0077   std::vector<bool> cleanedJets(jets.size(), false);
0078   std::transform(jets.begin(), jets.end(), cleanedJets.begin(), [this, &leptons](auto const& jet) -> bool {
0079     return cleanJet(jet, leptons);
0080   });
0081   // subtract jets out of sumPtUnclustered
0082   auto iCleaned = cleanedJets.begin();
0083   for (const auto& jet : jets) {
0084     // disambiguate jets and leptons
0085     if (!(*iCleaned++))
0086       continue;
0087     for (unsigned int n = 0; n < jet.numberOfSourceCandidatePtrs(); n++) {
0088       footprint.insert(jet.sourceCandidatePtr(n));
0089     }
0090   }
0091
0092   // calculate sumPtUnclustered
0093   for (size_t i = 0; i < pfCandidates.size(); ++i) {
0094     // check if candidate exists in a lepton or jet
0095     bool cleancand = true;
0096     if (footprint.find(pfCandidates.ptrAt(i)) == footprint.end()) {
0097       float weight = (weights != nullptr) ? (*weights)[pfCandidates.ptrAt(i)] : 1.0;
0098       //dP4 recovery
0099       for (const auto& it : footprint) {
0100         if (it.isNonnull() && it.isAvailable() && (reco::deltaR2(*it, pfCandidates[i]) < 0.00000025)) {
0101           cleancand = false;
0102           break;
0103         }
0104       }
0105       // if not, add to sumPtUnclustered
0106       if (cleancand) {
0107         sumPtUnclustered += pfCandidates[i].pt() * weight;
0108       }
0109     }
0110   }
0111
0112   // add jets to metsig covariance matrix and subtract them from sumPtUnclustered
0113   iCleaned = cleanedJets.begin();
0114   for (const auto& jet : jets) {
0115     // disambiguate jets and leptons
0116     if (!(*iCleaned++))
0117       continue;
0118
0119     double jpt = jet.pt();
0120
0121     // split into high-pt and low-pt sector
0122     if (jpt > jetThreshold_) {
0123       // high-pt jets enter into the covariance matrix via JER
0124
0125       double jeta = jet.eta();
0126       double feta = std::abs(jeta);
0127       double c = jet.px() / jpt;
0128       double s = jet.py() / jpt;
0129
0130       JME::JetParameters parameters;
0131       parameters.setJetPt(jpt).setJetEta(jeta).setRho(rho);
0132
0133       // jet energy resolutions
0134       double sigmapt = resPtObj.getResolution(parameters);
0135       double sigmaphi = resPhiObj.getResolution(parameters);
0136       // SF not needed since is already embedded in the sigma in the dataGlobalTag
0137       //      double sigmaSF = isRealData ? resSFObj.getScaleFactor(parameters) : 1.0;
0138
0139       double scale = 0;
0140       if (feta < jetEtas_[0])
0141         scale = jetParams_[0];
0142       else if (feta < jetEtas_[1])
0143         scale = jetParams_[1];
0144       else if (feta < jetEtas_[2])
0145         scale = jetParams_[2];
0146       else if (feta < jetEtas_[3])
0147         scale = jetParams_[3];
0148       else
0149         scale = jetParams_[4];
0150
0151       //         double dpt = sigmaSF*scale*jpt*sigmapt;
0152       double dpt = scale * jpt * sigmapt;
0153       double dph = jpt * sigmaphi;
0154
0155       cov_xx += dpt * dpt * c * c + dph * dph * s * s;
0156       cov_xy += (dpt * dpt - dph * dph) * c * s;
0157       cov_yy += dph * dph * c * c + dpt * dpt * s * s;
0158
0159     } else {
0160       // add the (corrected) jet to the sumPtUnclustered
0161       sumPtUnclustered += jpt;
0162     }
0163   }
0164
0165   //protection against unphysical events
0166   if (sumPtUnclustered < 0)
0167     sumPtUnclustered = 0;
0168
0169   // add pseudo-jet to metsig covariance matrix
0170   double pseudoJetCov = pjetParams_[0] * pjetParams_[0] + pjetParams_[1] * pjetParams_[1] * sumPtUnclustered;
0171   cov_xx += pseudoJetCov;
0172   cov_yy += pseudoJetCov;
0173
0174   reco::METCovMatrix cov;
0175   cov(0, 0) = cov_xx;
0176   cov(1, 0) = cov_xy;
0177   cov(0, 1) = cov_xy;
0178   cov(1, 1) = cov_yy;
0179
0180   return cov;
0181 }
0182
0183 double metsig::METSignificance::getSignificance(const reco::METCovMatrix& cov, const reco::MET& met) {
0184   // covariance matrix determinant
0185   double det = cov(0, 0) * cov(1, 1) - cov(0, 1) * cov(1, 0);
0186
0187   // invert matrix
0188   double ncov_xx = cov(1, 1) / det;
0189   double ncov_xy = -cov(0, 1) / det;
0190   double ncov_yy = cov(0, 0) / det;
0191
0192   // product of met and inverse of covariance
0193   double sig = met.px() * met.px() * ncov_xx + 2 * met.px() * met.py() * ncov_xy + met.py() * met.py() * ncov_yy;
0194
0195   return sig;
0196 }
0197
0198 bool metsig::METSignificance::cleanJet(const reco::Jet& jet,
0199                                        const std::vector<edm::Handle<reco::CandidateView> >& leptons) {
0200   for (const auto& lep_i : leptons) {
0201     for (const auto& lep : *lep_i) {
0202       if (reco::deltaR2(lep, jet) < dR2match_) {
0203         return false;
0204       }
0205     }
0206   }
0207   return true;
0208 }