CMS 3D CMS Logo

CMSTopTagger.h
Go to the documentation of this file.
1 // Original JHTopTagger code 2011 Matteo Cacciari, Gregory Soyez,
2 // and Gavin Salam.
3 // Modifications to make it CMSTopTagger 2011 Christopher
4 // Vermilion.
5 //
6 //----------------------------------------------------------------------
7 // This file is free software; you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation; either version 3 of the License, or
10 // (at your option) any later version.
11 //
12 // This file is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
16 //
17 // The GNU General Public License is available at
18 // http://www.gnu.org/licenses/gpl.html or you can write to the Free Software
19 // Foundation, Inc.:
20 // 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 //----------------------------------------------------------------------
22 
23 #ifndef __CMSTOPTAGGER_HH__
24 #define __CMSTOPTAGGER_HH__
25 
26 #include <fastjet/tools/JHTopTagger.hh>
27 #include <fastjet/CompositeJetStructure.hh>
28 #include <fastjet/LimitedWarning.hh>
29 #include <fastjet/Error.hh>
30 #include <fastjet/JetDefinition.hh>
31 #include <fastjet/ClusterSequence.hh>
32 
33 #include <sstream>
34 #include <limits>
35 
36 FASTJET_BEGIN_NAMESPACE
37 
52 
53 
54 class CMSTopTagger : public TopTaggerBase {
55 public:
66  CMSTopTagger(double delta_p=0.05, double delta_r=0.4, double A=0.0004);
67 
69  std::string description() const override;
70 
75  PseudoJet result(const PseudoJet & jet) const override;
76 
77  // the type of the associated structure
79 
80 protected:
82  std::vector<PseudoJet> _split_once(const PseudoJet & jet_to_split,
83  const PseudoJet & reference_jet) const;
84 
87  void _find_min_mass(const std::vector<PseudoJet>& subjets, int& i, int& j) const;
88 
89  double _delta_p, _delta_r, _A;
90  mutable LimitedWarning _warnings_nonca;
91 };
92 
95 class CMSTopTaggerStructure : public JHTopTaggerStructure {
96 public:
97  CMSTopTaggerStructure(const std::vector<PseudoJet>& pieces,
98  const JetDefinition::Recombiner *recombiner = nullptr)
99  : JHTopTaggerStructure(pieces, recombiner) {}
100 
101 protected:
102  friend class CMSTopTagger;
103 };
104 
105 
106 
107 
108 //----------------------------------------------------------------------
109 inline CMSTopTagger::CMSTopTagger(double delta_p, double delta_r, double A)
110  : _delta_p(delta_p), _delta_r(delta_r), _A(A) {}
111 
112 //------------------------------------------------------------------------
113 // description of the tagger
115  std::ostringstream oss;
116  oss << "CMSTopTagger with delta_p=" << _delta_p << ", delta_r=" << _delta_r
117  << ", and A=" << _A;
118  oss << description_of_selectors();
119  return oss.str();
120 }
121 
122 
123 //------------------------------------------------------------------------
124 // returns the tagged PseudoJet if successful, 0 otherwise
125 // - jet the PseudoJet to tag
126 inline PseudoJet CMSTopTagger::result(const PseudoJet & jet) const{
127  // make sure that there is a "regular" cluster sequence associated
128  // with the jet. Note that we also check it is valid (to avoid a
129  // more criptic error later on)
130  if (!jet.has_valid_cluster_sequence()){
131  throw Error("CMSTopTagger can only be applied on jets having an associated (and valid) ClusterSequence");
132  }
133 
134  // warn if the jet has not been clustered with a Cambridge/Aachen
135  // algorithm
136  if (jet.validated_cs()->jet_def().jet_algorithm() != cambridge_algorithm)
137  _warnings_nonca.warn("CMSTopTagger should only be applied on jets from a Cambridge/Aachen clustering; use it with other algorithms at your own risk.");
138 
139 
140  // do the first splitting
141  std::vector<PseudoJet> split0 = _split_once(jet, jet);
142  if (split0.empty()) return PseudoJet();
143 
144  // now try a second splitting on each of the resulting objects
145  std::vector<PseudoJet> subjets;
146  for (unsigned i = 0; i < 2; i++) {
147  std::vector<PseudoJet> split1 = _split_once(split0[i], jet);
148  if (!split1.empty()) {
149  subjets.push_back(split1[0]);
150  subjets.push_back(split1[1]);
151  } else {
152  subjets.push_back(split0[i]);
153  }
154  }
155 
156  // make sure things make sense
157  if (subjets.size() < 3) return PseudoJet();
158 
159  // now find the pair of subjets with minimum mass (only taking hardest three)
160  int ii=-1, jj=-1;
161  _find_min_mass(subjets, ii, jj);
162 
163  // order the subjets in the following order:
164  // - hardest of the W subjets
165  // - softest of the W subjets
166  // - hardest of the remaining subjets
167  // - softest of the remaining subjets (if any)
168  if (ii>0) std::swap(subjets[ii], subjets[0]);
169  if (jj>1) std::swap(subjets[jj], subjets[1]);
170  if (subjets[0].perp2() < subjets[1].perp2()) std::swap(subjets[0], subjets[1]);
171  if ((subjets.size()>3) && (subjets[2].perp2() < subjets[3].perp2()))
172  std::swap(subjets[2], subjets[3]);
173 
174  // create the result and its structure
175  const JetDefinition::Recombiner *rec
176  = jet.associated_cluster_sequence()->jet_def().recombiner();
177 
178  PseudoJet W = join(subjets[0], subjets[1], *rec);
179  PseudoJet non_W;
180  if (subjets.size()>3) {
181  non_W = join(subjets[2], subjets[3], *rec);
182  } else {
183  non_W = join(subjets[2], *rec);
184  }
185  PseudoJet result = join<CMSTopTaggerStructure>(W, non_W, *rec);
186  CMSTopTaggerStructure *s = (CMSTopTaggerStructure*) result.structure_non_const_ptr();
187  s->_cos_theta_w = _cos_theta_W(result);
188 
189  // Note that we could perhaps ensure this cut before constructing
190  // the result structure but this has the advantage that the top
191  // 4-vector is already available and does not have to de re-computed
192  if (!_top_selector.pass(result) || ! _W_selector.pass(W)) {
193  result *= 0.0;
194  }
195 
196  result = join(subjets); //Added by J. Pilot to combine the (up to 4) subjets identified in the decomposition instead of just the W and non_W components
197 
198  return result;
199 }
200 
201 
202 // runs the Johns Hopkins decomposition procedure
203 inline std::vector<PseudoJet> CMSTopTagger::_split_once(const PseudoJet & jet_to_split,
204  const PseudoJet & reference_jet) const{
205  PseudoJet this_jet = jet_to_split;
206  PseudoJet p1, p2;
207  std::vector<PseudoJet> result;
208  while (this_jet.has_parents(p1, p2)) {
209  if (p2.perp2() > p1.perp2()) std::swap(p1,p2); // order with hardness
210  if (p1.perp() < _delta_p * reference_jet.perp()) break; // harder is too soft wrt original jet
211  double DR = p1.delta_R(p2);
212  if (DR < _delta_r - _A * this_jet.perp()) break; // distance is too small
213  if (p2.perp() < _delta_p * reference_jet.perp()) {
214  this_jet = p1; // softer is too soft wrt original, so ignore it
215  continue;
216  }
217  //result.push_back(this_jet);
218  result.push_back(p1);
219  result.push_back(p2);
220  break;
221  }
222  return result;
223 }
224 
225 
226 // find the indices corresponding to the minimum mass pairing in subjets
227 inline void CMSTopTagger::_find_min_mass(const std::vector<PseudoJet>& subjets, int& i, int& j) const{
228  assert(subjets.size() > 1 && subjets.size() < 5);
229 
230  // if four subjets found, only consider three hardest
231  unsigned softest = 5; // invalid value
232  if (subjets.size() == 4) {
234  for (unsigned ii = 0; ii < subjets.size(); ++ii) {
235  if (subjets[ii].perp() < min_pt) {
236  min_pt = subjets[ii].perp();
237  softest = ii;
238  }
239  }
240  }
241 
242  double min_mass = std::numeric_limits<double>::max();
243  for (unsigned ii = 0; ii+1 < subjets.size(); ++ii) { // don't do size()-1: (unsigned(0)-1 != -1) !!
244  if (ii == softest) continue;
245  for (unsigned jj = ii + 1; jj < subjets.size(); ++jj) {
246  if (jj == softest) continue;
247  if ((subjets[ii]+subjets[jj]).m() < min_mass) {
248  min_mass = (subjets[ii]+subjets[jj]).m();
249  i = ii;
250  j = jj;
251  }
252  }
253  }
254 }
255 
256 
257 FASTJET_END_NAMESPACE
258 
259 #endif // __CMSTOPTAGGER_HH__
edm::ErrorSummaryEntry Error
LimitedWarning _warnings_nonca
Definition: CMSTopTagger.h:90
double _delta_p
Definition: CMSTopTagger.h:89
std::string description() const override
returns a textual description of the tagger
Definition: CMSTopTagger.h:114
void _find_min_mass(const std::vector< PseudoJet > &subjets, int &i, int &j) const
Definition: CMSTopTagger.h:227
CMSTopTaggerStructure StructureType
Definition: CMSTopTagger.h:78
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
CMSTopTagger(double delta_p=0.05, double delta_r=0.4, double A=0.0004)
Definition: CMSTopTagger.h:109
double delta_r(const Fourvec &a, const Fourvec &b)
Find the distance between two four-vectors in the two-dimensional space .
Definition: fourvec.cc:249
CMSTopTaggerStructure(const std::vector< PseudoJet > &pieces, const JetDefinition::Recombiner *recombiner=nullptr)
Definition: CMSTopTagger.h:97
double p2[4]
Definition: TauolaWrapper.h:90
ii
Definition: cuy.py:590
static std::string join(char **cmd)
Definition: RemoteFile.cc:18
T perp2() const
Squared magnitude of transverse component.
PseudoJet result(const PseudoJet &jet) const override
Definition: CMSTopTagger.h:126
T perp() const
Magnitude of transverse component.
double p1[4]
Definition: TauolaWrapper.h:89
std::vector< PseudoJet > _split_once(const PseudoJet &jet_to_split, const PseudoJet &reference_jet) const
runs the Johns Hopkins decomposition procedure
Definition: CMSTopTagger.h:203
double _delta_r
Definition: CMSTopTagger.h:89