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