Fork me on GitHub

source: git/external/fastjet/tools/Pruner.cc@ 4586f8f

Last change on this file since 4586f8f was b7b836a, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

update FastJet library to 3.3.1 and FastJet Contrib library to 1.036

  • Property mode set to 100644
File size: 12.9 KB
Line 
1//FJSTARTHEADER
2// $Id: Pruner.cc 4354 2018-04-22 07:12:37Z salam $
3//
4// Copyright (c) 2005-2018, 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/Pruner.hh"
32#include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
33#include "fastjet/Selector.hh"
34#include <cassert>
35#include <algorithm>
36#include <sstream>
37#include <typeinfo>
38
39using namespace std;
40
41
42FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
43
44
45//----------------------------------------------------------------------
46// class Pruner
47//----------------------------------------------------------------------
48
49//----------------------------------------------------------------------
50// alternative (dynamic) ctor
51// \param jet_def the jet definition for the internal clustering
52// \param zcut_dyn dynamic pt-fraction cut in the pruning
53// \param Rcut_dyn dynamic angular distance cut in the pruning
54Pruner::Pruner(const JetDefinition &jet_def,
55 const FunctionOfPseudoJet<double> *zcut_dyn,
56 const FunctionOfPseudoJet<double> *Rcut_dyn)
57 : _jet_def(jet_def), _zcut(0), _Rcut_factor(0),
58 _zcut_dyn(zcut_dyn), _Rcut_dyn(Rcut_dyn), _get_recombiner_from_jet(false) {
59 assert(_zcut_dyn != 0 && _Rcut_dyn != 0);
60}
61
62//----------------------------------------------------------------------
63// action on a single jet
64PseudoJet Pruner::result(const PseudoJet &jet) const{
65 // pruning can only be applied to jets that have constituents
66 if (!jet.has_constituents()){
67 throw Error("Pruner: trying to apply the Pruner transformer to a jet that has no constituents");
68 }
69
70 // if the jet has area support and there are explicit ghosts, we can
71 // transfer that support to the internal re-clustering
72 bool do_areas = jet.has_area() && _check_explicit_ghosts(jet);
73
74 // build the pruning plugin
75 double Rcut = (_Rcut_dyn) ? (*_Rcut_dyn)(jet) : _Rcut_factor * 2.0*jet.m()/jet.perp();
76 double zcut = (_zcut_dyn) ? (*_zcut_dyn)(jet) : _zcut;
77 PruningPlugin * pruning_plugin;
78
79 // for some constructors, we get the recombiner from the
80 // input jet -- some acrobatics are needed
81 if (_get_recombiner_from_jet) {
82 JetDefinition jet_def = _jet_def;
83
84 // if all the pieces have a shared recombiner, we'll use that
85 // one. Otherwise, use the one from _jet_def as a fallback.
86 JetDefinition jet_def_for_recombiner;
87 if (_check_common_recombiner(jet, jet_def_for_recombiner)){
88 jet_def.set_recombiner(jet_def_for_recombiner);
89 }
90 pruning_plugin = new PruningPlugin(jet_def, zcut, Rcut);
91 } else {
92 pruning_plugin = new PruningPlugin(_jet_def, zcut, Rcut);
93 }
94
95 // now recluster the constituents of the jet with that plugin
96 JetDefinition internal_jet_def(pruning_plugin);
97 // flag the plugin for automatic deletion _before_ we make
98 // copies (so that as long as the copies are also present
99 // it doesn't get deleted).
100 internal_jet_def.delete_plugin_when_unused();
101
102 ClusterSequence * cs;
103 if (do_areas){
104 vector<PseudoJet> particles, ghosts;
105 SelectorIsPureGhost().sift(jet.constituents(), ghosts, particles);
106 // determine the ghost area from the 1st ghost (if none, any value
107 // will do, as the area will be 0 and subtraction will have
108 // no effect!)
109 double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
110 cs = new ClusterSequenceActiveAreaExplicitGhosts(particles, internal_jet_def,
111 ghosts, ghost_area);
112 } else {
113 cs = new ClusterSequence(jet.constituents(), internal_jet_def);
114 }
115
116 PseudoJet result_local = SelectorNHardest(1)(cs->inclusive_jets())[0];
117 PrunerStructure * s = new PrunerStructure(result_local);
118 s->_Rcut = Rcut;
119 s->_zcut = zcut;
120 result_local.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(s));
121
122 // make sure things remain persistent -- i.e. tell the jet definition
123 // and the cluster sequence that it is their responsibility to clean
124 // up memory once the "result" reaches the end of its life in the user's
125 // code. (The CS deletes itself when the result goes out of scope and
126 // that also triggers deletion of the plugin)
127 cs->delete_self_when_unused();
128
129 return result_local;
130}
131
132// check if the jet has explicit_ghosts (knowing that there is an
133// area support)
134bool Pruner::_check_explicit_ghosts(const PseudoJet &jet) const{
135 // if the jet comes from a Clustering check explicit ghosts in that
136 // clustering
137 if (jet.has_associated_cluster_sequence())
138 return jet.validated_csab()->has_explicit_ghosts();
139
140 // if the jet has pieces, recurse in the pieces
141 if (jet.has_pieces()){
142 vector<PseudoJet> pieces = jet.pieces();
143 for (unsigned int i=0;i<pieces.size(); i++)
144 if (!_check_explicit_ghosts(pieces[i])) return false;
145 // never returned false, so we're OK.
146 return true;
147 }
148
149 // return false for any other (unknown) structure
150 return false;
151}
152
153// see if there is a common recombiner among the pieces; if there is
154// return true and set jet_def_for_recombiner so that the recombiner
155// can be taken from that JetDefinition. Otherwise, return
156// false. 'assigned' is initially false; when true, each time we meet
157// a new jet definition, we'll check it shares the same recombiner as
158// jet_def_for_recombiner.
159bool Pruner::_check_common_recombiner(const PseudoJet &jet,
160 JetDefinition &jet_def_for_recombiner,
161 bool assigned) const{
162 if (jet.has_associated_cluster_sequence()){
163 // if the jet def for recombination has already been assigned, check if we have the same
164 if (assigned)
165 return jet.validated_cs()->jet_def().has_same_recombiner(jet_def_for_recombiner);
166
167 // otherwise, assign it.
168 jet_def_for_recombiner = jet.validated_cs()->jet_def();
169 assigned = true;
170 return true;
171 }
172
173 // if the jet has pieces, recurse in the pieces
174 if (jet.has_pieces()){
175 vector<PseudoJet> pieces = jet.pieces();
176 if (pieces.size() == 0) return false;
177 for (unsigned int i=0;i<pieces.size(); i++)
178 if (!_check_common_recombiner(pieces[i], jet_def_for_recombiner, assigned)) return false;
179 // never returned false, so we're OK.
180 return true;
181 }
182
183 // return false for any other (unknown) structure
184 return false;
185}
186
187
188// transformer description
189std::string Pruner::description() const{
190 ostringstream oss;
191 oss << "Pruner with jet_definition = (" << _jet_def.description() << ")";
192 if (_zcut_dyn) {
193 oss << ", dynamic zcut (" << _zcut_dyn->description() << ")"
194 << ", dynamic Rcut (" << _Rcut_dyn->description() << ")";
195 } else {
196 oss << ", zcut = " << _zcut
197 << ", Rcut_factor = " << _Rcut_factor;
198 }
199 return oss.str();
200}
201
202
203
204//----------------------------------------------------------------------
205// class PrunerStructure
206//----------------------------------------------------------------------
207
208// return the other jets that may have been found along with the
209// result of the pruning
210// The resulting vector is sorted in pt
211vector<PseudoJet> PrunerStructure::extra_jets() const{
212 return sorted_by_pt((!SelectorNHardest(1))(validated_cs()->inclusive_jets()));;
213}
214
215
216//----------------------------------------------------------------------
217// class PruningRecombiner
218//----------------------------------------------------------------------
219
220// decide whether to recombine things or not
221void PruningRecombiner::recombine(const PseudoJet &pa,
222 const PseudoJet &pb,
223 PseudoJet &pab) const{
224 PseudoJet p;
225 _recombiner->recombine(pa, pb, p);
226
227 // if the 2 particles are close enough, do the recombination
228 if (pa.squared_distance(pb)<=_Rcut2){
229 pab=p; return;
230 }
231
232 double pt2a = pa.perp2();
233 double pt2b = pb.perp2();
234
235 // check which is the softest
236 if (pt2a < pt2b){
237 if (pt2a<_zcut2*p.perp2()){
238 pab = pb; _rejected.push_back(pa.cluster_hist_index());
239 } else {
240 pab = p;
241 }
242 } else {
243 if (pt2b<_zcut2*p.perp2()) {
244 pab = pa; _rejected.push_back(pb.cluster_hist_index());
245 } else {
246 pab = p;
247 }
248 }
249}
250
251// description
252string PruningRecombiner::description() const{
253 ostringstream oss;
254 oss << "Pruning recombiner with zcut = " << sqrt(_zcut2)
255 << ", Rcut = " << sqrt(_Rcut2)
256 << ", and underlying recombiner = " << _recombiner->description();
257 return oss.str();
258}
259
260
261
262
263//----------------------------------------------------------------------
264// class PruningPlugin
265//----------------------------------------------------------------------
266// the actual clustering work for the plugin
267void PruningPlugin::run_clustering(ClusterSequence &input_cs) const{
268 // declare a pruning recombiner
269 PruningRecombiner pruning_recombiner(_zcut, _Rcut, _jet_def.recombiner());
270 JetDefinition jet_def = _jet_def;
271 jet_def.set_recombiner(&pruning_recombiner);
272
273 // cluster the particles using that recombiner
274 ClusterSequence internal_cs(input_cs.jets(), jet_def);
275 const vector<ClusterSequence::history_element> & internal_hist = internal_cs.history();
276
277 // transfer the list of "childless" elements into a bool vector
278 vector<bool> kept(internal_hist.size(), true);
279 const vector<unsigned int> &pr_rej = pruning_recombiner.rejected();
280 for (unsigned int i=0;i<pr_rej.size(); i++) kept[pr_rej[i]]=false;
281
282 // browse the history, keeping only the elements that have not been
283 // vetoed.
284 //
285 // In the process we build a map for the history indices
286 vector<unsigned int> internal2input(internal_hist.size());
287 for (unsigned int i=0; i<input_cs.jets().size(); i++)
288 internal2input[i] = i;
289
290 for (unsigned int i=input_cs.jets().size(); i<internal_hist.size(); i++){
291 const ClusterSequence::history_element &he = internal_hist[i];
292
293 // deal with recombinations with the beam
294 if (he.parent2 == ClusterSequence::BeamJet){
295 int internal_jetp_index = internal_hist[he.parent1].jetp_index;
296 int internal_hist_index = internal_cs.jets()[internal_jetp_index].cluster_hist_index();
297
298 int input_jetp_index = input_cs.history()[internal2input[internal_hist_index]].jetp_index;
299
300 // cout << "Beam recomb for internal " << internal_hist_index
301 // << " (input jet index=" << input_jetp_index << endl;
302
303 input_cs.plugin_record_iB_recombination(input_jetp_index, he.dij);
304 continue;
305 }
306
307 // now, deal with two-body recombinations
308 if (!kept[he.parent1]){ // 1 is rejected, we keep only 2
309 internal2input[i]=internal2input[he.parent2];
310 // cout << "rejecting internal " << he.parent1
311 // << ", mapping internal " << i
312 // << " to internal " << he.parent2
313 // << " i.e. " << internal2input[i] << endl;
314 } else if (!kept[he.parent2]){ // 2 is rejected, we keep only 1
315 internal2input[i]=internal2input[he.parent1];
316 // cout << "rejecting internal " << he.parent2
317 // << ", mapping internal " << i
318 // << " to internal " << he.parent1
319 // << " i.e. " << internal2input[i] << endl;
320 } else { // do the recombination
321 int new_index;
322 input_cs.plugin_record_ij_recombination(input_cs.history()[internal2input[he.parent1]].jetp_index,
323 input_cs.history()[internal2input[he.parent2]].jetp_index,
324 he.dij, internal_cs.jets()[he.jetp_index], new_index);
325 internal2input[i]=input_cs.jets()[new_index].cluster_hist_index();
326 // cout << "merging " << internal2input[he.parent1] << " (int: " << he.parent1 << ")"
327 // << " and " << internal2input[he.parent2] << " (int: " << he.parent2 << ")"
328 // << " into " << internal2input[i] << " (int: " << i << ")" << endl;
329 }
330 }
331}
332
333// returns the plugin description
334string PruningPlugin::description() const{
335 ostringstream oss;
336 oss << "Pruning plugin with jet_definition = (" << _jet_def.description()
337 <<"), zcut = " << _zcut
338 << ", Rcut = " << _Rcut;
339 return oss.str();
340}
341
342
343FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.