FastJet  3.0.6
CASubJetTagger.cc
1 //STARTHEADER
2 // $Id: CASubJetTagger.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/CASubJetTagger.hh>
30 #include <fastjet/ClusterSequence.hh>
31 
32 #include <algorithm>
33 #include <cmath>
34 #include <sstream>
35 
36 using namespace std;
37 
38 FASTJET_BEGIN_NAMESPACE
39 
40 
41 LimitedWarning CASubJetTagger::_non_ca_warnings;
42 
43 // the tagger's description
44 //----------------------------------------------------------------------
45 string CASubJetTagger::description() const{
46  ostringstream oss;
47  oss << "CASubJetTagger with z_threshold=" << _z_threshold ;
48  if (_absolute_z_cut) oss << " (defined wrt original jet)";
49  oss << " and scale choice ";
50  switch (_scale_choice) {
51  case kt2_distance: oss << "kt2_distance"; break;
52  case jade_distance: oss << "jade_distance"; break;
53  case jade2_distance: oss << "jade2_distance"; break;
54  case plain_distance: oss << "plain_distance"; break;
55  case mass_drop_distance: oss << "mass_drop_distance"; break;
56  case dot_product_distance: oss << "dot_product_distance"; break;
57  default:
58  throw Error("unrecognized scale choice");
59  }
60 
61  return oss.str();
62 }
63 
64 // run the tagger on the given cs/jet
65 // returns the tagged PseudoJet if successful, 0 otherwise
66 //----------------------------------------------------------------------
67 PseudoJet CASubJetTagger::result(const PseudoJet & jet) const{
68  // make sure that the jet results from a Cambridge/Aachen clustering
70  _non_ca_warnings.warn("CASubJetTagger should only be applied on jets from a Cambridge/Aachen clustering; use it with other algorithms at your own risk");
71 
72  // recurse in the jet to find the max distance
73  JetAux aux;
74  aux.jet = PseudoJet();
75  aux.aux_distance = -numeric_limits<double>::max();
76  aux.delta_r = 0.0;
77  aux.z = 1.0;
78  _recurse_through_jet(jet, aux, jet); // last arg remains original jet
79 
80  // create the result and its associated structure
81  PseudoJet result_local = aux.jet;
82 
83  // the tagger is considered to have failed if aux has never been set
84  // (in which case it will not have parents).
85  if (result_local == PseudoJet()) return result_local;
86 
87  // otherwise sort out the structure
88  CASubJetTaggerStructure * s = new CASubJetTaggerStructure(result_local);
89 // s->_original_jet = jet;
90  s->_scale_choice = _scale_choice;
91  s->_distance = aux.aux_distance;
92  s->_absolute_z = _absolute_z_cut;
93  s->_z = aux.z;
94 
96 
97  return result_local;
98 }
99 
100 
101 ///----------------------------------------------------------------------
102 /// work through the jet, establishing a distance at each branching
103 inline void CASubJetTagger::_recurse_through_jet(const PseudoJet & jet, JetAux &aux, const PseudoJet & original_jet) const {
104 
105  PseudoJet parent1, parent2;
106  if (! jet.has_parents(parent1, parent2)) return;
107 
108  /// make sure the objects are not _too_ close together
109  if (parent1.squared_distance(parent2) < _dr2_min) return;
110 
111  // distance
112  double dist=0.0;
113  switch (_scale_choice) {
114  case kt2_distance:
115  // a standard (LI) kt distance
116  dist = parent1.kt_distance(parent2);
117  break;
118  case jade_distance:
119  // something a bit like a mass: pti ptj Delta R_ij^2
120  dist = parent1.perp()*parent2.perp()*parent1.squared_distance(parent2);
121  break;
122  case jade2_distance:
123  // something a bit like a mass*deltaR^2: pti ptj Delta R_ij^4
124  dist = parent1.perp()*parent2.perp()*pow(parent1.squared_distance(parent2),2);
125  break;
126  case plain_distance:
127  // Delta R_ij^2
128  dist = parent1.squared_distance(parent2);
129  break;
130  case mass_drop_distance:
131  // Delta R_ij^2
132  dist = jet.m() - std::max(parent1.m(),parent2.m());
133  break;
134  case dot_product_distance:
135  // parent1 . parent2
136  // ( = jet.m2() - parent1.m2() - parent2.m() in a
137  // 4-vector recombination scheme)
138  dist = dot_product(parent1, parent2);
139  break;
140  default:
141  throw Error("unrecognized scale choice");
142  }
143 
144  // check the z cut
145  bool zcut1 = true;
146  bool zcut2 = true;
147  double z2 = 0.0;
148 
149  // not very efficient -- sort out later
150  if (parent1.perp2() < parent2.perp2()) std::swap(parent1,parent2);
151 
152  if (_absolute_z_cut) {
153  z2 = parent2.perp() / original_jet.perp();
154  zcut1 = parent1.perp() / original_jet.perp() >= _z_threshold;
155  } else {
156  z2 = parent2.perp()/(parent1.perp()+parent2.perp());
157  }
158  zcut2 = z2 >= _z_threshold;
159 
160  if (zcut1 && zcut2){
161  if (dist > aux.aux_distance){
162  aux.jet = jet;
163  aux.aux_distance = dist;
164  aux.delta_r = sqrt(parent1.squared_distance(parent2));
165  aux.z = z2; // the softest
166  }
167  }
168 
169  if (zcut1) _recurse_through_jet(parent1, aux, original_jet);
170  if (zcut2) _recurse_through_jet(parent2, aux, original_jet);
171 }
172 
173 FASTJET_END_NAMESPACE
double _z
the transverse momentum fraction
JetAlgorithm jet_algorithm() const
return information about the definition...
the structure returned by a CASubJetTagger
class that contains the result internally
const ClusterSequence * validated_cs() const
shorthand for validated_cluster_sequence()
Definition: PseudoJet.cc:443
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:533
const JetDefinition & jet_def() const
return a reference to the jet definition
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:41
double squared_distance(const PseudoJet &other) const
returns squared cylinder (rap-phi) distance between this jet and another
Definition: PseudoJet.hh:184
double kt_distance(const PseudoJet &other) const
returns kt distance (R=1) between this jet and another
Definition: PseudoJet.cc:364
an implementation of C++0x shared pointers (or boost&#39;s)
Definition: SharedPtr.hh:114
void set_structure_shared_ptr(const SharedPtr< PseudoJetStructureBase > &structure)
set the associated structure
Definition: PseudoJet.cc:450
double perp() const
returns the scalar transverse momentum
Definition: PseudoJet.hh:141
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
double m() const
returns the invariant mass (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP) ...
Definition: PseudoJet.hh:909
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:65
bool _absolute_z
whether z is computed wrt to the original jet or not
double perp2() const
returns the squared transverse momentum
Definition: PseudoJet.hh:139
CASubJetTagger::ScaleChoice _scale_choice
the user scale choice
double _distance
the maximal distance associated with the result