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
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
|
//FJSTARTHEADER
// $Id$
//
// Copyright (c) 2005-2021, 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. They are described in the original FastJet paper,
// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
// FastJet as part of work towards a scientific publication, please
// quote the version you use and include a citation to the manual and
// optionally also to hep-ph/0512210.
//
// 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/>.
//----------------------------------------------------------------------
//FJENDHEADER
#include "fastjet/CDFJetCluPlugin.hh"
#include "fastjet/ClusterSequence.hh"
#include <sstream>
#include <cassert>
// CDF stuff
#include "JetCluAlgorithm.hh"
#include "PhysicsTower.hh"
#include "Cluster.hh"
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
using namespace std;
using namespace cdf;
thread_safety_helpers::FirstTimeTrue CDFJetCluPlugin::_first_time;
string CDFJetCluPlugin::description () const {
ostringstream desc;
desc << "CDF JetClu jet algorithm with "
<< "seed_threshold = " << seed_threshold () << ", "
<< "cone_radius = " << cone_radius () << ", "
<< "adjacency_cut = " << adjacency_cut () << ", "
<< "max_iterations = " << max_iterations () << ", "
<< "iratch = " << iratch () << ", "
<< "overlap_threshold = " << overlap_threshold () ;
return desc.str();
}
void CDFJetCluPlugin::run_clustering(ClusterSequence & clust_seq) const {
// print a banner if we run this for the first time
_print_banner(clust_seq.fastjet_banner_stream());
// create the physics towers needed by the CDF code
vector<PhysicsTower> towers;
towers.reserve(clust_seq.jets().size());
// create a map to identify jets (actually just the input particles)...
//map<double,int> jetmap;
for (unsigned i = 0; i < clust_seq.jets().size(); i++) {
PseudoJet particle(clust_seq.jets()[i]);
//_insert_unique(particle, jetmap);
LorentzVector fourvect(particle.px(), particle.py(),
particle.pz(), particle.E());
PhysicsTower tower(fourvect);
// add tracking information for later
tower.fjindex = i;
towers.push_back(tower);
}
// prepare the CDF algorithm
JetCluAlgorithm j(seed_threshold(), cone_radius(), adjacency_cut(),
max_iterations(), iratch(), overlap_threshold());
// run the CDF algorithm
std::vector<Cluster> jets;
j.run(towers,jets);
// now transfer the jets back into our own structure -- we will
// mimic the cone code with a sequential recombination sequence in
// which the jets are built up by adding one particle at a time
// NB: with g++-4.0, the reverse iterator code gave problems, so switch
// to indices instead
//for(vector<Cluster>::const_reverse_iterator jetIter = jets.rbegin();
// jetIter != jets.rend(); jetIter++) {
// const vector<PhysicsTower> & tower_list = jetIter->towerList;
// int jet_k = jetmap[tower_list[0].fourVector.E];
//
// int ntow = int(jetIter->towerList.size());
for(int iCDFjets = jets.size()-1; iCDFjets >= 0; iCDFjets--) {
const vector<PhysicsTower> & tower_list = jets[iCDFjets].towerList;
int ntow = int(tower_list.size());
// 2008-09-04: sort the towerList (according to fjindex) so as
// to have a consistent order for particles in jet
// (necessary because addition of ultra-soft particles
// sometimes often modifies the order, while maintaining
// the same overall set)
vector<int> jc_indices(ntow);
vector<double> fj_indices(ntow); // use double: benefit from existing routine
for (int itow = 0; itow < ntow; itow++) {
jc_indices[itow] = itow;
fj_indices[itow] = tower_list[itow].fjindex;
}
sort_indices(jc_indices, fj_indices);
int jet_k = tower_list[jc_indices[0]].fjindex;
for (int itow = 1; itow < ntow; itow++) {
if (tower_list[jc_indices[itow]].Et() > 1e-50) {
}
int jet_i = jet_k;
// retrieve our index for the jet
int jet_j;
jet_j = tower_list[jc_indices[itow]].fjindex;
// safety check
assert (jet_j >= 0 && jet_j < int(towers.size()));
// do a fake recombination step with dij=0
double dij = 0.0;
// JetClu does E-scheme recombination so we can stick with the
// simple option
clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k);
}
// NB: put a sensible looking d_iB just to be nice...
double d_iB = clust_seq.jets()[jet_k].perp2();
clust_seq.plugin_record_iB_recombination(jet_k, d_iB);
}
// following code is for testing only
//cout << endl;
//for(vector<Cluster>::const_iterator jetIter = jets.begin();
// jetIter != jets.end(); jetIter++) {
// cout << jetIter->fourVector.pt() << " " << jetIter->fourVector.y() << endl;
//}
//cout << "-----------------------------------------------------\n";
//vector<PseudoJet> ourjets(clust_seq.inclusive_jets());
//for (vector<PseudoJet>::const_iterator ourjet = ourjets.begin();
// ourjet != ourjets.end(); ourjet++) {
// cout << ourjet->perp() << " " << ourjet->rap() << endl;
//}
//cout << endl;
}
// print a banner for reference to the 3rd-party code
void CDFJetCluPlugin::_print_banner(ostream *ostr) const{
if (! _first_time()) return;
// make sure the user has not set the banner stream to NULL
if (!ostr) return;
(*ostr) << "#-------------------------------------------------------------------------" << endl;
(*ostr) << "# You are running the CDF JetClu plugin for FastJet " << endl;
(*ostr) << "# This is based on an implementation provided by Joey Huston. " << endl;
(*ostr) << "# If you use this plugin, please cite " << endl;
(*ostr) << "# F. Abe et al. [CDF Collaboration], Phys. Rev. D 45 (1992) 1448. " << endl;
(*ostr) << "# in addition to the usual FastJet reference. " << endl;
(*ostr) << "#-------------------------------------------------------------------------" << endl;
// make sure we really have the output done.
ostr->flush();
}
FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
|
