METPUSubtraction/src/NoPileUpMEtUtilities.cc

0001 #include "RecoMET/METPUSubtraction/interface/NoPileUpMEtUtilities.h"
0002
0003 #include "DataFormats/Math/interface/deltaR.h"
0004
0005 #include <algorithm>
0006 #include <cmath>
0007
0008 NoPileUpMEtUtilities::NoPileUpMEtUtilities() {
0009   minPtDef_ = -1;
0010   maxPtDef_ = 1000000;
0011 }
0012
0013 NoPileUpMEtUtilities::~NoPileUpMEtUtilities() {}
0014
0015 // namespace NoPileUpMEtUtilities
0016 // {
0017 //-------------------------------------------------------------------------------
0018 // general auxiliary functions
0019 void NoPileUpMEtUtilities::finalizeMEtData(CommonMETData& metData) {
0020   metData.met = sqrt(metData.mex * metData.mex + metData.mey * metData.mey);
0021   metData.mez = 0.;
0022   metData.phi = atan2(metData.mey, metData.mex);
0023 }
0024 //-------------------------------------------------------------------------------
0025
0026 //-------------------------------------------------------------------------------
0027 // auxiliary functions for jets
0028 reco::PUSubMETCandInfoCollection NoPileUpMEtUtilities::cleanJets(
0029     const reco::PUSubMETCandInfoCollection& jets,
0030     const std::vector<reco::Candidate::LorentzVector>& leptons,
0031     double dRoverlap,
0032     bool invert) {
0033   reco::PUSubMETCandInfoCollection retVal;
0034   for (reco::PUSubMETCandInfoCollection::const_iterator jet = jets.begin(); jet != jets.end(); ++jet) {
0035     bool isOverlap = false;
0036     for (std::vector<reco::Candidate::LorentzVector>::const_iterator lepton = leptons.begin(); lepton != leptons.end();
0037          ++lepton) {
0038       if (deltaR2(jet->p4(), *lepton) < dRoverlap * dRoverlap) {
0039         isOverlap = true;
0040         break;
0041       }
0042     }
0043     if ((!isOverlap && !invert) || (isOverlap && invert))
0044       retVal.push_back(*jet);
0045   }
0046   return retVal;
0047 }
0048
0049 CommonMETData NoPileUpMEtUtilities::computeCandidateSum(const reco::PUSubMETCandInfoCollection& cands,
0050                                                         bool neutralFracOnly,
0051                                                         double& sumAbsPx,
0052                                                         double& sumAbsPy) {
0053   CommonMETData retVal;
0054   retVal.mex = 0.;
0055   retVal.mey = 0.;
0056   retVal.sumet = 0.;
0057   double retVal_sumAbsPx = 0.;
0058   double retVal_sumAbsPy = 0.;
0059   double pFrac = 1;
0060   for (reco::PUSubMETCandInfoCollection::const_iterator cand = cands.begin(); cand != cands.end(); ++cand) {
0061     pFrac = 1;
0062     if (neutralFracOnly)
0063       pFrac = (1 - cand->chargedEnFrac());
0064
0065     retVal.mex += cand->p4().px() * pFrac;
0066     retVal.mey += cand->p4().py() * pFrac;
0067     retVal.sumet += cand->p4().pt() * pFrac;
0068     retVal_sumAbsPx += std::abs(cand->p4().px());
0069     retVal_sumAbsPy += std::abs(cand->p4().py());
0070   }
0071   finalizeMEtData(retVal);
0072   sumAbsPx = retVal_sumAbsPx;
0073   sumAbsPy = retVal_sumAbsPy;
0074   return retVal;
0075 }
0076
0077 //-------------------------------------------------------------------------------
0078
0079 //-------------------------------------------------------------------------------
0080 // auxiliary functions for PFCandidates
0081 reco::PUSubMETCandInfoCollection NoPileUpMEtUtilities::cleanPFCandidates(
0082     const reco::PUSubMETCandInfoCollection& pfCandidates,
0083     const std::vector<reco::Candidate::LorentzVector>& leptons,
0084     double dRoverlap,
0085     bool invert) {
0086   // invert: false = PFCandidates are required not to overlap with leptons
0087   //         true  = PFCandidates are required to overlap with leptons
0088
0089   reco::PUSubMETCandInfoCollection retVal;
0090   for (reco::PUSubMETCandInfoCollection::const_iterator pfCandidate = pfCandidates.begin();
0091        pfCandidate != pfCandidates.end();
0092        ++pfCandidate) {
0093     bool isOverlap = false;
0094     for (std::vector<reco::Candidate::LorentzVector>::const_iterator lepton = leptons.begin(); lepton != leptons.end();
0095          ++lepton) {
0096       if (deltaR2(pfCandidate->p4(), *lepton) < dRoverlap * dRoverlap) {
0097         isOverlap = true;
0098         break;
0099       }
0100     }
0101     if ((!isOverlap && !invert) || (isOverlap && invert))
0102       retVal.push_back(*pfCandidate);
0103   }
0104   return retVal;
0105 }
0106
0107 reco::PUSubMETCandInfoCollection NoPileUpMEtUtilities::selectCandidates(const reco::PUSubMETCandInfoCollection& cands,
0108                                                                         double minPt,
0109                                                                         double maxPt,
0110                                                                         int type,
0111                                                                         bool isCharged,
0112                                                                         int isWithinJet) {
0113   reco::PUSubMETCandInfoCollection retVal;
0114   for (reco::PUSubMETCandInfoCollection::const_iterator cand = cands.begin(); cand != cands.end(); ++cand) {
0115     if (isCharged && cand->charge() == 0)
0116       continue;
0117     double jetPt = cand->p4().pt();
0118     if (jetPt < minPt || jetPt > maxPt)
0119       continue;
0120     if (type != reco::PUSubMETCandInfo::kUndefined && cand->type() != type)
0121       continue;
0122
0123     //for pf candidates
0124     if (isWithinJet != NoPileUpMEtUtilities::kAll && (cand->isWithinJet() != isWithinJet))
0125       continue;
0126
0127     retVal.push_back(*cand);
0128   }
0129   return retVal;
0130 }
0131 //-------------------------------------------------------------------------------
0132
0133 //-------------------------------------------------------------------------------
0134 // auxiliary functions to compute different types of MEt/hadronic recoils
0135 //
0136 // NOTE: all pfCandidates and jets passed as function arguments
0137 //       need to be cleaned wrt. leptons
0138 //
0139
0140 void NoPileUpMEtUtilities::computeAllSums(const reco::PUSubMETCandInfoCollection& jets,
0141                                           const reco::PUSubMETCandInfoCollection& pfCandidates) {
0142   reco::PUSubMETCandInfoCollection pfcsCh = selectCandidates(
0143       pfCandidates, minPtDef_, maxPtDef_, reco::PUSubMETCandInfo::kUndefined, false, NoPileUpMEtUtilities::kAll);
0144
0145   reco::PUSubMETCandInfoCollection pfcsChHS = selectCandidates(
0146       pfCandidates, minPtDef_, maxPtDef_, reco::PUSubMETCandInfo::kChHS, true, NoPileUpMEtUtilities::kAll);
0147
0148   reco::PUSubMETCandInfoCollection pfcsChPU = selectCandidates(
0149       pfCandidates, minPtDef_, maxPtDef_, reco::PUSubMETCandInfo::kChPU, true, NoPileUpMEtUtilities::kAll);
0150
0151   reco::PUSubMETCandInfoCollection pfcsNUnclustered = selectCandidates(
0152       pfCandidates, minPtDef_, maxPtDef_, reco::PUSubMETCandInfo::kNeutral, false, NoPileUpMEtUtilities::kOutsideJet);
0153
0154   reco::PUSubMETCandInfoCollection jetsHS =
0155       selectCandidates(jets, 10.0, maxPtDef_, reco::PUSubMETCandInfo::kHS, false, NoPileUpMEtUtilities::kAll);
0156   reco::PUSubMETCandInfoCollection jetsPU =
0157       selectCandidates(jets, 10.0, maxPtDef_, reco::PUSubMETCandInfo::kPU, false, NoPileUpMEtUtilities::kAll);
0158
0159   //not used so far
0160   //_chPfcSum = computeCandidateSum(pfcsCh, false, &_chPfcSumAbsPx, &_chPfcSumAbsPy);
0161   chHSPfcSum_ = computeCandidateSum(pfcsChHS, false, chHSPfcSumAbsPx_, chHSPfcSumAbsPy_);
0162   chPUPfcSum_ = computeCandidateSum(pfcsChPU, false, chPUPfcSumAbsPx_, chPUPfcSumAbsPy_);
0163   nUncPfcSum_ = computeCandidateSum(pfcsNUnclustered, false, nUncPfcSumAbsPx_, nUncPfcSumAbsPy_);
0164
0165   nHSJetSum_ = computeCandidateSum(jetsHS, true, nHSJetSumAbsPx_, nHSJetSumAbsPy_);
0166   nPUJetSum_ = computeCandidateSum(jetsPU, true, nPUJetSumAbsPx_, nPUJetSumAbsPy_);
0167 }
0168
0169 CommonMETData NoPileUpMEtUtilities::computeRecoil(int metType, double& sumAbsPx, double& sumAbsPy) {
0170   CommonMETData retVal;
0171   double retSumAbsPx = 0.;
0172   double retSumAbsPy = 0.;
0173
0174   //never used....
0175   // if(metType==NoPileUpMEtUtilities::kChMET) {
0176   //   CommonMETData chPfcSum = computeCandidateSum(_pfcsCh, false, &trackSumAbsPx, &trackSumAbsPy);
0177   //   retVal.mex   = -chPfcSum.mex;
0178   //   retVal.mey   = -chPfcSum.mey;
0179   //   retVal.sumet = chPfcSum.sumet;
0180   //   retSumAbsPx = ;
0181   //   retSumAbsPy = ;
0182   // }
0183   if (metType == NoPileUpMEtUtilities::kChHSMET) {
0184     retVal.mex = chHSPfcSum_.mex;
0185     retVal.mey = chHSPfcSum_.mey;
0186     retVal.sumet = chHSPfcSum_.sumet;
0187     retSumAbsPx = chHSPfcSumAbsPx_;
0188     retSumAbsPy = chHSPfcSumAbsPy_;
0189   }
0190   if (metType == NoPileUpMEtUtilities::kChPUMET) {
0191     retVal.mex = chPUPfcSum_.mex;
0192     retVal.mey = chPUPfcSum_.mey;
0193     retVal.sumet = chPUPfcSum_.sumet;
0194     retSumAbsPx = chPUPfcSumAbsPx_;
0195     retSumAbsPy = chPUPfcSumAbsPy_;
0196   }
0197   if (metType == NoPileUpMEtUtilities::kNeutralUncMET) {
0198     retVal.mex = nUncPfcSum_.mex;
0199     retVal.mey = nUncPfcSum_.mey;
0200     retVal.sumet = nUncPfcSum_.sumet;
0201     retSumAbsPx = nUncPfcSumAbsPx_;
0202     retSumAbsPy = nUncPfcSumAbsPy_;
0203   }
0204   if (metType == NoPileUpMEtUtilities::kHadronicHSMET) {
0205     retVal.mex = chHSPfcSum_.mex + nHSJetSum_.mex;
0206     retVal.mey = chHSPfcSum_.mey + nHSJetSum_.mey;
0207     retVal.sumet = chHSPfcSum_.sumet + nHSJetSum_.sumet;
0208     retSumAbsPx = chHSPfcSumAbsPx_ + nHSJetSumAbsPx_;
0209     retSumAbsPy = chHSPfcSumAbsPy_ + nHSJetSumAbsPy_;
0210   }
0211   if (metType == NoPileUpMEtUtilities::kHadronicPUMET) {
0212     retVal.mex = chPUPfcSum_.mex + nHSJetSum_.mex;
0213     retVal.mey = chPUPfcSum_.mey + nHSJetSum_.mey;
0214     retVal.sumet = chPUPfcSum_.sumet + nHSJetSum_.sumet;
0215     retSumAbsPx = chPUPfcSumAbsPx_ + nHSJetSumAbsPx_;
0216     retSumAbsPy = chPUPfcSumAbsPy_ + nHSJetSumAbsPy_;
0217   }
0218
0219   sumAbsPx = retSumAbsPx;
0220   sumAbsPy = retSumAbsPy;
0221
0222   return retVal;
0223 }