CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PileupJetIdAlgo.cc
Go to the documentation of this file.
1 #include "RecoJets/JetProducers/interface/PileupJetIdAlgo.h"
2 
3 #include "DataFormats/JetReco/interface/PFJet.h"
4 #include "DataFormats/JetReco/interface/Jet.h"
5 #include "DataFormats/VertexReco/interface/Vertex.h"
6 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
7 #include "DataFormats/Math/interface/deltaR.h"
8 #include "FWCore/ParameterSet/interface/FileInPath.h"
9 #include "TMatrixDSym.h"
10 #include "TMatrixDSymEigen.h"
11 
12 // ------------------------------------------------------------------------------------------
13 const float large_val = std::numeric_limits<float>::max();
14 
15 // ------------------------------------------------------------------------------------------
16 PileupJetIdAlgo::PileupJetIdAlgo(const edm::ParameterSet & ps)
17 {
18  impactParTkThreshod_ = 1.;
19  cutBased_ = false;
20  //std::string label = ps.getParameter<std::string>("label");
21  cutBased_ = ps.getParameter<bool>("cutBased");
22  if(!cutBased_)
23  {
24  tmvaWeights_ = edm::FileInPath(ps.getParameter<std::string>("tmvaWeights")).fullPath();
25  tmvaMethod_ = ps.getParameter<std::string>("tmvaMethod");
26  tmvaVariables_ = ps.getParameter<std::vector<std::string> >("tmvaVariables");
27  tmvaSpectators_ = ps.getParameter<std::vector<std::string> >("tmvaSpectators");
28  version_ = ps.getParameter<int>("version");
29  }
30  reader_ = 0;
31  edm::ParameterSet jetConfig = ps.getParameter<edm::ParameterSet>("JetIdParams");
32  for(int i0 = 0; i0 < 3; i0++) {
33  std::string lCutType = "Tight";
34  if(i0 == PileupJetIdentifier::kMedium) lCutType = "Medium";
35  if(i0 == PileupJetIdentifier::kLoose) lCutType = "Loose";
36  int nCut = 1;
37  if(cutBased_) nCut++;
38  for(int i1 = 0; i1 < nCut; i1++) {
39  std::string lFullCutType = lCutType;
40  if(cutBased_ && i1 == 0) lFullCutType = "BetaStar"+ lCutType;
41  if(cutBased_ && i1 == 1) lFullCutType = "RMS" + lCutType;
42  std::vector<double> pt010 = jetConfig.getParameter<std::vector<double> >(("Pt010_" +lFullCutType).c_str());
43  std::vector<double> pt1020 = jetConfig.getParameter<std::vector<double> >(("Pt1020_"+lFullCutType).c_str());
44  std::vector<double> pt2030 = jetConfig.getParameter<std::vector<double> >(("Pt2030_"+lFullCutType).c_str());
45  std::vector<double> pt3050 = jetConfig.getParameter<std::vector<double> >(("Pt3050_"+lFullCutType).c_str());
46  if(!cutBased_) {
47  for(int i2 = 0; i2 < 4; i2++) mvacut_[i0][0][i2] = pt010 [i2];
48  for(int i2 = 0; i2 < 4; i2++) mvacut_[i0][1][i2] = pt1020[i2];
49  for(int i2 = 0; i2 < 4; i2++) mvacut_[i0][2][i2] = pt2030[i2];
50  for(int i2 = 0; i2 < 4; i2++) mvacut_[i0][3][i2] = pt3050[i2];
51  }
52  if(cutBased_ && i1 == 0) {
53  for(int i2 = 0; i2 < 4; i2++) betaStarCut_[i0][0][i2] = pt010 [i2];
54  for(int i2 = 0; i2 < 4; i2++) betaStarCut_[i0][1][i2] = pt1020[i2];
55  for(int i2 = 0; i2 < 4; i2++) betaStarCut_[i0][2][i2] = pt2030[i2];
56  for(int i2 = 0; i2 < 4; i2++) betaStarCut_[i0][3][i2] = pt3050[i2];
57  }
58  if(cutBased_ && i1 == 1) {
59  for(int i2 = 0; i2 < 4; i2++) rmsCut_[i0][0][i2] = pt010 [i2];
60  for(int i2 = 0; i2 < 4; i2++) rmsCut_[i0][1][i2] = pt1020[i2];
61  for(int i2 = 0; i2 < 4; i2++) rmsCut_[i0][2][i2] = pt2030[i2];
62  for(int i2 = 0; i2 < 4; i2++) rmsCut_[i0][3][i2] = pt3050[i2];
63  }
64  }
65  }
66  setup();
67 }
68 
69 
70 // ------------------------------------------------------------------------------------------
71 PileupJetIdAlgo::PileupJetIdAlgo(int version,
72  const std::string & tmvaWeights,
73  const std::string & tmvaMethod,
74  Float_t impactParTkThreshod,
75  const std::vector<std::string> & tmvaVariables
76  )
77 {
78  impactParTkThreshod_ = impactParTkThreshod;
79  tmvaWeights_ = tmvaWeights;
80  tmvaMethod_ = tmvaMethod;
81  tmvaVariables_ = tmvaVariables;
82  version_ = version;
83 
84  reader_ = 0;
85 
86  setup();
87 }
88 
89 // ------------------------------------------------------------------------------------------
90 void PileupJetIdAlgo::setup()
91 {
92  initVariables();
93 
94  if( version_ == PHILv0 ) {
95  tmvaVariables_.clear();
96  tmvaVariables_.push_back( "jspt_1" );
97  tmvaVariables_.push_back( "jseta_1" );
98  tmvaVariables_.push_back( "jsphi_1" );
99  tmvaVariables_.push_back( "jd0_1" );
100  tmvaVariables_.push_back( "jdZ_1" );
101  tmvaVariables_.push_back( "jm_1" );
102  tmvaVariables_.push_back( "npart_1" );
103  tmvaVariables_.push_back( "lpt_1" );
104  tmvaVariables_.push_back( "leta_1" );
105  tmvaVariables_.push_back( "lphi_1" );
106  tmvaVariables_.push_back( "spt_1" );
107  tmvaVariables_.push_back( "seta_1" );
108  tmvaVariables_.push_back( "sphi_1" );
109  tmvaVariables_.push_back( "lnept_1" );
110  tmvaVariables_.push_back( "lneeta_1");
111  tmvaVariables_.push_back( "lnephi_1");
112  tmvaVariables_.push_back( "lempt_1" );
113  tmvaVariables_.push_back( "lemeta_1");
114  tmvaVariables_.push_back( "lemphi_1");
115  tmvaVariables_.push_back( "lchpt_1" );
116  // tmvaVariables_.push_back( "lcheta_1"); FIXME missing in weights file
117  tmvaVariables_.push_back( "lchphi_1");
118  tmvaVariables_.push_back( "lLfr_1" );
119  tmvaVariables_.push_back( "drlc_1" );
120  tmvaVariables_.push_back( "drls_1" );
121  tmvaVariables_.push_back( "drm_1" );
122  tmvaVariables_.push_back( "drmne_1" );
123  tmvaVariables_.push_back( "drem_1" );
124  tmvaVariables_.push_back( "drch_1" );
125 
126  tmvaNames_["jspt_1"] = "jetPt";
127  tmvaNames_["jseta_1"] = "jetEta";
128  tmvaNames_["jsphi_1"] = "jetPhi";
129  tmvaNames_["jm_1"] = "jetM";
130  tmvaNames_["jd0_1"] = "d0";
131  tmvaNames_["jdZ_1"] = "dZ";
132  tmvaNames_["npart_1"] = "nParticles";
133 
134  tmvaNames_["lpt_1"] = "leadPt";
135  tmvaNames_["leta_1"] = "leadEta";
136  tmvaNames_["lphi_1"] = "leadPhi";
137  tmvaNames_["spt_1"] = "secondPt";
138  tmvaNames_["seta_1"] = "secondEta";
139  tmvaNames_["sphi_1"] = "secondPhi";
140  tmvaNames_["lnept_1"] = "leadNeutPt";
141  tmvaNames_["lneeta_1"] = "leadNeutEta";
142  tmvaNames_["lnephi_1"] = "leadNeutPhi";
143  tmvaNames_["lempt_1"] = "leadEmPt";
144  tmvaNames_["lemeta_1"] = "leadEmEta";
145  tmvaNames_["lemphi_1"] = "leadEmPhi";
146  tmvaNames_["lchpt_1"] = "leadChPt";
147  tmvaNames_["lcheta_1"] = "leadChEta";
148  tmvaNames_["lchphi_1"] = "leadChPhi";
149  tmvaNames_["lLfr_1"] = "leadFrac";
150 
151  tmvaNames_["drlc_1"] = "dRLeadCent";
152  tmvaNames_["drls_1"] = "dRLead2nd";
153  tmvaNames_["drm_1"] = "dRMean";
154  tmvaNames_["drmne_1"] = "dRMeanNeut";
155  tmvaNames_["drem_1"] = "dRMeanEm";
156  tmvaNames_["drch_1"] = "dRMeanCh";
157 
158  } else if( ! cutBased_ ){
159  assert( tmvaMethod_.empty() || (! tmvaVariables_.empty() && version_ == USER) );
160  }
161 }
162 
163 // ------------------------------------------------------------------------------------------
164 PileupJetIdAlgo::~PileupJetIdAlgo()
165 {
166  if( reader_ ) {
167  delete reader_;
168  }
169 }
170 
171 // ------------------------------------------------------------------------------------------
172 void assign(const std::vector<float> & vec, float & a, float & b, float & c, float & d )
173 {
174  size_t sz = vec.size();
175  a = ( sz > 0 ? vec[0] : 0. );
176  b = ( sz > 1 ? vec[1] : 0. );
177  c = ( sz > 2 ? vec[2] : 0. );
178  d = ( sz > 3 ? vec[3] : 0. );
179 }
180 // ------------------------------------------------------------------------------------------
181 void setPtEtaPhi(const reco::Candidate & p, float & pt, float & eta, float &phi )
182 {
183  pt = p.pt();
184  eta = p.eta();
185  phi = p.phi();
186 }
187 
188 // ------------------------------------------------------------------------------------------
189 void PileupJetIdAlgo::bookReader()
190 {
191  reader_ = new TMVA::Reader("!Color:!Silent");
192  assert( ! tmvaMethod_.empty() && ! tmvaWeights_.empty() );
193  for(std::vector<std::string>::iterator it=tmvaVariables_.begin(); it!=tmvaVariables_.end(); ++it) {
194  if( tmvaNames_[*it].empty() ) {
195  tmvaNames_[*it] = *it;
196  }
197  reader_->AddVariable( *it, variables_[ tmvaNames_[*it] ].first );
198  }
199  for(std::vector<std::string>::iterator it=tmvaSpectators_.begin(); it!=tmvaSpectators_.end(); ++it) {
200  if( tmvaNames_[*it].empty() ) {
201  tmvaNames_[*it] = *it;
202  }
203  reader_->AddSpectator( *it, variables_[ tmvaNames_[*it] ].first );
204  }
205  reader_->BookMVA(tmvaMethod_.c_str(), tmvaWeights_.c_str());
206 }
207 
208 // ------------------------------------------------------------------------------------------
209 void PileupJetIdAlgo::set(const PileupJetIdentifier & id)
210 {
211  internalId_ = id;
212 }
213 
214 // ------------------------------------------------------------------------------------------
215 void PileupJetIdAlgo::runMva()
216 {
217  if( cutBased_ ) {
218  internalId_.idFlag_ = computeCutIDflag(internalId_.betaStarClassic_,internalId_.dR2Mean_,internalId_.nvtx_,internalId_.jetPt_,internalId_.jetEta_);
219  } else {
220  if( ! reader_ ) { bookReader(); std::cerr << "Reader booked" << std::endl; }
221  if(fabs(internalId_.jetEta_) < 5.0) internalId_.mva_ = reader_->EvaluateMVA( tmvaMethod_.c_str() );
222  if(fabs(internalId_.jetEta_) >= 5.0) internalId_.mva_ = -2.;
223  internalId_.idFlag_ = computeIDflag(internalId_.mva_,internalId_.jetPt_,internalId_.jetEta_);
224  }
225 }
226 
227 // ------------------------------------------------------------------------------------------
228 std::pair<int,int> PileupJetIdAlgo::getJetIdKey(float jetPt, float jetEta)
229 {
230  int ptId = 0;
231  if(jetPt > 10 && jetPt < 20) ptId = 1;
232  if(jetPt > 20 && jetPt < 30) ptId = 2;
233  if(jetPt > 30 ) ptId = 3;
234 
235  int etaId = 0;
236  if(fabs(jetEta) > 2.5 && fabs(jetEta) < 2.75) etaId = 1;
237  if(fabs(jetEta) > 2.75 && fabs(jetEta) < 3.0 ) etaId = 2;
238  if(fabs(jetEta) > 3.0 && fabs(jetEta) < 5.0 ) etaId = 3;
239 
240  return std::pair<int,int>(ptId,etaId);
241 }
242 // ------------------------------------------------------------------------------------------
243 int PileupJetIdAlgo::computeCutIDflag(float betaStarClassic,float dR2Mean,float nvtx, float jetPt, float jetEta)
244 {
245  std::pair<int,int> jetIdKey = getJetIdKey(jetPt,jetEta);
246  float betaStarModified = betaStarClassic/log(nvtx-0.64);
247  int idFlag(0);
248  if(betaStarModified < betaStarCut_[PileupJetIdentifier::kTight ][jetIdKey.first][jetIdKey.second] &&
249  dR2Mean < rmsCut_ [PileupJetIdentifier::kTight ][jetIdKey.first][jetIdKey.second]
250  ) idFlag += 1 << PileupJetIdentifier::kTight;
251 
252  if(betaStarModified < betaStarCut_[PileupJetIdentifier::kMedium ][jetIdKey.first][jetIdKey.second] &&
253  dR2Mean < rmsCut_ [PileupJetIdentifier::kMedium ][jetIdKey.first][jetIdKey.second]
254  ) idFlag += 1 << PileupJetIdentifier::kMedium;
255 
256  if(betaStarModified < betaStarCut_[PileupJetIdentifier::kLoose ][jetIdKey.first][jetIdKey.second] &&
257  dR2Mean < rmsCut_ [PileupJetIdentifier::kLoose ][jetIdKey.first][jetIdKey.second]
258  ) idFlag += 1 << PileupJetIdentifier::kLoose;
259  return idFlag;
260 }
261 // ------------------------------------------------------------------------------------------
262 int PileupJetIdAlgo::computeIDflag(float mva, float jetPt, float jetEta)
263 {
264  std::pair<int,int> jetIdKey = getJetIdKey(jetPt,jetEta);
265  return computeIDflag(mva,jetIdKey.first,jetIdKey.second);
266 }
267 
268 // ------------------------------------------------------------------------------------------
269 int PileupJetIdAlgo::computeIDflag(float mva,int ptId,int etaId)
270 {
271  int idFlag(0);
272  if(mva > mvacut_[PileupJetIdentifier::kTight ][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kTight;
273  if(mva > mvacut_[PileupJetIdentifier::kMedium][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kMedium;
274  if(mva > mvacut_[PileupJetIdentifier::kLoose ][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kLoose;
275  return idFlag;
276 }
277 
278 
279 // ------------------------------------------------------------------------------------------
280 PileupJetIdentifier PileupJetIdAlgo::computeMva()
281 {
282  runMva();
283  return PileupJetIdentifier(internalId_);
284 }
285 
286 // ------------------------------------------------------------------------------------------
287 PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet * jet, float jec, const reco::Vertex * vtx,
288  const reco::VertexCollection & allvtx,
289  bool calculateMva)
290 {
291  static int printWarning = 10;
292  typedef std::vector <reco::PFCandidatePtr> constituents_type;
293  typedef std::vector <reco::PFCandidatePtr>::iterator constituents_iterator;
294 
295  // initialize all variables to 0
296  resetVariables();
297 
298  // loop over constituents, accumulate sums and find leading candidates
299  //const pat::Jet * patjet = dynamic_cast<const pat::Jet *>(jet);
300  const reco::PFJet * pfjet = dynamic_cast<const reco::PFJet *>(jet);
301  //assert( patjet != 0 || pfjet != 0 );
302  //if( patjet != 0 && jec == 0. ) { // if this is a pat jet and no jec has been passed take the jec from the object
303  //jec = patjet->pt()/patjet->correctedJet(0).pt();
304  //}
305  if( jec < 0. ) {
306  jec = 1.;
307  }
308  //constituents_type constituents = pfjet ? pfjet->getPFConstituents() : patjet->getPFConstituents();
309  constituents_type constituents = pfjet->getPFConstituents();
310 
311  reco::PFCandidatePtr lLead, lSecond, lLeadNeut, lLeadEm, lLeadCh, lTrail;
312  std::vector<float> frac, fracCh, fracEm, fracNeut;
313  float cones[] = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7 };
314  size_t ncones = sizeof(cones)/sizeof(float);
315  float * coneFracs[] = { &internalId_.frac01_, &internalId_.frac02_, &internalId_.frac03_, &internalId_.frac04_,
316  &internalId_.frac05_, &internalId_.frac06_, &internalId_.frac07_ };
317  float * coneEmFracs[] = { &internalId_.emFrac01_, &internalId_.emFrac02_, &internalId_.emFrac03_, &internalId_.emFrac04_,
318  &internalId_.emFrac05_, &internalId_.emFrac06_, &internalId_.emFrac07_ };
319  float * coneNeutFracs[] = { &internalId_.neutFrac01_, &internalId_.neutFrac02_, &internalId_.neutFrac03_, &internalId_.neutFrac04_,
320  &internalId_.neutFrac05_, &internalId_.neutFrac06_, &internalId_.neutFrac07_ };
321  float * coneChFracs[] = { &internalId_.chFrac01_, &internalId_.chFrac02_, &internalId_.chFrac03_, &internalId_.chFrac04_,
322  &internalId_.chFrac05_, &internalId_.chFrac06_, &internalId_.chFrac07_ };
323  TMatrixDSym covMatrix(2); covMatrix = 0.;
324 
325  reco::TrackRef impactTrack;
326  float jetPt = jet->pt() / jec; // use uncorrected pt for shape variables
327  float sumPt = 0., sumPt2 = 0., sumTkPt = 0.,sumPtCh=0,sumPtNe = 0;
328  setPtEtaPhi(*jet,internalId_.jetPt_,internalId_.jetEta_,internalId_.jetPhi_); // use corrected pt for jet kinematics
329  internalId_.jetM_ = jet->mass();
330  internalId_.nvtx_ = allvtx.size();
331 
332  for(constituents_iterator it=constituents.begin(); it!=constituents.end(); ++it) {
333  reco::PFCandidatePtr & icand = *it;
334  float candPt = icand->pt();
335  float candPtFrac = candPt/jetPt;
336  float candDr = reco::deltaR(**it,*jet);
337  float candDeta = fabs( (*it)->eta() - jet->eta() );
338  float candDphi = reco::deltaPhi(**it,*jet);
339  float candPtDr = candPt * candDr;
340  size_t icone = std::lower_bound(&cones[0],&cones[ncones],candDr) - &cones[0];
341 
342  // all particles
343  if( lLead.isNull() || candPt > lLead->pt() ) {
344  lSecond = lLead;
345  lLead = icand;
346  } else if( (lSecond.isNull() || candPt > lSecond->pt()) && (candPt < lLead->pt()) ) {
347  lSecond = icand;
348  }
349 
350  // average shapes
351  internalId_.dRMean_ += candPtDr;
352  internalId_.dR2Mean_ += candPtDr*candPtDr;
353  covMatrix(0,0) += candPt*candPt*candDeta*candDeta;
354  covMatrix(0,1) += candPt*candPt*candDeta*candDphi;
355  covMatrix(1,1) += candPt*candPt*candDphi*candDphi;
356  internalId_.ptD_ += candPt*candPt;
357  sumPt += candPt;
358  sumPt2 += candPt*candPt;
359 
360  // single most energetic candiates and jet shape profiles
361  frac.push_back(candPtFrac);
362  if( icone < ncones ) { *coneFracs[icone] += candPt; }
363 
364  // neutrals
365  if( icand->particleId() == reco::PFCandidate::h0 ) {
366  if (lLeadNeut.isNull() || candPt > lLeadNeut->pt()) { lLeadNeut = icand; }
367  internalId_.dRMeanNeut_ += candPtDr;
368  fracNeut.push_back(candPtFrac);
369  if( icone < ncones ) { *coneNeutFracs[icone] += candPt; }
370  internalId_.ptDNe_ += candPt*candPt;
371  sumPtNe += candPt;
372  }
373  // EM candidated
374  if( icand->particleId() == reco::PFCandidate::gamma ) {
375  if(lLeadEm.isNull() || candPt > lLeadEm->pt()) { lLeadEm = icand; }
376  internalId_.dRMeanEm_ += candPtDr;
377  fracEm.push_back(candPtFrac);
378  if( icone < ncones ) { *coneEmFracs[icone] += candPt; }
379  internalId_.ptDNe_ += candPt*candPt;
380  sumPtNe += candPt;
381  }
382  // Charged particles
383  if( icand->trackRef().isNonnull() && icand->trackRef().isAvailable() ) {
384  if (lLeadCh.isNull() || candPt > lLeadCh->pt()) { lLeadCh = icand; }
385  internalId_.dRMeanCh_ += candPtDr;
386  internalId_.ptDCh_ += candPt*candPt;
387  fracCh.push_back(candPtFrac);
388  if( icone < ncones ) { *coneChFracs[icone] += candPt; }
389  sumPtCh += candPt;
390  }
391  // beta and betastar
392  if( icand->trackRef().isNonnull() && icand->trackRef().isAvailable() ) {
393  try {
394  float tkpt = icand->trackRef()->pt();
395  sumTkPt += tkpt;
396  // 'classic' beta definition based on track-vertex association
397  bool inVtx0 = find( vtx->tracks_begin(), vtx->tracks_end(), reco::TrackBaseRef(icand->trackRef()) ) != vtx->tracks_end();
398  bool inAnyOther = false;
399  // alternative beta definition based on track-vertex distance of closest approach
400  double dZ0 = fabs(icand->trackRef()->dz(vtx->position()));
401  double dZ = dZ0;
402  for(reco::VertexCollection::const_iterator vi=allvtx.begin(); vi!=allvtx.end(); ++vi ) {
403  const reco::Vertex & iv = *vi;
404  if( iv.isFake() || iv.ndof() < 4 ) { continue; }
405  // the primary vertex may have been copied by the user: check identity by position
406  bool isVtx0 = (iv.position() - vtx->position()).r() < 0.02;
407  // 'classic' beta definition: check if the track is associated with any vertex other than the primary one
408  if( ! isVtx0 && ! inAnyOther ) {
409  inAnyOther = find( iv.tracks_begin(), iv.tracks_end(), reco::TrackBaseRef(icand->trackRef()) ) != iv.tracks_end();
410  }
411  // alternative beta: find closest vertex to the track
412  dZ = std::min(dZ,fabs(icand->trackRef()->dz(iv.position())));
413  }
414  // classic beta/betaStar
415  if( inVtx0 && ! inAnyOther ) {
416  internalId_.betaClassic_ += tkpt;
417  } else if( ! inVtx0 && inAnyOther ) {
418  internalId_.betaStarClassic_ += tkpt;
419  }
420  // alternative beta/betaStar
421  if( dZ0 < 0.2 ) {
422  internalId_.beta_ += tkpt;
423  } else if( dZ < 0.2 ) {
424  internalId_.betaStar_ += tkpt;
425  }
426  } catch (cms::Exception &e) {
427  if(printWarning-- > 0) { std::cerr << e << std::endl; }
428  }
429  }
430  // trailing candidate
431  if( lTrail.isNull() || candPt < lTrail->pt() ) {
432  lTrail = icand;
433  }
434  }
435 
436  // Finalize all variables
437  assert( lLead.isNonnull() );
438 
439  if ( lSecond.isNull() ) { lSecond = lTrail; }
440  if ( lLeadNeut.isNull() ) { lLeadNeut = lTrail; }
441  if ( lLeadEm.isNull() ) { lLeadEm = lTrail; }
442  if ( lLeadCh.isNull() ) { lLeadCh = lTrail; }
443  impactTrack = lLeadCh->trackRef();
444 
445  internalId_.nCharged_ = pfjet->chargedMultiplicity();
446  internalId_.nNeutrals_ = pfjet->neutralMultiplicity();
447  internalId_.chgEMfrac_ = pfjet->chargedEmEnergy() /jet->energy();
448  internalId_.neuEMfrac_ = pfjet->neutralEmEnergy() /jet->energy();
449  internalId_.chgHadrfrac_ = pfjet->chargedHadronEnergy()/jet->energy();
450  internalId_.neuHadrfrac_ = pfjet->neutralHadronEnergy()/jet->energy();
451 
452  if( impactTrack.isNonnull() && impactTrack.isAvailable() ) {
453  internalId_.d0_ = fabs(impactTrack->dxy(vtx->position()));
454  internalId_.dZ_ = fabs(impactTrack->dz(vtx->position()));
455  } else {
456  if(printWarning-- > 0) { std::cerr << "WARNING : did not find any valid track reference attached to the jet " << std::endl; }
457  }
458  internalId_.nParticles_ = constituents.size();
459 
460  setPtEtaPhi(*lLead,internalId_.leadPt_,internalId_.leadEta_,internalId_.leadPhi_);
461  setPtEtaPhi(*lSecond,internalId_.secondPt_,internalId_.secondEta_,internalId_.secondPhi_);
462  setPtEtaPhi(*lLeadNeut,internalId_.leadNeutPt_,internalId_.leadNeutEta_,internalId_.leadNeutPhi_);
463  setPtEtaPhi(*lLeadEm,internalId_.leadEmPt_,internalId_.leadEmEta_,internalId_.leadEmPhi_);
464  setPtEtaPhi(*lLeadCh,internalId_.leadChPt_,internalId_.leadChEta_,internalId_.leadChPhi_);
465 
466  std::sort(frac.begin(),frac.end(),std::greater<float>());
467  std::sort(fracCh.begin(),fracCh.end(),std::greater<float>());
468  std::sort(fracEm.begin(),fracEm.end(),std::greater<float>());
469  std::sort(fracNeut.begin(),fracNeut.end(),std::greater<float>());
470  assign(frac, internalId_.leadFrac_, internalId_.secondFrac_, internalId_.thirdFrac_, internalId_.fourthFrac_);
471  assign(fracCh, internalId_.leadChFrac_, internalId_.secondChFrac_, internalId_.thirdChFrac_, internalId_.fourthChFrac_);
472  assign(fracEm, internalId_.leadEmFrac_, internalId_.secondEmFrac_, internalId_.thirdEmFrac_, internalId_.fourthEmFrac_);
473  assign(fracNeut,internalId_.leadNeutFrac_,internalId_.secondNeutFrac_,internalId_.thirdNeutFrac_,internalId_.fourthNeutFrac_);
474 
475  covMatrix(0,0) /= sumPt2;
476  covMatrix(0,1) /= sumPt2;
477  covMatrix(1,1) /= sumPt2;
478  covMatrix(1,0) = covMatrix(0,1);
479  internalId_.etaW_ = sqrt(covMatrix(0,0));
480  internalId_.phiW_ = sqrt(covMatrix(1,1));
481  internalId_.jetW_ = 0.5*(internalId_.etaW_+internalId_.phiW_);
482  TVectorD eigVals(2); eigVals = TMatrixDSymEigen(covMatrix).GetEigenValues();
483  internalId_.majW_ = sqrt(fabs(eigVals(0)));
484  internalId_.minW_ = sqrt(fabs(eigVals(1)));
485  if( internalId_.majW_ < internalId_.minW_ ) { std::swap(internalId_.majW_,internalId_.minW_); }
486 
487  internalId_.dRLeadCent_ = reco::deltaR(*jet,*lLead);
488  if( lSecond.isNonnull() ) { internalId_.dRLead2nd_ = reco::deltaR(*jet,*lSecond); }
489  internalId_.dRMean_ /= jetPt;
490  internalId_.dRMeanNeut_ /= jetPt;
491  internalId_.dRMeanEm_ /= jetPt;
492  internalId_.dRMeanCh_ /= jetPt;
493  internalId_.dR2Mean_ /= sumPt2;
494 
495  for(size_t ic=0; ic<ncones; ++ic){
496  *coneFracs[ic] /= jetPt;
497  *coneEmFracs[ic] /= jetPt;
498  *coneNeutFracs[ic] /= jetPt;
499  *coneChFracs[ic] /= jetPt;
500  }
501  //http://jets.physics.harvard.edu/qvg/
502  double ptMean = sumPt/internalId_.nParticles_;
503  double ptRMS = 0;
504  for(unsigned int i0 = 0; i0 < frac.size(); i0++) {ptRMS+=(frac[i0]-ptMean)*(frac[i0]-ptMean);}
505  ptRMS/=internalId_.nParticles_;
506  ptRMS=sqrt(ptRMS);
507 
508  internalId_.ptMean_ = ptMean;
509  internalId_.ptRMS_ = ptRMS/jetPt;
510  internalId_.pt2A_ = sqrt( internalId_.ptD_ /internalId_.nParticles_)/jetPt;
511  internalId_.ptD_ = sqrt( internalId_.ptD_) / sumPt;
512  internalId_.ptDCh_ = sqrt( internalId_.ptDCh_) / sumPtCh;
513  internalId_.ptDNe_ = sqrt( internalId_.ptDNe_) / sumPtNe;
514  internalId_.sumPt_ = sumPt;
515  internalId_.sumChPt_ = sumPtCh;
516  internalId_.sumNePt_ = sumPtNe;
517 
518  if( sumTkPt != 0. ) {
519  internalId_.beta_ /= sumTkPt;
520  internalId_.betaStar_ /= sumTkPt;
521  internalId_.betaClassic_ /= sumTkPt;
522  internalId_.betaStarClassic_ /= sumTkPt;
523  } else {
524  assert( internalId_.beta_ == 0. && internalId_.betaStar_ == 0.&& internalId_.betaClassic_ == 0. && internalId_.betaStarClassic_ == 0. );
525  }
526 
527  if( calculateMva ) {
528  runMva();
529  }
530 
531  return PileupJetIdentifier(internalId_);
532 }
533 
534 
535 
536 // ------------------------------------------------------------------------------------------
537 std::string PileupJetIdAlgo::dumpVariables() const
538 {
539  std::stringstream out;
540  for(variables_list_t::const_iterator it=variables_.begin();
541  it!=variables_.end(); ++it ) {
542  out << std::setw(15) << it->first << std::setw(3) << "="
543  << std::setw(5) << *it->second.first
544  << " (" << std::setw(5) << it->second.second << ")" << std::endl;
545  }
546  return out.str();
547 }
548 
549 // ------------------------------------------------------------------------------------------
550 void PileupJetIdAlgo::resetVariables()
551 {
552  internalId_.idFlag_ = 0;
553  for(variables_list_t::iterator it=variables_.begin();
554  it!=variables_.end(); ++it ) {
555  *it->second.first = it->second.second;
556  }
557 }
558 
559 // ------------------------------------------------------------------------------------------
560 #define INIT_VARIABLE(NAME,TMVANAME,VAL) \
561  internalId_.NAME ## _ = VAL; \
562  variables_[ # NAME ] = std::make_pair(& internalId_.NAME ## _, VAL);
563 
564 // ------------------------------------------------------------------------------------------
565 void PileupJetIdAlgo::initVariables()
566 {
567  internalId_.idFlag_ = 0;
568  INIT_VARIABLE(mva , "", -100.);
569  INIT_VARIABLE(jetPt , "jspt_1", 0.);
570  INIT_VARIABLE(jetEta , "jseta_1", large_val);
571  INIT_VARIABLE(jetPhi , "jsphi_1", large_val);
572  INIT_VARIABLE(jetM , "jm_1", 0.);
573 
574  INIT_VARIABLE(nCharged , "", 0.);
575  INIT_VARIABLE(nNeutrals , "", 0.);
576 
577  INIT_VARIABLE(chgEMfrac , "", 0.);
578  INIT_VARIABLE(neuEMfrac , "", 0.);
579  INIT_VARIABLE(chgHadrfrac, "", 0.);
580  INIT_VARIABLE(neuHadrfrac, "", 0.);
581 
582  INIT_VARIABLE(d0 , "jd0_1" , -1000.);
583  INIT_VARIABLE(dZ , "jdZ_1" , -1000.);
584  INIT_VARIABLE(nParticles , "npart_1" , 0.);
585 
586  INIT_VARIABLE(leadPt , "lpt_1" , 0.);
587  INIT_VARIABLE(leadEta , "leta_1" , large_val);
588  INIT_VARIABLE(leadPhi , "lphi_1" , large_val);
589  INIT_VARIABLE(secondPt , "spt_1" , 0.);
590  INIT_VARIABLE(secondEta , "seta_1" , large_val);
591  INIT_VARIABLE(secondPhi , "sphi_1" , large_val);
592  INIT_VARIABLE(leadNeutPt , "lnept_1" , 0.);
593  INIT_VARIABLE(leadNeutEta, "lneeta_1" , large_val);
594  INIT_VARIABLE(leadNeutPhi, "lnephi_1" , large_val);
595  INIT_VARIABLE(leadEmPt , "lempt_1" , 0.);
596  INIT_VARIABLE(leadEmEta , "lemeta_1" , large_val);
597  INIT_VARIABLE(leadEmPhi , "lemphi_1" , large_val);
598  INIT_VARIABLE(leadChPt , "lchpt_1" , 0.);
599  INIT_VARIABLE(leadChEta , "lcheta_1" , large_val);
600  INIT_VARIABLE(leadChPhi , "lchphi_1" , large_val);
601  INIT_VARIABLE(leadFrac , "lLfr_1" , 0.);
602 
603  INIT_VARIABLE(dRLeadCent , "drlc_1" , 0.);
604  INIT_VARIABLE(dRLead2nd , "drls_1" , 0.);
605  INIT_VARIABLE(dRMean , "drm_1" , 0.);
606  INIT_VARIABLE(dRMeanNeut , "drmne_1" , 0.);
607  INIT_VARIABLE(dRMeanEm , "drem_1" , 0.);
608  INIT_VARIABLE(dRMeanCh , "drch_1" , 0.);
609  INIT_VARIABLE(dR2Mean , "" , 0.);
610 
611  INIT_VARIABLE(ptD , "", 0.);
612  INIT_VARIABLE(ptMean , "", 0.);
613  INIT_VARIABLE(ptRMS , "", 0.);
614  INIT_VARIABLE(pt2A , "", 0.);
615  INIT_VARIABLE(ptDCh , "", 0.);
616  INIT_VARIABLE(ptDNe , "", 0.);
617  INIT_VARIABLE(sumPt , "", 0.);
618  INIT_VARIABLE(sumChPt , "", 0.);
619  INIT_VARIABLE(sumNePt , "", 0.);
620 
621  INIT_VARIABLE(secondFrac ,"" ,0.);
622  INIT_VARIABLE(thirdFrac ,"" ,0.);
623  INIT_VARIABLE(fourthFrac ,"" ,0.);
624 
625  INIT_VARIABLE(leadChFrac ,"" ,0.);
626  INIT_VARIABLE(secondChFrac ,"" ,0.);
627  INIT_VARIABLE(thirdChFrac ,"" ,0.);
628  INIT_VARIABLE(fourthChFrac ,"" ,0.);
629 
630  INIT_VARIABLE(leadNeutFrac ,"" ,0.);
631  INIT_VARIABLE(secondNeutFrac ,"" ,0.);
632  INIT_VARIABLE(thirdNeutFrac ,"" ,0.);
633  INIT_VARIABLE(fourthNeutFrac ,"" ,0.);
634 
635  INIT_VARIABLE(leadEmFrac ,"" ,0.);
636  INIT_VARIABLE(secondEmFrac ,"" ,0.);
637  INIT_VARIABLE(thirdEmFrac ,"" ,0.);
638  INIT_VARIABLE(fourthEmFrac ,"" ,0.);
639 
640  INIT_VARIABLE(jetW ,"" ,1.);
641  INIT_VARIABLE(etaW ,"" ,1.);
642  INIT_VARIABLE(phiW ,"" ,1.);
643 
644  INIT_VARIABLE(majW ,"" ,1.);
645  INIT_VARIABLE(minW ,"" ,1.);
646 
647  INIT_VARIABLE(frac01 ,"" ,0.);
648  INIT_VARIABLE(frac02 ,"" ,0.);
649  INIT_VARIABLE(frac03 ,"" ,0.);
650  INIT_VARIABLE(frac04 ,"" ,0.);
651  INIT_VARIABLE(frac05 ,"" ,0.);
652  INIT_VARIABLE(frac06 ,"" ,0.);
653  INIT_VARIABLE(frac07 ,"" ,0.);
654 
655  INIT_VARIABLE(chFrac01 ,"" ,0.);
656  INIT_VARIABLE(chFrac02 ,"" ,0.);
657  INIT_VARIABLE(chFrac03 ,"" ,0.);
658  INIT_VARIABLE(chFrac04 ,"" ,0.);
659  INIT_VARIABLE(chFrac05 ,"" ,0.);
660  INIT_VARIABLE(chFrac06 ,"" ,0.);
661  INIT_VARIABLE(chFrac07 ,"" ,0.);
662 
663  INIT_VARIABLE(neutFrac01 ,"" ,0.);
664  INIT_VARIABLE(neutFrac02 ,"" ,0.);
665  INIT_VARIABLE(neutFrac03 ,"" ,0.);
666  INIT_VARIABLE(neutFrac04 ,"" ,0.);
667  INIT_VARIABLE(neutFrac05 ,"" ,0.);
668  INIT_VARIABLE(neutFrac06 ,"" ,0.);
669  INIT_VARIABLE(neutFrac07 ,"" ,0.);
670 
671  INIT_VARIABLE(emFrac01 ,"" ,0.);
672  INIT_VARIABLE(emFrac02 ,"" ,0.);
673  INIT_VARIABLE(emFrac03 ,"" ,0.);
674  INIT_VARIABLE(emFrac04 ,"" ,0.);
675  INIT_VARIABLE(emFrac05 ,"" ,0.);
676  INIT_VARIABLE(emFrac06 ,"" ,0.);
677  INIT_VARIABLE(emFrac07 ,"" ,0.);
678 
679  INIT_VARIABLE(beta ,"" ,0.);
680  INIT_VARIABLE(betaStar ,"" ,0.);
681  INIT_VARIABLE(betaClassic ,"" ,0.);
682  INIT_VARIABLE(betaStarClassic ,"" ,0.);
683 
684  INIT_VARIABLE(nvtx ,"" ,0.);
685 
686 }
687 
688 #undef INIT_VARIABLE
large_val
const float large_val
Definition: PileupJetIdAlgo.cc:13
beta
const double beta
Definition: CosmicMuonParameters.h:69
reco::LeafCandidate::energy
virtual double energy() const GCC11_FINAL
energy
Definition: LeafCandidate.h:134
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
PileupJetIdAlgo::set
void set(const PileupJetIdentifier &)
Definition: PileupJetIdAlgo.cc:209
PileupJetIdAlgo::impactParTkThreshod_
Float_t impactParTkThreshod_
Definition: PileupJetIdAlgo.h:73
PileupJetIdAlgo::tmvaVariables_
std::vector< std::string > tmvaVariables_
Definition: PileupJetIdAlgo.h:68
INIT_VARIABLE
#define INIT_VARIABLE(NAME, TMVANAME, VAL)
Definition: PileupJetIdAlgo.cc:560
PileupJetIdAlgo::mvacut_
Float_t mvacut_[3][4][4]
Definition: PileupJetIdAlgo.h:75
reco::PFJet::chargedEmEnergy
float chargedEmEnergy() const
chargedEmEnergy
Definition: PFJet.h:139
PileupJetIdAlgo::reader_
TMVA::Reader * reader_
Definition: PileupJetIdAlgo.h:66
reco::Vertex::tracks_end
trackRef_iterator tracks_end() const
last iterator over tracks
Definition: Vertex.cc:45
reco::Candidate::eta
virtual float eta() const =0
momentum pseudorapidity
reco::Jet
Base class for all types of Jets.
Definition: Jet.h:21
PileupJetIdAlgo::getJetIdKey
std::pair< int, int > getJetIdKey(float jetPt, float jetEta)
Definition: PileupJetIdAlgo.cc:228
create_public_lumi_plots.log
log
Definition: create_public_lumi_plots.py:1119
PileupJetIdAlgo::tmvaNames_
std::map< std::string, std::string > tmvaNames_
Definition: PileupJetIdAlgo.h:70
PileupJetIdAlgo::tmvaMethod_
std::string tmvaMethod_
Definition: PileupJetIdAlgo.h:67
reco::Candidate::phi
virtual float phi() const =0
momentum azimuthal angle
reco::VertexCollection
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
min
#define min(a, b)
Definition: mlp_lapack.h:161
PileupJetIdAlgo::computeCutIDflag
int computeCutIDflag(float betaStarClassic, float dR2Mean, float nvtx, float jetPt, float jetEta)
Definition: PileupJetIdAlgo.cc:243
edm::Ref::isAvailable
bool isAvailable() const
Definition: Ref.h:276
eta
T eta() const
Definition: Basic3DVectorLD.h:179
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
patTestJEC_cfi.jec
tuple jec
Definition: patTestJEC_cfi.py:5
PileupJetIdAlgo::betaStarCut_
Float_t betaStarCut_[3][4][4]
Definition: PileupJetIdAlgo.h:77
reco::PFJet::chargedMultiplicity
int chargedMultiplicity() const
chargedMultiplicity
Definition: PFJet.h:152
reco::Vertex::position
const Point & position() const
position
Definition: Vertex.h:93
reco::PFJet
Jets made from PFObjects.
Definition: PFJet.h:22
reco::PFJet::neutralEmEnergy
float neutralEmEnergy() const
neutralEmEnergy
Definition: PFJet.h:147
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:250
edm::Ptr::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ptr.h:152
reco::Candidate::pt
virtual float pt() const =0
transverse momentum
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
reco::TrackBaseRef
edm::RefToBase< reco::Track > TrackBaseRef
persistent reference to a Track, using views
Definition: TrackFwd.h:22
PileupJetIdAlgo::internalId_
PileupJetIdentifier internalId_
Definition: PileupJetIdAlgo.h:63
std::swap
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
Definition: DataFrameContainer.h:223
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
PileupJetIdAlgo::dumpVariables
std::string dumpVariables() const
Definition: PileupJetIdAlgo.cc:537
nParticles
#define nParticles
Definition: Histograms.cc:24
reco::deltaR
auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:30
assign
void assign(const std::vector< float > &vec, float &a, float &b, float &c, float &d)
Definition: PileupJetIdAlgo.cc:172
PileupJetIdAlgo::computeIdVariables
PileupJetIdentifier computeIdVariables(const reco::Jet *jet, float jec, const reco::Vertex *, const reco::VertexCollection &, bool calculateMva=false)
Definition: PileupJetIdAlgo.cc:287
first
bool first
Definition: L1TdeRCT.cc:94
PileupJetIdAlgo::computeIDflag
int computeIDflag(float mva, float jetPt, float jetEta)
Definition: PileupJetIdAlgo.cc:262
reco::Vertex::ndof
double ndof() const
Definition: Vertex.h:89
reco::PFJet::neutralMultiplicity
int neutralMultiplicity() const
neutralMultiplicity
Definition: PFJet.h:154
reco::LeafCandidate::eta
virtual float eta() const GCC11_FINAL
momentum pseudorapidity
Definition: LeafCandidate.h:159
dbtoconf.out
tuple out
Definition: dbtoconf.py:99
BeamSplash_cfg.version
tuple version
Definition: BeamSplash_cfg.py:45
setPtEtaPhi
void setPtEtaPhi(const reco::Candidate &p, float &pt, float &eta, float &phi)
Definition: PileupJetIdAlgo.cc:181
reco::deltaPhi
double deltaPhi(double phi1, double phi2)
Definition: deltaPhi.h:12
reco::LeafCandidate::mass
virtual float mass() const GCC11_FINAL
mass
Definition: LeafCandidate.h:138
reco::Vertex::isFake
bool isFake() const
Definition: Vertex.h:65
reco::tau::helpers::nCharged
unsigned int nCharged(const GenJet &jet)
Definition: RecoTauObjectEmbedderPlugin.cc:29
PileupJetIdAlgo::rmsCut_
Float_t rmsCut_[3][4][4]
Definition: PileupJetIdAlgo.h:76
PileupJetIdAlgo::PileupJetIdAlgo
PileupJetIdAlgo(int version=PHILv0, const std::string &tmvaWeight="", const std::string &tmvaMethod="", Float_t impactParTkThreshod_=1., const std::vector< std::string > &tmvaVariables=std::vector< std::string >())
Definition: PileupJetIdAlgo.cc:71
PileupJetIdAlgo::computeMva
PileupJetIdentifier computeMva()
Definition: PileupJetIdAlgo.cc:280
b
double b
Definition: hdecay.h:120
PileupJetIdAlgo::tmvaSpectators_
std::vector< std::string > tmvaSpectators_
Definition: PileupJetIdAlgo.h:69
dtNoiseDBValidation_cfg.cerr
tuple cerr
Definition: dtNoiseDBValidation_cfg.py:22
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
edm::Ptr::isNull
bool isNull() const
Checks for null.
Definition: Ptr.h:148
PileupJetIdAlgo::tmvaWeights_
std::string tmvaWeights_
Definition: PileupJetIdAlgo.h:67
a
double a
Definition: hdecay.h:121
PileupJetIdAlgo::variables_
variables_list_t variables_
Definition: PileupJetIdAlgo.h:64
reco::PFJet::getPFConstituents
virtual std::vector< reco::PFCandidatePtr > getPFConstituents() const
get all constituents
Definition: PFJet.cc:53
edm::Ptr::isAvailable
bool isAvailable() const
Definition: Ptr.h:158
reco::PFJet::neutralHadronEnergy
float neutralHadronEnergy() const
neutralHadronEnergy
Definition: PFJet.h:99
edm::FileInPath::fullPath
std::string fullPath() const
Definition: FileInPath.cc:171
reco::Vertex::tracks_begin
trackRef_iterator tracks_begin() const
first iterator over tracks
Definition: Vertex.cc:40
reco::LeafCandidate::pt
virtual float pt() const GCC11_FINAL
transverse momentum
Definition: LeafCandidate.h:153
reco::PFJet::chargedHadronEnergy
float chargedHadronEnergy() const
chargedHadronEnergy
Definition: PFJet.h:95
phi
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5