FastJet  3.0.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
JHTopTagger.cc
1 //STARTHEADER
2 // $Id: JHTopTagger.cc 2689 2011-11-14 14:51:06Z soyez $
3 //
4 // Copyright (c) 2005-2011, 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 and are described in hep-ph/0512210. If you use
16 // FastJet as part of work towards a scientific publication, please
17 // include a citation to the FastJet paper.
18 //
19 // FastJet is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 // GNU General Public License for more details.
23 //
24 // You should have received a copy of the GNU General Public License
25 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
26 //----------------------------------------------------------------------
27 //ENDHEADER
28 
29 #include <fastjet/tools/JHTopTagger.hh>
30 #include <fastjet/Error.hh>
31 #include <fastjet/JetDefinition.hh>
32 #include <fastjet/ClusterSequence.hh>
33 #include <sstream>
34 #include <limits>
35 
36 FASTJET_BEGIN_NAMESPACE
37 
38 using namespace std;
39 
40 //----------------------------------------------------------------------
41 // JHTopTagger class implementation
42 //----------------------------------------------------------------------
43 
44 LimitedWarning JHTopTagger::_warnings_nonca;
45 
46 //------------------------------------------------------------------------
47 // description of the tagger
48 string JHTopTagger::description() const{
49  ostringstream oss;
50  oss << "JHTopTagger with delta_p=" << _delta_p << ", delta_r=" << _delta_r
51  << ", cos_theta_W_max=" << _cos_theta_W_max
52  << " and mW = " << _mW;
53  oss << description_of_selectors();
54  return oss.str();
55 }
56 
57 //------------------------------------------------------------------------
58 // returns the tagged PseudoJet if successful, 0 otherwise
59 // - jet the PseudoJet to tag
60 PseudoJet JHTopTagger::result(const PseudoJet & jet) const{
61  // make sure that there is a "regular" cluster sequence associated
62  // with the jet. Note that we also check it is valid (to avoid a
63  // more criptic error later on)
64  if (!jet.has_valid_cluster_sequence()){
65  throw Error("JHTopTagger can only be applied on jets having an associated (and valid) ClusterSequence");
66  }
67 
68  // warn if the jet has not been clustered with a Cambridge/Aachen
69  // algorithm
71  _warnings_nonca.warn("JHTopTagger should only be applied on jets from a Cambridge/Aachen clustering; use it with other algorithms at your own risk.");
72 
73 
74  // do the first splitting
75  vector<PseudoJet> split0 = _split_once(jet, jet);
76  if (split0.size() == 0) return PseudoJet();
77 
78  // now try a second splitting on each of the resulting objects
79  vector<PseudoJet> subjets;
80  for (unsigned i = 0; i < 2; i++) {
81  vector<PseudoJet> split1 = _split_once(split0[i], jet);
82  if (split1.size() > 0) {
83  subjets.push_back(split1[0]);
84  subjets.push_back(split1[1]);
85  } else {
86  subjets.push_back(split0[i]);
87  }
88  }
89 
90  // make sure things make sense
91  if (subjets.size() < 3) return PseudoJet();
92 
93  // now find the pair of objects closest in mass to the W
94  double dmW_min = numeric_limits<double>::max();
95  int ii=-1, jj=-1;
96  for (unsigned i = 0 ; i < subjets.size()-1; i++) {
97  for (unsigned j = i+1 ; j < subjets.size(); j++) {
98  double dmW = abs(_mW - (subjets[i]+subjets[j]).m());
99  if (dmW < dmW_min) {
100  dmW_min = dmW; ii = i; jj = j;
101  }
102  }
103  }
104 
105  // order the subjets in the following order:
106  // - hardest of the W subjets
107  // - softest of the W subjets
108  // - hardest of the remaining subjets
109  // - softest of the remaining subjets (if any)
110  if (ii>0) std::swap(subjets[ii], subjets[0]);
111  if (jj>1) std::swap(subjets[jj], subjets[1]);
112  if (subjets[0].perp2() < subjets[1].perp2()) std::swap(subjets[0], subjets[1]);
113  if ((subjets.size()>3) && (subjets[2].perp2() < subjets[3].perp2()))
114  std::swap(subjets[2], subjets[3]);
115 
116  // create the result and its structure
117  const JetDefinition::Recombiner *rec
119 
120  PseudoJet W = join(subjets[0], subjets[1], *rec);
121  PseudoJet non_W;
122  if (subjets.size()>3) {
123  non_W = join(subjets[2], subjets[3], *rec);
124  } else {
125  non_W = join(subjets[2], *rec);
126  }
127  PseudoJet result_local = join<JHTopTaggerStructure>(W, non_W, *rec);
129  s->_cos_theta_w = _cos_theta_W(result_local);
130 
131  // if the polarisation angle does not pass the cut, consider that
132  // the tagging has failed
133  //
134  // Note that we could perhaps ensure this cut before constructing
135  // the result structure but this has the advantage that the top
136  // 4-vector is already available and does not have to de re-computed
137  if (s->_cos_theta_w >= _cos_theta_W_max ||
138  ! _top_selector.pass(result_local) || ! _W_selector.pass(W)
139  ) {
140  result_local *= 0.0;
141  }
142 
143  return result_local;
144 
145  // // old version
146  // PseudoJet result = join<JHTopTaggerStructure>(subjets, *rec);
147  // JHTopTaggerStructure *s = (JHTopTaggerStructure*) result.structure_non_const_ptr();
148  // // s->_original_jet = jet;
149  // s->_W = join(subjets[0], subjets[1], *rec);
150  // if (subjets.size()>3)
151  // s->_non_W = join(subjets[2], subjets[3], *rec);
152  // else
153  // s->_non_W = join(subjets[2], *rec);
154  // s->_cos_theta_w = _cos_theta_W(result);
155  //
156  // // if the polarisation angle does not pass the cut, consider that
157  // // the tagging has failed
158  // //
159  // // Note that we could perhaps ensure this cut before constructing
160  // // the result structure but this has the advantage that the top
161  // // 4-vector is already available and does not have to de re-computed
162  // if (s->_cos_theta_w >= _cos_theta_W_max)
163  // return PseudoJet();
164  //
165  // return result;
166 }
167 
168 // runs the Johns Hopkins decomposition procedure
169 vector<PseudoJet> JHTopTagger::_split_once(const PseudoJet & jet_to_split,
170  const PseudoJet & reference_jet) const{
171  PseudoJet this_jet = jet_to_split;
172  PseudoJet p1, p2;
173  vector<PseudoJet> result_local;
174  while (this_jet.has_parents(p1, p2)) {
175  if (p2.perp2() > p1.perp2()) std::swap(p1,p2); // order with hardness
176  if (p1.perp() < _delta_p * reference_jet.perp()) break; // harder is too soft wrt original jet
177  if ( (abs(p2.rap()-p1.rap()) + abs(p2.delta_phi_to(p1))) < _delta_r) break; // distance is too small
178  if (p2.perp() < _delta_p * reference_jet.perp()) {
179  this_jet = p1; // softer is too soft wrt original, so ignore it
180  continue;
181  }
182  //result.push_back(this_jet);
183  result_local.push_back(p1);
184  result_local.push_back(p2);
185  break;
186  }
187  return result_local;
188 }
189 
190 
191 
192 
193 FASTJET_END_NAMESPACE
194