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(etaBinnedWeights_){
114  tmvaNames_["DRweighted"] = "dR2Mean";
115  tmvaNames_["rho"] = "nvtx";
116  tmvaNames_["nTot"] = "nParticles";
117  tmvaNames_["nCh"] = "nCharged";
118  tmvaNames_["p4.fCoordinates.fPt"] ="jetPt";
119  tmvaNames_["p4.fCoordinates.fEta"] ="jetEta";
120  tmvaNames_["axisMajor"] = "majW";
121  tmvaNames_["axisMinor"] = "minW";
122  tmvaNames_["fRing0"] = "frac01";
123  tmvaNames_["fRing1"] = "frac02";
124  tmvaNames_["fRing2"] = "frac03";
125  tmvaNames_["fRing3"] = "frac04";
126  tmvaNames_["min(pull,0.1)"] = "dRMean";
127  tmvaNames_["dRMatch"] = "chFrac02";
128  }
129 
130 
131  if( version_ == PHILv0 ) {
132  tmvaVariables_.clear();
133  tmvaVariables_.push_back( "jspt_1" );
134  tmvaVariables_.push_back( "jseta_1" );
135  tmvaVariables_.push_back( "jsphi_1" );
136  tmvaVariables_.push_back( "jd0_1" );
137  tmvaVariables_.push_back( "jdZ_1" );
138  tmvaVariables_.push_back( "jm_1" );
139  tmvaVariables_.push_back( "npart_1" );
140  tmvaVariables_.push_back( "lpt_1" );
141  tmvaVariables_.push_back( "leta_1" );
142  tmvaVariables_.push_back( "lphi_1" );
143  tmvaVariables_.push_back( "spt_1" );
144  tmvaVariables_.push_back( "seta_1" );
145  tmvaVariables_.push_back( "sphi_1" );
146  tmvaVariables_.push_back( "lnept_1" );
147  tmvaVariables_.push_back( "lneeta_1");
148  tmvaVariables_.push_back( "lnephi_1");
149  tmvaVariables_.push_back( "lempt_1" );
150  tmvaVariables_.push_back( "lemeta_1");
151  tmvaVariables_.push_back( "lemphi_1");
152  tmvaVariables_.push_back( "lchpt_1" );
153  // tmvaVariables_.push_back( "lcheta_1"); FIXME missing in weights file
154  tmvaVariables_.push_back( "lchphi_1");
155  tmvaVariables_.push_back( "lLfr_1" );
156  tmvaVariables_.push_back( "drlc_1" );
157  tmvaVariables_.push_back( "drls_1" );
158  tmvaVariables_.push_back( "drm_1" );
159  tmvaVariables_.push_back( "drmne_1" );
160  tmvaVariables_.push_back( "drem_1" );
161  tmvaVariables_.push_back( "drch_1" );
162 
163  tmvaNames_["jspt_1"] = "jetPt";
164  tmvaNames_["jseta_1"] = "jetEta";
165  tmvaNames_["jsphi_1"] = "jetPhi";
166  tmvaNames_["jm_1"] = "jetM";
167  tmvaNames_["jd0_1"] = "d0";
168  tmvaNames_["jdZ_1"] = "dZ";
169  tmvaNames_["npart_1"] = "nParticles";
170 
171  tmvaNames_["lpt_1"] = "leadPt";
172  tmvaNames_["leta_1"] = "leadEta";
173  tmvaNames_["lphi_1"] = "leadPhi";
174  tmvaNames_["spt_1"] = "secondPt";
175  tmvaNames_["seta_1"] = "secondEta";
176  tmvaNames_["sphi_1"] = "secondPhi";
177  tmvaNames_["lnept_1"] = "leadNeutPt";
178  tmvaNames_["lneeta_1"] = "leadNeutEta";
179  tmvaNames_["lnephi_1"] = "leadNeutPhi";
180  tmvaNames_["lempt_1"] = "leadEmPt";
181  tmvaNames_["lemeta_1"] = "leadEmEta";
182  tmvaNames_["lemphi_1"] = "leadEmPhi";
183  tmvaNames_["lchpt_1"] = "leadChPt";
184  tmvaNames_["lcheta_1"] = "leadChEta";
185  tmvaNames_["lchphi_1"] = "leadChPhi";
186  tmvaNames_["lLfr_1"] = "leadFrac";
187 
188  tmvaNames_["drlc_1"] = "dRLeadCent";
189  tmvaNames_["drls_1"] = "dRLead2nd";
190  tmvaNames_["drm_1"] = "dRMean";
191  tmvaNames_["drmne_1"] = "dRMeanNeut";
192  tmvaNames_["drem_1"] = "dRMeanEm";
193  tmvaNames_["drch_1"] = "dRMeanCh";
194 
195  } else if( ! cutBased_ ){
196  assert( tmvaMethod_.empty() || ((! tmvaVariables_.empty() || ( !tmvaVariables_jteta_0_3_.empty() && !tmvaVariables_jteta_3_5_.empty() )) && version_ == USER) );
197  }
198  if(( ! cutBased_ ) && (runMvas_)) { bookReader();}
199 }
200 
201 // ------------------------------------------------------------------------------------------
203 {
204 }
205 
206 // ------------------------------------------------------------------------------------------
207 void assign(const std::vector<float> & vec, float & a, float & b, float & c, float & d )
208 {
209  size_t sz = vec.size();
210  a = ( sz > 0 ? vec[0] : 0. );
211  b = ( sz > 1 ? vec[1] : 0. );
212  c = ( sz > 2 ? vec[2] : 0. );
213  d = ( sz > 3 ? vec[3] : 0. );
214 }
215 // ------------------------------------------------------------------------------------------
216 void setPtEtaPhi(const reco::Candidate & p, float & pt, float & eta, float &phi )
217 {
218  pt = p.pt();
219  eta = p.eta();
220  phi = p.phi();
221 }
222 
223 // ------------------------------------------------------------------------------------------
225 {
226  if(etaBinnedWeights_){
227  reader_jteta_0_2_ = std::unique_ptr<TMVA::Reader>(new TMVA::Reader("!Color:Silent"));
228  reader_jteta_2_2p5_ = std::unique_ptr<TMVA::Reader>(new TMVA::Reader("!Color:Silent"));
229  reader_jteta_2p5_3_ = std::unique_ptr<TMVA::Reader>(new TMVA::Reader("!Color:Silent"));
230  reader_jteta_3_5_ = std::unique_ptr<TMVA::Reader>(new TMVA::Reader("!Color:Silent"));
231  } else {
232  reader_ = std::unique_ptr<TMVA::Reader>(new TMVA::Reader("!Color:Silent"));
233  }
234  if(etaBinnedWeights_){
235  for(std::vector<std::string>::iterator it=tmvaVariables_jteta_0_3_.begin(); it!=tmvaVariables_jteta_0_3_.end(); ++it) {
236  if( tmvaNames_[*it].empty() ) {
237  tmvaNames_[*it] = *it;
238  }
239  reader_jteta_0_2_->AddVariable( *it, variables_[ tmvaNames_[*it] ].first );
240  reader_jteta_2_2p5_->AddVariable( *it, variables_[ tmvaNames_[*it] ].first );
241  reader_jteta_2p5_3_->AddVariable( *it, variables_[ tmvaNames_[*it] ].first );
242  }
243  for(std::vector<std::string>::iterator it=tmvaVariables_jteta_3_5_.begin(); it!=tmvaVariables_jteta_3_5_.end(); ++it) {
244  if( tmvaNames_[*it].empty() ) {
245  tmvaNames_[*it] = *it;
246  }
247  reader_jteta_3_5_->AddVariable( *it, variables_[ tmvaNames_[*it] ].first );
248  }
249  } else {
250  for(std::vector<std::string>::iterator it=tmvaVariables_.begin(); it!=tmvaVariables_.end(); ++it) {
251  if( tmvaNames_[*it].empty() ) {
252  tmvaNames_[*it] = *it;
253  }
254  reader_->AddVariable( *it, variables_[ tmvaNames_[*it] ].first );
255  }
256  }
257  for(std::vector<std::string>::iterator it=tmvaSpectators_.begin(); it!=tmvaSpectators_.end(); ++it) {
258  if( tmvaNames_[*it].empty() ) {
259  tmvaNames_[*it] = *it;
260  }
261  if(etaBinnedWeights_){
262  reader_jteta_0_2_->AddSpectator( *it, variables_[ tmvaNames_[*it] ].first );
263  reader_jteta_2_2p5_->AddSpectator( *it, variables_[ tmvaNames_[*it] ].first );
264  reader_jteta_2p5_3_->AddSpectator( *it, variables_[ tmvaNames_[*it] ].first );
265  reader_jteta_3_5_->AddSpectator( *it, variables_[ tmvaNames_[*it] ].first );
266  } else {
267  reader_->AddSpectator( *it, variables_[ tmvaNames_[*it] ].first );
268  }
269  }
270  if(etaBinnedWeights_){
275  } else {
277  }
278 }
279 
280 // ------------------------------------------------------------------------------------------
282 {
283  internalId_ = id;
284 }
285 
286 // ------------------------------------------------------------------------------------------
288 {
289  if( cutBased_ ) {
290  internalId_.idFlag_ = computeCutIDflag(internalId_.betaStarClassic_,internalId_.dR2Mean_,internalId_.nvtx_,internalId_.jetPt_,internalId_.jetEta_);
291  } else {
292  if(std::abs(internalId_.jetEta_) >= 5.0) {
293  internalId_.mva_ = -2.;
294  } else {
295  if(etaBinnedWeights_){
296  if(std::abs(internalId_.jetEta_)<=2.) internalId_.mva_ = reader_jteta_0_2_->EvaluateMVA( tmvaMethod_.c_str() );
297  else if(std::abs(internalId_.jetEta_)<=2.5) internalId_.mva_ = reader_jteta_2_2p5_->EvaluateMVA( tmvaMethod_.c_str() );
298  else if(std::abs(internalId_.jetEta_)<=3.) internalId_.mva_ = reader_jteta_2p5_3_->EvaluateMVA( tmvaMethod_.c_str() );
299  else internalId_.mva_ = reader_jteta_3_5_->EvaluateMVA( tmvaMethod_.c_str() );
300  } else {
301  internalId_.mva_ = reader_->EvaluateMVA( tmvaMethod_.c_str() );
302  }
303  }
305  }
306 }
307 
308 // ------------------------------------------------------------------------------------------
309 std::pair<int,int> PileupJetIdAlgo::getJetIdKey(float jetPt, float jetEta)
310 {
311  int ptId = 0;
312  if(jetPt >= 10 && jetPt < 20) ptId = 1;
313  if(jetPt >= 20 && jetPt < 30) ptId = 2;
314  if(jetPt >= 30 ) ptId = 3;
315 
316  int etaId = 0;
317  if(std::abs(jetEta) >= 2.5 && std::abs(jetEta) < 2.75) etaId = 1;
318  if(std::abs(jetEta) >= 2.75 && std::abs(jetEta) < 3.0 ) etaId = 2;
319  if(std::abs(jetEta) >= 3.0 && std::abs(jetEta) < 5.0 ) etaId = 3;
320 
321  return std::pair<int,int>(ptId,etaId);
322 }
323 // ------------------------------------------------------------------------------------------
324 int PileupJetIdAlgo::computeCutIDflag(float betaStarClassic,float dR2Mean,float nvtx, float jetPt, float jetEta)
325 {
326  std::pair<int,int> jetIdKey = getJetIdKey(jetPt,jetEta);
327  float betaStarModified = betaStarClassic/log(nvtx-0.64);
328  int idFlag(0);
329  if(betaStarModified < betaStarCut_[PileupJetIdentifier::kTight ][jetIdKey.first][jetIdKey.second] &&
330  dR2Mean < rmsCut_ [PileupJetIdentifier::kTight ][jetIdKey.first][jetIdKey.second]
331  ) idFlag += 1 << PileupJetIdentifier::kTight;
332 
333  if(betaStarModified < betaStarCut_[PileupJetIdentifier::kMedium ][jetIdKey.first][jetIdKey.second] &&
334  dR2Mean < rmsCut_ [PileupJetIdentifier::kMedium ][jetIdKey.first][jetIdKey.second]
335  ) idFlag += 1 << PileupJetIdentifier::kMedium;
336 
337  if(betaStarModified < betaStarCut_[PileupJetIdentifier::kLoose ][jetIdKey.first][jetIdKey.second] &&
338  dR2Mean < rmsCut_ [PileupJetIdentifier::kLoose ][jetIdKey.first][jetIdKey.second]
339  ) idFlag += 1 << PileupJetIdentifier::kLoose;
340  return idFlag;
341 }
342 // ------------------------------------------------------------------------------------------
343 int PileupJetIdAlgo::computeIDflag(float mva, float jetPt, float jetEta)
344 {
345  std::pair<int,int> jetIdKey = getJetIdKey(jetPt,jetEta);
346  return computeIDflag(mva,jetIdKey.first,jetIdKey.second);
347 }
348 
349 // ------------------------------------------------------------------------------------------
350 int PileupJetIdAlgo::computeIDflag(float mva,int ptId,int etaId)
351 {
352  int idFlag(0);
353  if(mva > mvacut_[PileupJetIdentifier::kTight ][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kTight;
354  if(mva > mvacut_[PileupJetIdentifier::kMedium][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kMedium;
355  if(mva > mvacut_[PileupJetIdentifier::kLoose ][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kLoose;
356  return idFlag;
357 }
358 
359 
360 // ------------------------------------------------------------------------------------------
362 {
363  runMva();
365 }
366 
367 // ------------------------------------------------------------------------------------------
369  const reco::VertexCollection & allvtx, double rho)
370 {
371 
372  static std::atomic<int> printWarning{10};
373 
374  // initialize all variables to 0
375  resetVariables();
376 
377  // loop over constituents, accumulate sums and find leading candidates
378  const pat::Jet * patjet = dynamic_cast<const pat::Jet *>(jet);
379  const reco::PFJet * pfjet = dynamic_cast<const reco::PFJet *>(jet);
380  assert( patjet != nullptr || pfjet != nullptr );
381  if( patjet != nullptr && jec == 0. ) { // if this is a pat jet and no jec has been passed take the jec from the object
382  jec = patjet->pt()/patjet->correctedJet(0).pt();
383  }
384  if( jec <= 0. ) {
385  jec = 1.;
386  }
387 
388  const reco::Candidate* lLead = nullptr, *lSecond = nullptr, *lLeadNeut = nullptr, *lLeadEm = nullptr, *lLeadCh = nullptr, *lTrail = nullptr;
389  std::vector<float> frac, fracCh, fracEm, fracNeut;
390  float cones[] = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7 };
391  size_t ncones = sizeof(cones)/sizeof(float);
392  float * coneFracs[] = { &internalId_.frac01_, &internalId_.frac02_, &internalId_.frac03_, &internalId_.frac04_,
393  &internalId_.frac05_, &internalId_.frac06_, &internalId_.frac07_ };
400  TMatrixDSym covMatrix(2); covMatrix = 0.;
401  float jetPt = jet->pt() / jec; // use uncorrected pt for shape variables
402  float sumPt = 0., sumPt2 = 0., sumTkPt = 0.,sumPtCh=0,sumPtNe = 0;
403  setPtEtaPhi(*jet,internalId_.jetPt_,internalId_.jetEta_,internalId_.jetPhi_); // use corrected pt for jet kinematics
404  internalId_.jetM_ = jet->mass();
405  internalId_.nvtx_ = allvtx.size();
406  internalId_.rho_ = rho;
407 
408  float dRmin(1000);
409 
410  for ( unsigned i = 0; i < jet->numberOfSourceCandidatePtrs(); ++i ) {
411  reco::CandidatePtr pfJetConstituent = jet->sourceCandidatePtr(i);
412 
413  const reco::Candidate* icand = pfJetConstituent.get();
414  const pat::PackedCandidate* lPack = dynamic_cast<const pat::PackedCandidate *>( icand );
415  const reco::PFCandidate *lPF=dynamic_cast<const reco::PFCandidate*>( icand );
416  bool isPacked = true;
417  if (lPack == nullptr){
418  isPacked = false;
419  }
420 
421  float candPt = icand->pt();
422  float candPtFrac = candPt/jetPt;
423  float candDr = reco::deltaR(*icand,*jet);
424  float candDeta = std::abs( icand->eta() - jet->eta() );
425  float candDphi = reco::deltaPhi(*icand,*jet);
426  float candPtDr = candPt * candDr;
427  size_t icone = std::lower_bound(&cones[0],&cones[ncones],candDr) - &cones[0];
428 
429  if(candDr < dRmin) dRmin = candDr;
430 
431  // // all particles
432  if( lLead == nullptr || candPt > lLead->pt() ) {
433  lSecond = lLead;
434  lLead = icand;
435  } else if( (lSecond == nullptr || candPt > lSecond->pt()) && (candPt < lLead->pt()) ) {
436  lSecond = icand;
437  }
438 
439  // // average shapes
440  internalId_.dRMean_ += candPtDr;
441  internalId_.dR2Mean_ += candPtDr*candPtDr;
442  covMatrix(0,0) += candPt*candPt*candDeta*candDeta;
443  covMatrix(0,1) += candPt*candPt*candDeta*candDphi;
444  covMatrix(1,1) += candPt*candPt*candDphi*candDphi;
445  internalId_.ptD_ += candPt*candPt;
446  sumPt += candPt;
447  sumPt2 += candPt*candPt;
448 
449  // single most energetic candiates and jet shape profiles
450  frac.push_back(candPtFrac);
451 
452  if( icone < ncones ) { *coneFracs[icone] += candPt; }
453 
454  // neutrals
455  if( abs(icand->pdgId()) == 130) {
456  if (lLeadNeut == nullptr || candPt > lLeadNeut->pt()) { lLeadNeut = icand; }
457  internalId_.dRMeanNeut_ += candPtDr;
458  fracNeut.push_back(candPtFrac);
459  if( icone < ncones ) { *coneNeutFracs[icone] += candPt; }
460  internalId_.ptDNe_ += candPt*candPt;
461  sumPtNe += candPt;
462  }
463  // EM candidated
464  if( icand->pdgId() == 22 ) {
465  if(lLeadEm == nullptr || candPt > lLeadEm->pt()) { lLeadEm = icand; }
466  internalId_.dRMeanEm_ += candPtDr;
467  fracEm.push_back(candPtFrac);
468  if( icone < ncones ) { *coneEmFracs[icone] += candPt; }
469  internalId_.ptDNe_ += candPt*candPt;
470  sumPtNe += candPt;
471  }
472  // Charged particles
473  if( icand->charge() != 0 ) {
474  if (lLeadCh == nullptr || candPt > lLeadCh->pt()) {
475  lLeadCh = icand;
476 
477  const reco::Track* pfTrk = icand->bestTrack();
478  if (lPF && std::abs(icand->pdgId()) == 13 && pfTrk == nullptr){
479  reco::MuonRef lmuRef = lPF->muonRef();
480  if (lmuRef.isNonnull()){
481  const reco::Muon& lmu = *lmuRef.get();
482  pfTrk = lmu.bestTrack();
483  edm::LogWarning("BadMuon")<<"Found a PFCandidate muon without a trackRef: falling back to Muon::bestTrack ";
484  }
485  }
486  if(pfTrk==nullptr) { //protection against empty pointers for the miniAOD case
487  //To handle the electron case
488  if(lPF!=nullptr) {
489  pfTrk=(lPF->trackRef().get()==nullptr)?lPF->gsfTrackRef().get():lPF->trackRef().get();
490  }
491  const reco::Track& impactTrack = (lPack==nullptr)?(*pfTrk):(lPack->pseudoTrack());
492  internalId_.d0_ = std::abs(impactTrack.dxy(vtx->position()));
493  internalId_.dZ_ = std::abs(impactTrack.dz(vtx->position()));
494  }
495  else {
496  internalId_.d0_ = std::abs(pfTrk->dxy(vtx->position()));
497  internalId_.dZ_ = std::abs(pfTrk->dz(vtx->position()));
498  }
499  }
500  internalId_.dRMeanCh_ += candPtDr;
501  internalId_.ptDCh_ += candPt*candPt;
502  fracCh.push_back(candPtFrac);
503  if( icone < ncones ) { *coneChFracs[icone] += candPt; }
504  sumPtCh += candPt;
505  }
506  // // beta and betastar
507  if( icand->charge() != 0 ) {
508  if (!isPacked){
509  if(lPF->trackRef().isNonnull() ) {
510  float tkpt = candPt;
511  sumTkPt += tkpt;
512  // 'classic' beta definition based on track-vertex association
513  bool inVtx0 = vtx->trackWeight ( lPF->trackRef()) > 0 ;
514 
515  bool inAnyOther = false;
516  // alternative beta definition based on track-vertex distance of closest approach
517  double dZ0 = std::abs(lPF->trackRef()->dz(vtx->position()));
518  double dZ = dZ0;
519  for(reco::VertexCollection::const_iterator vi=allvtx.begin(); vi!=allvtx.end(); ++vi ) {
520  const reco::Vertex & iv = *vi;
521  if( iv.isFake() || iv.ndof() < 4 ) { continue; }
522  // the primary vertex may have been copied by the user: check identity by position
523  bool isVtx0 = (iv.position() - vtx->position()).r() < 0.02;
524  // 'classic' beta definition: check if the track is associated with
525  // any vertex other than the primary one
526  if( ! isVtx0 && ! inAnyOther ) {
527  inAnyOther = vtx->trackWeight ( lPF->trackRef()) <= 0 ;
528  }
529  // alternative beta: find closest vertex to the track
530  dZ = std::min(dZ,std::abs(lPF->trackRef()->dz(iv.position())));
531  }
532  // classic beta/betaStar
533  if( inVtx0 && ! inAnyOther ) {
534  internalId_.betaClassic_ += tkpt;
535  } else if( ! inVtx0 && inAnyOther ) {
536  internalId_.betaStarClassic_ += tkpt;
537  }
538  // alternative beta/betaStar
539  if( dZ0 < 0.2 ) {
540  internalId_.beta_ += tkpt;
541  } else if( dZ < 0.2 ) {
542  internalId_.betaStar_ += tkpt;
543  }
544  }
545  }
546  else{
547  // setting classic and alternative to be the same for now
548  float tkpt = candPt;
549  sumTkPt += tkpt;
550  bool inVtx0 = false;
551  bool inVtxOther = false;
552  if (lPack->fromPV() == pat::PackedCandidate::PVUsedInFit) inVtx0 = true;
553  if (lPack->fromPV() == 0) inVtxOther = true;
554  if (inVtx0){
555  internalId_.betaClassic_ += tkpt;
556  internalId_.beta_ += tkpt;
557  }
558  if (inVtxOther){
559  internalId_.betaStarClassic_ += tkpt;
560  internalId_.betaStar_ += tkpt;
561  }
562  }
563  }
564  // trailing candidate
565  if( lTrail == nullptr || candPt < lTrail->pt() ) {
566  lTrail = icand;
567  }
568 
569  }
570 
571  // // Finalize all variables
572  assert( !(lLead == nullptr) );
573 
574  if ( lSecond == nullptr ) { lSecond = lTrail; }
575  if ( lLeadNeut == nullptr ) { lLeadNeut = lTrail; }
576  if ( lLeadEm == nullptr ) { lLeadEm = lTrail; }
577  if ( lLeadCh == nullptr ) { lLeadCh = lTrail; }
578 
579  internalId_.nCharged_ = pfjet->chargedMultiplicity();
580  internalId_.nNeutrals_ = pfjet->neutralMultiplicity();
581  internalId_.chgEMfrac_ = pfjet->chargedEmEnergy() /jet->energy();
582  internalId_.neuEMfrac_ = pfjet->neutralEmEnergy() /jet->energy();
585  internalId_.nParticles_ = jet->numberOfDaughters();
586 
592 
593  std::sort(frac.begin(),frac.end(),std::greater<float>());
594  std::sort(fracCh.begin(),fracCh.end(),std::greater<float>());
595  std::sort(fracEm.begin(),fracEm.end(),std::greater<float>());
596  std::sort(fracNeut.begin(),fracNeut.end(),std::greater<float>());
601 
602  covMatrix(0,0) /= sumPt2;
603  covMatrix(0,1) /= sumPt2;
604  covMatrix(1,1) /= sumPt2;
605  covMatrix(1,0) = covMatrix(0,1);
606  internalId_.etaW_ = sqrt(covMatrix(0,0));
607  internalId_.phiW_ = sqrt(covMatrix(1,1));
609  TVectorD eigVals(2); eigVals = TMatrixDSymEigen(covMatrix).GetEigenValues();
610  internalId_.majW_ = sqrt(std::abs(eigVals(0)));
611  internalId_.minW_ = sqrt(std::abs(eigVals(1)));
613 
614  internalId_.dRLeadCent_ = reco::deltaR(*jet,*lLead);
615  if( lSecond == nullptr ) { internalId_.dRLead2nd_ = reco::deltaR(*jet,*lSecond); }
616  internalId_.dRMean_ /= jetPt;
620  internalId_.dR2Mean_ /= sumPt2;
621 
622  for(size_t ic=0; ic<ncones; ++ic){
623  *coneFracs[ic] /= jetPt;
624  *coneEmFracs[ic] /= jetPt;
625  *coneNeutFracs[ic] /= jetPt;
626  *coneChFracs[ic] /= jetPt;
627  }
628  //http://jets.physics.harvard.edu/qvg/
629  double ptMean = sumPt/internalId_.nParticles_;
630  double ptRMS = 0;
631  for(unsigned int i0 = 0; i0 < frac.size(); i0++) {ptRMS+=(frac[i0]-ptMean)*(frac[i0]-ptMean);}
632  ptRMS/=internalId_.nParticles_;
633  ptRMS=sqrt(ptRMS);
634 
635  internalId_.ptMean_ = ptMean;
636  internalId_.ptRMS_ = ptRMS/jetPt;
637  internalId_.pt2A_ = sqrt( internalId_.ptD_ /internalId_.nParticles_)/jetPt;
638  internalId_.ptD_ = sqrt( internalId_.ptD_) / sumPt;
639  internalId_.ptDCh_ = sqrt( internalId_.ptDCh_) / sumPtCh;
640  internalId_.ptDNe_ = sqrt( internalId_.ptDNe_) / sumPtNe;
641  internalId_.sumPt_ = sumPt;
642  internalId_.sumChPt_ = sumPtCh;
643  internalId_.sumNePt_ = sumPtNe;
644 
645  internalId_.jetR_ = lLead->pt()/sumPt;
646  internalId_.jetRchg_ = lLeadEm->pt()/sumPt;
647  internalId_.dRMatch_ = dRmin;
648 
649  if( sumTkPt != 0. ) {
650  internalId_.beta_ /= sumTkPt;
651  internalId_.betaStar_ /= sumTkPt;
652  internalId_.betaClassic_ /= sumTkPt;
653  internalId_.betaStarClassic_ /= sumTkPt;
654  } else {
655  assert( internalId_.beta_ == 0. && internalId_.betaStar_ == 0.&& internalId_.betaClassic_ == 0. && internalId_.betaStarClassic_ == 0. );
656  }
657 
658  if( runMvas_ ) {
659  runMva();
660  }
661 
663 }
664 
665 
666 
667 // ------------------------------------------------------------------------------------------
669 {
670  std::stringstream out;
671  for(variables_list_t::const_iterator it=variables_.begin();
672  it!=variables_.end(); ++it ) {
673  out << std::setw(15) << it->first << std::setw(3) << "="
674  << std::setw(5) << *it->second.first
675  << " (" << std::setw(5) << it->second.second << ")" << std::endl;
676  }
677  return out.str();
678 }
679 
680 // ------------------------------------------------------------------------------------------
682 {
683  internalId_.idFlag_ = 0;
684  for(variables_list_t::iterator it=variables_.begin();
685  it!=variables_.end(); ++it ) {
686  *it->second.first = it->second.second;
687  }
688 }
689 
690 // ------------------------------------------------------------------------------------------
691 #define INIT_VARIABLE(NAME,TMVANAME,VAL) \
692  internalId_.NAME ## _ = VAL; \
693  variables_[ # NAME ] = std::make_pair(& internalId_.NAME ## _, VAL);
694 
695 // ------------------------------------------------------------------------------------------
697 {
698  internalId_.idFlag_ = 0;
699  INIT_VARIABLE(mva , "", -100.);
700  //INIT_VARIABLE(jetPt , "jspt_1", 0.);
701  //INIT_VARIABLE(jetEta , "jseta_1", large_val);
702  INIT_VARIABLE(jetPt , "p4.fCoordinates.fPt", 0.);
703  INIT_VARIABLE(jetEta , "p4.fCoordinates.fEta", large_val);
704  INIT_VARIABLE(jetPhi , "jsphi_1", large_val);
705  INIT_VARIABLE(jetM , "jm_1", 0.);
706 
707  //INIT_VARIABLE(nCharged , "", 0.);
708  INIT_VARIABLE(nCharged , "nCh" , 0.);
709  INIT_VARIABLE(nNeutrals , "", 0.);
710 
711  INIT_VARIABLE(chgEMfrac , "", 0.);
712  INIT_VARIABLE(neuEMfrac , "", 0.);
713  INIT_VARIABLE(chgHadrfrac, "", 0.);
714  INIT_VARIABLE(neuHadrfrac, "", 0.);
715 
716  INIT_VARIABLE(d0 , "jd0_1" , -1000.);
717  INIT_VARIABLE(dZ , "jdZ_1" , -1000.);
718  //INIT_VARIABLE(nParticles , "npart_1" , 0.);
719  INIT_VARIABLE(nParticles , "nTot" , 0.);
720 
721  INIT_VARIABLE(leadPt , "lpt_1" , 0.);
722  INIT_VARIABLE(leadEta , "leta_1" , large_val);
723  INIT_VARIABLE(leadPhi , "lphi_1" , large_val);
724  INIT_VARIABLE(secondPt , "spt_1" , 0.);
725  INIT_VARIABLE(secondEta , "seta_1" , large_val);
726  INIT_VARIABLE(secondPhi , "sphi_1" , large_val);
727  INIT_VARIABLE(leadNeutPt , "lnept_1" , 0.);
728  INIT_VARIABLE(leadNeutEta, "lneeta_1" , large_val);
729  INIT_VARIABLE(leadNeutPhi, "lnephi_1" , large_val);
730  INIT_VARIABLE(leadEmPt , "lempt_1" , 0.);
731  INIT_VARIABLE(leadEmEta , "lemeta_1" , large_val);
732  INIT_VARIABLE(leadEmPhi , "lemphi_1" , large_val);
733  INIT_VARIABLE(leadChPt , "lchpt_1" , 0.);
734  INIT_VARIABLE(leadChEta , "lcheta_1" , large_val);
735  INIT_VARIABLE(leadChPhi , "lchphi_1" , large_val);
736  INIT_VARIABLE(leadFrac , "lLfr_1" , 0.);
737 
738  INIT_VARIABLE(dRLeadCent , "drlc_1" , 0.);
739  INIT_VARIABLE(dRLead2nd , "drls_1" , 0.);
740  //INIT_VARIABLE(dRMean , "drm_1" , 0.);
741  INIT_VARIABLE(dRMean , "min(pull,0.1)" , 0.);
742  INIT_VARIABLE(dRMeanNeut , "drmne_1" , 0.);
743  INIT_VARIABLE(dRMeanEm , "drem_1" , 0.);
744  INIT_VARIABLE(dRMeanCh , "drch_1" , 0.);
745  //INIT_VARIABLE(dR2Mean , "" , 0.);
746  INIT_VARIABLE(dR2Mean , "DRweighted" , 0.);
747 
748  INIT_VARIABLE(ptD , "", 0.);
749  INIT_VARIABLE(ptMean , "", 0.);
750  INIT_VARIABLE(ptRMS , "", 0.);
751  INIT_VARIABLE(pt2A , "", 0.);
752  INIT_VARIABLE(ptDCh , "", 0.);
753  INIT_VARIABLE(ptDNe , "", 0.);
754  INIT_VARIABLE(sumPt , "", 0.);
755  INIT_VARIABLE(sumChPt , "", 0.);
756  INIT_VARIABLE(sumNePt , "", 0.);
757 
758  INIT_VARIABLE(secondFrac ,"" ,0.);
759  INIT_VARIABLE(thirdFrac ,"" ,0.);
760  INIT_VARIABLE(fourthFrac ,"" ,0.);
761 
762  INIT_VARIABLE(leadChFrac ,"" ,0.);
763  INIT_VARIABLE(secondChFrac ,"" ,0.);
764  INIT_VARIABLE(thirdChFrac ,"" ,0.);
765  INIT_VARIABLE(fourthChFrac ,"" ,0.);
766 
767  INIT_VARIABLE(leadNeutFrac ,"" ,0.);
768  INIT_VARIABLE(secondNeutFrac ,"" ,0.);
769  INIT_VARIABLE(thirdNeutFrac ,"" ,0.);
770  INIT_VARIABLE(fourthNeutFrac ,"" ,0.);
771 
772  INIT_VARIABLE(leadEmFrac ,"" ,0.);
773  INIT_VARIABLE(secondEmFrac ,"" ,0.);
774  INIT_VARIABLE(thirdEmFrac ,"" ,0.);
775  INIT_VARIABLE(fourthEmFrac ,"" ,0.);
776 
777  INIT_VARIABLE(jetW ,"" ,1.);
778  INIT_VARIABLE(etaW ,"" ,1.);
779  INIT_VARIABLE(phiW ,"" ,1.);
780 
781  //INIT_VARIABLE(majW ,"" ,1.);
782  //INIT_VARIABLE(minW ,"" ,1.);
783  INIT_VARIABLE(majW ,"axisMajor" ,1.);
784  INIT_VARIABLE(minW ,"axisMinor" ,1.);
785 
786  //INIT_VARIABLE(frac01 ,"" ,0.);
787  //INIT_VARIABLE(frac02 ,"" ,0.);
788  //INIT_VARIABLE(frac03 ,"" ,0.);
789  //INIT_VARIABLE(frac04 ,"" ,0.);
790  INIT_VARIABLE(frac05 ,"" ,0.);
791  INIT_VARIABLE(frac06 ,"" ,0.);
792  INIT_VARIABLE(frac07 ,"" ,0.);
793  INIT_VARIABLE(frac01 ,"fRing0" ,0.);
794  INIT_VARIABLE(frac02 ,"fRing1" ,0.);
795  INIT_VARIABLE(frac03 ,"fRing2" ,0.);
796  INIT_VARIABLE(frac04 ,"fRing3" ,0.);
797 
798  INIT_VARIABLE(chFrac01 ,"" ,0.);
799  INIT_VARIABLE(chFrac02 ,"" ,0.);
800  INIT_VARIABLE(chFrac03 ,"" ,0.);
801  INIT_VARIABLE(chFrac04 ,"" ,0.);
802  INIT_VARIABLE(chFrac05 ,"" ,0.);
803  INIT_VARIABLE(chFrac06 ,"" ,0.);
804  INIT_VARIABLE(chFrac07 ,"" ,0.);
805 
806  INIT_VARIABLE(neutFrac01 ,"" ,0.);
807  INIT_VARIABLE(neutFrac02 ,"" ,0.);
808  INIT_VARIABLE(neutFrac03 ,"" ,0.);
809  INIT_VARIABLE(neutFrac04 ,"" ,0.);
810  INIT_VARIABLE(neutFrac05 ,"" ,0.);
811  INIT_VARIABLE(neutFrac06 ,"" ,0.);
812  INIT_VARIABLE(neutFrac07 ,"" ,0.);
813 
814  INIT_VARIABLE(emFrac01 ,"" ,0.);
815  INIT_VARIABLE(emFrac02 ,"" ,0.);
816  INIT_VARIABLE(emFrac03 ,"" ,0.);
817  INIT_VARIABLE(emFrac04 ,"" ,0.);
818  INIT_VARIABLE(emFrac05 ,"" ,0.);
819  INIT_VARIABLE(emFrac06 ,"" ,0.);
820  INIT_VARIABLE(emFrac07 ,"" ,0.);
821 
822  INIT_VARIABLE(beta ,"" ,0.);
823  INIT_VARIABLE(betaStar ,"" ,0.);
824  INIT_VARIABLE(betaClassic ,"" ,0.);
825  INIT_VARIABLE(betaStarClassic ,"" ,0.);
826 
827  INIT_VARIABLE(nvtx ,"" ,0.);
828  INIT_VARIABLE(rho ,"" ,0.);
829  INIT_VARIABLE(nTrueInt ,"" ,0.);
830 
831  INIT_VARIABLE(jetR , "", 0.);
832  INIT_VARIABLE(jetRchg , "", 0.);
833  INIT_VARIABLE(dRMatch , "", 0.);
834 
835 }
836 
837 #undef INIT_VARIABLE
const float large_val
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:250
#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 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:143
Base class for all types of Jets.
Definition: Jet.h:20
std::unique_ptr< TMVA::Reader > reader_jteta_2p5_3_
Definition: DDAxes.h:10
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
double deltaR(const T1 &t1, const T2 &t2)
Definition: deltaR.h:48
int computeCutIDflag(float betaStarClassic, float dR2Mean, float nvtx, float jetPt, float jetEta)
T eta() const
Float_t betaStarCut_[3][4][4]
int chargedMultiplicity() const
chargedMultiplicity
Definition: PFJet.h:155
const Point & position() const
position
Definition: Vertex.h:106
std::unique_ptr< TMVA::Reader > reader_jteta_3_5_
Jets made from PFObjects.
Definition: PFJet.h:21
virtual double eta() const
momentum pseudorapidity
virtual double pt() const
transverse momentum
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 double mass() const
mass
virtual double energy() const
energy
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)
float trackWeight(const TrackBaseRef &r) const
returns the weight with which a Track has contributed to the vertex-fit.
const PVAssoc fromPV(size_t ipv=0) const
T sqrt(T t)
Definition: SSEVec.h:48
std::string dumpVariables() const
std::vector< std::string > tmvaVariables_jteta_3_5_
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
virtual const reco::Track & pseudoTrack() const
Return reference to a pseudo track made with candidate kinematics, parameterized error for eta...
void assign(const std::vector< float > &vec, float &a, float &b, float &c, float &d)
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:242
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)
int computeIDflag(float mva, float jetPt, float jetEta)
double ndof() const
Definition: Vertex.h:102
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:562
tuple out
Definition: dbtoconf.py:99
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:12
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]
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)
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_
std::string fullPath() const
Definition: FileInPath.cc:165
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:544
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
tuple log
Definition: cmsBatch.py:347
Definition: DDAxes.h:10