CMS 3D CMS Logo
MvaMEtUtilities.cc
Go to the documentation of this file.
1 #include "RecoMET/METPUSubtraction/interface/MvaMEtUtilities.h"
2 
3 #include "DataFormats/Math/interface/deltaR.h"
4 
5 #include <algorithm>
6 #include <cmath>
7 
8 MvaMEtUtilities::MvaMEtUtilities(const edm::ParameterSet& cfg) {
9  // jet id
10  /* ===> this code parses an xml it uses parameter set
11  edm::ParameterSet jetConfig = iConfig.getParameter<edm::ParameterSet>("JetIdParams");
12  for(int i0 = 0; i0 < 3; i0++) {
13  std::string lCutType = "Tight";
14  if(i0 == PileupJetIdentifier::kMedium) lCutType = "Medium";
15  if(i0 == PileupJetIdentifier::kLoose) lCutType = "Loose";
16  std::vector<double> pt010 = jetConfig.getParameter<std::vector<double> >(("Pt010_" +lCutType).c_str());
17  std::vector<double> pt1020 = jetConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
18  std::vector<double> pt2030 = jetConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
19  std::vector<double> pt3050 = jetConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
20  for(int i1 = 0; i1 < 4; i1++) mvacut_[i0][0][i1] = pt010 [i1];
21  for(int i1 = 0; i1 < 4; i1++) mvacut_[i0][1][i1] = pt1020[i1];
22  for(int i1 = 0; i1 < 4; i1++) mvacut_[i0][2][i1] = pt2030[i1];
23  for(int i1 = 0; i1 < 4; i1++) mvacut_[i0][3][i1] = pt3050[i1];
24  }
25  */
26  //Tight Id => not used
27  mvaCut_[0][0][0] = 0.5;
28  mvaCut_[0][0][1] = 0.6;
29  mvaCut_[0][0][2] = 0.6;
30  mvaCut_[0][0][3] = 0.9;
31  mvaCut_[0][1][0] = -0.2;
32  mvaCut_[0][1][1] = 0.2;
33  mvaCut_[0][1][2] = 0.2;
34  mvaCut_[0][1][3] = 0.6;
35  mvaCut_[0][2][0] = 0.3;
36  mvaCut_[0][2][1] = 0.4;
37  mvaCut_[0][2][2] = 0.7;
38  mvaCut_[0][2][3] = 0.8;
39  mvaCut_[0][3][0] = 0.5;
40  mvaCut_[0][3][1] = 0.4;
41  mvaCut_[0][3][2] = 0.8;
42  mvaCut_[0][3][3] = 0.9;
43  //Medium id => not used
44  mvaCut_[1][0][0] = 0.2;
45  mvaCut_[1][0][1] = 0.4;
46  mvaCut_[1][0][2] = 0.2;
47  mvaCut_[1][0][3] = 0.6;
48  mvaCut_[1][1][0] = -0.3;
49  mvaCut_[1][1][1] = 0.;
50  mvaCut_[1][1][2] = 0.;
51  mvaCut_[1][1][3] = 0.5;
52  mvaCut_[1][2][0] = 0.2;
53  mvaCut_[1][2][1] = 0.2;
54  mvaCut_[1][2][2] = 0.5;
55  mvaCut_[1][2][3] = 0.7;
56  mvaCut_[1][3][0] = 0.3;
57  mvaCut_[1][3][1] = 0.2;
58  mvaCut_[1][3][2] = 0.7;
59  mvaCut_[1][3][3] = 0.8;
60  //Met Id => used
61  mvaCut_[2][0][0] = -0.2;
62  mvaCut_[2][0][1] = -0.3;
63  mvaCut_[2][0][2] = -0.5;
64  mvaCut_[2][0][3] = -0.5;
65  mvaCut_[2][1][0] = -0.2;
66  mvaCut_[2][1][1] = -0.2;
67  mvaCut_[2][1][2] = -0.5;
68  mvaCut_[2][1][3] = -0.3;
69  mvaCut_[2][2][0] = -0.2;
70  mvaCut_[2][2][1] = -0.2;
71  mvaCut_[2][2][2] = -0.2;
72  mvaCut_[2][2][3] = 0.1;
73  mvaCut_[2][3][0] = -0.2;
74  mvaCut_[2][3][1] = -0.2;
75  mvaCut_[2][3][2] = 0.;
76  mvaCut_[2][3][3] = 0.2;
77 
78  dzCut_ = cfg.getParameter<double>("dZcut");
79  ptThreshold_ = (cfg.exists("ptThreshold")) ? cfg.getParameter<int>("ptThreshold") : -1000;
80 }
81 
82 MvaMEtUtilities::~MvaMEtUtilities() {
83  // nothing to be done yet...
84 }
85 
86 bool MvaMEtUtilities::passesMVA(const reco::Candidate::LorentzVector& jetP4, double mvaJetId) {
87  int ptBin = 0;
88  if (jetP4.pt() >= 10. && jetP4.pt() < 20.)
89  ptBin = 1;
90  if (jetP4.pt() >= 20. && jetP4.pt() < 30.)
91  ptBin = 2;
92  if (jetP4.pt() >= 30.)
93  ptBin = 3;
94 
95  int etaBin = 0;
96  if (std::abs(jetP4.eta()) >= 2.5 && std::abs(jetP4.eta()) < 2.75)
97  etaBin = 1;
98  if (std::abs(jetP4.eta()) >= 2.75 && std::abs(jetP4.eta()) < 3.0)
99  etaBin = 2;
100  if (std::abs(jetP4.eta()) >= 3.0 && std::abs(jetP4.eta()) < 5.0)
101  etaBin = 3;
102 
103  return (mvaJetId > mvaCut_[2][ptBin][etaBin]);
104 }
105 
106 reco::Candidate::LorentzVector MvaMEtUtilities::leadJetP4(const std::vector<reco::PUSubMETCandInfo>& jets) {
107  return jetP4(jets, 0);
108 }
109 
110 reco::Candidate::LorentzVector MvaMEtUtilities::subleadJetP4(const std::vector<reco::PUSubMETCandInfo>& jets) {
111  return jetP4(jets, 1);
112 }
113 
114 reco::Candidate::LorentzVector MvaMEtUtilities::jetP4(const std::vector<reco::PUSubMETCandInfo>& jets, unsigned idx) {
115  reco::Candidate::LorentzVector retVal(0., 0., 0., 0.);
116  if (idx < jets.size()) {
117  std::vector<reco::PUSubMETCandInfo> jets_sorted = jets;
118  std::sort(jets_sorted.rbegin(), jets_sorted.rend());
119  retVal = jets_sorted[idx].p4();
120  }
121  return retVal;
122 }
123 unsigned MvaMEtUtilities::numJetsAboveThreshold(const std::vector<reco::PUSubMETCandInfo>& jets, double ptThreshold) {
124  unsigned retVal = 0;
125  for (std::vector<reco::PUSubMETCandInfo>::const_iterator jet = jets.begin(); jet != jets.end(); ++jet) {
126  if (jet->p4().pt() > ptThreshold)
127  ++retVal;
128  }
129  return retVal;
130 }
131 std::vector<reco::PUSubMETCandInfo> MvaMEtUtilities::cleanJets(const std::vector<reco::PUSubMETCandInfo>& jets,
132  const std::vector<reco::PUSubMETCandInfo>& leptons,
133  double ptThreshold,
134  double dRmatch) {
135  double dR2match = dRmatch * dRmatch;
136  std::vector<reco::PUSubMETCandInfo> retVal;
137  for (std::vector<reco::PUSubMETCandInfo>::const_iterator jet = jets.begin(); jet != jets.end(); ++jet) {
138  bool isOverlap = false;
139  for (std::vector<reco::PUSubMETCandInfo>::const_iterator lepton = leptons.begin(); lepton != leptons.end();
140  ++lepton) {
141  if (deltaR2(jet->p4(), lepton->p4()) < dR2match)
142  isOverlap = true;
143  }
144  if (jet->p4().pt() > ptThreshold && !isOverlap)
145  retVal.push_back(*jet);
146  }
147  return retVal;
148 }
149 
150 std::vector<reco::PUSubMETCandInfo> MvaMEtUtilities::cleanPFCands(
151  const std::vector<reco::PUSubMETCandInfo>& pfCandidates,
152  const std::vector<reco::PUSubMETCandInfo>& leptons,
153  double dRmatch,
154  bool invert) {
155  double dR2match = dRmatch * dRmatch;
156  std::vector<reco::PUSubMETCandInfo> retVal;
157  for (std::vector<reco::PUSubMETCandInfo>::const_iterator pfCandidate = pfCandidates.begin();
158  pfCandidate != pfCandidates.end();
159  ++pfCandidate) {
160  bool isOverlap = false;
161  for (std::vector<reco::PUSubMETCandInfo>::const_iterator lepton = leptons.begin(); lepton != leptons.end();
162  ++lepton) {
163  if (deltaR2(pfCandidate->p4(), lepton->p4()) < dR2match)
164  isOverlap = true;
165  }
166  if ((!isOverlap && !invert) || (isOverlap && invert))
167  retVal.push_back(*pfCandidate);
168  }
169  return retVal;
170 }
171 void MvaMEtUtilities::finalize(CommonMETData& metData) {
172  metData.met = sqrt(metData.mex * metData.mex + metData.mey * metData.mey);
173  metData.mez = 0.;
174  metData.phi = atan2(metData.mey, metData.mex);
175 }
176 
177 CommonMETData MvaMEtUtilities::computeCandSum(int compKey,
178  double dZmax,
179  int dZflag,
180  bool iCharged,
181  bool mvaPassFlag,
182  const std::vector<reco::PUSubMETCandInfo>& objects) {
183  CommonMETData retVal;
184  retVal.mex = 0.;
185  retVal.mey = 0.;
186  retVal.sumet = 0.;
187 
188  for (std::vector<reco::PUSubMETCandInfo>::const_iterator object = objects.begin(); object != objects.end();
189  ++object) {
190  double pFrac = 1;
191 
192  //pf candidates
193  // dZcut
194  // maximum distance within which tracks are
195  //considered to be associated to hard scatter vertex
196  // dZflag
197  // 0 : select charged PFCandidates originating from hard scatter vertex
198  // 1 : select charged PFCandidates originating from pile-up vertices
199  // 2 : select all PFCandidates
200  if (compKey == MvaMEtUtilities::kPFCands) {
201  if (object->dZ() < 0. && dZflag != 2)
202  continue;
203  if (object->dZ() > dZmax && dZflag == 0)
204  continue;
205  if (object->dZ() < dZmax && dZflag == 1)
206  continue;
207  }
208 
209  //leptons
210  if (compKey == MvaMEtUtilities::kLeptons) {
211  if (iCharged)
212  pFrac = object->chargedEnFrac();
213  }
214 
215  //jets
216  if (compKey == MvaMEtUtilities::kJets) {
217  bool passesMVAjetId = passesMVA(object->p4(), object->mva());
218 
219  if (passesMVAjetId && !mvaPassFlag)
220  continue;
221  if (!passesMVAjetId && mvaPassFlag)
222  continue;
223 
224  pFrac = 1 - object->chargedEnFrac(); //neutral energy fraction
225  }
226 
227  retVal.mex += object->p4().px() * pFrac;
228  retVal.mey += object->p4().py() * pFrac;
229  retVal.sumet += object->p4().pt() * pFrac;
230  }
231 
232  finalize(retVal);
233  return retVal;
234 }
235 
236 CommonMETData MvaMEtUtilities::computeRecoil(int metType) {
237  CommonMETData retVal;
238 
239  if (metType == MvaMEtUtilities::kPF) {
240  //MET = pfMET = - all candidates
241  // MET (1) in JME-13-003
242  retVal.mex = leptonsSum_.mex - pfCandSum_.mex;
243  retVal.mey = leptonsSum_.mey - pfCandSum_.mey;
244  retVal.sumet = pfCandSum_.sumet - leptonsSum_.sumet;
245  }
246  if (metType == MvaMEtUtilities::kChHS) {
247  //MET = - charged HS
248  // MET (2) in JME-13-003
249  retVal.mex = leptonsChSum_.mex - pfCandChHSSum_.mex;
250  retVal.mey = leptonsChSum_.mey - pfCandChHSSum_.mey;
251  retVal.sumet = pfCandChHSSum_.sumet - leptonsChSum_.sumet;
252  }
253  if (metType == MvaMEtUtilities::kHS) {
254  //MET = - charged HS - neutral HS in jets
255  // MET (3) in JME-13-003
256  retVal.mex = leptonsChSum_.mex - (pfCandChHSSum_.mex + neutralJetHSSum_.mex);
257  retVal.mey = leptonsChSum_.mey - (pfCandChHSSum_.mey + neutralJetHSSum_.mey);
258  retVal.sumet = pfCandChHSSum_.sumet + neutralJetHSSum_.sumet - leptonsChSum_.sumet;
259  }
260  if (metType == MvaMEtUtilities::kPU) {
261  //MET = - charged PU - neutral PU in jets
262  //MET = -recoil in that particular case, - sign not useful for the MVA and then discarded
263  //motivated as PU IS its own recoil
264  // MET (4) in JME-13-003
265  retVal.mex = -(pfCandChPUSum_.mex + neutralJetPUSum_.mex);
266  retVal.mey = -(pfCandChPUSum_.mey + neutralJetPUSum_.mey);
267  retVal.sumet = pfCandChPUSum_.sumet + neutralJetPUSum_.sumet;
268  }
269  if (metType == MvaMEtUtilities::kHSMinusNeutralPU) {
270  //MET = all candidates - charged PU - neutral PU in jets
271  // = all charged HS + all neutrals - neutral PU in jets
272  // MET (5) in JME-13-003
273  retVal.mex = leptonsSum_.mex - (pfCandSum_.mex - pfCandChPUSum_.mex - neutralJetPUSum_.mex);
274  retVal.mey = leptonsSum_.mey - (pfCandSum_.mey - pfCandChPUSum_.mey - neutralJetPUSum_.mey);
275  retVal.sumet = (pfCandSum_.sumet - pfCandChPUSum_.sumet - neutralJetPUSum_.sumet) - leptonsSum_.sumet;
276  }
277 
278  finalize(retVal);
279  return retVal;
280 }
281 
282 void MvaMEtUtilities::computeAllSums(const std::vector<reco::PUSubMETCandInfo>& jets,
283  const std::vector<reco::PUSubMETCandInfo>& leptons,
284  const std::vector<reco::PUSubMETCandInfo>& pfCandidates) {
285  cleanedJets_ = cleanJets(jets, leptons, ptThreshold_, 0.5);
286 
287  leptonsSum_ = computeCandSum(kLeptons, 0., 0, false, false, leptons);
288  leptonsChSum_ = computeCandSum(kLeptons, 0., 0, true, false, leptons);
289  pfCandSum_ = computeCandSum(kPFCands, dzCut_, 2, false, false, pfCandidates);
290  pfCandChHSSum_ = computeCandSum(kPFCands, dzCut_, 0, false, false, pfCandidates);
291  pfCandChPUSum_ = computeCandSum(kPFCands, dzCut_, 1, false, false, pfCandidates);
292  neutralJetHSSum_ = computeCandSum(kJets, 0., 0, false, true, jets);
293  neutralJetPUSum_ = computeCandSum(kJets, 0., 0, false, false, jets);
294 }
295 
296 double MvaMEtUtilities::getLeptonsSumMEX() const { return leptonsSum_.mex; }
297 
298 double MvaMEtUtilities::getLeptonsSumMEY() const { return leptonsSum_.mey; }
299 
300 double MvaMEtUtilities::getLeptonsChSumMEX() const { return leptonsChSum_.mex; }
301 
302 double MvaMEtUtilities::getLeptonsChSumMEY() const { return leptonsChSum_.mey; }
303 
304 const std::vector<reco::PUSubMETCandInfo>& MvaMEtUtilities::getCleanedJets() const { return cleanedJets_; }
ewkMuLumiMonitorDQM_cfi.ptThreshold
ptThreshold
Definition: ewkMuLumiMonitorDQM_cfi.py:13
CommonMETData::met
double met
Definition: CommonMETData.h:14
zmumugammaAnalyzer_cfi.pfCandidates
pfCandidates
Definition: zmumugammaAnalyzer_cfi.py:11
resolutioncreator_cfi.object
object
Definition: resolutioncreator_cfi.py:4
MvaMEtUtilities::leadJetP4
reco::Candidate::LorentzVector leadJetP4(const std::vector< reco::PUSubMETCandInfo > &)
Definition: MvaMEtUtilities.cc:106
etaBin
int etaBin(const l1t::HGCalMulticluster *cl)
Definition: L1EGammaEEProducer.cc:19
HLT_FULL_cff.leptons
leptons
Definition: HLT_FULL_cff.py:26277
MvaMEtUtilities::leptonsSum_
CommonMETData leptonsSum_
Definition: MvaMEtUtilities.h:19
MvaMEtUtilities::cleanedJets_
std::vector< reco::PUSubMETCandInfo > cleanedJets_
Definition: MvaMEtUtilities.h:27
objects
Definition: __init__.py:1
MvaMEtUtilities::ptThreshold_
double ptThreshold_
Definition: MvaMEtUtilities.h:30
MvaMEtUtilities::pfCandSum_
CommonMETData pfCandSum_
Definition: MvaMEtUtilities.h:21
MvaMEtUtilities::subleadJetP4
reco::Candidate::LorentzVector subleadJetP4(const std::vector< reco::PUSubMETCandInfo > &)
Definition: MvaMEtUtilities.cc:110
MvaMEtUtilities::getLeptonsSumMEX
double getLeptonsSumMEX() const
Definition: MvaMEtUtilities.cc:296
CommonMETData::mez
double mez
Definition: CommonMETData.h:17
singleTopDQM_cfi.jets
jets
Definition: singleTopDQM_cfi.py:42
MvaMEtUtilities::pfCandChHSSum_
CommonMETData pfCandChHSSum_
Definition: MvaMEtUtilities.h:22
CommonMETData
Structure containing data common to all types of MET.
Definition: CommonMETData.h:12
objects.autophobj.metType
metType
Definition: autophobj.py:214
MvaMEtUtilities::passesMVA
bool passesMVA(const reco::Candidate::LorentzVector &, double)
Definition: MvaMEtUtilities.cc:86
muonRecoAnalyzer_cfi.ptBin
ptBin
Definition: muonRecoAnalyzer_cfi.py:45
heavyIonCSV_trainingSettings.idx
idx
Definition: heavyIonCSV_trainingSettings.py:5
CommonMETData::mex
double mex
Definition: CommonMETData.h:15
MvaMEtUtilities::jetP4
reco::Candidate::LorentzVector jetP4(const std::vector< reco::PUSubMETCandInfo > &, unsigned)
Definition: MvaMEtUtilities.cc:114
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
MvaMEtUtilities::getCleanedJets
const std::vector< reco::PUSubMETCandInfo > & getCleanedJets() const
Definition: MvaMEtUtilities.cc:304
MvaMEtUtilities::kPF
Definition: MvaMEtUtilities.h:16
MvaMEtUtilities::pfCandChPUSum_
CommonMETData pfCandChPUSum_
Definition: MvaMEtUtilities.h:23
CommonMETData::sumet
double sumet
Definition: CommonMETData.h:18
MvaMEtUtilities::numJetsAboveThreshold
unsigned numJetsAboveThreshold(const std::vector< reco::PUSubMETCandInfo > &, double)
Definition: MvaMEtUtilities.cc:123
CommonMETData::phi
double phi
Definition: CommonMETData.h:19
MvaMEtUtilities::computeRecoil
CommonMETData computeRecoil(int metType)
Definition: MvaMEtUtilities.cc:236
edm::ParameterSet
Definition: ParameterSet.h:47
MvaMEtUtilities::kHS
Definition: MvaMEtUtilities.h:16
deltaR.h
jetUpdater_cfi.sort
sort
Definition: jetUpdater_cfi.py:29
MvaMEtUtilities::mvaCut_
double mvaCut_[3][4][4]
Definition: MvaMEtUtilities.h:60
MvaMEtUtilities::MvaMEtUtilities
MvaMEtUtilities(const edm::ParameterSet &cfg)
Definition: MvaMEtUtilities.cc:8
MvaMEtUtilities.h
MvaMEtUtilities::kChHS
Definition: MvaMEtUtilities.h:16
looper.cfg
cfg
Definition: looper.py:297
MvaMEtUtilities::computeCandSum
CommonMETData computeCandSum(int compKey, double dZmax, int dZflag, bool iCharged, bool mvaPassFlag, const std::vector< reco::PUSubMETCandInfo > &objects)
Definition: MvaMEtUtilities.cc:177
MvaMEtUtilities::getLeptonsChSumMEX
double getLeptonsChSumMEX() const
Definition: MvaMEtUtilities.cc:300
MvaMEtUtilities::computeAllSums
void computeAllSums(const std::vector< reco::PUSubMETCandInfo > &jets, const std::vector< reco::PUSubMETCandInfo > &leptons, const std::vector< reco::PUSubMETCandInfo > &pfCandidates)
Definition: MvaMEtUtilities.cc:282
MvaMEtUtilities::neutralJetHSSum_
CommonMETData neutralJetHSSum_
Definition: MvaMEtUtilities.h:24
MvaMEtUtilities::kLeptons
Definition: MvaMEtUtilities.h:15
HLTMuonOfflineAnalyzer_cfi.deltaR2
deltaR2
Definition: HLTMuonOfflineAnalyzer_cfi.py:105
l1tstage2_dqm_sourceclient-live_cfg.invert
invert
Definition: l1tstage2_dqm_sourceclient-live_cfg.py:85
MvaMEtUtilities::kPFCands
Definition: MvaMEtUtilities.h:15
metsig::jet
Definition: SignAlgoResolutions.h:47
MvaMEtUtilities::finalize
void finalize(CommonMETData &metData)
Definition: MvaMEtUtilities.cc:171
MvaMEtUtilities::kHSMinusNeutralPU
Definition: MvaMEtUtilities.h:16
MvaMEtUtilities::~MvaMEtUtilities
virtual ~MvaMEtUtilities()
Definition: MvaMEtUtilities.cc:82
MvaMEtUtilities::cleanPFCands
std::vector< reco::PUSubMETCandInfo > cleanPFCands(const std::vector< reco::PUSubMETCandInfo > &, const std::vector< reco::PUSubMETCandInfo > &, double, bool)
Definition: MvaMEtUtilities.cc:150
MvaMEtUtilities::neutralJetPUSum_
CommonMETData neutralJetPUSum_
Definition: MvaMEtUtilities.h:25
MvaMEtUtilities::dzCut_
double dzCut_
Definition: MvaMEtUtilities.h:29
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
MvaMEtUtilities::getLeptonsChSumMEY
double getLeptonsChSumMEY() const
Definition: MvaMEtUtilities.cc:302
MvaMEtUtilities::kJets
Definition: MvaMEtUtilities.h:15
MvaMEtUtilities::leptonsChSum_
CommonMETData leptonsChSum_
Definition: MvaMEtUtilities.h:20
reco::Candidate::LorentzVector
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Candidate.h:36
CommonMETData::mey
double mey
Definition: CommonMETData.h:16
MvaMEtUtilities::kPU
Definition: MvaMEtUtilities.h:16
MvaMEtUtilities::getLeptonsSumMEY
double getLeptonsSumMEY() const
Definition: MvaMEtUtilities.cc:298
MvaMEtUtilities::cleanJets
std::vector< reco::PUSubMETCandInfo > cleanJets(const std::vector< reco::PUSubMETCandInfo > &, const std::vector< reco::PUSubMETCandInfo > &, double, double)
Definition: MvaMEtUtilities.cc:131
Generated for CMSSW Reference Manual by doxygen 1.8.16