File indexing completed on 2023-03-17 11:18:22
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include "RecoJets/JetAlgorithms/interface/HEPTopTaggerWrapperV2.h"
0021
0022 #include "fastjet/Error.hh"
0023 #include "fastjet/JetDefinition.hh"
0024 #include "fastjet/ClusterSequence.hh"
0025 #include "fastjet/PseudoJet.hh"
0026 #include "fastjet/tools/Pruner.hh"
0027 #include "fastjet/tools/Filter.hh"
0028
0029 #include <cmath>
0030 #include <limits>
0031 #include <cassert>
0032 using namespace std;
0033
0034 #include "RecoJets/JetAlgorithms/interface/HEPTopTaggerV2.h"
0035
0036 FASTJET_BEGIN_NAMESPACE
0037
0038
0039
0040
0041
0042
0043 double R_min_expected_function(double x) {
0044 if (x > 700)
0045 x = 700;
0046
0047 double A = -9.42052;
0048 double B = 0.202773;
0049 double C = 4498.45;
0050 double D = -1.05737e+06;
0051 double E = 9.95494e+07;
0052
0053 return A + B * sqrt(x) + C / x + D / (x * x) + E / (x * x * x);
0054 }
0055
0056
0057
0058
0059 PseudoJet HEPTopTaggerV2::result(const PseudoJet &jet) const {
0060
0061
0062
0063 if (!jet.has_valid_cluster_sequence()) {
0064 throw Error("HEPTopTagger can only be applied on jets having an associated (and valid) ClusterSequence");
0065 }
0066
0067 external::HEPTopTaggerV2 tagger(jet);
0068
0069 external::HEPTopTaggerV2 best_tagger;
0070
0071
0072
0073 double mw_over_mt = 80.4 / 172.3;
0074 double ratio_min = mw_over_mt * (100. - massRatioWidth_) / 100.;
0075 double ratio_max = mw_over_mt * (100. + massRatioWidth_) / 100.;
0076
0077
0078 tagger.set_max_subjet_mass(subjetMass_);
0079 tagger.set_mass_drop_threshold(muCut_);
0080 tagger.set_filtering_R(filtR_);
0081 tagger.set_filtering_n(filtN_);
0082 tagger.set_filtering_minpt_subjet(minSubjetPt_);
0083
0084
0085 tagger.do_optimalR(DoOptimalR_);
0086 tagger.set_optimalR_reject_minimum(optRrejectMin_);
0087
0088
0089 tagger.set_mode((external::Mode)mode_);
0090
0091
0092 tagger.set_top_minpt(minCandPt_);
0093 tagger.set_top_mass_range(minCandMass_, maxCandMass_);
0094 tagger.set_mass_ratio_cut(minM23Cut_, minM13Cut_, maxM13Cut_);
0095 tagger.set_mass_ratio_range(ratio_min, ratio_max);
0096
0097
0098 tagger.set_optimalR_calc_fun(R_min_expected_function);
0099
0100 double qweight = -1;
0101 double qepsilon = -1;
0102 double qsigmaM = -1;
0103
0104 if (DoQjets_) {
0105 int niter(100);
0106 double q_zcut(0.1);
0107 double q_dcut_fctr(0.5);
0108 double q_exp_min(0.);
0109 double q_exp_max(0.);
0110 double q_rigidity(0.1);
0111 double q_truncation_fctr(0.0);
0112
0113 double weight_q1 = -1.;
0114 double m_sum = 0.;
0115 double m2_sum = 0.;
0116 int qtags = 0;
0117
0118 tagger.set_qjets(q_zcut, q_dcut_fctr, q_exp_min, q_exp_max, q_rigidity, q_truncation_fctr);
0119 tagger.set_qjets_rng(engine_);
0120 tagger.do_qjets(true);
0121 tagger.run();
0122
0123 for (int iq = 0; iq < niter; iq++) {
0124 tagger.run();
0125 if (tagger.is_tagged()) {
0126 qtags++;
0127 m_sum += tagger.t().m();
0128 m2_sum += tagger.t().m() * tagger.t().m();
0129 if (tagger.q_weight() > weight_q1) {
0130 best_tagger = tagger;
0131 weight_q1 = tagger.q_weight();
0132 }
0133 }
0134 }
0135
0136 tagger = best_tagger;
0137 qweight = weight_q1;
0138 qepsilon = float(qtags) / float(niter);
0139
0140
0141 if (qtags > 0) {
0142 double mean_m = m_sum / qtags;
0143 double mean_m2 = m2_sum / qtags;
0144 qsigmaM = sqrt(mean_m2 - mean_m * mean_m);
0145 }
0146 } else {
0147 tagger.run();
0148 }
0149
0150
0151
0152
0153
0154
0155 if (!tagger.is_tagged())
0156 return PseudoJet();
0157
0158
0159 const JetDefinition::Recombiner *rec = jet.associated_cluster_sequence()->jet_def().recombiner();
0160
0161 const vector<PseudoJet> &subjets = tagger.top_subjets();
0162 assert(subjets.size() == 3);
0163
0164 PseudoJet non_W = subjets[0];
0165 PseudoJet W1 = subjets[1];
0166 PseudoJet W2 = subjets[2];
0167 PseudoJet W = join(subjets[1], subjets[2], *rec);
0168
0169 PseudoJet result = join<HEPTopTaggerV2Structure>(W1, W2, non_W, *rec);
0170 HEPTopTaggerV2Structure *s = (HEPTopTaggerV2Structure *)result.structure_non_const_ptr();
0171
0172 s->_fj_mass = jet.m();
0173 s->_fj_pt = jet.perp();
0174 s->_fj_eta = jet.eta();
0175 s->_fj_phi = jet.phi();
0176
0177 s->_top_mass = tagger.t().m();
0178 s->_pruned_mass = tagger.pruned_mass();
0179 s->_unfiltered_mass = tagger.unfiltered_mass();
0180 s->_fRec = tagger.f_rec();
0181 s->_mass_ratio_passed = tagger.is_masscut_passed();
0182
0183 if (DoOptimalR_) {
0184 s->_tau1Unfiltered = tagger.nsub_unfiltered(1);
0185 s->_tau2Unfiltered = tagger.nsub_unfiltered(2);
0186 s->_tau3Unfiltered = tagger.nsub_unfiltered(3);
0187 s->_tau1Filtered = tagger.nsub_filtered(1);
0188 s->_tau2Filtered = tagger.nsub_filtered(2);
0189 s->_tau3Filtered = tagger.nsub_filtered(3);
0190 }
0191
0192 s->_qweight = qweight;
0193 s->_qepsilon = qepsilon;
0194 s->_qsigmaM = qsigmaM;
0195
0196 if (DoOptimalR_) {
0197 s->_ropt = tagger.Ropt();
0198 s->_roptCalc = tagger.Ropt_calc();
0199 s->_ptForRoptCalc = tagger.pt_for_Ropt_calc();
0200 } else {
0201 s->_ropt = -1;
0202 s->_roptCalc = -1;
0203 s->_ptForRoptCalc = -1;
0204 }
0205
0206
0207
0208 return result;
0209 }
0210
0211 FASTJET_END_NAMESPACE