CMS 3D CMS Logo

HEPTopTaggerV2.h
Go to the documentation of this file.
1 #ifndef RecoJets_JetAlgorithms_interface_HEPTopTaggerV2_h
2 #define RecoJets_JetAlgorithms_interface_HEPTopTaggerV2_h
3 
4 #include <vector>
5 #include <algorithm>
6 #include <cmath>
7 #include "fastjet/PseudoJet.hh"
8 #include "fastjet/ClusterSequence.hh"
9 #include "fastjet/tools/Pruner.hh"
10 #include "fastjet/tools/Filter.hh"
11 #include "fastjet/contrib/Njettiness.hh"
12 #include "fastjet/contrib/Nsubjettiness.hh"
13 #include "QjetsPlugin.h"
14 #include "CLHEP/Random/RandomEngine.h"
15 
17 
18 // Allow putting evertything into a separate namepsace
19 // Do not change next line, it's needed by the sed-code that makes the tagger CMSSW-compatible.
20 namespace external {
21 
22 using namespace std;
23 using namespace fastjet;
24 
30 
31 
33 public:
34  typedef fastjet::ClusterSequence ClusterSequence;
35  typedef fastjet::JetAlgorithm JetAlgorithm;
36  typedef fastjet::JetDefinition JetDefinition;
37  typedef fastjet::PseudoJet PseudoJet;
38 
40 
41  HEPTopTaggerV2_fixed_R(fastjet::PseudoJet jet);
42 
43  HEPTopTaggerV2_fixed_R(fastjet::PseudoJet jet,
44  double mtmass, double mwmass);
45 
46  //run tagger
47  void run();
48 
49  //settings
50  void do_qjets(bool qjets) {_do_qjets = qjets;}
51 
52  void set_mass_drop_threshold(double x) {_mass_drop_threshold = x;}
53  void set_max_subjet_mass(double x) {_max_subjet_mass = x;}
54 
55  void set_filtering_n(unsigned nfilt) {_nfilt = nfilt;}
56  void set_filtering_R(double Rfilt) {_Rfilt = Rfilt;}
57  void set_filtering_minpt_subjet(double x) {_minpt_subjet = x;}
58  void set_filtering_jetalgorithm(JetAlgorithm jet_algorithm) {_jet_algorithm_filter = jet_algorithm;}
59 
60  void set_reclustering_jetalgorithm(JetAlgorithm jet_algorithm) {_jet_algorithm_recluster = jet_algorithm;}
61 
62  void set_mode(enum Mode mode) {_mode = mode;}
63  void set_mt(double x) {_mtmass = x;}
64  void set_mw(double x) {_mwmass = x;}
65  void set_top_mass_range(double xmin, double xmax) {_mtmin = xmin; _mtmax = xmax;}
66  void set_fw(double fw) {_rmin = (1.-fw)*_mwmass/_mtmass; _rmax=(1.+fw)*_mwmass/_mtmass;}
67  void set_mass_ratio_range(double rmin, double rmax) {_rmin = rmin; _rmax = rmax;}
68  void set_mass_ratio_cut(double m23cut, double m13cutmin,double m13cutmax) {_m23cut = m23cut; _m13cutmin = m13cutmin; _m13cutmax = m13cutmax;}
69  void set_top_minpt(double x) {_minpt_tag = x;}
70 
71  void set_pruning_zcut(double zcut) {_zcut = zcut;}
72  void set_pruning_rcut_factor(double rcut_factor) {_rcut_factor = rcut_factor;}
73 
74  void set_debug(bool debug) {_debug = debug;}
75  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) {
76  _q_zcut = q_zcut; _q_dcut_fctr = q_dcut_fctr; _q_exp_min = q_exp_min; _q_exp_max = q_exp_max; _q_rigidity = q_rigidity; _q_truncation_fctr = q_truncation_fctr;
77  }
78  void set_qjets_rng(CLHEP::HepRandomEngine* engine){ _rnEngine = engine;}
79 
80  //get information
81  bool is_maybe_top() const {return _is_maybe_top;}
82  bool is_masscut_passed() const {return _is_masscut_passed;}
83  bool is_minptcut_passed() const {return _is_ptmincut_passed;}
84  bool is_tagged() const {return (_is_masscut_passed && _is_ptmincut_passed);}
85 
86  double delta_top() const {return _delta_top;}
87  double djsum() const {return _djsum;}
88  double pruned_mass() const {return _pruned_mass;}
89  double unfiltered_mass() const {return _unfiltered_mass;}
90 
91  double f_rec();
92  const PseudoJet & t() const {return _top_candidate;}
93  const PseudoJet & b() const {return _top_subjets[0];}
94  const PseudoJet & W() const {return _W;}
95  const PseudoJet & W1() const {return _top_subjets[1];}
96  const PseudoJet & W2() const {return _top_subjets[2];}
97  const std::vector<PseudoJet> & top_subjets() const {return _top_subjets;}
98  const PseudoJet & j1() const {return _top_subs[0];}
99  const PseudoJet & j2() const {return _top_subs[1];}
100  const PseudoJet & j3() const {return _top_subs[2];}
101  const std::vector<PseudoJet> & top_hadrons() const {return _top_hadrons;}
102  const std::vector<PseudoJet> & hardparts() const {return _top_parts;}
103  const PseudoJet & fat_initial() {return _fat;}
104 
105  void get_setting() const;
106  void get_info() const;
107 
108  double nsub(fastjet::PseudoJet jet, int order, fastjet::contrib::Njettiness::AxesMode axes = fastjet::contrib::Njettiness::kt_axes, double beta = 1., double R0 = 1.);
109  double q_weight() {return _qweight;}
110 
111 private:
112  bool _do_qjets;
113 
114  PseudoJet _jet;
115  PseudoJet _initial_jet;
116 
119 
121  double _mtmass, _mwmass;
122  double _mtmin, _mtmax;
123  double _rmin, _rmax;
124  double _m23cut, _m13cutmin, _m13cutmax;
125  double _minpt_tag;
126 
127  unsigned _nfilt;
128  double _Rfilt;
129  JetAlgorithm _jet_algorithm_filter;
131 
133 
134  double _zcut;
135  double _rcut_factor;
136 
137  double _q_zcut, _q_dcut_fctr, _q_exp_min, _q_exp_max, _q_rigidity, _q_truncation_fctr;
138  JetDefinition _qjet_def;
139 
140  PseudoJet _fat;
141 
142  CLHEP::HepRandomEngine* _rnEngine;
143 
144  bool _debug;
145 
149 
150  double _delta_top;
151  double _djsum;
152 
153  double _pruned_mass;
155 
156  double _fw;
157  unsigned _parts_size;
158 
159  PseudoJet _top_candidate;
160  PseudoJet _W;
161  std::vector<PseudoJet> _top_subs;
162  std::vector<PseudoJet> _top_subjets;
163  std::vector<PseudoJet> _top_hadrons;
164  std::vector<PseudoJet> _top_parts;
165 
167  double _qweight;
168 
169  //internal functions
170  void FindHardSubst(const PseudoJet& jet, std::vector<fastjet::PseudoJet>& t_parts);
171  std::vector<PseudoJet> Filtering(const std::vector <PseudoJet> & top_constits, const JetDefinition & filtering_def);
172  void store_topsubjets(const std::vector<PseudoJet>& top_subs);
173  bool check_mass_criteria(const std::vector<fastjet::PseudoJet> & top_subs) const;
174  double perp(const PseudoJet & vec, const fastjet::PseudoJet & ref);
175  double djademod (const fastjet::PseudoJet & subjet_i, const fastjet::PseudoJet & subjet_j, const fastjet::PseudoJet & ref);
176 
177  void print_banner();
178 
179 };
180 
182 public:
183  HEPTopTaggerV2();
184 
185  HEPTopTaggerV2(const fastjet::PseudoJet & jet);
186 
187 
188  HEPTopTaggerV2(const fastjet::PseudoJet & jet,
189  double mtmass, double mwmass
190  );
191 
192  ~HEPTopTaggerV2();
193 
194  //run tagger
195  void run();
196 
197  //get information
198  bool is_maybe_top() const {return _HEPTopTaggerV2_opt.is_maybe_top();}
199  bool is_masscut_passed() const {return _HEPTopTaggerV2_opt.is_masscut_passed();}
200  bool is_minptcut_passed() const {return _HEPTopTaggerV2_opt.is_minptcut_passed();}
201  bool is_tagged() const {return _HEPTopTaggerV2_opt.is_tagged();}
202 
203  double delta_top() const {return _HEPTopTaggerV2_opt.delta_top();}
204  double djsum() const {return _HEPTopTaggerV2_opt.djsum();}
205  double pruned_mass() const {return _HEPTopTaggerV2_opt.pruned_mass();}
206  double unfiltered_mass() const {return _HEPTopTaggerV2_opt.unfiltered_mass();}
207 
208  double f_rec() {return _HEPTopTaggerV2_opt.f_rec();}
209  const PseudoJet & t() const {return _HEPTopTaggerV2_opt.t();}
210  const PseudoJet & b() const {return _HEPTopTaggerV2_opt.b();}
211  const PseudoJet & W() const {return _HEPTopTaggerV2_opt.W();}
212  const PseudoJet & W1() const {return _HEPTopTaggerV2_opt.W1();}
213  const PseudoJet & W2() const {return _HEPTopTaggerV2_opt.W2();}
214  const std::vector<PseudoJet> & top_subjets() const {return _HEPTopTaggerV2_opt.top_subjets();}
215  const PseudoJet & j1() const {return _HEPTopTaggerV2_opt.j1();}
216  const PseudoJet & j2() const {return _HEPTopTaggerV2_opt.j2();}
217  const PseudoJet & j3() const {return _HEPTopTaggerV2_opt.j3();}
218  const std::vector<PseudoJet> & top_hadrons() const {return _HEPTopTaggerV2_opt.top_hadrons();}
219  const std::vector<PseudoJet> & hardparts() const {return _HEPTopTaggerV2_opt.hardparts();}
220  const PseudoJet & fat_initial() {return _fat;}
221  const PseudoJet & fat_Ropt() {return _HEPTopTaggerV2_opt.fat_initial();}
222  //HEPTopTaggerV2_fixed_R cand_Ropt(){return _HEPTopTaggerV2[_Ropt];}
223  HEPTopTaggerV2_fixed_R HEPTopTaggerV2agger(int i) {return _HEPTopTaggerV2[i];}
224 
225  double Ropt() const {return _Ropt/10.;}
226  double Ropt_calc() const {return _R_opt_calc;}
227  double pt_for_Ropt_calc() const {return _pt_for_R_opt_calc;}
228 
229  int optimalR_type();
230  double nsub_unfiltered(int order, fastjet::contrib::Njettiness::AxesMode axes = fastjet::contrib::Njettiness::kt_axes, double beta = 1., double R0 = 1.);
231  double nsub_filtered(int order, fastjet::contrib::Njettiness::AxesMode axes = fastjet::contrib::Njettiness::kt_axes, double beta = 1., double R0 = 1.);
232 
233  void get_setting() const {return _HEPTopTaggerV2_opt.get_setting();};
234  void get_info() const {return _HEPTopTaggerV2_opt.get_info();};
235 
236  double q_weight() {return _qweight;}
237 
238  //settings
239  void do_optimalR(bool optimalR) {_do_optimalR = optimalR;}
240 
241  void set_mass_drop_threshold(double x) {_mass_drop_threshold = x;}
242  void set_max_subjet_mass(double x) {_max_subjet_mass = x;}
243 
244  void set_filtering_n(unsigned nfilt) {_nfilt = nfilt;}
245  void set_filtering_R(double Rfilt) {_Rfilt = Rfilt;}
246  void set_filtering_minpt_subjet(double x) {_minpt_subjet = x;}
247  void set_filtering_jetalgorithm(JetAlgorithm jet_algorithm) {_jet_algorithm_filter = jet_algorithm;}
248 
249  void set_reclustering_jetalgorithm(JetAlgorithm jet_algorithm) {_jet_algorithm_recluster = jet_algorithm;}
250 
251  void set_mode(enum Mode mode) {_mode = mode;}
252  void set_mt(double x) {_mtmass = x;}
253  void set_mw(double x) {_mwmass = x;}
254  void set_top_mass_range(double xmin, double xmax) {_mtmin = xmin; _mtmax = xmax;}
255  void set_fw(double fw) {_rmin = (1.-fw)*_mwmass/_mtmass; _rmax=(1.+fw)*_mwmass/_mtmass;}
256  void set_mass_ratio_range(double rmin, double rmax) {_rmin = rmin; _rmax = rmax;}
257  void set_mass_ratio_cut(double m23cut, double m13cutmin,double m13cutmax) {_m23cut = m23cut; _m13cutmin = m13cutmin; _m13cutmax = m13cutmax;}
258  void set_top_minpt(double x) {_minpt_tag = x;}
259 
260  void set_optimalR_max(double x) {_max_fatjet_R = x;}
261  void set_optimalR_min(double x) {_min_fatjet_R = x;}
262  void set_optimalR_step(double x) {_step_R = x;}
263  void set_optimalR_threshold(double x) {_optimalR_threshold = x;}
264 
265  void set_filtering_optimalR_calc_R(double x) {_R_filt_optimalR_calc = x;}
266  void set_filtering_optimalR_calc_n(unsigned x) {_N_filt_optimalR_calc = x;}
267  void set_optimalR_calc_fun(double (*f)(double)) {_r_min_exp_function = f;}
268 
269  void set_optimalR_type_top_mass_range(double x, double y) {_optimalR_mmin = x; _optimalR_mmax = y;}
270  void set_optimalR_type_fw(double x) {_optimalR_fw = x;}
271  void set_optimalR_type_max_diff(double x) {_R_opt_diff = x;}
272 
273  void set_optimalR_reject_minimum(bool x) {_R_opt_reject_min = x;}
274 
275  void set_filtering_optimalR_pass_R(double x) {_R_filt_optimalR_pass = x;}
276  void set_filtering_optimalR_pass_n(unsigned x) {_N_filt_optimalR_pass = x;}
277  void set_filtering_optimalR_fail_R(double x) {_R_filt_optimalR_fail = x;}
278  void set_filtering_optimalR_fail_n(unsigned x) {_N_filt_optimalR_fail = x;}
279 
280  void set_pruning_zcut(double zcut) {_zcut = zcut;}
281  void set_pruning_rcut_factor(double rcut_factor) {_rcut_factor = rcut_factor;}
282 
283  void set_debug(bool debug) {_debug = debug;}
284  void do_qjets(bool qjets) {_do_qjets = qjets;}
285  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) {
286  _q_zcut = q_zcut; _q_dcut_fctr = q_dcut_fctr; _q_exp_min = q_exp_min; _q_exp_max = q_exp_max; _q_rigidity = q_rigidity; _q_truncation_fctr = q_truncation_fctr;
287  }
288  void set_qjets_rng(CLHEP::HepRandomEngine* engine){ _rnEngine = engine;}
289 
290 
291 private:
292  bool _do_optimalR, _do_qjets;
293 
294  PseudoJet _jet;
295  PseudoJet _initial_jet;
296 
299 
301  double _mtmass, _mwmass;
302  double _mtmin, _mtmax;
303  double _rmin, _rmax;
304  double _m23cut, _m13cutmin, _m13cutmax;
305  double _minpt_tag;
306 
307  unsigned _nfilt;
308  double _Rfilt;
309  fastjet::JetAlgorithm _jet_algorithm_filter;
311 
312  fastjet::JetAlgorithm _jet_algorithm_recluster;
313 
314  double _zcut;
315  double _rcut_factor;
316 
317  double _max_fatjet_R, _min_fatjet_R, _step_R, _optimalR_threshold;
318 
319  double _R_filt_optimalR_calc, _N_filt_optimalR_calc;
320  double (*_r_min_exp_function)(double);
321 
322  double _optimalR_mmin, _optimalR_mmax, _optimalR_fw, _R_opt_calc, _pt_for_R_opt_calc, _R_opt_diff;
324  double _R_filt_optimalR_pass, _N_filt_optimalR_pass, _R_filt_optimalR_fail, _N_filt_optimalR_fail;
325 
326  double _q_zcut, _q_dcut_fctr, _q_exp_min, _q_exp_max, _q_rigidity, _q_truncation_fctr;
327  JetDefinition _qjet_def;
328 
329  PseudoJet _fat, _filt_fat;
330  map<int,int> _n_small_fatjets;
331  map<int,HEPTopTaggerV2_fixed_R> _HEPTopTaggerV2;
333 
334  int _Ropt;
335 
336  CLHEP::HepRandomEngine* _rnEngine;
337 
338  bool _debug;
339  double _qweight;
340 
341  void UnclusterFatjets(const vector<fastjet::PseudoJet> & big_fatjets, vector<fastjet::PseudoJet> & small_fatjets, const ClusterSequence & cs, const double small_radius);
342 
343 };
344 //--------------------------------------------------------------------
345 // Do not change next line, it's needed by the sed-code that makes the tagger CMSSW-compatible.
346 };
347 
348 #endif // __HEPTOPTAGGERV2_HH__
349 
const PseudoJet & W1() const
void set_top_mass_range(double xmin, double xmax)
void do_qjets(bool qjets)
const PseudoJet & j3() const
void set_pruning_rcut_factor(double rcut_factor)
std::vector< PseudoJet > _top_subs
const PseudoJet & j3() const
void set_top_mass_range(double xmin, double xmax)
const PseudoJet & b() const
void set_filtering_n(unsigned nfilt)
static void get_info(const dqmstorepb::ROOTFilePB::Histo &h, std::string &dirname, std::string &objname, TObject **obj)
Definition: fastHadd.cc:189
void set_optimalR_type_max_diff(double x)
std::vector< PseudoJet > _top_subjets
const PseudoJet & j2() const
void set_optimalR_reject_minimum(bool x)
fastjet::JetAlgorithm _jet_algorithm_filter
unique_ptr< ClusterSequence > cs
void set_mass_ratio_range(double rmin, double rmax)
const PseudoJet & t() const
void set_qjets_rng(CLHEP::HepRandomEngine *engine)
void set_filtering_minpt_subjet(double x)
const PseudoJet & W() const
fastjet::JetAlgorithm JetAlgorithm
void set_top_minpt(double x)
void set_filtering_jetalgorithm(JetAlgorithm jet_algorithm)
void set_filtering_optimalR_pass_R(double x)
void set_filtering_optimalR_calc_n(unsigned x)
void set_mass_ratio_range(double rmin, double rmax)
const std::vector< PseudoJet > & hardparts() const
void set_mode(enum Mode mode)
const std::vector< PseudoJet > & top_subjets() const
void set_optimalR_threshold(double x)
bool is_minptcut_passed() const
CLHEP::HepRandomEngine * _rnEngine
std::vector< PseudoJet > _top_hadrons
void set_pruning_rcut_factor(double rcut_factor)
const PseudoJet & b() const
void do_optimalR(bool optimalR)
const PseudoJet & j1() const
double delta_top() const
void set_filtering_optimalR_fail_n(unsigned x)
void set_filtering_optimalR_calc_R(double x)
const std::vector< PseudoJet > & top_hadrons() const
HEPTopTaggerV2_fixed_R HEPTopTaggerV2agger(int i)
void set_optimalR_type_top_mass_range(double x, double y)
void set_filtering_optimalR_fail_R(double x)
void set_optimalR_type_fw(double x)
const PseudoJet & W2() const
double pt_for_Ropt_calc() const
void set_filtering_minpt_subjet(double x)
void set_pruning_zcut(double zcut)
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)
double unfiltered_mass() const
map< int, int > _n_small_fatjets
void set_filtering_R(double Rfilt)
void set_pruning_zcut(double zcut)
const PseudoJet & fat_Ropt()
double f[11][100]
void set_filtering_n(unsigned nfilt)
void set_debug(bool debug)
double Ropt_calc() const
void set_mode(enum Mode mode)
const std::vector< PseudoJet > & top_hadrons() const
const std::vector< PseudoJet > & top_subjets() const
void set_mass_ratio_cut(double m23cut, double m13cutmin, double m13cutmax)
void set_mass_drop_threshold(double x)
void set_max_subjet_mass(double x)
void set_filtering_optimalR_pass_n(unsigned x)
fastjet::ClusterSequence ClusterSequence
void set_filtering_jetalgorithm(JetAlgorithm jet_algorithm)
#define debug
Definition: HDRShower.cc:19
map< int, HEPTopTaggerV2_fixed_R > _HEPTopTaggerV2
const PseudoJet & W() const
fastjet::JetAlgorithm _jet_algorithm_recluster
std::vector< PseudoJet > _top_parts
void set_qjets_rng(CLHEP::HepRandomEngine *engine)
void set_optimalR_step(double x)
fastjet::JetDefinition JetDefinition
const std::vector< PseudoJet > & hardparts() const
T perp() const
Magnitude of transverse component.
bool is_masscut_passed() const
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
const PseudoJet & W1() const
double pruned_mass() const
void set_mass_ratio_cut(double m23cut, double m13cutmin, double m13cutmax)
const PseudoJet & j1() const
const PseudoJet & fat_initial()
void set_optimalR_min(double x)
const PseudoJet & W2() const
const PseudoJet & j2() const
void set_reclustering_jetalgorithm(JetAlgorithm jet_algorithm)
void set_optimalR_max(double x)
HEPTopTaggerV2_fixed_R _HEPTopTaggerV2_opt
void set_filtering_R(double Rfilt)
CLHEP::HepRandomEngine * _rnEngine
void set_optimalR_calc_fun(double(*f)(double))
void set_reclustering_jetalgorithm(JetAlgorithm jet_algorithm)