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  //std::string label = ps.getParameter<std::string>("label");
23  cutBased_ = ps.getParameter<bool>("cutBased");
24  if(!cutBased_)
25  {
27  tmvaMethod_ = ps.getParameter<std::string>("tmvaMethod");
28  tmvaVariables_ = ps.getParameter<std::vector<std::string> >("tmvaVariables");
29  tmvaSpectators_ = ps.getParameter<std::vector<std::string> >("tmvaSpectators");
30  version_ = ps.getParameter<int>("version");
31  }
32  else version_ = USER;
33  reader_ = nullptr;
34  edm::ParameterSet jetConfig = ps.getParameter<edm::ParameterSet>("JetIdParams");
35  for(int i0 = 0; i0 < 3; i0++) {
36  std::string lCutType = "Tight";
37  if(i0 == PileupJetIdentifier::kMedium) lCutType = "Medium";
38  if(i0 == PileupJetIdentifier::kLoose) lCutType = "Loose";
39  int nCut = 1;
40  if(cutBased_) nCut++;
41  for(int i1 = 0; i1 < nCut; i1++) {
42  std::string lFullCutType = lCutType;
43  if(cutBased_ && i1 == 0) lFullCutType = "BetaStar"+ lCutType;
44  if(cutBased_ && i1 == 1) lFullCutType = "RMS" + lCutType;
45  std::vector<double> pt010 = jetConfig.getParameter<std::vector<double> >(("Pt010_" +lFullCutType).c_str());
46  std::vector<double> pt1020 = jetConfig.getParameter<std::vector<double> >(("Pt1020_"+lFullCutType).c_str());
47  std::vector<double> pt2030 = jetConfig.getParameter<std::vector<double> >(("Pt2030_"+lFullCutType).c_str());
48  std::vector<double> pt3050 = jetConfig.getParameter<std::vector<double> >(("Pt3050_"+lFullCutType).c_str());
49  if(!cutBased_) {
50  for(int i2 = 0; i2 < 4; i2++) mvacut_[i0][0][i2] = pt010 [i2];
51  for(int i2 = 0; i2 < 4; i2++) mvacut_[i0][1][i2] = pt1020[i2];
52  for(int i2 = 0; i2 < 4; i2++) mvacut_[i0][2][i2] = pt2030[i2];
53  for(int i2 = 0; i2 < 4; i2++) mvacut_[i0][3][i2] = pt3050[i2];
54  }
55  if(cutBased_ && i1 == 0) {
56  for(int i2 = 0; i2 < 4; i2++) betaStarCut_[i0][0][i2] = pt010 [i2];
57  for(int i2 = 0; i2 < 4; i2++) betaStarCut_[i0][1][i2] = pt1020[i2];
58  for(int i2 = 0; i2 < 4; i2++) betaStarCut_[i0][2][i2] = pt2030[i2];
59  for(int i2 = 0; i2 < 4; i2++) betaStarCut_[i0][3][i2] = pt3050[i2];
60  }
61  if(cutBased_ && i1 == 1) {
62  for(int i2 = 0; i2 < 4; i2++) rmsCut_[i0][0][i2] = pt010 [i2];
63  for(int i2 = 0; i2 < 4; i2++) rmsCut_[i0][1][i2] = pt1020[i2];
64  for(int i2 = 0; i2 < 4; i2++) rmsCut_[i0][2][i2] = pt2030[i2];
65  for(int i2 = 0; i2 < 4; i2++) rmsCut_[i0][3][i2] = pt3050[i2];
66  }
67  }
68  }
69  setup();
70 }
71 
72 
73 // ------------------------------------------------------------------------------------------
75  const std::string & tmvaWeights,
76  const std::string & tmvaMethod,
77  Float_t impactParTkThreshod,
78  const std::vector<std::string> & tmvaVariables
79  )
80 {
81  impactParTkThreshod_ = impactParTkThreshod;
85  version_ = version;
86 
87  reader_ = nullptr;
88 
89  setup();
90 }
91 
92 // ------------------------------------------------------------------------------------------
94 {
95  initVariables();
96 
97  if( version_ == PHILv0 ) {
98  tmvaVariables_.clear();
99  tmvaVariables_.push_back( "jspt_1" );
100  tmvaVariables_.push_back( "jseta_1" );
101  tmvaVariables_.push_back( "jsphi_1" );
102  tmvaVariables_.push_back( "jd0_1" );
103  tmvaVariables_.push_back( "jdZ_1" );
104  tmvaVariables_.push_back( "jm_1" );
105  tmvaVariables_.push_back( "npart_1" );
106  tmvaVariables_.push_back( "lpt_1" );
107  tmvaVariables_.push_back( "leta_1" );
108  tmvaVariables_.push_back( "lphi_1" );
109  tmvaVariables_.push_back( "spt_1" );
110  tmvaVariables_.push_back( "seta_1" );
111  tmvaVariables_.push_back( "sphi_1" );
112  tmvaVariables_.push_back( "lnept_1" );
113  tmvaVariables_.push_back( "lneeta_1");
114  tmvaVariables_.push_back( "lnephi_1");
115  tmvaVariables_.push_back( "lempt_1" );
116  tmvaVariables_.push_back( "lemeta_1");
117  tmvaVariables_.push_back( "lemphi_1");
118  tmvaVariables_.push_back( "lchpt_1" );
119  // tmvaVariables_.push_back( "lcheta_1"); FIXME missing in weights file
120  tmvaVariables_.push_back( "lchphi_1");
121  tmvaVariables_.push_back( "lLfr_1" );
122  tmvaVariables_.push_back( "drlc_1" );
123  tmvaVariables_.push_back( "drls_1" );
124  tmvaVariables_.push_back( "drm_1" );
125  tmvaVariables_.push_back( "drmne_1" );
126  tmvaVariables_.push_back( "drem_1" );
127  tmvaVariables_.push_back( "drch_1" );
128 
129  tmvaNames_["jspt_1"] = "jetPt";
130  tmvaNames_["jseta_1"] = "jetEta";
131  tmvaNames_["jsphi_1"] = "jetPhi";
132  tmvaNames_["jm_1"] = "jetM";
133  tmvaNames_["jd0_1"] = "d0";
134  tmvaNames_["jdZ_1"] = "dZ";
135  tmvaNames_["npart_1"] = "nParticles";
136 
137  tmvaNames_["lpt_1"] = "leadPt";
138  tmvaNames_["leta_1"] = "leadEta";
139  tmvaNames_["lphi_1"] = "leadPhi";
140  tmvaNames_["spt_1"] = "secondPt";
141  tmvaNames_["seta_1"] = "secondEta";
142  tmvaNames_["sphi_1"] = "secondPhi";
143  tmvaNames_["lnept_1"] = "leadNeutPt";
144  tmvaNames_["lneeta_1"] = "leadNeutEta";
145  tmvaNames_["lnephi_1"] = "leadNeutPhi";
146  tmvaNames_["lempt_1"] = "leadEmPt";
147  tmvaNames_["lemeta_1"] = "leadEmEta";
148  tmvaNames_["lemphi_1"] = "leadEmPhi";
149  tmvaNames_["lchpt_1"] = "leadChPt";
150  tmvaNames_["lcheta_1"] = "leadChEta";
151  tmvaNames_["lchphi_1"] = "leadChPhi";
152  tmvaNames_["lLfr_1"] = "leadFrac";
153 
154  tmvaNames_["drlc_1"] = "dRLeadCent";
155  tmvaNames_["drls_1"] = "dRLead2nd";
156  tmvaNames_["drm_1"] = "dRMean";
157  tmvaNames_["drmne_1"] = "dRMeanNeut";
158  tmvaNames_["drem_1"] = "dRMeanEm";
159  tmvaNames_["drch_1"] = "dRMeanCh";
160 
161  } else if( ! cutBased_ ){
162  assert( tmvaMethod_.empty() || (! tmvaVariables_.empty() && version_ == USER) );
163  }
164 }
165 
166 // ------------------------------------------------------------------------------------------
168 {
169  if( reader_ ) {
170  delete reader_;
171  }
172 }
173 
174 // ------------------------------------------------------------------------------------------
175 void assign(const std::vector<float> & vec, float & a, float & b, float & c, float & d )
176 {
177  size_t sz = vec.size();
178  a = ( sz > 0 ? vec[0] : 0. );
179  b = ( sz > 1 ? vec[1] : 0. );
180  c = ( sz > 2 ? vec[2] : 0. );
181  d = ( sz > 3 ? vec[3] : 0. );
182 }
183 // ------------------------------------------------------------------------------------------
184 void setPtEtaPhi(const reco::Candidate & p, float & pt, float & eta, float &phi )
185 {
186  pt = p.pt();
187  eta = p.eta();
188  phi = p.phi();
189 }
190 
191 // ------------------------------------------------------------------------------------------
193 {
194  reader_ = new TMVA::Reader("!Color:Silent");
195  assert( ! tmvaMethod_.empty() && ! tmvaWeights_.empty() );
196  for(std::vector<std::string>::iterator it=tmvaVariables_.begin(); it!=tmvaVariables_.end(); ++it) {
197  if( tmvaNames_[*it].empty() ) {
198  tmvaNames_[*it] = *it;
199  }
200  reader_->AddVariable( *it, variables_[ tmvaNames_[*it] ].first );
201  }
202  for(std::vector<std::string>::iterator it=tmvaSpectators_.begin(); it!=tmvaSpectators_.end(); ++it) {
203  if( tmvaNames_[*it].empty() ) {
204  tmvaNames_[*it] = *it;
205  }
206  reader_->AddSpectator( *it, variables_[ tmvaNames_[*it] ].first );
207  }
209 }
210 
211 // ------------------------------------------------------------------------------------------
213 {
214  internalId_ = id;
215 }
216 
217 // ------------------------------------------------------------------------------------------
219 {
220  if( cutBased_ ) {
221  internalId_.idFlag_ = computeCutIDflag(internalId_.betaStarClassic_,internalId_.dR2Mean_,internalId_.nvtx_,internalId_.jetPt_,internalId_.jetEta_);
222  } else {
223  if( ! reader_ ) { bookReader();}
224  if(std::abs(internalId_.jetEta_) < 5.0) internalId_.mva_ = reader_->EvaluateMVA( tmvaMethod_.c_str() );
225  if(std::abs(internalId_.jetEta_) >= 5.0) internalId_.mva_ = -2.;
227  }
228 }
229 
230 // ------------------------------------------------------------------------------------------
231 std::pair<int,int> PileupJetIdAlgo::getJetIdKey(float jetPt, float jetEta)
232 {
233  int ptId = 0;
234  if(jetPt >= 10 && jetPt < 20) ptId = 1;
235  if(jetPt >= 20 && jetPt < 30) ptId = 2;
236  if(jetPt >= 30 ) ptId = 3;
237 
238  int etaId = 0;
239  if(std::abs(jetEta) >= 2.5 && std::abs(jetEta) < 2.75) etaId = 1;
240  if(std::abs(jetEta) >= 2.75 && std::abs(jetEta) < 3.0 ) etaId = 2;
241  if(std::abs(jetEta) >= 3.0 && std::abs(jetEta) < 5.0 ) etaId = 3;
242 
243  return std::pair<int,int>(ptId,etaId);
244 }
245 // ------------------------------------------------------------------------------------------
246 int PileupJetIdAlgo::computeCutIDflag(float betaStarClassic,float dR2Mean,float nvtx, float jetPt, float jetEta)
247 {
248  std::pair<int,int> jetIdKey = getJetIdKey(jetPt,jetEta);
249  float betaStarModified = betaStarClassic/log(nvtx-0.64);
250  int idFlag(0);
251  if(betaStarModified < betaStarCut_[PileupJetIdentifier::kTight ][jetIdKey.first][jetIdKey.second] &&
252  dR2Mean < rmsCut_ [PileupJetIdentifier::kTight ][jetIdKey.first][jetIdKey.second]
253  ) idFlag += 1 << PileupJetIdentifier::kTight;
254 
255  if(betaStarModified < betaStarCut_[PileupJetIdentifier::kMedium ][jetIdKey.first][jetIdKey.second] &&
256  dR2Mean < rmsCut_ [PileupJetIdentifier::kMedium ][jetIdKey.first][jetIdKey.second]
257  ) idFlag += 1 << PileupJetIdentifier::kMedium;
258 
259  if(betaStarModified < betaStarCut_[PileupJetIdentifier::kLoose ][jetIdKey.first][jetIdKey.second] &&
260  dR2Mean < rmsCut_ [PileupJetIdentifier::kLoose ][jetIdKey.first][jetIdKey.second]
261  ) idFlag += 1 << PileupJetIdentifier::kLoose;
262  return idFlag;
263 }
264 // ------------------------------------------------------------------------------------------
265 int PileupJetIdAlgo::computeIDflag(float mva, float jetPt, float jetEta)
266 {
267  std::pair<int,int> jetIdKey = getJetIdKey(jetPt,jetEta);
268  return computeIDflag(mva,jetIdKey.first,jetIdKey.second);
269 }
270 
271 // ------------------------------------------------------------------------------------------
272 int PileupJetIdAlgo::computeIDflag(float mva,int ptId,int etaId)
273 {
274  int idFlag(0);
275  if(mva > mvacut_[PileupJetIdentifier::kTight ][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kTight;
276  if(mva > mvacut_[PileupJetIdentifier::kMedium][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kMedium;
277  if(mva > mvacut_[PileupJetIdentifier::kLoose ][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kLoose;
278  return idFlag;
279 }
280 
281 
282 // ------------------------------------------------------------------------------------------
284 {
285  runMva();
287 }
288 
289 // ------------------------------------------------------------------------------------------
291  const reco::VertexCollection & allvtx,
292  bool calculateMva)
293 {
294 
295  static std::atomic<int> printWarning{10};
296 
297  // initialize all variables to 0
298  resetVariables();
299 
300  // loop over constituents, accumulate sums and find leading candidates
301  const pat::Jet * patjet = dynamic_cast<const pat::Jet *>(jet);
302  const reco::PFJet * pfjet = dynamic_cast<const reco::PFJet *>(jet);
303  assert( patjet != nullptr || pfjet != nullptr );
304  if( patjet != nullptr && jec == 0. ) { // if this is a pat jet and no jec has been passed take the jec from the object
305  jec = patjet->pt()/patjet->correctedJet(0).pt();
306  }
307  if( jec <= 0. ) {
308  jec = 1.;
309  }
310 
311  const reco::Candidate* lLead = nullptr, *lSecond = nullptr, *lLeadNeut = nullptr, *lLeadEm = nullptr, *lLeadCh = nullptr, *lTrail = nullptr;
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_ };
323  TMatrixDSym covMatrix(2); covMatrix = 0.;
324  float jetPt = jet->pt() / jec; // use uncorrected pt for shape variables
325  float sumPt = 0., sumPt2 = 0., sumTkPt = 0.,sumPtCh=0,sumPtNe = 0;
326  setPtEtaPhi(*jet,internalId_.jetPt_,internalId_.jetEta_,internalId_.jetPhi_); // use corrected pt for jet kinematics
327  internalId_.jetM_ = jet->mass();
328  internalId_.nvtx_ = allvtx.size();
329 
330  for ( unsigned i = 0; i < jet->numberOfSourceCandidatePtrs(); ++i ) {
331  reco::CandidatePtr pfJetConstituent = jet->sourceCandidatePtr(i);
332 
333  const reco::Candidate* icand = pfJetConstituent.get();
334  const pat::PackedCandidate* lPack = dynamic_cast<const pat::PackedCandidate *>( icand );
335  const reco::PFCandidate *lPF=dynamic_cast<const reco::PFCandidate*>( icand );
336  bool isPacked = true;
337  if (lPack == nullptr){
338  isPacked = false;
339  }
340 
341  float candPt = icand->pt();
342  float candPtFrac = candPt/jetPt;
343  float candDr = reco::deltaR(*icand,*jet);
344  float candDeta = std::abs( icand->eta() - jet->eta() );
345  float candDphi = reco::deltaPhi(*icand,*jet);
346  float candPtDr = candPt * candDr;
347  size_t icone = std::lower_bound(&cones[0],&cones[ncones],candDr) - &cones[0];
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 
510 
511  std::sort(frac.begin(),frac.end(),std::greater<float>());
512  std::sort(fracCh.begin(),fracCh.end(),std::greater<float>());
513  std::sort(fracEm.begin(),fracEm.end(),std::greater<float>());
514  std::sort(fracNeut.begin(),fracNeut.end(),std::greater<float>());
519 
520  covMatrix(0,0) /= sumPt2;
521  covMatrix(0,1) /= sumPt2;
522  covMatrix(1,1) /= sumPt2;
523  covMatrix(1,0) = covMatrix(0,1);
524  internalId_.etaW_ = sqrt(covMatrix(0,0));
525  internalId_.phiW_ = sqrt(covMatrix(1,1));
527  TVectorD eigVals(2); eigVals = TMatrixDSymEigen(covMatrix).GetEigenValues();
528  internalId_.majW_ = sqrt(std::abs(eigVals(0)));
529  internalId_.minW_ = sqrt(std::abs(eigVals(1)));
531 
532  internalId_.dRLeadCent_ = reco::deltaR(*jet,*lLead);
533  if( lSecond == nullptr ) { internalId_.dRLead2nd_ = reco::deltaR(*jet,*lSecond); }
534  internalId_.dRMean_ /= jetPt;
538  internalId_.dR2Mean_ /= sumPt2;
539 
540  for(size_t ic=0; ic<ncones; ++ic){
541  *coneFracs[ic] /= jetPt;
542  *coneEmFracs[ic] /= jetPt;
543  *coneNeutFracs[ic] /= jetPt;
544  *coneChFracs[ic] /= jetPt;
545  }
546  //http://jets.physics.harvard.edu/qvg/
547  double ptMean = sumPt/internalId_.nParticles_;
548  double ptRMS = 0;
549  for(unsigned int i0 = 0; i0 < frac.size(); i0++) {ptRMS+=(frac[i0]-ptMean)*(frac[i0]-ptMean);}
550  ptRMS/=internalId_.nParticles_;
551  ptRMS=sqrt(ptRMS);
552 
553  internalId_.ptMean_ = ptMean;
554  internalId_.ptRMS_ = ptRMS/jetPt;
555  internalId_.pt2A_ = sqrt( internalId_.ptD_ /internalId_.nParticles_)/jetPt;
556  internalId_.ptD_ = sqrt( internalId_.ptD_) / sumPt;
557  internalId_.ptDCh_ = sqrt( internalId_.ptDCh_) / sumPtCh;
558  internalId_.ptDNe_ = sqrt( internalId_.ptDNe_) / sumPtNe;
559  internalId_.sumPt_ = sumPt;
560  internalId_.sumChPt_ = sumPtCh;
561  internalId_.sumNePt_ = sumPtNe;
562 
563  if( sumTkPt != 0. ) {
564  internalId_.beta_ /= sumTkPt;
565  internalId_.betaStar_ /= sumTkPt;
566  internalId_.betaClassic_ /= sumTkPt;
567  internalId_.betaStarClassic_ /= sumTkPt;
568  } else {
569  assert( internalId_.beta_ == 0. && internalId_.betaStar_ == 0.&& internalId_.betaClassic_ == 0. && internalId_.betaStarClassic_ == 0. );
570  }
571 
572  if( calculateMva ) {
573  runMva();
574  }
575 
577 }
578 
579 
580 
581 // ------------------------------------------------------------------------------------------
583 {
584  std::stringstream out;
585  for(variables_list_t::const_iterator it=variables_.begin();
586  it!=variables_.end(); ++it ) {
587  out << std::setw(15) << it->first << std::setw(3) << "="
588  << std::setw(5) << *it->second.first
589  << " (" << std::setw(5) << it->second.second << ")" << std::endl;
590  }
591  return out.str();
592 }
593 
594 // ------------------------------------------------------------------------------------------
596 {
597  internalId_.idFlag_ = 0;
598  for(variables_list_t::iterator it=variables_.begin();
599  it!=variables_.end(); ++it ) {
600  *it->second.first = it->second.second;
601  }
602 }
603 
604 // ------------------------------------------------------------------------------------------
605 #define INIT_VARIABLE(NAME,TMVANAME,VAL) \
606  internalId_.NAME ## _ = VAL; \
607  variables_[ # NAME ] = std::make_pair(& internalId_.NAME ## _, VAL);
608 
609 // ------------------------------------------------------------------------------------------
611 {
612  internalId_.idFlag_ = 0;
613  INIT_VARIABLE(mva , "", -100.);
614  INIT_VARIABLE(jetPt , "jspt_1", 0.);
615  INIT_VARIABLE(jetEta , "jseta_1", large_val);
616  INIT_VARIABLE(jetPhi , "jsphi_1", large_val);
617  INIT_VARIABLE(jetM , "jm_1", 0.);
618 
619  INIT_VARIABLE(nCharged , "", 0.);
620  INIT_VARIABLE(nNeutrals , "", 0.);
621 
622  INIT_VARIABLE(chgEMfrac , "", 0.);
623  INIT_VARIABLE(neuEMfrac , "", 0.);
624  INIT_VARIABLE(chgHadrfrac, "", 0.);
625  INIT_VARIABLE(neuHadrfrac, "", 0.);
626 
627  INIT_VARIABLE(d0 , "jd0_1" , -1000.);
628  INIT_VARIABLE(dZ , "jdZ_1" , -1000.);
629  INIT_VARIABLE(nParticles , "npart_1" , 0.);
630 
631  INIT_VARIABLE(leadPt , "lpt_1" , 0.);
632  INIT_VARIABLE(leadEta , "leta_1" , large_val);
633  INIT_VARIABLE(leadPhi , "lphi_1" , large_val);
634  INIT_VARIABLE(secondPt , "spt_1" , 0.);
635  INIT_VARIABLE(secondEta , "seta_1" , large_val);
636  INIT_VARIABLE(secondPhi , "sphi_1" , large_val);
637  INIT_VARIABLE(leadNeutPt , "lnept_1" , 0.);
638  INIT_VARIABLE(leadNeutEta, "lneeta_1" , large_val);
639  INIT_VARIABLE(leadNeutPhi, "lnephi_1" , large_val);
640  INIT_VARIABLE(leadEmPt , "lempt_1" , 0.);
641  INIT_VARIABLE(leadEmEta , "lemeta_1" , large_val);
642  INIT_VARIABLE(leadEmPhi , "lemphi_1" , large_val);
643  INIT_VARIABLE(leadChPt , "lchpt_1" , 0.);
644  INIT_VARIABLE(leadChEta , "lcheta_1" , large_val);
645  INIT_VARIABLE(leadChPhi , "lchphi_1" , large_val);
646  INIT_VARIABLE(leadFrac , "lLfr_1" , 0.);
647 
648  INIT_VARIABLE(dRLeadCent , "drlc_1" , 0.);
649  INIT_VARIABLE(dRLead2nd , "drls_1" , 0.);
650  INIT_VARIABLE(dRMean , "drm_1" , 0.);
651  INIT_VARIABLE(dRMeanNeut , "drmne_1" , 0.);
652  INIT_VARIABLE(dRMeanEm , "drem_1" , 0.);
653  INIT_VARIABLE(dRMeanCh , "drch_1" , 0.);
654  INIT_VARIABLE(dR2Mean , "" , 0.);
655 
656  INIT_VARIABLE(ptD , "", 0.);
657  INIT_VARIABLE(ptMean , "", 0.);
658  INIT_VARIABLE(ptRMS , "", 0.);
659  INIT_VARIABLE(pt2A , "", 0.);
660  INIT_VARIABLE(ptDCh , "", 0.);
661  INIT_VARIABLE(ptDNe , "", 0.);
662  INIT_VARIABLE(sumPt , "", 0.);
663  INIT_VARIABLE(sumChPt , "", 0.);
664  INIT_VARIABLE(sumNePt , "", 0.);
665 
666  INIT_VARIABLE(secondFrac ,"" ,0.);
667  INIT_VARIABLE(thirdFrac ,"" ,0.);
668  INIT_VARIABLE(fourthFrac ,"" ,0.);
669 
670  INIT_VARIABLE(leadChFrac ,"" ,0.);
671  INIT_VARIABLE(secondChFrac ,"" ,0.);
672  INIT_VARIABLE(thirdChFrac ,"" ,0.);
673  INIT_VARIABLE(fourthChFrac ,"" ,0.);
674 
675  INIT_VARIABLE(leadNeutFrac ,"" ,0.);
676  INIT_VARIABLE(secondNeutFrac ,"" ,0.);
677  INIT_VARIABLE(thirdNeutFrac ,"" ,0.);
678  INIT_VARIABLE(fourthNeutFrac ,"" ,0.);
679 
680  INIT_VARIABLE(leadEmFrac ,"" ,0.);
681  INIT_VARIABLE(secondEmFrac ,"" ,0.);
682  INIT_VARIABLE(thirdEmFrac ,"" ,0.);
683  INIT_VARIABLE(fourthEmFrac ,"" ,0.);
684 
685  INIT_VARIABLE(jetW ,"" ,1.);
686  INIT_VARIABLE(etaW ,"" ,1.);
687  INIT_VARIABLE(phiW ,"" ,1.);
688 
689  INIT_VARIABLE(majW ,"" ,1.);
690  INIT_VARIABLE(minW ,"" ,1.);
691 
692  INIT_VARIABLE(frac01 ,"" ,0.);
693  INIT_VARIABLE(frac02 ,"" ,0.);
694  INIT_VARIABLE(frac03 ,"" ,0.);
695  INIT_VARIABLE(frac04 ,"" ,0.);
696  INIT_VARIABLE(frac05 ,"" ,0.);
697  INIT_VARIABLE(frac06 ,"" ,0.);
698  INIT_VARIABLE(frac07 ,"" ,0.);
699 
700  INIT_VARIABLE(chFrac01 ,"" ,0.);
701  INIT_VARIABLE(chFrac02 ,"" ,0.);
702  INIT_VARIABLE(chFrac03 ,"" ,0.);
703  INIT_VARIABLE(chFrac04 ,"" ,0.);
704  INIT_VARIABLE(chFrac05 ,"" ,0.);
705  INIT_VARIABLE(chFrac06 ,"" ,0.);
706  INIT_VARIABLE(chFrac07 ,"" ,0.);
707 
708  INIT_VARIABLE(neutFrac01 ,"" ,0.);
709  INIT_VARIABLE(neutFrac02 ,"" ,0.);
710  INIT_VARIABLE(neutFrac03 ,"" ,0.);
711  INIT_VARIABLE(neutFrac04 ,"" ,0.);
712  INIT_VARIABLE(neutFrac05 ,"" ,0.);
713  INIT_VARIABLE(neutFrac06 ,"" ,0.);
714  INIT_VARIABLE(neutFrac07 ,"" ,0.);
715 
716  INIT_VARIABLE(emFrac01 ,"" ,0.);
717  INIT_VARIABLE(emFrac02 ,"" ,0.);
718  INIT_VARIABLE(emFrac03 ,"" ,0.);
719  INIT_VARIABLE(emFrac04 ,"" ,0.);
720  INIT_VARIABLE(emFrac05 ,"" ,0.);
721  INIT_VARIABLE(emFrac06 ,"" ,0.);
722  INIT_VARIABLE(emFrac07 ,"" ,0.);
723 
724  INIT_VARIABLE(beta ,"" ,0.);
725  INIT_VARIABLE(betaStar ,"" ,0.);
726  INIT_VARIABLE(betaClassic ,"" ,0.);
727  INIT_VARIABLE(betaStarClassic ,"" ,0.);
728 
729  INIT_VARIABLE(nvtx ,"" ,0.);
730 
731 }
732 
733 #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
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]
float chargedEmEnergy() const
chargedEmEnergy
Definition: PFJet.h:142
TMVA::Reader * reader_
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
assert(m_qm.get())
std::pair< int, int > getJetIdKey(float jetPt, float jetEta)
std::map< std::string, std::string > tmvaNames_
std::string tmvaMethod_
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
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
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
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
PileupJetIdentifier computeIdVariables(const reco::Jet *jet, float jec, const reco::Vertex *, const reco::VertexCollection &, bool calculateMva=false)
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
bool isFake() const
Definition: Vertex.h:64
unsigned int nCharged(const GenJet &jet)
Float_t rmsCut_[3][4][4]
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 >())
PileupJetIdentifier computeMva()
double b
Definition: hdecay.h:120
Analysis-level calorimeter jet class.
Definition: Jet.h:77
std::vector< std::string > tmvaSpectators_
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:39
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::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
virtual double phi() const =0
momentum azimuthal angle
virtual double eta() const =0
momentum pseudorapidity
void loadTMVAWeights(TMVA::Reader *reader, const std::string &method, const std::string &weightFile, bool verbose=false)
tuple log
Definition: cmsBatch.py:347
Definition: DDAxes.h:10