FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
CASubJetTagger.cc
1 //FJSTARTHEADER
2 // $Id: CASubJetTagger.cc 3433 2014-07-23 08:17:03Z salam $
3 //
4 // Copyright (c) 2005-2014, 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/CASubJetTagger.hh>
32 #include <fastjet/ClusterSequence.hh>
33 
34 #include <algorithm>
35 #include <cmath>
36 #include <sstream>
37 
38 using namespace std;
39 
40 FASTJET_BEGIN_NAMESPACE
41 
42 
43 LimitedWarning CASubJetTagger::_non_ca_warnings;
44 
45 // the tagger's description
46 //----------------------------------------------------------------------
47 string CASubJetTagger::description() const{
48  ostringstream oss;
49  oss << "CASubJetTagger with z_threshold=" << _z_threshold ;
50  if (_absolute_z_cut) oss << " (defined wrt original jet)";
51  oss << " and scale choice ";
52  switch (_scale_choice) {
53  case kt2_distance: oss << "kt2_distance"; break;
54  case jade_distance: oss << "jade_distance"; break;
55  case jade2_distance: oss << "jade2_distance"; break;
56  case plain_distance: oss << "plain_distance"; break;
57  case mass_drop_distance: oss << "mass_drop_distance"; break;
58  case dot_product_distance: oss << "dot_product_distance"; break;
59  default:
60  throw Error("unrecognized scale choice");
61  }
62 
63  return oss.str();
64 }
65 
66 // run the tagger on the given cs/jet
67 // returns the tagged PseudoJet if successful, 0 otherwise
68 //----------------------------------------------------------------------
69 PseudoJet CASubJetTagger::result(const PseudoJet & jet) const{
70  // make sure that the jet results from a Cambridge/Aachen clustering
72  _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");
73 
74  // recurse in the jet to find the max distance
75  JetAux aux;
76  aux.jet = PseudoJet();
77  aux.aux_distance = -numeric_limits<double>::max();
78  aux.delta_r = 0.0;
79  aux.z = 1.0;
80  _recurse_through_jet(jet, aux, jet); // last arg remains original jet
81 
82  // create the result and its associated structure
83  PseudoJet result_local = aux.jet;
84 
85  // the tagger is considered to have failed if aux has never been set
86  // (in which case it will not have parents).
87  if (result_local == PseudoJet()) return result_local;
88 
89  // otherwise sort out the structure
90  CASubJetTaggerStructure * s = new CASubJetTaggerStructure(result_local);
91 // s->_original_jet = jet;
92  s->_scale_choice = _scale_choice;
93  s->_distance = aux.aux_distance;
94  s->_absolute_z = _absolute_z_cut;
95  s->_z = aux.z;
96 
98 
99  return result_local;
100 }
101 
102 
103 ///----------------------------------------------------------------------
104 /// work through the jet, establishing a distance at each branching
105 inline void CASubJetTagger::_recurse_through_jet(const PseudoJet & jet, JetAux &aux, const PseudoJet & original_jet) const {
106 
107  PseudoJet parent1, parent2;
108  if (! jet.has_parents(parent1, parent2)) return;
109 
110  /// make sure the objects are not _too_ close together
111  if (parent1.squared_distance(parent2) < _dr2_min) return;
112 
113  // distance
114  double dist=0.0;
115  switch (_scale_choice) {
116  case kt2_distance:
117  // a standard (LI) kt distance
118  dist = parent1.kt_distance(parent2);
119  break;
120  case jade_distance:
121  // something a bit like a mass: pti ptj Delta R_ij^2
122  dist = parent1.perp()*parent2.perp()*parent1.squared_distance(parent2);
123  break;
124  case jade2_distance:
125  // something a bit like a mass*deltaR^2: pti ptj Delta R_ij^4
126  dist = parent1.perp()*parent2.perp()*pow(parent1.squared_distance(parent2),2);
127  break;
128  case plain_distance:
129  // Delta R_ij^2
130  dist = parent1.squared_distance(parent2);
131  break;
132  case mass_drop_distance:
133  // Delta R_ij^2
134  dist = jet.m() - std::max(parent1.m(),parent2.m());
135  break;
136  case dot_product_distance:
137  // parent1 . parent2
138  // ( = jet.m2() - parent1.m2() - parent2.m() in a
139  // 4-vector recombination scheme)
140  dist = dot_product(parent1, parent2);
141  break;
142  default:
143  throw Error("unrecognized scale choice");
144  }
145 
146  // check the z cut
147  bool zcut1 = true;
148  bool zcut2 = true;
149  double z2 = 0.0;
150 
151  // not very efficient -- sort out later
152  if (parent1.perp2() < parent2.perp2()) std::swap(parent1,parent2);
153 
154  if (_absolute_z_cut) {
155  z2 = parent2.perp() / original_jet.perp();
156  zcut1 = parent1.perp() / original_jet.perp() >= _z_threshold;
157  } else {
158  z2 = parent2.perp()/(parent1.perp()+parent2.perp());
159  }
160  zcut2 = z2 >= _z_threshold;
161 
162  if (zcut1 && zcut2){
163  if (dist > aux.aux_distance){
164  aux.jet = jet;
165  aux.aux_distance = dist;
166  aux.delta_r = sqrt(parent1.squared_distance(parent2));
167  aux.z = z2; // the softest
168  }
169  }
170 
171  if (zcut1) _recurse_through_jet(parent1, aux, original_jet);
172  if (zcut2) _recurse_through_jet(parent2, aux, original_jet);
173 }
174 
175 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:457
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 'parent...
Definition: PseudoJet.cc:547
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:47
double squared_distance(const PseudoJet &other) const
returns squared cylinder (rap-phi) distance between this jet and another
Definition: PseudoJet.hh:186
double kt_distance(const PseudoJet &other) const
returns kt distance (R=1) between this jet and another
Definition: PseudoJet.cc:378
an implementation of C++0x shared pointers (or boost's)
Definition: SharedPtr.hh:116
void set_structure_shared_ptr(const SharedPtr< PseudoJetStructureBase > &structure)
set the associated structure
Definition: PseudoJet.cc:464
double perp() const
returns the scalar transverse momentum
Definition: PseudoJet.hh:143
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:920
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
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:141
CASubJetTagger::ScaleChoice _scale_choice
the user scale choice
double _distance
the maximal distance associated with the result