16 edm::LogInfo(
"HEPTopTaggerV2") <<
"#--------------------------------------------------------------------------\n";
17 edm::LogInfo(
"HEPTopTaggerV2") <<
"# HEPTopTaggerV2 - under construction \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";
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();
30 per_per -= vec_ref * vec_ref / ref_ref;
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.);
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);
66 if (parent1.m() < parent2.m())
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);
84 if (dm23 <= dm12 && dm23 <= dm13) {
88 }
else if (dm13 <= dm12 && dm13 < dm23) {
92 }
else if (dm12 < dm23 && dm12 < dm13) {
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;
113 double m13m12_square_p1 = (1 + (m13/m12) * (m13/m12));
114 double m12m13_square_p1 = (1 + (m12/m13) * (m12/m13));
117 && (m23_over_m123 > _rmin && _rmax > m23_over_m123))
119 ((m23_over_m123_square < 1 - rmin_square * m13m12_square_p1)
121 (m23_over_m123_square > 1 - rmax_square * m13m12_square_p1)
125 ((m23_over_m123_square < 1 - rmin_square * m12m13_square_p1)
127 (m23_over_m123_square > 1 - rmax_square * m12m13_square_p1)
137 fastjet::contrib::Nsubjettiness
nsub(order, axes, beta, R0);
138 return nsub.result(jet);
202 edm::LogError(
"HEPTopTaggerV2") <<
"ERROR: UNKNOWN MODE" << std::endl;
208 _qjet_def = fastjet::JetDefinition(&_qjet_plugin);
210 vector<fastjet::PseudoJet> _q_constits;
216 _qjet = sorted_by_pt(_qjet_seq->inclusive_jets())[0];
217 _qjet_seq->delete_self_when_unused();
243 if (
_debug) {
edm::LogInfo(
"HEPTopTaggerV2") <<
"< 3 hard substructures " << std::endl;}
254 for (
unsigned ll =
rr + 1; ll <
_top_parts.size(); ll++) {
280 if (topcandidate.m() <
_mtmin ||
_mtmax < topcandidate.m())
continue;
283 if (topcandidate.pieces().size() < 3)
289 std::vector <PseudoJet> top_subs = sorted_by_pt(cs_top_sub->exclusive_jets(3));
290 cs_top_sub->delete_self_when_unused();
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);
326 edm::LogError(
"HEPTopTaggerV2") <<
"ERROR: UNKNOWN MODE (IN DISTANCE MEASURE SELECTION)" << std::endl;
339 double _Rprun =
_initial_jet.validated_cluster_sequence()->jet_def().R();
340 JetDefinition jet_def_prune(fastjet::cambridge_algorithm, _Rprun);
364 edm::LogInfo(
"HEPTopTaggerV2") <<
"#--------------------------------------------------------------------------\n";
365 edm::LogInfo(
"HEPTopTaggerV2") <<
"# HEPTopTaggerV2 Result" << std::endl;
370 edm::LogInfo(
"HEPTopTaggerV2") <<
"# top candidate (pt, eta, phi): (" 379 edm::LogInfo(
"HEPTopTaggerV2") <<
"#--------------------------------------------------------------------------\n";
384 edm::LogInfo(
"HEPTopTaggerV2") <<
"#--------------------------------------------------------------------------\n";
385 edm::LogInfo(
"HEPTopTaggerV2") <<
"# HEPTopTaggerV2 Settings" << 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;
397 edm::LogInfo(
"HEPTopTaggerV2") <<
"# mass plane cuts: (m23cut, m13min, m13max) = (" 405 edm::LogInfo(
"HEPTopTaggerV2") <<
"# internal jet algorithms (0 = kt, 1 = C/A, 2 = anti-kt): " << std::endl;
408 edm::LogInfo(
"HEPTopTaggerV2") <<
"#--------------------------------------------------------------------------\n";
415 vector<fastjet::PseudoJet> & small_fatjets,
417 const double small_radius) {
418 for (
unsigned i=0;
i < big_fatjets.size();
i++) {
420 PseudoJet parent1(0, 0, 0, 0), parent2(0, 0, 0, 0);
421 bool test = cseq.has_parents(this_jet, parent1, parent2);
424 if(test) dR =
sqrt(parent1.squared_distance(parent2));
426 if (!test || dR<small_radius) {
427 small_fatjets.push_back(this_jet);
429 vector<fastjet::PseudoJet>
parents;
430 parents.push_back(parent1);
431 parents.push_back(parent2);
432 UnclusterFatjets(parents, small_fatjets, cseq, small_radius);
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),
471 double mtmass,
double mwmass
525 _qjet_def = fastjet::JetDefinition(&_qjet_plugin);
526 vector<fastjet::PseudoJet> _q_constits;
527 ClusterSequence* _qjet_seq;
529 const ClusterSequence*
_seq;
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();
545 vector<fastjet::PseudoJet> big_fatjets;
546 vector<fastjet::PseudoJet> small_fatjets;
548 big_fatjets.push_back(
_jet);
551 for (
int R = maxR;
R >= minR;
R -= stepR) {
554 if (
_debug) {
edm::LogInfo(
"HEPTopTaggerV2") <<
"R = " <<
R <<
" -> n_small_fatjets = " << small_fatjets.size() << endl;}
559 double dummy = -99999;
561 for (
unsigned i = 0;
i < small_fatjets.size();
i++) {
584 if (htt.
t().perp() > dummy) {
585 dummy = htt.
t().perp();
591 if (
_Ropt == 0 &&
R < maxR) {
598 big_fatjets = small_fatjets;
599 small_fatjets.clear();
645 fastjet::contrib::Nsubjettiness nsub(order, axes, beta, R0);
646 return nsub.result(
_fat);
650 fastjet::contrib::Nsubjettiness nsub(order, axes, beta, R0);
double _R_filt_optimalR_fail
std::vector< PseudoJet > _top_subs
double _N_filt_optimalR_fail
void SetRNEngine(CLHEP::HepRandomEngine *rnEngine)
void set_top_mass_range(double xmin, double xmax)
void set_debug(bool debug)
void set_filtering_n(unsigned nfilt)
std::vector< PseudoJet > _top_subjets
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)
void set_filtering_jetalgorithm(JetAlgorithm jet_algorithm)
void set_mode(enum Mode mode)
double perp(const PseudoJet &vec, const fastjet::PseudoJet &ref)
double _R_filt_optimalR_calc
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)
double _q_truncation_fctr
void set_top_minpt(double x)
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
void set_filtering_minpt_subjet(double x)
map< int, int > _n_small_fatjets
void set_mass_drop_threshold(double x)
double _mass_drop_threshold
double _R_filt_optimalR_pass
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)
double _mass_drop_threshold
fastjet::PseudoJet PseudoJet
double delta_R(double eta1, double phi1, double eta2, double phi2)
void set_max_subjet_mass(double x)
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)
map< int, HEPTopTaggerV2_fixed_R > _HEPTopTaggerV2
fastjet::JetAlgorithm _jet_algorithm_recluster
std::vector< PseudoJet > _top_parts
double _optimalR_threshold
JetAlgorithm _jet_algorithm_filter
double _N_filt_optimalR_calc
fastjet::JetDefinition JetDefinition
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
JetAlgorithm _jet_algorithm_recluster
double _q_truncation_fctr
double _pt_for_R_opt_calc
double(* _r_min_exp_function)(double)
Power< A, B >::type pow(const A &a, const B &b)
double _N_filt_optimalR_pass
HEPTopTaggerV2_fixed_R _HEPTopTaggerV2_opt
void set_filtering_R(double Rfilt)
CLHEP::HepRandomEngine * _rnEngine
void set_reclustering_jetalgorithm(JetAlgorithm jet_algorithm)