CMS 3D CMS Logo

HEPTopTaggerV2.cc
Go to the documentation of this file.
2 
3 // Do not change next line, it's needed by the sed-code that makes the tagger CMSSW-compatible.
4 namespace external {
5 
6 //optimal_R fit
7 double R_opt_calc_funct(double pt_filt) {
8  return 327./pt_filt;
9 }
10 
12  if (!_first_time) {return;}
13  _first_time = false;
14 
15 
16  edm::LogInfo("HEPTopTaggerV2") << "#--------------------------------------------------------------------------\n";
17  edm::LogInfo("HEPTopTaggerV2") << "# HEPTopTaggerV2 - under construction \n";
18  edm::LogInfo("HEPTopTaggerV2") << "# \n";
19  edm::LogInfo("HEPTopTaggerV2") << "# Please cite JHEP 1010 (2010) 078 [arXiv:1006.2833 [hep-ph]] \n";
20  edm::LogInfo("HEPTopTaggerV2") << "# and Phys.Rev. D89 (2014) 074047 [arXiv:1312.1504 [hep-ph]] \n";
21  edm::LogInfo("HEPTopTaggerV2") << "#--------------------------------------------------------------------------\n";
22 }
23 
24 //pt wrt a reference vector
25 double HEPTopTaggerV2_fixed_R::perp(const fastjet::PseudoJet & vec, const fastjet::PseudoJet & ref) {
26  double ref_ref = ref.px() * ref.px() + ref.py() * ref.py() + ref.pz() * ref.pz();
27  double vec_ref = vec.px() * ref.px() + vec.py() * ref.py() + vec.pz() * ref.pz();
28  double per_per = vec.px() * vec.px() + vec.py() * vec.py() + vec.pz() * vec.pz();
29  if (ref_ref > 0.)
30  per_per -= vec_ref * vec_ref / ref_ref;
31  if (per_per < 0.)
32  per_per = 0.;
33  return sqrt(per_per);
34 }
35 
36 //modified Jade distance
37 double HEPTopTaggerV2_fixed_R::djademod (const fastjet::PseudoJet& subjet_i, const fastjet::PseudoJet& subjet_j, const fastjet::PseudoJet& ref) {
38  double dj = -1.0;
39  double delta_phi = subjet_i.delta_phi_to(subjet_j);
40  double delta_eta = subjet_i.eta() - subjet_j.eta();
41  double delta_R = sqrt(delta_eta * delta_eta + delta_phi * delta_phi);
42  dj = perp(subjet_i, ref) * perp(subjet_j, ref) * pow(delta_R, 4.);
43  return dj;
44 }
45 
46 //minimal |(m_ij / m_123) / (m_w/ m_t) - 1|
48  double m12 = (_top_subs[0] + _top_subs[1]).m();
49  double m13 = (_top_subs[0] + _top_subs[2]).m();
50  double m23 = (_top_subs[1] + _top_subs[2]).m();
51  double m123 = (_top_subs[0] + _top_subs[1] + _top_subs[2]).m();
52 
53  double fw12 = fabs( (m12/m123) / (_mwmass/_mtmass) - 1);
54  double fw13 = fabs( (m13/m123) / (_mwmass/_mtmass) - 1);
55  double fw23 = fabs( (m23/m123) / (_mwmass/_mtmass) - 1);
56 
57  return std::min(fw12, std::min(fw13, fw23));
58 }
59 
60 //find hard substructures
61 void HEPTopTaggerV2_fixed_R::FindHardSubst(const PseudoJet & this_jet, std::vector<fastjet::PseudoJet> & t_parts) {
62  PseudoJet parent1(0, 0, 0, 0), parent2(0, 0, 0, 0);
63  if (this_jet.m() < _max_subjet_mass || !this_jet.validated_cs()->has_parents(this_jet, parent1, parent2)) {
64  t_parts.push_back(this_jet);
65  } else {
66  if (parent1.m() < parent2.m())
67  std::swap(parent1, parent2);
68  FindHardSubst(parent1, t_parts);
69  if (parent1.m() < _mass_drop_threshold * this_jet.m())
70  FindHardSubst(parent2, t_parts);
71  }
72 }
73 
74 //store subjets as vector<PseudoJet> with [0]->b [1]->W-jet 1 [2]->W-jet 2
75 void HEPTopTaggerV2_fixed_R::store_topsubjets(const std::vector<PseudoJet>& top_subs) {
76  _top_subjets.resize(0);
77  double m12 = (top_subs[0] + top_subs[1]).m();
78  double m13 = (top_subs[0] + top_subs[2]).m();
79  double m23 = (top_subs[1] + top_subs[2]).m();
80  double dm12 = fabs(m12 - _mwmass);
81  double dm13 = fabs(m13 - _mwmass);
82  double dm23 = fabs(m23 - _mwmass);
83 
84  if (dm23 <= dm12 && dm23 <= dm13) {
85  _top_subjets.push_back(top_subs[0]);
86  _top_subjets.push_back(top_subs[1]);
87  _top_subjets.push_back(top_subs[2]);
88  } else if (dm13 <= dm12 && dm13 < dm23) {
89  _top_subjets.push_back(top_subs[1]);
90  _top_subjets.push_back(top_subs[0]);
91  _top_subjets.push_back(top_subs[2]);
92  } else if (dm12 < dm23 && dm12 < dm13) {
93  _top_subjets.push_back(top_subs[2]);
94  _top_subjets.push_back(top_subs[0]);
95  _top_subjets.push_back(top_subs[1]);
96  }
97  _W = _top_subjets[1] + _top_subjets[2];
98  return;
99 }
100 
101 //check mass plane cuts
102 bool HEPTopTaggerV2_fixed_R::check_mass_criteria(const std::vector<PseudoJet> & top_subs) const {
103  bool is_passed = false;
104  double m12 = (top_subs[0] + top_subs[1]).m();
105  double m13 = (top_subs[0] + top_subs[2]).m();
106  double m23 = (top_subs[1] + top_subs[2]).m();
107  double m123 = (top_subs[0] + top_subs[1] + top_subs[2]).m();
108  double atan1312 = atan(m13/m12);
109  double m23_over_m123 = m23/m123;
110  double m23_over_m123_square = m23_over_m123 * m23_over_m123;
111  double rmin_square = _rmin * _rmin;
112  double rmax_square = _rmax * _rmax;
113  double m13m12_square_p1 = (1 + (m13/m12) * (m13/m12));
114  double m12m13_square_p1 = (1 + (m12/m13) * (m12/m13));
115  if (
116  (atan1312 > _m13cutmin && _m13cutmax > atan1312
117  && (m23_over_m123 > _rmin && _rmax > m23_over_m123))
118  ||
119  ((m23_over_m123_square < 1 - rmin_square * m13m12_square_p1)
120  &&
121  (m23_over_m123_square > 1 - rmax_square * m13m12_square_p1)
122  &&
123  (m23_over_m123 > _m23cut))
124  ||
125  ((m23_over_m123_square < 1 - rmin_square * m12m13_square_p1)
126  &&
127  (m23_over_m123_square > 1 - rmax_square * m12m13_square_p1)
128  &&
129  (m23_over_m123 > _m23cut))
130  ) {
131  is_passed = true;
132  }
133  return is_passed;
134 }
135 
136 double HEPTopTaggerV2_fixed_R::nsub(fastjet::PseudoJet jet, int order, fastjet::contrib::Njettiness::AxesMode axes, double beta, double R0) {
137  fastjet::contrib::Nsubjettiness nsub(order, axes, beta, R0);
138  return nsub.result(jet);
139 }
140 
143  _mode(Mode(0)), _mtmass(172.3), _mwmass(80.4), _mtmin(150.), _mtmax(200.), _rmin(0.85*80.4/172.3), _rmax(1.15*80.4/172.3),
144  _m23cut(0.35), _m13cutmin(0.2), _m13cutmax(1.3), _minpt_tag(200.),
145  _nfilt(5), _Rfilt(0.3), _jet_algorithm_filter(fastjet::cambridge_algorithm), _minpt_subjet(0.),
146  _jet_algorithm_recluster(fastjet::cambridge_algorithm),
147  _zcut(0.1), _rcut_factor(0.5),
149  //_qjet_plugin(_q_zcut, _q_dcut_fctr, _q_exp_min, _q_exp_max, _q_rigidity, _q_truncation_fctr),
151 {
152  _djsum = 0.;
153  _delta_top = 1000000000000.0;
154  _pruned_mass = 0.;
155  _unfiltered_mass = 0.;
156  _top_candidate.reset(0., 0., 0., 0.);
157  _parts_size = 0;
159  _top_subs.clear();
160  _top_subjets.clear();
161  _top_hadrons.clear();
162  _top_parts.clear();
163  _qweight= -1.;
164 
165 }
166 
168  _jet(jet), _initial_jet(jet),
170  _mode(Mode(0)), _mtmass(172.3), _mwmass(80.4), _mtmin(150.), _mtmax(200.), _rmin(0.85*80.4/172.3), _rmax(1.15*80.4/172.3),
171  _m23cut(0.35), _m13cutmin(0.2), _m13cutmax(1.3), _minpt_tag(200.),
172  _nfilt(5), _Rfilt(0.3), _jet_algorithm_filter(fastjet::cambridge_algorithm), _minpt_subjet(0.),
173  _jet_algorithm_recluster(fastjet::cambridge_algorithm),
174  _zcut(0.1), _rcut_factor(0.5),
175  _q_zcut(0.1), _q_dcut_fctr(0.5), _q_exp_min(0.), _q_exp_max(0.), _q_rigidity(0.1), _q_truncation_fctr(0.0),
176  _fat(jet), _rnEngine(nullptr),
178 {}
179 
181  double mtmass, double mwmass) : _do_qjets(false),
182  _jet(jet), _initial_jet(jet),
184  _mode(Mode(0)), _mtmass(mtmass), _mwmass(mwmass), _rmin(0.85*80.4/172.3), _rmax(1.15*80.4/172.3),
185  _m23cut(0.35), _m13cutmin(0.2), _m13cutmax(1.3), _minpt_tag(200.),
186  _nfilt(5), _Rfilt(0.3), _jet_algorithm_filter(fastjet::cambridge_algorithm), _minpt_subjet(0.),
187  _jet_algorithm_recluster(fastjet::cambridge_algorithm),
188  _zcut(0.1), _rcut_factor(0.5),
189  _q_zcut(0.1), _q_dcut_fctr(0.5), _q_exp_min(0.), _q_exp_max(0.), _q_rigidity(0.1), _q_truncation_fctr(0.0),
190  _fat(jet), _rnEngine(nullptr),
192 {}
193 
195  print_banner();
196 
201  && (_mode != TWO_STEP_FILTER) ) {
202  edm::LogError("HEPTopTaggerV2") << "ERROR: UNKNOWN MODE" << std::endl;
203  return;
204  }
205 
206  //Qjets
208  _qjet_def = fastjet::JetDefinition(&_qjet_plugin);
209  _qweight=-1;
210  vector<fastjet::PseudoJet> _q_constits;
211  ClusterSequence* _qjet_seq;
212  PseudoJet _qjet;
213  if (_do_qjets){
214  _q_constits = _initial_jet.associated_cluster_sequence()->constituents(_initial_jet);
215  _qjet_seq = new ClusterSequence(_q_constits, _qjet_def);
216  _qjet = sorted_by_pt(_qjet_seq->inclusive_jets())[0];
217  _qjet_seq->delete_self_when_unused();
218  const QjetsBaseExtras* ext =
219  dynamic_cast<const QjetsBaseExtras*>(_qjet_seq->extras());
220  _qweight=ext->weight();
221  _jet = _qjet;
222  _fat = _qjet;
223  _qjet_plugin.SetRNEngine(_rnEngine);
224  }
225 
226  //initialization
227  _djsum = 0.;
228  _delta_top = 1000000000000.0;
229  _pruned_mass = 0.;
230  _unfiltered_mass = 0.;
231  _top_candidate.reset(0., 0., 0., 0.);
232  _parts_size = 0;
234  _top_subs.clear();
235  _top_subjets.clear();
236  _top_hadrons.clear();
237  _top_parts.clear();
238 
239  //find hard substructures
241 
242  if (_top_parts.size() < 3) {
243  if (_debug) {edm::LogInfo("HEPTopTaggerV2") << "< 3 hard substructures " << std::endl;}
244  return; //such events are not interesting
245  }
246 
247  // Sort subjets-after-unclustering by pT.
248  // Necessary so that two-step-filtering can use the leading-three.
249  _top_parts=sorted_by_pt(_top_parts);
250 
251  // loop over triples
252  _top_parts = sorted_by_pt(_top_parts);
253  for (unsigned rr = 0; rr < _top_parts.size(); rr++) {
254  for (unsigned ll = rr + 1; ll < _top_parts.size(); ll++) {
255  for (unsigned kk = ll + 1; kk < _top_parts.size(); kk++) {
256 
257  // two-step filtering
258  // This means that we only look at the triplet formed by the
259  // three leading-in-pT subjets-after-unclustering.
260  if((_mode==TWO_STEP_FILTER) && rr>0)
261  continue;
262  if((_mode==TWO_STEP_FILTER) && ll>1)
263  continue;
264  if((_mode==TWO_STEP_FILTER) && kk>2)
265  continue;
266 
267  //pick triple
268  PseudoJet triple = join(_top_parts[rr], _top_parts[ll], _top_parts[kk]);
269 
270  //filtering
271  double filt_top_R
272  = std::min(_Rfilt, 0.5*sqrt(std::min(_top_parts[kk].squared_distance(_top_parts[ll]),
273  std::min(_top_parts[rr].squared_distance(_top_parts[ll]),
274  _top_parts[kk].squared_distance(_top_parts[rr])))));
275  JetDefinition filtering_def(_jet_algorithm_filter, filt_top_R);
276  fastjet::Filter filter(filtering_def, fastjet::SelectorNHardest(_nfilt) * fastjet::SelectorPtMin(_minpt_subjet));
277  PseudoJet topcandidate = filter(triple);
278 
279  //mass window cut
280  if (topcandidate.m() < _mtmin || _mtmax < topcandidate.m()) continue;
281 
282  // Sanity cut: can't recluster less than 3 objects into three subjets
283  if (topcandidate.pieces().size() < 3)
284  continue;
285 
286  // Recluster to 3 subjets and apply mass plane cuts
287  JetDefinition reclustering(_jet_algorithm_recluster, 3.14);
288  ClusterSequence * cs_top_sub = new ClusterSequence(topcandidate.constituents(), reclustering);
289  std::vector <PseudoJet> top_subs = sorted_by_pt(cs_top_sub->exclusive_jets(3));
290  cs_top_sub->delete_self_when_unused();
291 
292  // Require the third subjet to be above the pT threshold
293  if (top_subs[2].perp() < _minpt_subjet)
294  continue;
295 
296  // Modes with early 2d-massplane cuts
297  if (_mode == EARLY_MASSRATIO_SORT_MASS && !check_mass_criteria(top_subs)) {continue;}
298  if (_mode == EARLY_MASSRATIO_SORT_MODDJADE && !check_mass_criteria(top_subs)) {continue;}
299 
300  //is this candidate better than the other? -> update
301  double deltatop = fabs(topcandidate.m() - _mtmass);
302  double djsum = djademod(top_subs[0], top_subs[1], topcandidate)
303  + djademod(top_subs[0], top_subs[2], topcandidate)
304  + djademod(top_subs[1], top_subs[2], topcandidate);
305  bool better = false;
306 
307  // Modes 0 and 1 sort by top mass
310  if (deltatop < _delta_top)
311  better = true;
312  }
313  // Modes 2 and 3 sort by modified jade distance
314  else if ( (_mode == EARLY_MASSRATIO_SORT_MODDJADE)
316  if (djsum > _djsum)
317  better = true;
318  }
319  // Mode 4 is the two-step filtering. No sorting necessary as
320  // we just look at the triplet of highest pT objects after
321  // unclustering
322  else if (_mode == TWO_STEP_FILTER) {
323  better = true;
324  }
325  else {
326  edm::LogError("HEPTopTaggerV2") << "ERROR: UNKNOWN MODE (IN DISTANCE MEASURE SELECTION)" << std::endl;
327  return;
328  }
329 
330  if (better) {
331  _djsum = djsum;
332  _delta_top = deltatop;
333  _is_maybe_top = true;
334  _top_candidate = topcandidate;
335  _top_subs = top_subs;
336  store_topsubjets(top_subs);
337  _top_hadrons = topcandidate.constituents();
338  // Pruning
339  double _Rprun = _initial_jet.validated_cluster_sequence()->jet_def().R();
340  JetDefinition jet_def_prune(fastjet::cambridge_algorithm, _Rprun);
341  fastjet::Pruner pruner(jet_def_prune, _zcut, _rcut_factor);
342  PseudoJet prunedjet = pruner(triple);
343  _pruned_mass = prunedjet.m();
344  _unfiltered_mass = triple.m();
345 
346  //are all criteria fulfilled?
347  _is_masscut_passed = false;
348  if (check_mass_criteria(top_subs)) {
349  _is_masscut_passed = true;
350  }
351  _is_ptmincut_passed = false;
352  if (_top_candidate.pt() > _minpt_tag) {
353  _is_ptmincut_passed = true;
354  }
355  }//end better
356  }//end kk
357  }//end ll
358  }//end rr
359 
360  return;
361 }
362 
364  edm::LogInfo("HEPTopTaggerV2") << "#--------------------------------------------------------------------------\n";
365  edm::LogInfo("HEPTopTaggerV2") << "# HEPTopTaggerV2 Result" << std::endl;
366  edm::LogInfo("HEPTopTaggerV2") << "#" << std::endl;
367  edm::LogInfo("HEPTopTaggerV2") << "# is top candidate: " << _is_maybe_top << std::endl;
368  edm::LogInfo("HEPTopTaggerV2") << "# mass plane cuts passed: " << _is_masscut_passed << std::endl;
369  edm::LogInfo("HEPTopTaggerV2") << "# top candidate mass: " << _top_candidate.m() << std::endl;
370  edm::LogInfo("HEPTopTaggerV2") << "# top candidate (pt, eta, phi): ("
371  << _top_candidate.perp() << ", "
372  << _top_candidate.eta() << ", "
373  << _top_candidate.phi_std() << ")" << std::endl;
374  edm::LogInfo("HEPTopTaggerV2") << "# top hadrons: " << _top_hadrons.size() << std::endl;
375  edm::LogInfo("HEPTopTaggerV2") << "# hard substructures: " << _parts_size << std::endl;
376  edm::LogInfo("HEPTopTaggerV2") << "# |m - mtop| : " << _delta_top << std::endl;
377  edm::LogInfo("HEPTopTaggerV2") << "# djsum : " << _djsum << std::endl;
378  edm::LogInfo("HEPTopTaggerV2") << "# is consistency cut passed: " << _is_ptmincut_passed << std::endl;
379  edm::LogInfo("HEPTopTaggerV2") << "#--------------------------------------------------------------------------\n";
380  return;
381 }
382 
384  edm::LogInfo("HEPTopTaggerV2") << "#--------------------------------------------------------------------------\n";
385  edm::LogInfo("HEPTopTaggerV2") << "# HEPTopTaggerV2 Settings" << std::endl;
386  edm::LogInfo("HEPTopTaggerV2") << "#" << std::endl;
387  edm::LogInfo("HEPTopTaggerV2") << "# mode: " << _mode << " (0 = EARLY_MASSRATIO_SORT_MASS) " << std::endl;
388  edm::LogInfo("HEPTopTaggerV2") << "# " << " (1 = LATE_MASSRATIO_SORT_MASS) " << std::endl;
389  edm::LogInfo("HEPTopTaggerV2") << "# " << " (2 = EARLY_MASSRATIO_SORT_MODDJADE) " << std::endl;
390  edm::LogInfo("HEPTopTaggerV2") << "# " << " (3 = LATE_MASSRATIO_SORT_MODDJADE) " << std::endl;
391  edm::LogInfo("HEPTopTaggerV2") << "# " << " (4 = TWO_STEP_FILTER) " << std::endl;
392  edm::LogInfo("HEPTopTaggerV2") << "# top mass: " << _mtmass << " ";
393  edm::LogInfo("HEPTopTaggerV2") << "W mass: " << _mwmass << std::endl;
394  edm::LogInfo("HEPTopTaggerV2") << "# top mass window: [" << _mtmin << ", " << _mtmax << "]" << std::endl;
395  edm::LogInfo("HEPTopTaggerV2") << "# W mass ratio: [" << _rmin << ", " << _rmax << "] (["
396  <<_rmin*_mtmass/_mwmass<< "%, "<< _rmax*_mtmass/_mwmass << "%])"<< std::endl;
397  edm::LogInfo("HEPTopTaggerV2") << "# mass plane cuts: (m23cut, m13min, m13max) = ("
398  << _m23cut << ", " << _m13cutmin << ", " << _m13cutmax << ")" << std::endl;
399  edm::LogInfo("HEPTopTaggerV2") << "# mass_drop_threshold: " << _mass_drop_threshold << " ";
400  edm::LogInfo("HEPTopTaggerV2") << "max_subjet_mass: " << _max_subjet_mass << std::endl;
401  edm::LogInfo("HEPTopTaggerV2") << "# R_filt: " << _Rfilt << " ";
402  edm::LogInfo("HEPTopTaggerV2") << "n_filt: " << _nfilt << std::endl;
403  edm::LogInfo("HEPTopTaggerV2") << "# minimal subjet pt: " << _minpt_subjet << std::endl;
404  edm::LogInfo("HEPTopTaggerV2") << "# minimal reconstructed pt: " << _minpt_tag << std::endl;
405  edm::LogInfo("HEPTopTaggerV2") << "# internal jet algorithms (0 = kt, 1 = C/A, 2 = anti-kt): " << std::endl;
406  edm::LogInfo("HEPTopTaggerV2") << "# filtering: "<< _jet_algorithm_filter << std::endl;
407  edm::LogInfo("HEPTopTaggerV2") << "# reclustering: "<< _jet_algorithm_recluster << std::endl;
408  edm::LogInfo("HEPTopTaggerV2") << "#--------------------------------------------------------------------------\n";
409 
410  return;
411 }
412 
413 //uncluster a fat jet to subjets of given cone size
414 void HEPTopTaggerV2::UnclusterFatjets(const vector<fastjet::PseudoJet> & big_fatjets,
415  vector<fastjet::PseudoJet> & small_fatjets,
416  const ClusterSequence & cseq,
417  const double small_radius) {
418  for (unsigned i=0; i < big_fatjets.size(); i++) {
419  PseudoJet this_jet = big_fatjets[i];
420  PseudoJet parent1(0, 0, 0, 0), parent2(0, 0, 0, 0);
421  bool test = cseq.has_parents(this_jet, parent1, parent2);
422  double dR = 100;
423 
424  if(test) dR = sqrt(parent1.squared_distance(parent2));
425 
426  if (!test || dR<small_radius) {
427  small_fatjets.push_back(this_jet);
428  } else {
429  vector<fastjet::PseudoJet> parents;
430  parents.push_back(parent1);
431  parents.push_back(parent2);
432  UnclusterFatjets(parents, small_fatjets, cseq, small_radius);
433  }
434  }
435 }
436 
439  _mode(Mode(0)), _mtmass(172.3), _mwmass(80.4), _mtmin(150.), _mtmax(200.), _rmin(0.85*80.4/172.3), _rmax(1.15*80.4/172.3),
440  _m23cut(0.35), _m13cutmin(0.2), _m13cutmax(1.3), _minpt_tag(200.),
441  _nfilt(5), _Rfilt(0.3), _jet_algorithm_filter(fastjet::cambridge_algorithm), _minpt_subjet(0.),
442  _jet_algorithm_recluster(fastjet::cambridge_algorithm),
443  _zcut(0.1), _rcut_factor(0.5),
444  _max_fatjet_R(1.8), _min_fatjet_R(0.5), _step_R(0.1), _optimalR_threshold(0.2),
445  _R_filt_optimalR_calc(0.2), _N_filt_optimalR_calc(10), _r_min_exp_function(&R_opt_calc_funct),
446  _optimalR_mmin(150.), _optimalR_mmax(200.), _optimalR_fw(0.175), _R_opt_diff(0.3), _R_opt_reject_min(false),
447  _R_filt_optimalR_pass(0.2), _N_filt_optimalR_pass(5), _R_filt_optimalR_fail(0.3), _N_filt_optimalR_fail(3),
449  _debug(false)
450 {}
451 
452 HEPTopTaggerV2::HEPTopTaggerV2(const fastjet::PseudoJet & jet
454  _jet(jet), _initial_jet(jet),
456  _mode(Mode(0)), _mtmass(172.3), _mwmass(80.4), _mtmin(150.), _mtmax(200.), _rmin(0.85*80.4/172.3), _rmax(1.15*80.4/172.3),
457  _m23cut(0.35), _m13cutmin(0.2), _m13cutmax(1.3), _minpt_tag(200.),
458  _nfilt(5), _Rfilt(0.3), _jet_algorithm_filter(fastjet::cambridge_algorithm), _minpt_subjet(0.),
459  _jet_algorithm_recluster(fastjet::cambridge_algorithm),
460  _zcut(0.1), _rcut_factor(0.5),
461  _max_fatjet_R(jet.validated_cluster_sequence()->jet_def().R()), _min_fatjet_R(0.5), _step_R(0.1), _optimalR_threshold(0.2),
465  _q_zcut(0.1), _q_dcut_fctr(0.5), _q_exp_min(0.), _q_exp_max(0.), _q_rigidity(0.1), _q_truncation_fctr(0.0),
466  _fat(jet),_rnEngine(nullptr),
467  _debug(false)
468 {}
469 
470 HEPTopTaggerV2::HEPTopTaggerV2(const fastjet::PseudoJet & jet,
471  double mtmass, double mwmass
473  _jet(jet), _initial_jet(jet),
475  _mode(Mode(0)), _mtmass(mtmass), _mwmass(mwmass), _mtmin(150.), _mtmax(200.), _rmin(0.85*80.4/172.3), _rmax(1.15*80.4/172.3),
476  _m23cut(0.35), _m13cutmin(0.2), _m13cutmax(1.3), _minpt_tag(200.),
477  _nfilt(5), _Rfilt(0.3), _jet_algorithm_filter(fastjet::cambridge_algorithm), _minpt_subjet(0.),
478  _jet_algorithm_recluster(fastjet::cambridge_algorithm),
479  _zcut(0.1), _rcut_factor(0.5),
480  _max_fatjet_R(jet.validated_cluster_sequence()->jet_def().R()), _min_fatjet_R(0.5), _step_R(0.1), _optimalR_threshold(0.2),
484  _q_zcut(0.1), _q_dcut_fctr(0.5), _q_exp_min(0.), _q_exp_max(0.), _q_rigidity(0.1), _q_truncation_fctr(0.0),
485  _fat(jet), _rnEngine(nullptr),
486  _debug(false)
487 {}
488 
490 
492  int maxR = int(_max_fatjet_R * 10);
493  int minR = int(_min_fatjet_R * 10);
494  int stepR = int(_step_R * 10);
495  _qweight=-1;
496  _qjet_plugin.SetRNEngine(_rnEngine);
497 
498  if (!_do_optimalR) {
502  htt.set_filtering_n(_nfilt);
503  htt.set_filtering_R(_Rfilt);
507  htt.set_mode(_mode );
508  htt.set_mt(_mtmass);
509  htt.set_mw(_mwmass);
514  htt.set_pruning_zcut(_zcut);
516  htt.set_debug(_debug);
518  htt.run();
519 
520  _HEPTopTaggerV2[maxR] = htt;
521  _Ropt = maxR;
522  _qweight = htt.q_weight();
524  } else {
525  _qjet_def = fastjet::JetDefinition(&_qjet_plugin);
526  vector<fastjet::PseudoJet> _q_constits;
527  ClusterSequence* _qjet_seq;
528  PseudoJet _qjet;
529  const ClusterSequence* _seq;
530  _seq = _initial_jet.validated_cluster_sequence();
531  if (_do_qjets){
532  _q_constits = _initial_jet.associated_cluster_sequence()->constituents(_initial_jet);
533  _qjet_seq = new ClusterSequence(_q_constits, _qjet_def);
534  _qjet = sorted_by_pt(_qjet_seq->inclusive_jets())[0];
535  _qjet_seq->delete_self_when_unused();
536  const QjetsBaseExtras* ext =
537  dynamic_cast<const QjetsBaseExtras*>(_qjet_seq->extras());
538  _qweight=ext->weight();
539  _jet = _qjet;
540  _seq = _qjet_seq;
541  _fat = _qjet;
542  }
543 
544  // Do MultiR procedure
545  vector<fastjet::PseudoJet> big_fatjets;
546  vector<fastjet::PseudoJet> small_fatjets;
547 
548  big_fatjets.push_back(_jet);
549  _Ropt = 0;
550 
551  for (int R = maxR; R >= minR; R -= stepR) {
552  UnclusterFatjets(big_fatjets, small_fatjets, *_seq, R / 10.);
553 
554  if (_debug) {edm::LogInfo("HEPTopTaggerV2") << "R = " << R << " -> n_small_fatjets = " << small_fatjets.size() << endl;}
555 
556  _n_small_fatjets[R] = small_fatjets.size();
557 
558  // We are sorting by pt - so start with a negative dummy
559  double dummy = -99999;
560 
561  for (unsigned i = 0; i < small_fatjets.size(); i++) {
562  HEPTopTaggerV2_fixed_R htt(small_fatjets[i]);
565  htt.set_filtering_n(_nfilt);
566  htt.set_filtering_R(_Rfilt);
570  htt.set_mode(_mode );
571  htt.set_mt(_mtmass);
572  htt.set_mw(_mwmass);
577  htt.set_pruning_zcut(_zcut);
579  htt.set_debug(_debug);
581 
582  htt.run();
583 
584  if (htt.t().perp() > dummy) {
585  dummy = htt.t().perp();
586  _HEPTopTaggerV2[R] = htt;
587  }
588  } //End of loop over small_fatjets
589 
590  // Only check if we have not found Ropt yet
591  if (_Ropt == 0 && R < maxR) {
592  // If the new mass is OUTSIDE the window ..
593  if (_HEPTopTaggerV2[R].t().m() < (1-_optimalR_threshold)*_HEPTopTaggerV2[maxR].t().m())
594  // .. set _Ropt to the previous mass
595  _Ropt = R + stepR;
596  }
597 
598  big_fatjets = small_fatjets;
599  small_fatjets.clear();
600  }//End of loop over R
601 
602  // if we did not find Ropt in the loop:
603  if (_Ropt == 0 && _HEPTopTaggerV2[maxR].t().m() > 0){
604  // either pick the last value (_R_opt_reject_min=false)
605  // or leace Ropt at zero (_R_opt_reject_min=true)
606  if (_R_opt_reject_min==false)
607  _Ropt = minR;
608  }
609 
610  //for the case that there is no tag at all (< 3 hard substructures)
611  if (_Ropt == 0 && _HEPTopTaggerV2[maxR].t().m() == 0)
612  _Ropt = maxR;
613 
615 
616  Filter filter_optimalR_calc(_R_filt_optimalR_calc, SelectorNHardest(_N_filt_optimalR_calc));
617  _R_opt_calc = _r_min_exp_function(filter_optimalR_calc(_fat).pt());
618  _pt_for_R_opt_calc = filter_optimalR_calc(_fat).pt();
619 
620  Filter filter_optimalR_pass(_R_filt_optimalR_pass, SelectorNHardest(_N_filt_optimalR_pass));
621  Filter filter_optimalR_fail(_R_filt_optimalR_fail, SelectorNHardest(_N_filt_optimalR_fail));
622  if(optimalR_type() == 1) {
623  _filt_fat = filter_optimalR_pass(_fat);
624  } else {
625  _filt_fat = filter_optimalR_fail(_fat);
626  }
627  }
628 }
629 
630 //optimal_R type
633  return 0;
634  }
636  return 0;
637  }
638  if(_Ropt/10. - _R_opt_calc > _R_opt_diff) {
639  return 0;
640  }
641  return 1;
642 }
643 
644 double HEPTopTaggerV2::nsub_unfiltered(int order, fastjet::contrib::Njettiness::AxesMode axes, double beta, double R0) {
645  fastjet::contrib::Nsubjettiness nsub(order, axes, beta, R0);
646  return nsub.result(_fat);
647 }
648 
649 double HEPTopTaggerV2::nsub_filtered(int order, fastjet::contrib::Njettiness::AxesMode axes, double beta, double R0) {
650  fastjet::contrib::Nsubjettiness nsub(order, axes, beta, R0);
651  return nsub.result(_filt_fat);
652 }
653 
655 // Do not change next line, it's needed by the sed-code that makes the tagger CMSSW-compatible.
656 };
std::vector< PseudoJet > _top_subs
void SetRNEngine(CLHEP::HepRandomEngine *rnEngine)
Definition: QjetsPlugin.h:26
TPRegexp parents
Definition: eve_filter.cc:21
void set_top_mass_range(double xmin, double xmax)
void set_filtering_n(unsigned nfilt)
std::vector< PseudoJet > _top_subjets
JetDefinition jet_def
fastjet::JetAlgorithm _jet_algorithm_filter
void set_mass_ratio_range(double rmin, double rmax)
const PseudoJet & t() const
double R_opt_calc_funct(double pt_filt)
#define nullptr
void set_filtering_jetalgorithm(JetAlgorithm jet_algorithm)
void set_mode(enum Mode mode)
double perp(const PseudoJet &vec, const fastjet::PseudoJet &ref)
CLHEP::HepRandomEngine * _rnEngine
std::vector< PseudoJet > _top_hadrons
void set_pruning_rcut_factor(double rcut_factor)
bool check_mass_criteria(const std::vector< fastjet::PseudoJet > &top_subs) const
double nsub_filtered(int order, fastjet::contrib::Njettiness::AxesMode axes=fastjet::contrib::Njettiness::kt_axes, double beta=1., double R0=1.)
void store_topsubjets(const std::vector< PseudoJet > &top_subs)
void FindHardSubst(const PseudoJet &jet, std::vector< fastjet::PseudoJet > &t_parts)
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
void set_filtering_minpt_subjet(double x)
T sqrt(T t)
Definition: SSEVec.h:18
map< int, int > _n_small_fatjets
void UnclusterFatjets(const vector< fastjet::PseudoJet > &big_fatjets, vector< fastjet::PseudoJet > &small_fatjets, const ClusterSequence &cs, const double small_radius)
double nsub_unfiltered(int order, fastjet::contrib::Njettiness::AxesMode axes=fastjet::contrib::Njettiness::kt_axes, double beta=1., double R0=1.)
double djademod(const fastjet::PseudoJet &subjet_i, const fastjet::PseudoJet &subjet_j, const fastjet::PseudoJet &ref)
void set_pruning_zcut(double zcut)
T min(T a, T b)
Definition: MathUtil.h:58
double delta_R(double eta1, double phi1, double eta2, double phi2)
Definition: AnglesUtil.h:106
void set_mass_ratio_cut(double m23cut, double m13cutmin, double m13cutmax)
double nsub(fastjet::PseudoJet jet, int order, fastjet::contrib::Njettiness::AxesMode axes=fastjet::contrib::Njettiness::kt_axes, double beta=1., double R0=1.)
fastjet::ClusterSequence ClusterSequence
static std::string join(char **cmd)
Definition: RemoteFile.cc:18
map< int, HEPTopTaggerV2_fixed_R > _HEPTopTaggerV2
cms::Filter Filter
fastjet::JetAlgorithm _jet_algorithm_recluster
std::vector< PseudoJet > _top_parts
fastjet::JetDefinition JetDefinition
Definition: Filter.py:1
void set_qjets(double q_zcut, double q_dcut_fctr, double q_exp_min, double q_exp_max, double q_rigidity, double q_truncation_fctr)
const PseudoJet & t() const
Definition: memstream.h:15
double(* _r_min_exp_function)(double)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
virtual double weight() const
Definition: Qjets.h:73
HEPTopTaggerV2_fixed_R _HEPTopTaggerV2_opt
void set_filtering_R(double Rfilt)
CLHEP::HepRandomEngine * _rnEngine
void set_reclustering_jetalgorithm(JetAlgorithm jet_algorithm)