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.
2 
10 
11 #include "TMatrixDSym.h"
12 #include "TMatrixDSymEigen.h"
13 
14 // ------------------------------------------------------------------------------------------
16 
17 // ------------------------------------------------------------------------------------------
19 {
21  cutBased_ = false;
22  etaBinnedWeights_ = false;
23  runMvas_=runMvas;
24  //std::string label = ps.getParameter<std::string>("label");
25  cutBased_ = ps.getParameter<bool>("cutBased");
26  if(!cutBased_)
27  {
28  etaBinnedWeights_ = ps.getParameter<bool>("etaBinnedWeights");
31  tmvaWeights_jteta_2_2p5_ = edm::FileInPath(ps.getParameter<std::string>("tmvaWeights_jteta_2_2p5")).fullPath();
32  tmvaWeights_jteta_2p5_3_ = edm::FileInPath(ps.getParameter<std::string>("tmvaWeights_jteta_2p5_3")).fullPath();
34  }
35  else{
37  }
38  tmvaMethod_ = ps.getParameter<std::string>("tmvaMethod");
40  tmvaVariables_jteta_0_3_ = ps.getParameter<std::vector<std::string> >("tmvaVariables_jteta_0_3");
41  tmvaVariables_jteta_3_5_ = ps.getParameter<std::vector<std::string> >("tmvaVariables_jteta_3_5");
42  } else {
43  tmvaVariables_ = ps.getParameter<std::vector<std::string> >("tmvaVariables");
44  }
45  tmvaSpectators_ = ps.getParameter<std::vector<std::string> >("tmvaSpectators");
46  version_ = ps.getParameter<int>("version");
47  }
48  else version_ = USER;
49  edm::ParameterSet jetConfig = ps.getParameter<edm::ParameterSet>("JetIdParams");
50  for(int i0 = 0; i0 < 3; i0++) {
51  std::string lCutType = "Tight";
52  if(i0 == PileupJetIdentifier::kMedium) lCutType = "Medium";
53  if(i0 == PileupJetIdentifier::kLoose) lCutType = "Loose";
54  int nCut = 1;
55  if(cutBased_) nCut++;
56  for(int i1 = 0; i1 < nCut; i1++) {
57  std::string lFullCutType = lCutType;
58  if(cutBased_ && i1 == 0) lFullCutType = "BetaStar"+ lCutType;
59  if(cutBased_ && i1 == 1) lFullCutType = "RMS" + lCutType;
60  std::vector<double> pt010 = jetConfig.getParameter<std::vector<double> >(("Pt010_" +lFullCutType).c_str());
61  std::vector<double> pt1020 = jetConfig.getParameter<std::vector<double> >(("Pt1020_"+lFullCutType).c_str());
62  std::vector<double> pt2030 = jetConfig.getParameter<std::vector<double> >(("Pt2030_"+lFullCutType).c_str());
63  std::vector<double> pt3050 = jetConfig.getParameter<std::vector<double> >(("Pt3050_"+lFullCutType).c_str());
64  if(!cutBased_) {
65  for(int i2 = 0; i2 < 4; i2++) mvacut_[i0][0][i2] = pt010 [i2];
66  for(int i2 = 0; i2 < 4; i2++) mvacut_[i0][1][i2] = pt1020[i2];
67  for(int i2 = 0; i2 < 4; i2++) mvacut_[i0][2][i2] = pt2030[i2];
68  for(int i2 = 0; i2 < 4; i2++) mvacut_[i0][3][i2] = pt3050[i2];
69  }
70  if(cutBased_ && i1 == 0) {
71  for(int i2 = 0; i2 < 4; i2++) betaStarCut_[i0][0][i2] = pt010 [i2];
72  for(int i2 = 0; i2 < 4; i2++) betaStarCut_[i0][1][i2] = pt1020[i2];
73  for(int i2 = 0; i2 < 4; i2++) betaStarCut_[i0][2][i2] = pt2030[i2];
74  for(int i2 = 0; i2 < 4; i2++) betaStarCut_[i0][3][i2] = pt3050[i2];
75  }
76  if(cutBased_ && i1 == 1) {
77  for(int i2 = 0; i2 < 4; i2++) rmsCut_[i0][0][i2] = pt010 [i2];
78  for(int i2 = 0; i2 < 4; i2++) rmsCut_[i0][1][i2] = pt1020[i2];
79  for(int i2 = 0; i2 < 4; i2++) rmsCut_[i0][2][i2] = pt2030[i2];
80  for(int i2 = 0; i2 < 4; i2++) rmsCut_[i0][3][i2] = pt3050[i2];
81  }
82  }
83  }
84  setup();
85 }
86 
87 
88 // ------------------------------------------------------------------------------------------
90  const std::string & tmvaWeights,
91  const std::string & tmvaMethod,
92  Float_t impactParTkThreshod,
93  const std::vector<std::string> & tmvaVariables,
94  bool runMvas
95  )
96 {
97  impactParTkThreshod_ = impactParTkThreshod;
101  version_ = version;
102 
103  runMvas_=runMvas;
104 
105  setup();
106 }
107 
108 // ------------------------------------------------------------------------------------------
110 {
111  initVariables();
112 
113  if( ! cutBased_ ){
114  assert( tmvaMethod_.empty() || ((! tmvaVariables_.empty() || ( !tmvaVariables_jteta_0_3_.empty() && !tmvaVariables_jteta_3_5_.empty() )) && version_ == USER) );
115  }
116  if(( ! cutBased_ ) && (runMvas_)) { bookReader();}
117 }
118 
119 // ------------------------------------------------------------------------------------------
121 {
122 }
123 
124 // ------------------------------------------------------------------------------------------
125 void assign(const std::vector<float> & vec, float & a, float & b, float & c, float & d )
126 {
127  size_t sz = vec.size();
128  a = ( sz > 0 ? vec[0] : 0. );
129  b = ( sz > 1 ? vec[1] : 0. );
130  c = ( sz > 2 ? vec[2] : 0. );
131  d = ( sz > 3 ? vec[3] : 0. );
132 }
133 // ------------------------------------------------------------------------------------------
134 void setPtEtaPhi(const reco::Candidate & p, float & pt, float & eta, float &phi )
135 {
136  pt = p.pt();
137  eta = p.eta();
138  phi = p.phi();
139 }
140 
141 // ------------------------------------------------------------------------------------------
143 {
144  if(etaBinnedWeights_){
145  reader_jteta_0_2_ = std::unique_ptr<TMVA::Reader>(new TMVA::Reader("!Color:Silent"));
146  reader_jteta_2_2p5_ = std::unique_ptr<TMVA::Reader>(new TMVA::Reader("!Color:Silent"));
147  reader_jteta_2p5_3_ = std::unique_ptr<TMVA::Reader>(new TMVA::Reader("!Color:Silent"));
148  reader_jteta_3_5_ = std::unique_ptr<TMVA::Reader>(new TMVA::Reader("!Color:Silent"));
149  } else {
150  reader_ = std::unique_ptr<TMVA::Reader>(new TMVA::Reader("!Color:Silent"));
151  }
152  if(etaBinnedWeights_){
153  for(std::vector<std::string>::iterator it=tmvaVariables_jteta_0_3_.begin(); it!=tmvaVariables_jteta_0_3_.end(); ++it) {
154  if( tmvaNames_[*it].empty() ) {
155  tmvaNames_[*it] = *it;
156  }
157  reader_jteta_0_2_->AddVariable( *it, variables_[ tmvaNames_[*it] ].first );
158  reader_jteta_2_2p5_->AddVariable( *it, variables_[ tmvaNames_[*it] ].first );
159  reader_jteta_2p5_3_->AddVariable( *it, variables_[ tmvaNames_[*it] ].first );
160  }
161  for(std::vector<std::string>::iterator it=tmvaVariables_jteta_3_5_.begin(); it!=tmvaVariables_jteta_3_5_.end(); ++it) {
162  if( tmvaNames_[*it].empty() ) {
163  tmvaNames_[*it] = *it;
164  }
165  reader_jteta_3_5_->AddVariable( *it, variables_[ tmvaNames_[*it] ].first );
166  }
167  } else {
168  for(std::vector<std::string>::iterator it=tmvaVariables_.begin(); it!=tmvaVariables_.end(); ++it) {
169  if( tmvaNames_[*it].empty() ) {
170  tmvaNames_[*it] = *it;
171  }
172  reader_->AddVariable( *it, variables_[ tmvaNames_[*it] ].first );
173  }
174  }
175  for(std::vector<std::string>::iterator it=tmvaSpectators_.begin(); it!=tmvaSpectators_.end(); ++it) {
176  if( tmvaNames_[*it].empty() ) {
177  tmvaNames_[*it] = *it;
178  }
179  if(etaBinnedWeights_){
180  reader_jteta_0_2_->AddSpectator( *it, variables_[ tmvaNames_[*it] ].first );
181  reader_jteta_2_2p5_->AddSpectator( *it, variables_[ tmvaNames_[*it] ].first );
182  reader_jteta_2p5_3_->AddSpectator( *it, variables_[ tmvaNames_[*it] ].first );
183  reader_jteta_3_5_->AddSpectator( *it, variables_[ tmvaNames_[*it] ].first );
184  } else {
185  reader_->AddSpectator( *it, variables_[ tmvaNames_[*it] ].first );
186  }
187  }
188  if(etaBinnedWeights_){
193  } else {
195  }
196 }
197 
198 // ------------------------------------------------------------------------------------------
200 {
201  internalId_ = id;
202 }
203 
204 // ------------------------------------------------------------------------------------------
206 {
207  if( cutBased_ ) {
208  internalId_.idFlag_ = computeCutIDflag(internalId_.betaStarClassic_,internalId_.dR2Mean_,internalId_.nvtx_,internalId_.jetPt_,internalId_.jetEta_);
209  } else {
210  if(std::abs(internalId_.jetEta_) >= 5.0) {
211  internalId_.mva_ = -2.;
212  } else {
213  if(etaBinnedWeights_){
214  if(std::abs(internalId_.jetEta_)<=2.) internalId_.mva_ = reader_jteta_0_2_->EvaluateMVA( tmvaMethod_.c_str() );
215  else if(std::abs(internalId_.jetEta_)<=2.5) internalId_.mva_ = reader_jteta_2_2p5_->EvaluateMVA( tmvaMethod_.c_str() );
216  else if(std::abs(internalId_.jetEta_)<=3.) internalId_.mva_ = reader_jteta_2p5_3_->EvaluateMVA( tmvaMethod_.c_str() );
217  else internalId_.mva_ = reader_jteta_3_5_->EvaluateMVA( tmvaMethod_.c_str() );
218  } else {
219  internalId_.mva_ = reader_->EvaluateMVA( tmvaMethod_.c_str() );
220  }
221  }
223  }
224 }
225 
226 // ------------------------------------------------------------------------------------------
227 std::pair<int,int> PileupJetIdAlgo::getJetIdKey(float jetPt, float jetEta)
228 {
229  int ptId = 0;
230  if(jetPt >= 10 && jetPt < 20) ptId = 1;
231  if(jetPt >= 20 && jetPt < 30) ptId = 2;
232  if(jetPt >= 30 ) ptId = 3;
233 
234  int etaId = 0;
235  if(std::abs(jetEta) >= 2.5 && std::abs(jetEta) < 2.75) etaId = 1;
236  if(std::abs(jetEta) >= 2.75 && std::abs(jetEta) < 3.0 ) etaId = 2;
237  if(std::abs(jetEta) >= 3.0 && std::abs(jetEta) < 5.0 ) etaId = 3;
238 
239  return std::pair<int,int>(ptId,etaId);
240 }
241 // ------------------------------------------------------------------------------------------
242 int PileupJetIdAlgo::computeCutIDflag(float betaStarClassic,float dR2Mean,float nvtx, float jetPt, float jetEta)
243 {
244  std::pair<int,int> jetIdKey = getJetIdKey(jetPt,jetEta);
245  float betaStarModified = betaStarClassic/log(nvtx-0.64);
246  int idFlag(0);
247  if(betaStarModified < betaStarCut_[PileupJetIdentifier::kTight ][jetIdKey.first][jetIdKey.second] &&
248  dR2Mean < rmsCut_ [PileupJetIdentifier::kTight ][jetIdKey.first][jetIdKey.second]
249  ) idFlag += 1 << PileupJetIdentifier::kTight;
250 
251  if(betaStarModified < betaStarCut_[PileupJetIdentifier::kMedium ][jetIdKey.first][jetIdKey.second] &&
252  dR2Mean < rmsCut_ [PileupJetIdentifier::kMedium ][jetIdKey.first][jetIdKey.second]
253  ) idFlag += 1 << PileupJetIdentifier::kMedium;
254 
255  if(betaStarModified < betaStarCut_[PileupJetIdentifier::kLoose ][jetIdKey.first][jetIdKey.second] &&
256  dR2Mean < rmsCut_ [PileupJetIdentifier::kLoose ][jetIdKey.first][jetIdKey.second]
257  ) idFlag += 1 << PileupJetIdentifier::kLoose;
258  return idFlag;
259 }
260 // ------------------------------------------------------------------------------------------
262 {
263  std::pair<int,int> jetIdKey = getJetIdKey(jetPt,jetEta);
264  return computeIDflag(mva,jetIdKey.first,jetIdKey.second);
265 }
266 
267 // ------------------------------------------------------------------------------------------
268 int PileupJetIdAlgo::computeIDflag(float mva,int ptId,int etaId)
269 {
270  int idFlag(0);
271  if(mva > mvacut_[PileupJetIdentifier::kTight ][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kTight;
272  if(mva > mvacut_[PileupJetIdentifier::kMedium][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kMedium;
273  if(mva > mvacut_[PileupJetIdentifier::kLoose ][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kLoose;
274  return idFlag;
275 }
276 
277 
278 // ------------------------------------------------------------------------------------------
280 {
281  runMva();
283 }
284 
285 // ------------------------------------------------------------------------------------------
287  const reco::VertexCollection & allvtx, double rho)
288 {
289 
290  static std::atomic<int> printWarning{10};
291 
292  // initialize all variables to 0
293  resetVariables();
294 
295  // loop over constituents, accumulate sums and find leading candidates
296  const pat::Jet * patjet = dynamic_cast<const pat::Jet *>(jet);
297  const reco::PFJet * pfjet = dynamic_cast<const reco::PFJet *>(jet);
298  assert( patjet != nullptr || pfjet != nullptr );
299  if( patjet != nullptr && jec == 0. ) { // if this is a pat jet and no jec has been passed take the jec from the object
300  jec = patjet->pt()/patjet->correctedJet(0).pt();
301  }
302  if( jec <= 0. ) {
303  jec = 1.;
304  }
305 
306  const reco::Candidate* lLead = nullptr, *lSecond = nullptr, *lLeadNeut = nullptr, *lLeadEm = nullptr, *lLeadCh = nullptr, *lTrail = nullptr;
307  std::vector<float> frac, fracCh, fracEm, fracNeut;
308  float cones[] = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7 };
309  size_t ncones = sizeof(cones)/sizeof(float);
310  float * coneFracs[] = { &internalId_.frac01_, &internalId_.frac02_, &internalId_.frac03_, &internalId_.frac04_,
311  &internalId_.frac05_, &internalId_.frac06_, &internalId_.frac07_ };
318  TMatrixDSym covMatrix(2); covMatrix = 0.;
319  float jetPt = jet->pt() / jec; // use uncorrected pt for shape variables
320  float sumPt = 0., sumPt2 = 0., sumTkPt = 0.,sumPtCh=0,sumPtNe = 0;
321  setPtEtaPhi(*jet,internalId_.jetPt_,internalId_.jetEta_,internalId_.jetPhi_); // use corrected pt for jet kinematics
322  internalId_.jetM_ = jet->mass();
323  internalId_.nvtx_ = allvtx.size();
324  internalId_.rho_ = rho;
325 
326  float dRmin(1000);
327 
328  for ( unsigned i = 0; i < jet->numberOfSourceCandidatePtrs(); ++i ) {
329  reco::CandidatePtr pfJetConstituent = jet->sourceCandidatePtr(i);
330 
331  const reco::Candidate* icand = pfJetConstituent.get();
332  const pat::PackedCandidate* lPack = dynamic_cast<const pat::PackedCandidate *>( icand );
333  const reco::PFCandidate *lPF=dynamic_cast<const reco::PFCandidate*>( icand );
334  bool isPacked = true;
335  if (lPack == nullptr){
336  isPacked = false;
337  }
338 
339  float candPt = icand->pt();
340  float candPtFrac = candPt/jetPt;
341  float candDr = reco::deltaR(*icand,*jet);
342  float candDeta = std::abs( icand->eta() - jet->eta() );
343  float candDphi = reco::deltaPhi(*icand,*jet);
344  float candPtDr = candPt * candDr;
345  size_t icone = std::lower_bound(&cones[0],&cones[ncones],candDr) - &cones[0];
346 
347  if(candDr < dRmin) dRmin = candDr;
348 
349  // // all particles
350  if( lLead == nullptr || candPt > lLead->pt() ) {
351  lSecond = lLead;
352  lLead = icand;
353  } else if( (lSecond == nullptr || candPt > lSecond->pt()) && (candPt < lLead->pt()) ) {
354  lSecond = icand;
355  }
356 
357  // // average shapes
358  internalId_.dRMean_ += candPtDr;
359  internalId_.dR2Mean_ += candPtDr*candPtDr;
360  covMatrix(0,0) += candPt*candPt*candDeta*candDeta;
361  covMatrix(0,1) += candPt*candPt*candDeta*candDphi;
362  covMatrix(1,1) += candPt*candPt*candDphi*candDphi;
363  internalId_.ptD_ += candPt*candPt;
364  sumPt += candPt;
365  sumPt2 += candPt*candPt;
366 
367  // single most energetic candiates and jet shape profiles
368  frac.push_back(candPtFrac);
369 
370  if( icone < ncones ) { *coneFracs[icone] += candPt; }
371 
372  // neutrals
373  if( abs(icand->pdgId()) == 130) {
374  if (lLeadNeut == nullptr || candPt > lLeadNeut->pt()) { lLeadNeut = icand; }
375  internalId_.dRMeanNeut_ += candPtDr;
376  fracNeut.push_back(candPtFrac);
377  if( icone < ncones ) { *coneNeutFracs[icone] += candPt; }
378  internalId_.ptDNe_ += candPt*candPt;
379  sumPtNe += candPt;
380  }
381  // EM candidated
382  if( icand->pdgId() == 22 ) {
383  if(lLeadEm == nullptr || candPt > lLeadEm->pt()) { lLeadEm = icand; }
384  internalId_.dRMeanEm_ += candPtDr;
385  fracEm.push_back(candPtFrac);
386  if( icone < ncones ) { *coneEmFracs[icone] += candPt; }
387  internalId_.ptDNe_ += candPt*candPt;
388  sumPtNe += candPt;
389  }
390  // Charged particles
391  if( icand->charge() != 0 ) {
392  if (lLeadCh == nullptr || candPt > lLeadCh->pt()) {
393  lLeadCh = icand;
394 
395  const reco::Track* pfTrk = icand->bestTrack();
396  if (lPF && std::abs(icand->pdgId()) == 13 && pfTrk == nullptr){
397  reco::MuonRef lmuRef = lPF->muonRef();
398  if (lmuRef.isNonnull()){
399  const reco::Muon& lmu = *lmuRef.get();
400  pfTrk = lmu.bestTrack();
401  edm::LogWarning("BadMuon")<<"Found a PFCandidate muon without a trackRef: falling back to Muon::bestTrack ";
402  }
403  }
404  if(pfTrk==nullptr) { //protection against empty pointers for the miniAOD case
405  //To handle the electron case
406  if(lPF!=nullptr) {
407  pfTrk=(lPF->trackRef().get()==nullptr)?lPF->gsfTrackRef().get():lPF->trackRef().get();
408  }
409  const reco::Track& impactTrack = (lPack==nullptr)?(*pfTrk):(lPack->pseudoTrack());
410  internalId_.d0_ = std::abs(impactTrack.dxy(vtx->position()));
411  internalId_.dZ_ = std::abs(impactTrack.dz(vtx->position()));
412  }
413  else {
414  internalId_.d0_ = std::abs(pfTrk->dxy(vtx->position()));
415  internalId_.dZ_ = std::abs(pfTrk->dz(vtx->position()));
416  }
417  }
418  internalId_.dRMeanCh_ += candPtDr;
419  internalId_.ptDCh_ += candPt*candPt;
420  fracCh.push_back(candPtFrac);
421  if( icone < ncones ) { *coneChFracs[icone] += candPt; }
422  sumPtCh += candPt;
423  }
424  // // beta and betastar
425  if( icand->charge() != 0 ) {
426  if (!isPacked){
427  if(lPF->trackRef().isNonnull() ) {
428  float tkpt = candPt;
429  sumTkPt += tkpt;
430  // 'classic' beta definition based on track-vertex association
431  bool inVtx0 = vtx->trackWeight ( lPF->trackRef()) > 0 ;
432 
433  bool inAnyOther = false;
434  // alternative beta definition based on track-vertex distance of closest approach
435  double dZ0 = std::abs(lPF->trackRef()->dz(vtx->position()));
436  double dZ = dZ0;
437  for(reco::VertexCollection::const_iterator vi=allvtx.begin(); vi!=allvtx.end(); ++vi ) {
438  const reco::Vertex & iv = *vi;
439  if( iv.isFake() || iv.ndof() < 4 ) { continue; }
440  // the primary vertex may have been copied by the user: check identity by position
441  bool isVtx0 = (iv.position() - vtx->position()).r() < 0.02;
442  // 'classic' beta definition: check if the track is associated with
443  // any vertex other than the primary one
444  if( ! isVtx0 && ! inAnyOther ) {
445  inAnyOther = vtx->trackWeight ( lPF->trackRef()) <= 0 ;
446  }
447  // alternative beta: find closest vertex to the track
448  dZ = std::min(dZ,std::abs(lPF->trackRef()->dz(iv.position())));
449  }
450  // classic beta/betaStar
451  if( inVtx0 && ! inAnyOther ) {
452  internalId_.betaClassic_ += tkpt;
453  } else if( ! inVtx0 && inAnyOther ) {
454  internalId_.betaStarClassic_ += tkpt;
455  }
456  // alternative beta/betaStar
457  if( dZ0 < 0.2 ) {
458  internalId_.beta_ += tkpt;
459  } else if( dZ < 0.2 ) {
460  internalId_.betaStar_ += tkpt;
461  }
462  }
463  }
464  else{
465  // setting classic and alternative to be the same for now
466  float tkpt = candPt;
467  sumTkPt += tkpt;
468  bool inVtx0 = false;
469  bool inVtxOther = false;
470  if (lPack->fromPV() == pat::PackedCandidate::PVUsedInFit) inVtx0 = true;
471  if (lPack->fromPV() == 0) inVtxOther = true;
472  if (inVtx0){
473  internalId_.betaClassic_ += tkpt;
474  internalId_.beta_ += tkpt;
475  }
476  if (inVtxOther){
477  internalId_.betaStarClassic_ += tkpt;
478  internalId_.betaStar_ += tkpt;
479  }
480  }
481  }
482  // trailing candidate
483  if( lTrail == nullptr || candPt < lTrail->pt() ) {
484  lTrail = icand;
485  }
486 
487  }
488 
489  // // Finalize all variables
490  assert( !(lLead == nullptr) );
491 
492  if ( lSecond == nullptr ) { lSecond = lTrail; }
493  if ( lLeadNeut == nullptr ) { lLeadNeut = lTrail; }
494  if ( lLeadEm == nullptr ) { lLeadEm = lTrail; }
495  if ( lLeadCh == nullptr ) { lLeadCh = lTrail; }
496 
497  internalId_.nCharged_ = pfjet->chargedMultiplicity();
498  internalId_.nNeutrals_ = pfjet->neutralMultiplicity();
499  internalId_.chgEMfrac_ = pfjet->chargedEmEnergy() /jet->energy();
500  internalId_.neuEMfrac_ = pfjet->neutralEmEnergy() /jet->energy();
503  internalId_.nParticles_ = jet->numberOfDaughters();
504 
505 
507  float sumW(0.0), sumW2(0.0);
508  float sum_deta(0.0),sum_dphi(0.0);
509  float ave_deta(0.0), ave_dphi(0.0);
510  for (size_t j = 0; j < jet->numberOfDaughters(); j++) {
511  const auto& part = jet->daughterPtr(j);
512  if (!(part.isAvailable() && part.isNonnull()) ){
513  continue;
514  }
515 
516  float weight = part->pt();
517  float weight2 = weight * weight;
518  sumW2 += weight2;
519  float deta = part->eta() - jet->eta();
520  float dphi = reco::deltaPhi(*part, *jet);
521  sum_deta += deta*weight2;
522  sum_dphi += dphi*weight2;
523  if (sumW > 0) {
524  ave_deta = sum_deta/sumW2;
525  ave_dphi = sum_dphi/sumW2;
526  }
527  }
528 
529  float ddetaR_sum(0.0), ddphiR_sum(0.0), pull_tmp(0.0);
530  for (size_t i = 0; i < jet->numberOfDaughters(); i++) {
531  const auto& part = jet->daughterPtr(i);
532  if (!(part.isAvailable() && part.isNonnull()) ){
533  continue;
534  }
535  float weight =part->pt()*part->pt();
536  float deta = part->eta() - jet->eta();
537  float dphi = reco::deltaPhi(*part, *jet);
538  float ddeta, ddphi, ddR;
539  ddeta = deta - ave_deta ;
540  ddphi = dphi-ave_dphi;
541  ddR = sqrt(ddeta*ddeta + ddphi*ddphi);
542  ddetaR_sum += ddR*ddeta*weight;
543  ddphiR_sum += ddR*ddphi*weight;
544  }
545  if (sumW2 > 0) {
546  float ddetaR_ave = ddetaR_sum/sumW2;
547  float ddphiR_ave = ddphiR_sum/sumW2;
548  pull_tmp = sqrt(ddetaR_ave*ddetaR_ave+ddphiR_ave*ddphiR_ave);
549  }
550  internalId_.pull_ = pull_tmp;
552 
553 
559 
560  std::sort(frac.begin(),frac.end(),std::greater<float>());
561  std::sort(fracCh.begin(),fracCh.end(),std::greater<float>());
562  std::sort(fracEm.begin(),fracEm.end(),std::greater<float>());
563  std::sort(fracNeut.begin(),fracNeut.end(),std::greater<float>());
568 
569  covMatrix(0,0) /= sumPt2;
570  covMatrix(0,1) /= sumPt2;
571  covMatrix(1,1) /= sumPt2;
572  covMatrix(1,0) = covMatrix(0,1);
573  internalId_.etaW_ = sqrt(covMatrix(0,0));
574  internalId_.phiW_ = sqrt(covMatrix(1,1));
576  TVectorD eigVals(2); eigVals = TMatrixDSymEigen(covMatrix).GetEigenValues();
577  internalId_.majW_ = sqrt(std::abs(eigVals(0)));
578  internalId_.minW_ = sqrt(std::abs(eigVals(1)));
579  if( internalId_.majW_ < internalId_.minW_ ) { std::swap(internalId_.majW_,internalId_.minW_); }
580 
581  internalId_.dRLeadCent_ = reco::deltaR(*jet,*lLead);
582  if( lSecond == nullptr ) { internalId_.dRLead2nd_ = reco::deltaR(*jet,*lSecond); }
583  internalId_.dRMean_ /= jetPt;
587  internalId_.dR2Mean_ /= sumPt2;
588 
589  for(size_t ic=0; ic<ncones; ++ic){
590  *coneFracs[ic] /= jetPt;
591  *coneEmFracs[ic] /= jetPt;
592  *coneNeutFracs[ic] /= jetPt;
593  *coneChFracs[ic] /= jetPt;
594  }
595  //http://jets.physics.harvard.edu/qvg/
596  double ptMean = sumPt/internalId_.nParticles_;
597  double ptRMS = 0;
598  for(unsigned int i0 = 0; i0 < frac.size(); i0++) {ptRMS+=(frac[i0]-ptMean)*(frac[i0]-ptMean);}
599  ptRMS/=internalId_.nParticles_;
600  ptRMS=sqrt(ptRMS);
601 
602  internalId_.ptMean_ = ptMean;
603  internalId_.ptRMS_ = ptRMS/jetPt;
604  internalId_.pt2A_ = sqrt( internalId_.ptD_ /internalId_.nParticles_)/jetPt;
605  internalId_.ptD_ = sqrt( internalId_.ptD_) / sumPt;
606  internalId_.ptDCh_ = sqrt( internalId_.ptDCh_) / sumPtCh;
607  internalId_.ptDNe_ = sqrt( internalId_.ptDNe_) / sumPtNe;
608  internalId_.sumPt_ = sumPt;
609  internalId_.sumChPt_ = sumPtCh;
610  internalId_.sumNePt_ = sumPtNe;
611 
612  internalId_.jetR_ = lLead->pt()/sumPt;
613  internalId_.jetRchg_ = lLeadEm->pt()/sumPt;
614  internalId_.dRMatch_ = dRmin;
615 
616  if( sumTkPt != 0. ) {
617  internalId_.beta_ /= sumTkPt;
618  internalId_.betaStar_ /= sumTkPt;
619  internalId_.betaClassic_ /= sumTkPt;
620  internalId_.betaStarClassic_ /= sumTkPt;
621  } else {
622  assert( internalId_.beta_ == 0. && internalId_.betaStar_ == 0.&& internalId_.betaClassic_ == 0. && internalId_.betaStarClassic_ == 0. );
623  }
624 
625  if( runMvas_ ) {
626  runMva();
627  }
628 
630 }
631 
632 
633 
634 // ------------------------------------------------------------------------------------------
636 {
637  std::stringstream out;
638  for(variables_list_t::const_iterator it=variables_.begin();
639  it!=variables_.end(); ++it ) {
640  out << std::setw(15) << it->first << std::setw(3) << "="
641  << std::setw(5) << *it->second.first
642  << " (" << std::setw(5) << it->second.second << ")" << std::endl;
643  }
644  return out.str();
645 }
646 
647 // ------------------------------------------------------------------------------------------
649 {
650  internalId_.idFlag_ = 0;
651  for(variables_list_t::iterator it=variables_.begin();
652  it!=variables_.end(); ++it ) {
653  *it->second.first = it->second.second;
654  }
655 }
656 
657 // ------------------------------------------------------------------------------------------
658 #define INIT_VARIABLE(NAME,TMVANAME,VAL) \
659  internalId_.NAME ## _ = VAL; \
660  variables_[ # NAME ] = std::make_pair(& internalId_.NAME ## _, VAL);
661 
662 // ------------------------------------------------------------------------------------------
664 {
665  internalId_.idFlag_ = 0;
666  INIT_VARIABLE(mva , "", -100.);
667  //INIT_VARIABLE(jetPt , "jspt_1", 0.);
668  //INIT_VARIABLE(jetEta , "jseta_1", large_val);
669  INIT_VARIABLE(jetPt , "", 0.);
671  INIT_VARIABLE(jetPhi , "jsphi_1", large_val);
672  INIT_VARIABLE(jetM , "jm_1", 0.);
673 
674  INIT_VARIABLE(nCharged , "", 0.);
675  INIT_VARIABLE(nNeutrals , "", 0.);
676 
677  INIT_VARIABLE(chgEMfrac , "", 0.);
678  INIT_VARIABLE(neuEMfrac , "", 0.);
679  INIT_VARIABLE(chgHadrfrac, "", 0.);
680  INIT_VARIABLE(neuHadrfrac, "", 0.);
681 
682  INIT_VARIABLE(d0 , "jd0_1" , -1000.);
683  INIT_VARIABLE(dZ , "jdZ_1" , -1000.);
684  //INIT_VARIABLE(nParticles , "npart_1" , 0.);
685  INIT_VARIABLE(nParticles , "" , 0.);
686 
687  INIT_VARIABLE(leadPt , "lpt_1" , 0.);
688  INIT_VARIABLE(leadEta , "leta_1" , large_val);
689  INIT_VARIABLE(leadPhi , "lphi_1" , large_val);
690  INIT_VARIABLE(secondPt , "spt_1" , 0.);
691  INIT_VARIABLE(secondEta , "seta_1" , large_val);
692  INIT_VARIABLE(secondPhi , "sphi_1" , large_val);
693  INIT_VARIABLE(leadNeutPt , "lnept_1" , 0.);
694  INIT_VARIABLE(leadNeutEta, "lneeta_1" , large_val);
695  INIT_VARIABLE(leadNeutPhi, "lnephi_1" , large_val);
696  INIT_VARIABLE(leadEmPt , "lempt_1" , 0.);
697  INIT_VARIABLE(leadEmEta , "lemeta_1" , large_val);
698  INIT_VARIABLE(leadEmPhi , "lemphi_1" , large_val);
699  INIT_VARIABLE(leadChPt , "lchpt_1" , 0.);
700  INIT_VARIABLE(leadChEta , "lcheta_1" , large_val);
701  INIT_VARIABLE(leadChPhi , "lchphi_1" , large_val);
702  INIT_VARIABLE(leadFrac , "lLfr_1" , 0.);
703 
704  INIT_VARIABLE(dRLeadCent , "drlc_1" , 0.);
705  INIT_VARIABLE(dRLead2nd , "drls_1" , 0.);
706  INIT_VARIABLE(dRMean , "drm_1" , 0.);
707  INIT_VARIABLE(dRMean , "" , 0.);
708  INIT_VARIABLE(pull , "" , 0.);
709  INIT_VARIABLE(dRMeanNeut , "drmne_1" , 0.);
710  INIT_VARIABLE(dRMeanEm , "drem_1" , 0.);
711  INIT_VARIABLE(dRMeanCh , "drch_1" , 0.);
712  INIT_VARIABLE(dR2Mean , "" , 0.);
713 
714  INIT_VARIABLE(ptD , "", 0.);
715  INIT_VARIABLE(ptMean , "", 0.);
716  INIT_VARIABLE(ptRMS , "", 0.);
717  INIT_VARIABLE(pt2A , "", 0.);
718  INIT_VARIABLE(ptDCh , "", 0.);
719  INIT_VARIABLE(ptDNe , "", 0.);
720  INIT_VARIABLE(sumPt , "", 0.);
721  INIT_VARIABLE(sumChPt , "", 0.);
722  INIT_VARIABLE(sumNePt , "", 0.);
723 
724  INIT_VARIABLE(secondFrac ,"" ,0.);
725  INIT_VARIABLE(thirdFrac ,"" ,0.);
726  INIT_VARIABLE(fourthFrac ,"" ,0.);
727 
728  INIT_VARIABLE(leadChFrac ,"" ,0.);
729  INIT_VARIABLE(secondChFrac ,"" ,0.);
730  INIT_VARIABLE(thirdChFrac ,"" ,0.);
731  INIT_VARIABLE(fourthChFrac ,"" ,0.);
732 
733  INIT_VARIABLE(leadNeutFrac ,"" ,0.);
734  INIT_VARIABLE(secondNeutFrac ,"" ,0.);
735  INIT_VARIABLE(thirdNeutFrac ,"" ,0.);
736  INIT_VARIABLE(fourthNeutFrac ,"" ,0.);
737 
738  INIT_VARIABLE(leadEmFrac ,"" ,0.);
739  INIT_VARIABLE(secondEmFrac ,"" ,0.);
740  INIT_VARIABLE(thirdEmFrac ,"" ,0.);
741  INIT_VARIABLE(fourthEmFrac ,"" ,0.);
742 
743  INIT_VARIABLE(jetW ,"" ,1.);
744  INIT_VARIABLE(etaW ,"" ,1.);
745  INIT_VARIABLE(phiW ,"" ,1.);
746 
747  INIT_VARIABLE(majW ,"" ,1.);
748  INIT_VARIABLE(minW ,"" ,1.);
749 
750  INIT_VARIABLE(frac01 ,"" ,0.);
751  INIT_VARIABLE(frac02 ,"" ,0.);
752  INIT_VARIABLE(frac03 ,"" ,0.);
753  INIT_VARIABLE(frac04 ,"" ,0.);
754  INIT_VARIABLE(frac05 ,"" ,0.);
755  INIT_VARIABLE(frac06 ,"" ,0.);
756  INIT_VARIABLE(frac07 ,"" ,0.);
757 
758  INIT_VARIABLE(chFrac01 ,"" ,0.);
759  INIT_VARIABLE(chFrac02 ,"" ,0.);
760  INIT_VARIABLE(chFrac03 ,"" ,0.);
761  INIT_VARIABLE(chFrac04 ,"" ,0.);
762  INIT_VARIABLE(chFrac05 ,"" ,0.);
763  INIT_VARIABLE(chFrac06 ,"" ,0.);
764  INIT_VARIABLE(chFrac07 ,"" ,0.);
765 
766  INIT_VARIABLE(neutFrac01 ,"" ,0.);
767  INIT_VARIABLE(neutFrac02 ,"" ,0.);
768  INIT_VARIABLE(neutFrac03 ,"" ,0.);
769  INIT_VARIABLE(neutFrac04 ,"" ,0.);
770  INIT_VARIABLE(neutFrac05 ,"" ,0.);
771  INIT_VARIABLE(neutFrac06 ,"" ,0.);
772  INIT_VARIABLE(neutFrac07 ,"" ,0.);
773 
774  INIT_VARIABLE(emFrac01 ,"" ,0.);
775  INIT_VARIABLE(emFrac02 ,"" ,0.);
776  INIT_VARIABLE(emFrac03 ,"" ,0.);
777  INIT_VARIABLE(emFrac04 ,"" ,0.);
778  INIT_VARIABLE(emFrac05 ,"" ,0.);
779  INIT_VARIABLE(emFrac06 ,"" ,0.);
780  INIT_VARIABLE(emFrac07 ,"" ,0.);
781 
782  INIT_VARIABLE(beta ,"" ,0.);
783  INIT_VARIABLE(betaStar ,"" ,0.);
784  INIT_VARIABLE(betaClassic ,"" ,0.);
785  INIT_VARIABLE(betaStarClassic ,"" ,0.);
786 
787  INIT_VARIABLE(nvtx ,"" ,0.);
788  INIT_VARIABLE(rho ,"" ,0.);
789  INIT_VARIABLE(nTrueInt ,"" ,0.);
790 
791  INIT_VARIABLE(jetR , "", 0.);
792  INIT_VARIABLE(jetRchg , "", 0.);
793  INIT_VARIABLE(dRMatch , "", 0.);
794 
795 }
796 
797 #undef INIT_VARIABLE
const double beta
T getParameter(std::string const &) const
void set(const PileupJetIdentifier &)
int i
Definition: DBlmapReader.cc:9
std::string tmvaWeights_jteta_2p5_3_
Float_t impactParTkThreshod_
std::vector< std::string > tmvaVariables_
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:252
#define INIT_VARIABLE(NAME, TMVANAME, VAL)
Float_t mvacut_[3][4][4]
std::string tmvaWeights_jteta_2_2p5_
float chargedEmEnergy() const
chargedEmEnergy
Definition: PFJet.h:142
virtual double energy() const final
energy
CandidatePtr daughterPtr(size_type i) const
reference to daughter at given position
virtual const Track * bestTrack() const
Definition: Candidate.h:247
virtual double pt() const =0
transverse momentum
T const * get() const
Returns C++ pointer to the item.
Definition: Ptr.h:160
Base class for all types of Jets.
Definition: Jet.h:20
std::unique_ptr< TMVA::Reader > reader_jteta_2p5_3_
assert(m_qm.get())
std::pair< int, int > getJetIdKey(float jetPt, float jetEta)
std::map< std::string, std::string > tmvaNames_
std::string tmvaMethod_
std::unique_ptr< TMVA::Reader > reader_jteta_0_2_
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
int computeCutIDflag(float betaStarClassic, float dR2Mean, float nvtx, float jetPt, float jetEta)
Float_t betaStarCut_[3][4][4]
int chargedMultiplicity() const
chargedMultiplicity
Definition: PFJet.h:155
const Point & position() const
position
Definition: Vertex.h:99
std::unique_ptr< TMVA::Reader > reader_jteta_3_5_
Jets made from PFObjects.
Definition: PFJet.h:21
float neutralEmEnergy() const
neutralEmEnergy
Definition: PFJet.h:150
std::vector< std::string > tmvaVariables_jteta_0_3_
reco::TrackRef trackRef() const
Definition: PFCandidate.cc:433
tuple d
Definition: ztail.py:151
virtual const Track * bestTrack() const
best track pointer
Definition: Muon.h:61
virtual size_t numberOfDaughters() const
number of daughters
virtual CandidatePtr sourceCandidatePtr(size_type i) const
PileupJetIdentifier internalId_
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
const PVAssoc fromPV(size_t ipv=0) const
T sqrt(T t)
Definition: SSEVec.h:18
std::string dumpVariables() const
std::vector< std::string > tmvaVariables_jteta_3_5_
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
virtual const reco::Track & pseudoTrack() const
Return reference to a pseudo track made with candidate kinematics, parameterized error for eta...
auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:28
void assign(const std::vector< float > &vec, float &a, float &b, float &c, float &d)
float trackWeight(const TREF &r) const
returns the weight with which a Track has contributed to the vertex-fit.
Definition: Vertex.h:73
virtual int charge() const =0
electric charge
virtual size_type numberOfSourceCandidatePtrs() const
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:244
T min(T a, T b)
Definition: MathUtil.h:58
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 >(), bool runMvas=true)
virtual double mass() const final
mass
int computeIDflag(float mva, float jetPt, float jetEta)
double ndof() const
Definition: Vertex.h:95
int neutralMultiplicity() const
neutralMultiplicity
Definition: PFJet.h:157
reco::MuonRef muonRef() const
Definition: PFCandidate.cc:450
double dz() const
dz parameter (= dsz/cos(lambda)). This is the track z0 w.r.t (0,0,0) only if the refPoint is close to...
Definition: TrackBase.h:604
virtual int pdgId() const =0
PDG identifier.
void setPtEtaPhi(const reco::Candidate &p, float &pt, float &eta, float &phi)
double deltaPhi(double phi1, double phi2)
Definition: deltaPhi.h:22
std::string tmvaWeights_jteta_3_5_
bool isFake() const
Definition: Vertex.h:64
PileupJetIdentifier computeIdVariables(const reco::Jet *jet, float jec, const reco::Vertex *, const reco::VertexCollection &, double rho)
unsigned int nCharged(const GenJet &jet)
Float_t rmsCut_[3][4][4]
part
Definition: HCALResponse.h:20
const float large_val
Definition: MVAJetPuId.cc:14
PileupJetIdentifier computeMva()
double b
Definition: hdecay.h:120
Analysis-level calorimeter jet class.
Definition: Jet.h:77
std::vector< std::string > tmvaSpectators_
TMVA::IMethod * loadTMVAWeights(TMVA::Reader *reader, const std::string &method, const std::string &weightFile, bool verbose=false)
Geom::Phi< T > phi() const
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:39
std::unique_ptr< TMVA::Reader > reader_
std::string tmvaWeights_
double a
Definition: hdecay.h:121
reco::GsfTrackRef gsfTrackRef() const
Definition: PFCandidate.cc:471
variables_list_t variables_
float neutralHadronEnergy() const
neutralHadronEnergy
Definition: PFJet.h:102
std::unique_ptr< TMVA::Reader > reader_jteta_2_2p5_
virtual double eta() const final
momentum pseudorapidity
std::string fullPath() const
Definition: FileInPath.cc:184
int weight
Definition: histoStyle.py:50
double dxy() const
dxy parameter. (This is the transverse impact parameter w.r.t. to (0,0,0) ONLY if refPoint is close t...
Definition: TrackBase.h:586
Jet correctedJet(const std::string &level, const std::string &flavor="none", const std::string &set="") const
float chargedHadronEnergy() const
chargedHadronEnergy
Definition: PFJet.h:98
std::string tmvaWeights_jteta_0_2_
virtual double phi() const =0
momentum azimuthal angle
virtual double eta() const =0
momentum pseudorapidity
virtual double pt() const final
transverse momentum