1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
|
//STARTHEADER
// $Id: CASubJetTagger.cc 2689 2011-11-14 14:51:06Z soyez $
//
// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
// FastJet is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// The algorithms that underlie FastJet have required considerable
// development and are described in hep-ph/0512210. If you use
// FastJet as part of work towards a scientific publication, please
// include a citation to the FastJet paper.
//
// FastJet is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
//ENDHEADER
#include <fastjet/tools/CASubJetTagger.hh>
#include <fastjet/ClusterSequence.hh>
#include <algorithm>
#include <cmath>
#include <sstream>
using namespace std;
FASTJET_BEGIN_NAMESPACE
LimitedWarning CASubJetTagger::_non_ca_warnings;
// the tagger's description
//----------------------------------------------------------------------
string CASubJetTagger::description() const{
ostringstream oss;
oss << "CASubJetTagger with z_threshold=" << _z_threshold ;
if (_absolute_z_cut) oss << " (defined wrt original jet)";
oss << " and scale choice ";
switch (_scale_choice) {
case kt2_distance: oss << "kt2_distance"; break;
case jade_distance: oss << "jade_distance"; break;
case jade2_distance: oss << "jade2_distance"; break;
case plain_distance: oss << "plain_distance"; break;
case mass_drop_distance: oss << "mass_drop_distance"; break;
case dot_product_distance: oss << "dot_product_distance"; break;
default:
throw Error("unrecognized scale choice");
}
return oss.str();
}
// run the tagger on the given cs/jet
// returns the tagged PseudoJet if successful, 0 otherwise
//----------------------------------------------------------------------
PseudoJet CASubJetTagger::result(const PseudoJet & jet) const{
// make sure that the jet results from a Cambridge/Aachen clustering
if (jet.validated_cs()->jet_def().jet_algorithm() != cambridge_algorithm)
_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");
// recurse in the jet to find the max distance
JetAux aux;
aux.jet = PseudoJet();
aux.aux_distance = -numeric_limits<double>::max();
aux.delta_r = 0.0;
aux.z = 1.0;
_recurse_through_jet(jet, aux, jet); // last arg remains original jet
// create the result and its associated structure
PseudoJet result_local = aux.jet;
// the tagger is considered to have failed if aux has never been set
// (in which case it will not have parents).
if (result_local == PseudoJet()) return result_local;
// otherwise sort out the structure
CASubJetTaggerStructure * s = new CASubJetTaggerStructure(result_local);
// s->_original_jet = jet;
s->_scale_choice = _scale_choice;
s->_distance = aux.aux_distance;
s->_absolute_z = _absolute_z_cut;
s->_z = aux.z;
result_local.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(s));
return result_local;
}
///----------------------------------------------------------------------
/// work through the jet, establishing a distance at each branching
inline void CASubJetTagger::_recurse_through_jet(const PseudoJet & jet, JetAux &aux, const PseudoJet & original_jet) const {
PseudoJet parent1, parent2;
if (! jet.has_parents(parent1, parent2)) return;
/// make sure the objects are not _too_ close together
if (parent1.squared_distance(parent2) < _dr2_min) return;
// distance
double dist=0.0;
switch (_scale_choice) {
case kt2_distance:
// a standard (LI) kt distance
dist = parent1.kt_distance(parent2);
break;
case jade_distance:
// something a bit like a mass: pti ptj Delta R_ij^2
dist = parent1.perp()*parent2.perp()*parent1.squared_distance(parent2);
break;
case jade2_distance:
// something a bit like a mass*deltaR^2: pti ptj Delta R_ij^4
dist = parent1.perp()*parent2.perp()*pow(parent1.squared_distance(parent2),2);
break;
case plain_distance:
// Delta R_ij^2
dist = parent1.squared_distance(parent2);
break;
case mass_drop_distance:
// Delta R_ij^2
dist = jet.m() - std::max(parent1.m(),parent2.m());
break;
case dot_product_distance:
// parent1 . parent2
// ( = jet.m2() - parent1.m2() - parent2.m() in a
// 4-vector recombination scheme)
dist = dot_product(parent1, parent2);
break;
default:
throw Error("unrecognized scale choice");
}
// check the z cut
bool zcut1 = true;
bool zcut2 = true;
double z2 = 0.0;
// not very efficient -- sort out later
if (parent1.perp2() < parent2.perp2()) std::swap(parent1,parent2);
if (_absolute_z_cut) {
z2 = parent2.perp() / original_jet.perp();
zcut1 = parent1.perp() / original_jet.perp() >= _z_threshold;
} else {
z2 = parent2.perp()/(parent1.perp()+parent2.perp());
}
zcut2 = z2 >= _z_threshold;
if (zcut1 && zcut2){
if (dist > aux.aux_distance){
aux.jet = jet;
aux.aux_distance = dist;
aux.delta_r = sqrt(parent1.squared_distance(parent2));
aux.z = z2; // the softest
}
}
if (zcut1) _recurse_through_jet(parent1, aux, original_jet);
if (zcut2) _recurse_through_jet(parent2, aux, original_jet);
}
FASTJET_END_NAMESPACE
|
