FastJet  3.3.3
JHTopTagger.cc
1 //FJSTARTHEADER
2 // $Id: JHTopTagger.cc 4420 2019-11-29 09:28:20Z soyez $
3 //
4 // Copyright (c) 2005-2019, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
5 //
6 //----------------------------------------------------------------------
7 // This file is part of FastJet.
8 //
9 // FastJet is free software; you can redistribute it and/or modify
10 // it under the terms of the GNU General Public License as published by
11 // the Free Software Foundation; either version 2 of the License, or
12 // (at your option) any later version.
13 //
14 // The algorithms that underlie FastJet have required considerable
15 // development. They are described in the original FastJet paper,
16 // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
17 // FastJet as part of work towards a scientific publication, please
18 // quote the version you use and include a citation to the manual and
19 // optionally also to hep-ph/0512210.
20 //
21 // FastJet is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
28 //----------------------------------------------------------------------
29 //FJENDHEADER
30 
31 #include <fastjet/tools/JHTopTagger.hh>
32 #include <fastjet/Error.hh>
33 #include <fastjet/JetDefinition.hh>
34 #include <fastjet/ClusterSequence.hh>
35 #include <sstream>
36 #include <limits>
37 
38 FASTJET_BEGIN_NAMESPACE
39 
40 using namespace std;
41 
42 //----------------------------------------------------------------------
43 // JHTopTagger class implementation
44 //----------------------------------------------------------------------
45 
46 LimitedWarning JHTopTagger::_warnings_nonca;
47 
48 //------------------------------------------------------------------------
49 // description of the tagger
50 string JHTopTagger::description() const{
51  ostringstream oss;
52  oss << "JHTopTagger with delta_p=" << _delta_p << ", delta_r=" << _delta_r
53  << ", cos_theta_W_max=" << _cos_theta_W_max
54  << " and mW = " << _mW;
55  oss << description_of_selectors();
56  return oss.str();
57 }
58 
59 //------------------------------------------------------------------------
60 // returns the tagged PseudoJet if successful, 0 otherwise
61 // - jet the PseudoJet to tag
62 PseudoJet JHTopTagger::result(const PseudoJet & jet) const{
63  // make sure that there is a "regular" cluster sequence associated
64  // with the jet. Note that we also check it is valid (to avoid a
65  // more criptic error later on)
66  if (!jet.has_valid_cluster_sequence()){
67  throw Error("JHTopTagger can only be applied on jets having an associated (and valid) ClusterSequence");
68  }
69 
70  // warn if the jet has not been clustered with a Cambridge/Aachen
71  // algorithm
73  _warnings_nonca.warn("JHTopTagger should only be applied on jets from a Cambridge/Aachen clustering; use it with other algorithms at your own risk.");
74 
75 
76  // do the first splitting
77  vector<PseudoJet> split0 = _split_once(jet, jet);
78  if (split0.size() == 0) return PseudoJet();
79 
80  // now try a second splitting on each of the resulting objects
81  vector<PseudoJet> subjets;
82  for (unsigned i = 0; i < 2; i++) {
83  vector<PseudoJet> split1 = _split_once(split0[i], jet);
84  if (split1.size() > 0) {
85  subjets.push_back(split1[0]);
86  subjets.push_back(split1[1]);
87  } else {
88  subjets.push_back(split0[i]);
89  }
90  }
91 
92  // make sure things make sense
93  if (subjets.size() < 3) return PseudoJet();
94 
95  // now find the pair of objects closest in mass to the W
96  double dmW_min = numeric_limits<double>::max();
97  int ii=-1, jj=-1;
98  for (unsigned i = 0 ; i < subjets.size()-1; i++) {
99  for (unsigned j = i+1 ; j < subjets.size(); j++) {
100  double dmW = abs(_mW - (subjets[i]+subjets[j]).m());
101  if (dmW < dmW_min) {
102  dmW_min = dmW; ii = i; jj = j;
103  }
104  }
105  }
106 
107  // order the subjets in the following order:
108  // - hardest of the W subjets
109  // - softest of the W subjets
110  // - hardest of the remaining subjets
111  // - softest of the remaining subjets (if any)
112  if (ii>0) std::swap(subjets[ii], subjets[0]);
113  if (jj>1) std::swap(subjets[jj], subjets[1]);
114  if (subjets[0].perp2() < subjets[1].perp2()) std::swap(subjets[0], subjets[1]);
115  if ((subjets.size()>3) && (subjets[2].perp2() < subjets[3].perp2()))
116  std::swap(subjets[2], subjets[3]);
117 
118  // create the result and its structure
119  const JetDefinition::Recombiner *rec
121 
122  PseudoJet W = join(subjets[0], subjets[1], *rec);
123  PseudoJet non_W;
124  if (subjets.size()>3) {
125  non_W = join(subjets[2], subjets[3], *rec);
126  } else {
127  non_W = join(subjets[2], *rec);
128  }
129  PseudoJet result_local = join<JHTopTaggerStructure>(W, non_W, *rec);
131  s->_cos_theta_w = _cos_theta_W(result_local);
132 
133  // if the polarisation angle does not pass the cut, consider that
134  // the tagging has failed
135  //
136  // Note that we could perhaps ensure this cut before constructing
137  // the result structure but this has the advantage that the top
138  // 4-vector is already available and does not have to de re-computed
139  if (s->_cos_theta_w >= _cos_theta_W_max ||
140  ! _top_selector.pass(result_local) || ! _W_selector.pass(W)
141  ) {
142  result_local *= 0.0;
143  }
144 
145  return result_local;
146 
147  // // old version
148  // PseudoJet result = join<JHTopTaggerStructure>(subjets, *rec);
149  // JHTopTaggerStructure *s = (JHTopTaggerStructure*) result.structure_non_const_ptr();
150  // // s->_original_jet = jet;
151  // s->_W = join(subjets[0], subjets[1], *rec);
152  // if (subjets.size()>3)
153  // s->_non_W = join(subjets[2], subjets[3], *rec);
154  // else
155  // s->_non_W = join(subjets[2], *rec);
156  // s->_cos_theta_w = _cos_theta_W(result);
157  //
158  // // if the polarisation angle does not pass the cut, consider that
159  // // the tagging has failed
160  // //
161  // // Note that we could perhaps ensure this cut before constructing
162  // // the result structure but this has the advantage that the top
163  // // 4-vector is already available and does not have to de re-computed
164  // if (s->_cos_theta_w >= _cos_theta_W_max)
165  // return PseudoJet();
166  //
167  // return result;
168 }
169 
170 // runs the Johns Hopkins decomposition procedure
171 vector<PseudoJet> JHTopTagger::_split_once(const PseudoJet & jet_to_split,
172  const PseudoJet & reference_jet) const{
173  PseudoJet this_jet = jet_to_split;
174  PseudoJet p1, p2;
175  vector<PseudoJet> result_local;
176  while (this_jet.has_parents(p1, p2)) {
177  if (p2.perp2() > p1.perp2()) std::swap(p1,p2); // order with hardness
178  if (p1.perp() < _delta_p * reference_jet.perp()) break; // harder is too soft wrt original jet
179  if ( (abs(p2.rap()-p1.rap()) + abs(p2.delta_phi_to(p1))) < _delta_r) break; // distance is too small
180  if (p2.perp() < _delta_p * reference_jet.perp()) {
181  this_jet = p1; // softer is too soft wrt original, so ignore it
182  continue;
183  }
184  //result.push_back(this_jet);
185  result_local.push_back(p1);
186  result_local.push_back(p2);
187  break;
188  }
189  return result_local;
190 }
191 
192 
193 
194 
195 FASTJET_END_NAMESPACE
196 
const JetDefinition & jet_def() const
return a reference to the jet definition
JetAlgorithm jet_algorithm() const
return information about the definition...
virtual bool has_parents(PseudoJet &parent1, PseudoJet &parent2) const
check if it is the product of a recombination, in which case return the 2 parents through the &#39;parent...
Definition: PseudoJet.cc:551
const ClusterSequence * validated_cs() const
shorthand for validated_cluster_sequence()
Definition: PseudoJet.cc:461
PseudoJetStructureBase * structure_non_const_ptr()
return a non-const pointer to the structure (of type PseudoJetStructureBase*) associated with this Ps...
Definition: PseudoJet.cc:498
double perp() const
returns the scalar transverse momentum
Definition: PseudoJet.hh:145
const Recombiner * recombiner() const
returns a pointer to the currently defined recombiner.
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:47
double perp2() const
returns the squared transverse momentum
Definition: PseudoJet.hh:143
const ClusterSequence * associated_cluster_sequence() const
get a (const) pointer to the parent ClusterSequence (NULL if inexistent)
Definition: PseudoJet.cc:441
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:125
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
An abstract base class that will provide the recombination scheme facilities and/or allow a user to e...
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
bool has_valid_cluster_sequence() const
returns true if this PseudoJet has an associated and still valid(ated) ClusterSequence.
Definition: PseudoJet.cc:451
the structure returned by the JHTopTagger transformer.
Definition: JHTopTagger.hh:169
double delta_phi_to(const PseudoJet &other) const
returns other.phi() - this.phi(), constrained to be in range -pi .
Definition: PseudoJet.cc:405