File indexing completed on 2023-03-17 11:18:21
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023 #ifndef __CMSTOPTAGGER_HH__
0024 #define __CMSTOPTAGGER_HH__
0025
0026 #include <fastjet/tools/JHTopTagger.hh>
0027 #include <fastjet/CompositeJetStructure.hh>
0028 #include <fastjet/LimitedWarning.hh>
0029 #include <fastjet/Error.hh>
0030 #include <fastjet/JetDefinition.hh>
0031 #include <fastjet/ClusterSequence.hh>
0032
0033 #include <sstream>
0034 #include <limits>
0035
0036 FASTJET_BEGIN_NAMESPACE
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051 class CMSTopTaggerStructure;
0052
0053 class CMSTopTagger : public TopTaggerBase {
0054 public:
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065 CMSTopTagger(double delta_p = 0.05, double delta_r = 0.4, double A = 0.0004);
0066
0067
0068 std::string description() const override;
0069
0070
0071
0072
0073
0074 PseudoJet result(const PseudoJet& jet) const override;
0075
0076
0077 typedef CMSTopTaggerStructure StructureType;
0078
0079 protected:
0080
0081 std::vector<PseudoJet> _split_once(const PseudoJet& jet_to_split, const PseudoJet& reference_jet) const;
0082
0083
0084
0085 void _find_min_mass(const std::vector<PseudoJet>& subjets, int& i, int& j) const;
0086
0087 double _delta_p, _delta_r, _A;
0088 mutable LimitedWarning _warnings_nonca;
0089 };
0090
0091
0092
0093 class CMSTopTaggerStructure : public JHTopTaggerStructure {
0094 public:
0095 CMSTopTaggerStructure(const std::vector<PseudoJet>& pieces, const JetDefinition::Recombiner* recombiner = nullptr)
0096 : JHTopTaggerStructure(pieces, recombiner) {}
0097
0098 protected:
0099 friend class CMSTopTagger;
0100 };
0101
0102
0103 inline CMSTopTagger::CMSTopTagger(double delta_p, double delta_r, double A)
0104 : _delta_p(delta_p), _delta_r(delta_r), _A(A) {}
0105
0106
0107
0108 inline std::string CMSTopTagger::description() const {
0109 std::ostringstream oss;
0110 oss << "CMSTopTagger with delta_p=" << _delta_p << ", delta_r=" << _delta_r << ", and A=" << _A;
0111 oss << description_of_selectors();
0112 return oss.str();
0113 }
0114
0115
0116
0117
0118 inline PseudoJet CMSTopTagger::result(const PseudoJet& jet) const {
0119
0120
0121
0122 if (!jet.has_valid_cluster_sequence()) {
0123 throw Error("CMSTopTagger can only be applied on jets having an associated (and valid) ClusterSequence");
0124 }
0125
0126
0127
0128 if (jet.validated_cs()->jet_def().jet_algorithm() != cambridge_algorithm)
0129 _warnings_nonca.warn(
0130 "CMSTopTagger should only be applied on jets from a Cambridge/Aachen clustering; use it with other algorithms "
0131 "at your own risk.");
0132
0133
0134 std::vector<PseudoJet> split0 = _split_once(jet, jet);
0135 if (split0.empty())
0136 return PseudoJet();
0137
0138
0139 std::vector<PseudoJet> subjets;
0140 for (unsigned i = 0; i < 2; i++) {
0141 std::vector<PseudoJet> split1 = _split_once(split0[i], jet);
0142 if (!split1.empty()) {
0143 subjets.push_back(split1[0]);
0144 subjets.push_back(split1[1]);
0145 } else {
0146 subjets.push_back(split0[i]);
0147 }
0148 }
0149
0150
0151 if (subjets.size() < 3)
0152 return PseudoJet();
0153
0154
0155 int ii = -1, jj = -1;
0156 _find_min_mass(subjets, ii, jj);
0157
0158
0159
0160
0161
0162
0163 if (ii > 0)
0164 std::swap(subjets[ii], subjets[0]);
0165 if (jj > 1)
0166 std::swap(subjets[jj], subjets[1]);
0167 if (subjets[0].perp2() < subjets[1].perp2())
0168 std::swap(subjets[0], subjets[1]);
0169 if ((subjets.size() > 3) && (subjets[2].perp2() < subjets[3].perp2()))
0170 std::swap(subjets[2], subjets[3]);
0171
0172
0173 const JetDefinition::Recombiner* rec = jet.associated_cluster_sequence()->jet_def().recombiner();
0174
0175 PseudoJet W = join(subjets[0], subjets[1], *rec);
0176 PseudoJet non_W;
0177 if (subjets.size() > 3) {
0178 non_W = join(subjets[2], subjets[3], *rec);
0179 } else {
0180 non_W = join(subjets[2], *rec);
0181 }
0182 PseudoJet result = join<CMSTopTaggerStructure>(W, non_W, *rec);
0183 CMSTopTaggerStructure* s = (CMSTopTaggerStructure*)result.structure_non_const_ptr();
0184 s->_cos_theta_w = _cos_theta_W(result);
0185
0186
0187
0188
0189 if (!_top_selector.pass(result) || !_W_selector.pass(W)) {
0190 result *= 0.0;
0191 }
0192
0193 result = join(
0194 subjets);
0195
0196 return result;
0197 }
0198
0199
0200 inline std::vector<PseudoJet> CMSTopTagger::_split_once(const PseudoJet& jet_to_split,
0201 const PseudoJet& reference_jet) const {
0202 PseudoJet this_jet = jet_to_split;
0203 PseudoJet p1, p2;
0204 std::vector<PseudoJet> result;
0205 while (this_jet.has_parents(p1, p2)) {
0206 if (p2.perp2() > p1.perp2())
0207 std::swap(p1, p2);
0208 if (p1.perp() < _delta_p * reference_jet.perp())
0209 break;
0210 double DR = p1.delta_R(p2);
0211 if (DR < _delta_r - _A * this_jet.perp())
0212 break;
0213 if (p2.perp() < _delta_p * reference_jet.perp()) {
0214 this_jet = p1;
0215 continue;
0216 }
0217
0218 result.push_back(p1);
0219 result.push_back(p2);
0220 break;
0221 }
0222 return result;
0223 }
0224
0225
0226 inline void CMSTopTagger::_find_min_mass(const std::vector<PseudoJet>& subjets, int& i, int& j) const {
0227 assert(subjets.size() > 1 && subjets.size() < 5);
0228
0229
0230 unsigned softest = 5;
0231 if (subjets.size() == 4) {
0232 double min_pt = std::numeric_limits<double>::max();
0233 ;
0234 for (unsigned ii = 0; ii < subjets.size(); ++ii) {
0235 if (subjets[ii].perp() < min_pt) {
0236 min_pt = subjets[ii].perp();
0237 softest = ii;
0238 }
0239 }
0240 }
0241
0242 double min_mass = std::numeric_limits<double>::max();
0243 for (unsigned ii = 0; ii + 1 < subjets.size(); ++ii) {
0244 if (ii == softest)
0245 continue;
0246 for (unsigned jj = ii + 1; jj < subjets.size(); ++jj) {
0247 if (jj == softest)
0248 continue;
0249 if ((subjets[ii] + subjets[jj]).m() < min_mass) {
0250 min_mass = (subjets[ii] + subjets[jj]).m();
0251 i = ii;
0252 j = jj;
0253 }
0254 }
0255 }
0256 }
0257
0258 FASTJET_END_NAMESPACE
0259
0260 #endif