CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
METSignificance.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: METAlgorithms
4 // Class: METSignificance
5 //
11 //
12 // Original Author: Nathan Mirman (Cornell University)
13 // Created: Thu May 30 16:39:52 CDT 2013
14 //
15 //
16 
18 
19 
21 
22  edm::ParameterSet cfgParams = iConfig.getParameter<edm::ParameterSet>("parameters");
23 
24 
25  std::string ptResFileName = cfgParams.getParameter<std::string>("ptResFile");
26  std::string phiResFileName = cfgParams.getParameter<std::string>("phiResFile");
27 
28  double dRmatch = cfgParams.getParameter<double>("dRMatch");
29  dR2match_ = dRmatch*dRmatch;
30 
31  jetThreshold_ = cfgParams.getParameter<double>("jetThreshold");
32  jetEtas_ = cfgParams.getParameter<std::vector<double> >("jeta");
33  jetParams_ = cfgParams.getParameter<std::vector<double> >("jpar");
34  pjetParams_ = cfgParams.getParameter<std::vector<double> >("pjpar");
35 
36 
37  edm::FileInPath fpt("CondFormats/JetMETObjects/data/"+ptResFileName);
38  edm::FileInPath fphi("CondFormats/JetMETObjects/data/"+phiResFileName);
39 
40  ptRes_ = new JetResolution(fpt.fullPath().c_str(),false);
41  phiRes_ = new JetResolution(fphi.fullPath().c_str(),false);
42 
43 }
44 
46  delete ptRes_;
47  delete phiRes_;
48 }
49 
50 
53  const std::vector< edm::Handle<reco::CandidateView> >& leptons,
55 
56  // metsig covariance
57  double cov_xx = 0;
58  double cov_xy = 0;
59  double cov_yy = 0;
60 
61  // for lepton and jet subtraction
62  std::vector<reco::CandidatePtr> footprint;
63 
64  // subtract leptons out of sumPt
65  for ( std::vector< edm::Handle<reco::CandidateView> >::const_iterator lep_i = leptons.begin();
66  lep_i != leptons.end(); ++lep_i ) {
67  for( reco::CandidateView::const_iterator lep = (*lep_i)->begin(); lep != (*lep_i)->end(); lep++ ){
68  if( lep->pt() > 10 ){
69  for( unsigned int n=0; n < lep->numberOfSourceCandidatePtrs(); n++ ){
70  if( lep->sourceCandidatePtr(n).isNonnull() and lep->sourceCandidatePtr(n).isAvailable() ){
71  footprint.push_back(lep->sourceCandidatePtr(n));
72  }
73  }
74  }
75  }
76  }
77  // subtract jets out of sumPt
78  for(edm::View<reco::Jet>::const_iterator jet = jets.begin(); jet != jets.end(); ++jet) {
79 
80  // disambiguate jets and leptons
81  if(!cleanJet(*jet, leptons) ) continue;
82 
83  for( unsigned int n=0; n < jet->numberOfSourceCandidatePtrs(); n++){
84  if( jet->sourceCandidatePtr(n).isNonnull() and jet->sourceCandidatePtr(n).isAvailable() ){
85  footprint.push_back(jet->sourceCandidatePtr(n));
86  }
87  }
88 
89  }
90 
91  // calculate sumPt
92  double sumPt = 0;
93  for( edm::View<reco::Candidate>::const_iterator cand = pfCandidates.begin();
94  cand != pfCandidates.end(); ++cand){
95 
96  // check if candidate exists in a lepton or jet
97  bool cleancand = true;
98  for(unsigned int i=0; i < footprint.size(); i++){
99  if( footprint[i]->p4() == cand->p4() ){
100  cleancand = false;
101  }
102  }
103  // if not, add to sumPt
104  if( cleancand ){
105  sumPt += cand->pt();
106  }
107 
108  }
109 
110  // add jets to metsig covariance matrix and subtract them from sumPt
111  for(edm::View<reco::Jet>::const_iterator jet = jets.begin(); jet != jets.end(); ++jet) {
112 
113  // disambiguate jets and leptons
114  if(!cleanJet(*jet, leptons) ) continue;
115 
116  double jpt = jet->pt();
117  double jeta = jet->eta();
118  double feta = std::abs(jeta);
119  double c = jet->px()/jet->pt();
120  double s = jet->py()/jet->pt();
121 
122  // jet energy resolutions
123  double jeta_res = (std::abs(jeta) < 9.9) ? jeta : 9.89; // JetResolutions defined for |eta|<9.9
124  double sigmapt = ptRes_->parameterEtaEval("sigma",jeta_res,jpt);
125  double sigmaphi = phiRes_->parameterEtaEval("sigma",jeta_res,jpt);
126 
127  // split into high-pt and low-pt sector
128  if( jpt > jetThreshold_ ){
129  // high-pt jets enter into the covariance matrix via JER
130 
131  double scale = 0;
132  if(feta<jetEtas_[0]) scale = jetParams_[0];
133  else if(feta<jetEtas_[1]) scale = jetParams_[1];
134  else if(feta<jetEtas_[2]) scale = jetParams_[2];
135  else if(feta<jetEtas_[3]) scale = jetParams_[3];
136  else scale = jetParams_[4];
137 
138  double dpt = scale*jpt*sigmapt;
139  double dph = jpt*sigmaphi;
140 
141  cov_xx += dpt*dpt*c*c + dph*dph*s*s;
142  cov_xy += (dpt*dpt-dph*dph)*c*s;
143  cov_yy += dph*dph*c*c + dpt*dpt*s*s;
144 
145  } else {
146 
147  // add the (corrected) jet to the sumPt
148  sumPt += jpt;
149 
150  }
151 
152  }
153 
154 
155  //protection against unphysical events
156  if(sumPt<0) sumPt=0;
157 
158  // add pseudo-jet to metsig covariance matrix
159  cov_xx += pjetParams_[0]*pjetParams_[0] + pjetParams_[1]*pjetParams_[1]*sumPt;
160  cov_yy += pjetParams_[0]*pjetParams_[0] + pjetParams_[1]*pjetParams_[1]*sumPt;
161 
162  reco::METCovMatrix cov;
163  cov(0,0) = cov_xx;
164  cov(1,0) = cov_xy;
165  cov(0,1) = cov_xy;
166  cov(1,1) = cov_yy;
167 
168  return cov;
169 }
170 
171 double
173 
174  // covariance matrix determinant
175  double det = cov(0,0)*cov(1,1) - cov(0,1)*cov(1,0);
176 
177  // invert matrix
178  double ncov_xx = cov(1,1) / det;
179  double ncov_xy = -cov(0,1) / det;
180  double ncov_yy = cov(0,0) / det;
181 
182  // product of met and inverse of covariance
183  double sig = met.px()*met.px()*ncov_xx + 2*met.px()*met.py()*ncov_xy + met.py()*met.py()*ncov_yy;
184 
185  return sig;
186 }
187 
188 bool
190  const std::vector< edm::Handle<reco::CandidateView> >& leptons ){
191 
192  for ( std::vector< edm::Handle<reco::CandidateView> >::const_iterator lep_i = leptons.begin();
193  lep_i != leptons.end(); ++lep_i ) {
194  for( reco::CandidateView::const_iterator lep = (*lep_i)->begin(); lep != (*lep_i)->end(); lep++ ){
195  if ( reco::deltaR2(*lep, jet) < dR2match_ ) {
196  return false;
197  }
198  }
199  }
200  return true;
201 }
T getParameter(std::string const &) const
int i
Definition: DBlmapReader.cc:9
bool cleanJet(const reco::Jet &jet, const std::vector< edm::Handle< reco::CandidateView > > &leptons)
Base class for all types of Jets.
Definition: Jet.h:20
ROOT::Math::SMatrix< double, 2 > METCovMatrix
Definition: MET.h:40
std::vector< double > jetEtas_
std::vector< PFCandidatePtr > pfCandidates(const PFJet &jet, int particleId, bool sort=true)
const_iterator begin() const
Definition: MET.h:42
double p4[4]
Definition: TauolaWrapper.h:92
vector< PseudoJet > jets
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double deltaR2(const T1 &t1, const T2 &t2)
Definition: deltaR.h:36
std::vector< double > pjetParams_
JetResolution * phiRes_
std::vector< double > jetParams_
METSignificance(const edm::ParameterSet &iConfig)
reco::METCovMatrix getCovariance(const edm::View< reco::Jet > &jets, const std::vector< edm::Handle< reco::CandidateView > > &leptons, const edm::View< reco::Candidate > &pfCandidates)
virtual double px() const
x coordinate of momentum vector
JetResolution * ptRes_
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:85
const_iterator end() const
double getSignificance(const reco::METCovMatrix &cov, const reco::MET &met) const
std::string fullPath() const
Definition: FileInPath.cc:165
virtual double py() const
y coordinate of momentum vector