FastJet  3.3.3
Pruner.cc
1 //FJSTARTHEADER
2 // $Id: Pruner.cc 4420 2019-11-29 09:28:20Z soyez $
3 //
4 // Copyright (c) 2005-2019, 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 
39 using namespace std;
40 
41 
42 FASTJET_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
54 Pruner::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
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)
128 
129  return result_local;
130 }
131 
132 // check if the jet has explicit_ghosts (knowing that there is an
133 // area support)
134 bool Pruner::_check_explicit_ghosts(const PseudoJet &jet) const{
135  // if the jet comes from a Clustering check explicit ghosts in that
136  // clustering
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.
159 bool Pruner::_check_common_recombiner(const PseudoJet &jet,
160  JetDefinition &jet_def_for_recombiner,
161  bool assigned) const{
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
189 std::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
211 vector<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
221 void 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
252 string 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
267 void 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
334 string 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 
343 FASTJET_END_NAMESPACE
const JetDefinition & jet_def() const
return a reference to the jet definition
const ClusterSequenceAreaBase * validated_csab() const
shorthand for validated_cluster_sequence_area_base()
Definition: PseudoJet.cc:703
virtual PseudoJet result(const PseudoJet &jet) const
action on a single jet
Definition: Pruner.cc:64
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition: PseudoJet.cc:774
virtual bool has_area() const
check if it has a defined area
Definition: PseudoJet.cc:712
deals with clustering
Selector SelectorNHardest(unsigned int n)
select the n hardest objects
Definition: Selector.cc:1074
const std::vector< history_element > & history() const
allow the user to access the raw internal history.
bool has_associated_cluster_sequence() const
returns true if this PseudoJet has an associated ClusterSequence.
Definition: PseudoJet.cc:434
std::vector< PseudoJet > extra_jets() const
return the other jets that may have been found along with the result of the pruning The resulting vec...
Definition: Pruner.cc:211
void delete_plugin_when_unused()
calling this causes the JetDefinition to handle the deletion of the plugin when it is no longer used ...
const ClusterSequence * validated_cs() const
shorthand for validated_cluster_sequence()
Definition: PseudoJet.cc:461
double squared_distance(const PseudoJet &other) const
returns squared cylinder (rap-phi) distance between this jet and another
Definition: PseudoJet.hh:196
std::vector< PseudoJet > inclusive_jets(const double ptmin=0.0) const
return a vector of all jets (in the sense of the inclusive algorithm) with pt >= ptmin.
int jetp_index
index in _history where the current jet is recombined with another jet to form its child...
void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, int &newjet_k)
record the fact that there has been a recombination between jets()[jet_i] and jets()[jet_k], with the specified dij, and return the index (newjet_k) allocated to the new jet, whose momentum is assumed to be the 4-vector sum of that of jet_i and jet_j
virtual bool has_pieces() const
returns true if a jet has pieces
Definition: PseudoJet.cc:675
int parent2
index in _history where first parent of this jet was created (InexistentParent if this jet is an orig...
double dij
index in the _jets vector where we will find the
void plugin_record_iB_recombination(int jet_i, double diB)
record the fact that there has been a recombination between jets()[jet_i] and the beam...
virtual std::vector< PseudoJet > constituents() const
retrieve the constituents.
Definition: PseudoJet.cc:584
double perp() const
returns the scalar transverse momentum
Definition: PseudoJet.hh:145
const Recombiner * recombiner() const
returns a pointer to the currently defined recombiner.
The structure associated with a PseudoJet thas has gone through a Pruner transformer.
Definition: Pruner.hh:189
virtual std::string description() const
description
Definition: Pruner.cc:189
virtual std::vector< PseudoJet > pieces() const
retrieve the pieces that make up the jet.
Definition: PseudoJet.cc:684
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:47
double m() const
returns the invariant mass (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP) ...
Definition: PseudoJet.hh:970
Like ClusterSequence with computation of the active jet area with the addition of explicit ghosts...
int cluster_hist_index() const
return the cluster_hist_index, intended to be used by clustering routines.
Definition: PseudoJet.hh:784
double perp2() const
returns the squared transverse momentum
Definition: PseudoJet.hh:143
an implementation of C++0x shared pointers (or boost&#39;s)
Definition: SharedPtr.hh:121
virtual bool has_constituents() const
returns true if the PseudoJet has constituents
Definition: PseudoJet.cc:578
virtual bool has_explicit_ghosts() const
returns true if ghosts are explicitly included within jets for this ClusterSequence; ...
void sift(const std::vector< PseudoJet > &jets, std::vector< PseudoJet > &jets_that_pass, std::vector< PseudoJet > &jets_that_fail) const
sift the input jets into two vectors – those that pass the selector and those that do not ...
Definition: Selector.cc:153
a single element in the clustering history
void set_recombiner(const Recombiner *recomb)
set the recombiner class to the one provided
bool has_same_recombiner(const JetDefinition &other_jd) const
returns true if the current jet definitions shares the same recombiner as the one passed as an argume...
std::string description() const
return a string describing what kind of PseudoJet we are dealing with
Definition: PseudoJet.cc:413
Selector SelectorIsPureGhost()
select objects that are (or are only made of) ghosts.
Definition: Selector.cc:1420
std::string description() const
return a textual description of the current jet definition
void delete_self_when_unused()
by calling this routine you tell the ClusterSequence to delete itself when all the Pseudojets associa...
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
virtual std::string description() const
returns a description of the function (an empty string by default)
class that is intended to hold a full definition of the jet clusterer