JetAlgorithms/interface/HEPTopTaggerV2.h

0001 #ifndef RecoJets_JetAlgorithms_interface_HEPTopTaggerV2_h
0002 #define RecoJets_JetAlgorithms_interface_HEPTopTaggerV2_h
0003
0004 #include <vector>
0005 #include <algorithm>
0006 #include <cmath>
0007 #include "fastjet/PseudoJet.hh"
0008 #include "fastjet/ClusterSequence.hh"
0009 #include "fastjet/tools/Pruner.hh"
0010 #include "fastjet/tools/Filter.hh"
0011 #include "fastjet/contrib/Njettiness.hh"
0012 #include "fastjet/contrib/Nsubjettiness.hh"
0013 #include "QjetsPlugin.h"
0014 #include "CLHEP/Random/RandomEngine.h"
0015
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017
0018 // Allow putting evertything into a separate namepsace
0019 // Do not change next line, it's needed by the sed-code that makes the tagger CMSSW-compatible.
0020 namespace external {
0021
0022   using namespace std;
0023   using namespace fastjet;
0024
0025   enum Mode {
0026     EARLY_MASSRATIO_SORT_MASS,
0027     LATE_MASSRATIO_SORT_MASS,
0028     EARLY_MASSRATIO_SORT_MODDJADE,
0029     LATE_MASSRATIO_SORT_MODDJADE,
0030     TWO_STEP_FILTER
0031   };
0032
0033   class HEPTopTaggerV2_fixed_R {
0034   public:
0035     typedef fastjet::ClusterSequence ClusterSequence;
0036     typedef fastjet::JetAlgorithm JetAlgorithm;
0037     typedef fastjet::JetDefinition JetDefinition;
0038     typedef fastjet::PseudoJet PseudoJet;
0039
0040     HEPTopTaggerV2_fixed_R();
0041
0042     HEPTopTaggerV2_fixed_R(fastjet::PseudoJet jet);
0043
0044     HEPTopTaggerV2_fixed_R(fastjet::PseudoJet jet, double mtmass, double mwmass);
0045
0046     //run tagger
0047     void run();
0048
0049     //settings
0050     void do_qjets(bool qjets) { _do_qjets = qjets; }
0051
0052     void set_mass_drop_threshold(double x) { _mass_drop_threshold = x; }
0053     void set_max_subjet_mass(double x) { _max_subjet_mass = x; }
0054
0055     void set_filtering_n(unsigned nfilt) { _nfilt = nfilt; }
0056     void set_filtering_R(double Rfilt) { _Rfilt = Rfilt; }
0057     void set_filtering_minpt_subjet(double x) { _minpt_subjet = x; }
0058     void set_filtering_jetalgorithm(JetAlgorithm jet_algorithm) { _jet_algorithm_filter = jet_algorithm; }
0059
0060     void set_reclustering_jetalgorithm(JetAlgorithm jet_algorithm) { _jet_algorithm_recluster = jet_algorithm; }
0061
0062     void set_mode(enum Mode mode) { _mode = mode; }
0063     void set_mt(double x) { _mtmass = x; }
0064     void set_mw(double x) { _mwmass = x; }
0065     void set_top_mass_range(double xmin, double xmax) {
0066       _mtmin = xmin;
0067       _mtmax = xmax;
0068     }
0069     void set_fw(double fw) {
0070       _rmin = (1. - fw) * _mwmass / _mtmass;
0071       _rmax = (1. + fw) * _mwmass / _mtmass;
0072     }
0073     void set_mass_ratio_range(double rmin, double rmax) {
0074       _rmin = rmin;
0075       _rmax = rmax;
0076     }
0077     void set_mass_ratio_cut(double m23cut, double m13cutmin, double m13cutmax) {
0078       _m23cut = m23cut;
0079       _m13cutmin = m13cutmin;
0080       _m13cutmax = m13cutmax;
0081     }
0082     void set_top_minpt(double x) { _minpt_tag = x; }
0083
0084     void set_pruning_zcut(double zcut) { _zcut = zcut; }
0085     void set_pruning_rcut_factor(double rcut_factor) { _rcut_factor = rcut_factor; }
0086
0087     void set_debug(bool debug) { _debug = debug; }
0088     void set_qjets(double q_zcut,
0089                    double q_dcut_fctr,
0090                    double q_exp_min,
0091                    double q_exp_max,
0092                    double q_rigidity,
0093                    double q_truncation_fctr) {
0094       _q_zcut = q_zcut;
0095       _q_dcut_fctr = q_dcut_fctr;
0096       _q_exp_min = q_exp_min;
0097       _q_exp_max = q_exp_max;
0098       _q_rigidity = q_rigidity;
0099       _q_truncation_fctr = q_truncation_fctr;
0100     }
0101     void set_qjets_rng(CLHEP::HepRandomEngine* engine) { _rnEngine = engine; }
0102
0103     //get information
0104     bool is_maybe_top() const { return _is_maybe_top; }
0105     bool is_masscut_passed() const { return _is_masscut_passed; }
0106     bool is_minptcut_passed() const { return _is_ptmincut_passed; }
0107     bool is_tagged() const { return (_is_masscut_passed && _is_ptmincut_passed); }
0108
0109     double delta_top() const { return _delta_top; }
0110     double djsum() const { return _djsum; }
0111     double pruned_mass() const { return _pruned_mass; }
0112     double unfiltered_mass() const { return _unfiltered_mass; }
0113
0114     double f_rec();
0115     const PseudoJet& t() const { return _top_candidate; }
0116     const PseudoJet& b() const { return _top_subjets[0]; }
0117     const PseudoJet& W() const { return _W; }
0118     const PseudoJet& W1() const { return _top_subjets[1]; }
0119     const PseudoJet& W2() const { return _top_subjets[2]; }
0120     const std::vector<PseudoJet>& top_subjets() const { return _top_subjets; }
0121     const PseudoJet& j1() const { return _top_subs[0]; }
0122     const PseudoJet& j2() const { return _top_subs[1]; }
0123     const PseudoJet& j3() const { return _top_subs[2]; }
0124     const std::vector<PseudoJet>& top_hadrons() const { return _top_hadrons; }
0125     const std::vector<PseudoJet>& hardparts() const { return _top_parts; }
0126     const PseudoJet& fat_initial() { return _fat; }
0127
0128     void get_setting() const;
0129     void get_info() const;
0130
0131     double nsub(fastjet::PseudoJet jet,
0132                 int order,
0133                 fastjet::contrib::Njettiness::AxesMode axes = fastjet::contrib::Njettiness::kt_axes,
0134                 double beta = 1.,
0135                 double R0 = 1.);
0136     double q_weight() { return _qweight; }
0137
0138   private:
0139     bool _do_qjets;
0140
0141     PseudoJet _jet;
0142     PseudoJet _initial_jet;
0143
0144     double _mass_drop_threshold;
0145     double _max_subjet_mass;
0146
0147     Mode _mode;
0148     double _mtmass, _mwmass;
0149     double _mtmin, _mtmax;
0150     double _rmin, _rmax;
0151     double _m23cut, _m13cutmin, _m13cutmax;
0152     double _minpt_tag;
0153
0154     unsigned _nfilt;
0155     double _Rfilt;
0156     JetAlgorithm _jet_algorithm_filter;
0157     double _minpt_subjet;
0158
0159     JetAlgorithm _jet_algorithm_recluster;
0160
0161     double _zcut;
0162     double _rcut_factor;
0163
0164     double _q_zcut, _q_dcut_fctr, _q_exp_min, _q_exp_max, _q_rigidity, _q_truncation_fctr;
0165     JetDefinition _qjet_def;
0166
0167     PseudoJet _fat;
0168
0169     CLHEP::HepRandomEngine* _rnEngine;
0170
0171     bool _debug;
0172
0173     bool _is_masscut_passed;
0174     bool _is_ptmincut_passed;
0175     bool _is_maybe_top;
0176
0177     double _delta_top;
0178     double _djsum;
0179
0180     double _pruned_mass;
0181     double _unfiltered_mass;
0182
0183     double _fw;
0184     unsigned _parts_size;
0185
0186     PseudoJet _top_candidate;
0187     PseudoJet _W;
0188     std::vector<PseudoJet> _top_subs;
0189     std::vector<PseudoJet> _top_subjets;
0190     std::vector<PseudoJet> _top_hadrons;
0191     std::vector<PseudoJet> _top_parts;
0192
0193     bool _first_time;
0194     double _qweight;
0195
0196     //internal functions
0197     void FindHardSubst(const PseudoJet& jet, std::vector<fastjet::PseudoJet>& t_parts);
0198     std::vector<PseudoJet> Filtering(const std::vector<PseudoJet>& top_constits, const JetDefinition& filtering_def);
0199     void store_topsubjets(const std::vector<PseudoJet>& top_subs);
0200     bool check_mass_criteria(const std::vector<fastjet::PseudoJet>& top_subs) const;
0201     double perp(const PseudoJet& vec, const fastjet::PseudoJet& ref);
0202     double djademod(const fastjet::PseudoJet& subjet_i,
0203                     const fastjet::PseudoJet& subjet_j,
0204                     const fastjet::PseudoJet& ref);
0205
0206     void print_banner();
0207   };
0208
0209   class HEPTopTaggerV2 {
0210   public:
0211     HEPTopTaggerV2();
0212
0213     HEPTopTaggerV2(const fastjet::PseudoJet& jet);
0214
0215     HEPTopTaggerV2(const fastjet::PseudoJet& jet, double mtmass, double mwmass);
0216
0217     ~HEPTopTaggerV2();
0218
0219     //run tagger
0220     void run();
0221
0222     //get information
0223     bool is_maybe_top() const { return _HEPTopTaggerV2_opt.is_maybe_top(); }
0224     bool is_masscut_passed() const { return _HEPTopTaggerV2_opt.is_masscut_passed(); }
0225     bool is_minptcut_passed() const { return _HEPTopTaggerV2_opt.is_minptcut_passed(); }
0226     bool is_tagged() const { return _HEPTopTaggerV2_opt.is_tagged(); }
0227
0228     double delta_top() const { return _HEPTopTaggerV2_opt.delta_top(); }
0229     double djsum() const { return _HEPTopTaggerV2_opt.djsum(); }
0230     double pruned_mass() const { return _HEPTopTaggerV2_opt.pruned_mass(); }
0231     double unfiltered_mass() const { return _HEPTopTaggerV2_opt.unfiltered_mass(); }
0232
0233     double f_rec() { return _HEPTopTaggerV2_opt.f_rec(); }
0234     const PseudoJet& t() const { return _HEPTopTaggerV2_opt.t(); }
0235     const PseudoJet& b() const { return _HEPTopTaggerV2_opt.b(); }
0236     const PseudoJet& W() const { return _HEPTopTaggerV2_opt.W(); }
0237     const PseudoJet& W1() const { return _HEPTopTaggerV2_opt.W1(); }
0238     const PseudoJet& W2() const { return _HEPTopTaggerV2_opt.W2(); }
0239     const std::vector<PseudoJet>& top_subjets() const { return _HEPTopTaggerV2_opt.top_subjets(); }
0240     const PseudoJet& j1() const { return _HEPTopTaggerV2_opt.j1(); }
0241     const PseudoJet& j2() const { return _HEPTopTaggerV2_opt.j2(); }
0242     const PseudoJet& j3() const { return _HEPTopTaggerV2_opt.j3(); }
0243     const std::vector<PseudoJet>& top_hadrons() const { return _HEPTopTaggerV2_opt.top_hadrons(); }
0244     const std::vector<PseudoJet>& hardparts() const { return _HEPTopTaggerV2_opt.hardparts(); }
0245     const PseudoJet& fat_initial() { return _fat; }
0246     const PseudoJet& fat_Ropt() { return _HEPTopTaggerV2_opt.fat_initial(); }
0247     //HEPTopTaggerV2_fixed_R cand_Ropt(){return _HEPTopTaggerV2[_Ropt];}
0248     HEPTopTaggerV2_fixed_R HEPTopTaggerV2agger(int i) { return _HEPTopTaggerV2[i]; }
0249
0250     double Ropt() const { return _Ropt / 10.; }
0251     double Ropt_calc() const { return _R_opt_calc; }
0252     double pt_for_Ropt_calc() const { return _pt_for_R_opt_calc; }
0253
0254     int optimalR_type();
0255     double nsub_unfiltered(int order,
0256                            fastjet::contrib::Njettiness::AxesMode axes = fastjet::contrib::Njettiness::kt_axes,
0257                            double beta = 1.,
0258                            double R0 = 1.);
0259     double nsub_filtered(int order,
0260                          fastjet::contrib::Njettiness::AxesMode axes = fastjet::contrib::Njettiness::kt_axes,
0261                          double beta = 1.,
0262                          double R0 = 1.);
0263
0264     void get_setting() const { return _HEPTopTaggerV2_opt.get_setting(); };
0265     void get_info() const { return _HEPTopTaggerV2_opt.get_info(); };
0266
0267     double q_weight() { return _qweight; }
0268
0269     //settings
0270     void do_optimalR(bool optimalR) { _do_optimalR = optimalR; }
0271
0272     void set_mass_drop_threshold(double x) { _mass_drop_threshold = x; }
0273     void set_max_subjet_mass(double x) { _max_subjet_mass = x; }
0274
0275     void set_filtering_n(unsigned nfilt) { _nfilt = nfilt; }
0276     void set_filtering_R(double Rfilt) { _Rfilt = Rfilt; }
0277     void set_filtering_minpt_subjet(double x) { _minpt_subjet = x; }
0278     void set_filtering_jetalgorithm(JetAlgorithm jet_algorithm) { _jet_algorithm_filter = jet_algorithm; }
0279
0280     void set_reclustering_jetalgorithm(JetAlgorithm jet_algorithm) { _jet_algorithm_recluster = jet_algorithm; }
0281
0282     void set_mode(enum Mode mode) { _mode = mode; }
0283     void set_mt(double x) { _mtmass = x; }
0284     void set_mw(double x) { _mwmass = x; }
0285     void set_top_mass_range(double xmin, double xmax) {
0286       _mtmin = xmin;
0287       _mtmax = xmax;
0288     }
0289     void set_fw(double fw) {
0290       _rmin = (1. - fw) * _mwmass / _mtmass;
0291       _rmax = (1. + fw) * _mwmass / _mtmass;
0292     }
0293     void set_mass_ratio_range(double rmin, double rmax) {
0294       _rmin = rmin;
0295       _rmax = rmax;
0296     }
0297     void set_mass_ratio_cut(double m23cut, double m13cutmin, double m13cutmax) {
0298       _m23cut = m23cut;
0299       _m13cutmin = m13cutmin;
0300       _m13cutmax = m13cutmax;
0301     }
0302     void set_top_minpt(double x) { _minpt_tag = x; }
0303
0304     void set_optimalR_max(double x) { _max_fatjet_R = x; }
0305     void set_optimalR_min(double x) { _min_fatjet_R = x; }
0306     void set_optimalR_step(double x) { _step_R = x; }
0307     void set_optimalR_threshold(double x) { _optimalR_threshold = x; }
0308
0309     void set_filtering_optimalR_calc_R(double x) { _R_filt_optimalR_calc = x; }
0310     void set_filtering_optimalR_calc_n(unsigned x) { _N_filt_optimalR_calc = x; }
0311     void set_optimalR_calc_fun(double (*f)(double)) { _r_min_exp_function = f; }
0312
0313     void set_optimalR_type_top_mass_range(double x, double y) {
0314       _optimalR_mmin = x;
0315       _optimalR_mmax = y;
0316     }
0317     void set_optimalR_type_fw(double x) { _optimalR_fw = x; }
0318     void set_optimalR_type_max_diff(double x) { _R_opt_diff = x; }
0319
0320     void set_optimalR_reject_minimum(bool x) { _R_opt_reject_min = x; }
0321
0322     void set_filtering_optimalR_pass_R(double x) { _R_filt_optimalR_pass = x; }
0323     void set_filtering_optimalR_pass_n(unsigned x) { _N_filt_optimalR_pass = x; }
0324     void set_filtering_optimalR_fail_R(double x) { _R_filt_optimalR_fail = x; }
0325     void set_filtering_optimalR_fail_n(unsigned x) { _N_filt_optimalR_fail = x; }
0326
0327     void set_pruning_zcut(double zcut) { _zcut = zcut; }
0328     void set_pruning_rcut_factor(double rcut_factor) { _rcut_factor = rcut_factor; }
0329
0330     void set_debug(bool debug) { _debug = debug; }
0331     void do_qjets(bool qjets) { _do_qjets = qjets; }
0332     void set_qjets(double q_zcut,
0333                    double q_dcut_fctr,
0334                    double q_exp_min,
0335                    double q_exp_max,
0336                    double q_rigidity,
0337                    double q_truncation_fctr) {
0338       _q_zcut = q_zcut;
0339       _q_dcut_fctr = q_dcut_fctr;
0340       _q_exp_min = q_exp_min;
0341       _q_exp_max = q_exp_max;
0342       _q_rigidity = q_rigidity;
0343       _q_truncation_fctr = q_truncation_fctr;
0344     }
0345     void set_qjets_rng(CLHEP::HepRandomEngine* engine) { _rnEngine = engine; }
0346
0347   private:
0348     bool _do_optimalR, _do_qjets;
0349
0350     PseudoJet _jet;
0351     PseudoJet _initial_jet;
0352
0353     double _mass_drop_threshold;
0354     double _max_subjet_mass;
0355
0356     Mode _mode;
0357     double _mtmass, _mwmass;
0358     double _mtmin, _mtmax;
0359     double _rmin, _rmax;
0360     double _m23cut, _m13cutmin, _m13cutmax;
0361     double _minpt_tag;
0362
0363     unsigned _nfilt;
0364     double _Rfilt;
0365     fastjet::JetAlgorithm _jet_algorithm_filter;
0366     double _minpt_subjet;
0367
0368     fastjet::JetAlgorithm _jet_algorithm_recluster;
0369
0370     double _zcut;
0371     double _rcut_factor;
0372
0373     double _max_fatjet_R, _min_fatjet_R, _step_R, _optimalR_threshold;
0374
0375     double _R_filt_optimalR_calc, _N_filt_optimalR_calc;
0376     double (*_r_min_exp_function)(double);
0377
0378     double _optimalR_mmin, _optimalR_mmax, _optimalR_fw, _R_opt_calc, _pt_for_R_opt_calc, _R_opt_diff;
0379     bool _R_opt_reject_min;
0380     double _R_filt_optimalR_pass, _N_filt_optimalR_pass, _R_filt_optimalR_fail, _N_filt_optimalR_fail;
0381
0382     double _q_zcut, _q_dcut_fctr, _q_exp_min, _q_exp_max, _q_rigidity, _q_truncation_fctr;
0383     JetDefinition _qjet_def;
0384
0385     PseudoJet _fat, _filt_fat;
0386     map<int, int> _n_small_fatjets;
0387     map<int, HEPTopTaggerV2_fixed_R> _HEPTopTaggerV2;
0388     HEPTopTaggerV2_fixed_R _HEPTopTaggerV2_opt;
0389
0390     int _Ropt;
0391
0392     CLHEP::HepRandomEngine* _rnEngine;
0393
0394     bool _debug;
0395     double _qweight;
0396
0397     void UnclusterFatjets(const vector<fastjet::PseudoJet>& big_fatjets,
0398                           vector<fastjet::PseudoJet>& small_fatjets,
0399                           const ClusterSequence& cs,
0400                           const double small_radius);
0401   };
0402   //--------------------------------------------------------------------
0403   // Do not change next line, it's needed by the sed-code that makes the tagger CMSSW-compatible.
0404 };  // namespace external
0405
0406 #endif  // __HEPTOPTAGGERV2_HH__