FastJet  3.1.0-beta.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Recluster.cc
1 // $Id: Recluster.cc 3629 2014-08-14 17:21:15Z salam $
2 //
3 // Copyright (c) 2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
4 //
5 //----------------------------------------------------------------------
6 // This file is part of FastJet contrib.
7 //
8 // It is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU General Public License as published by the
10 // Free Software Foundation; either version 2 of the License, or (at
11 // your option) any later version.
12 //
13 // It is distributed in the hope that it will be useful, but WITHOUT
14 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 // or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
16 // License for more details.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with this code. If not, see <http://www.gnu.org/licenses/>.
20 //----------------------------------------------------------------------
21 
22 #include "fastjet/tools/Recluster.hh"
23 #include "fastjet/CompositeJetStructure.hh"
24 #include <fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh>
25 #include <sstream>
26 #include <typeinfo>
27 
28 using namespace std;
29 
30 // Comments:
31 //
32 // - If the jet comes from a C/A clustering (or is a composite jet
33 // made of C/A clusterings) and we recluster it with a C/A
34 // algorithm, we just use exclusive jets instead of doing the
35 // clustering explicitly. In this specific case, we need to make
36 // sure that all the pieces share the same cluster sequence.
37 //
38 // - If the recombiner has to be taken from the original jet and that
39 // jet is composite, we need to check that all the pieces were
40 // obtained with the same recombiner.
41 //
42 // - Note that a preliminary version of this code has been
43 // distributed in the RecursiveTools fastjet-contrib. The present
44 // version has a few minor fixed
45 
46 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
47 
48 LimitedWarning Recluster::_explicit_ghost_warning;
49 
50 // ctor
51 Recluster::Recluster(JetAlgorithm new_jet_alg, double new_jet_radius, Keep keep_in)
52  : _new_jet_def(JetDefinition(new_jet_alg, new_jet_radius)),
53  _acquire_recombiner(true), _keep(keep_in), _cambridge_optimisation_enabled(true){}
54 
55 Recluster::Recluster(JetAlgorithm new_jet_alg, Keep keep_in)
56  : _acquire_recombiner(true), _keep(keep_in), _cambridge_optimisation_enabled(true){
57  switch (JetDefinition::n_parameters_for_algorithm(new_jet_alg)){
58  case 0: _new_jet_def = JetDefinition(new_jet_alg); break;
59  case 1: _new_jet_def = JetDefinition(new_jet_alg, JetDefinition::max_allowable_R); break;
60  default:
61  throw Error("Recluster(): tried to construct specifying only a jet algorithm ("+JetDefinition::algorithm_description(new_jet_alg)+") which takes more than 1 parameter");
62  };
63 }
64 
65 
66 // class description
67 string Recluster::description() const {
68  ostringstream ostr;
69  ostr << "Recluster with new_jet_def = ";
70  if (_acquire_recombiner){
71  ostr << _new_jet_def.description_no_recombiner();
72  ostr << ", using a recombiner obtained from the jet being reclustered";
73  } else {
74  ostr << _new_jet_def.description();
75  }
76 
77  if (_keep == keep_only_hardest)
78  ostr << " and keeping the hardest inclusive jet";
79  else
80  ostr << " and joining all inclusive jets into a composite jet";
81 
82  return ostr.str();
83 }
84 
85 
86 // the main piece of code that performs the reclustering
88  // get the incljets and the exact jet definition that has been used
89  // to get them
90  vector<PseudoJet> incljets;
91  bool ca_optimised = get_new_jets_and_def(jet, incljets);
92 
93  return generate_output_jet(incljets, ca_optimised);
94 }
95 
96 
97 // a lower-level method that does the actual work of reclustering the
98 // input jet. The resulting incljets are stored in output_jets and the
99 // jet definition that has been used can be deduced from their
100 // associated ClusterSequence
101 //
102 // - input_jet the (input) jet that one wants to recluster
103 // - output_jets incljets resulting from the new clustering
104 //
105 // returns true if the C/A optimisation has been used (this means
106 // that generate_output_jet will watch out for non-explicit-ghost
107 // areas that might be leftover)
109  vector<PseudoJet> & output_jets) const{
110  // generic sanity checks
111  //-------------------------------------------------------------------
112  // make sure that the jet has constituents
113  if (! input_jet.has_constituents())
114  throw Error("Recluster can only be applied on jets having constituents");
115 
116  // tests particular to certain configurations
117  //-------------------------------------------------------------------
118 
119  // for a whole variety of tests, we shall need the "recursive"
120  // pieces of the jet (the jet itself or recursing down to its most
121  // fundamental pieces). So we do compute these once and for all.
122  //
123  // Note that the pieces are always needed (either for C/A or for the
124  // area checks)
125  vector<PseudoJet> all_pieces;
126  if ((!_get_all_pieces(input_jet, all_pieces)) || (all_pieces.size()==0)){
127  throw Error("Recluster: failed to retrieve all the pieces composing the jet.");
128  }
129 
130  // decide which jet definition to use
131  //-------------------------------------------------------------------
132  JetDefinition new_jet_def = _new_jet_def;
133  if (_acquire_recombiner){
134  _acquire_recombiner_from_pieces(all_pieces, new_jet_def);
135  }
136 
137  // the vector that will ultimately hold the incljets
138  output_jets.clear();
139 
140  // check if we can apply the simplification for C/A jets reclustered
141  // with C/A
142  //
143  // we apply C/A clustering iff
144  // - the requested new_jet_def is C/A
145  // - the jet is either directly coming from C/A or if it is a
146  // superposition of C/A jets from the same cluster sequence
147  // - the pieces agree with the recombination scheme of new_jet_def
148  //
149  // In this case area support will be automatically inherited so we
150  // can only worry about this later
151  // -------------------------------------------------------------------
152  if (_check_ca(all_pieces, new_jet_def)){
153  _recluster_ca(all_pieces, output_jets, new_jet_def.R());
154  output_jets = sorted_by_pt(output_jets);
155  return true;
156  }
157 
158  // decide if area support has to be kept
159  //-------------------------------------------------------------------
160  bool include_area_support = input_jet.has_area();
161  if ((include_area_support) && (!_check_explicit_ghosts(all_pieces))){
162  _explicit_ghost_warning.warn("Recluster: the original cluster sequence is lacking explicit ghosts; area support will no longer be available after re-clustering");
163  include_area_support = false;
164  }
165 
166  // extract the incljets
167  //-------------------------------------------------------------------
168  _recluster_generic(input_jet, output_jets, new_jet_def, include_area_support);
169  output_jets = sorted_by_pt(output_jets);
170 
171  return false;
172 }
173 
174 // given a set of incljets and a jet definition used, create the
175 // resulting PseudoJet
176 PseudoJet Recluster::generate_output_jet(std::vector<PseudoJet> & incljets,
177  bool ca_optimisation_used) const{
178  // first handle the case where we only need to keep the hardest incljet
179  if (_keep == keep_only_hardest) {
180  if (incljets.size() > 0) {
181  return incljets[0];
182  } else {
183  return PseudoJet();
184  }
185  }
186 
187  // now the case where all incljets have to be joined
188 
189  // safekeeper
190  if (incljets.size()==0) return join(incljets);
191 
192  PseudoJet reclustered = join(incljets,
193  *(incljets[0].associated_cluster_sequence()->jet_def().recombiner()));
194 
195  // if we've used C/A optimisation, we need to get rid of the area
196  // information if it comes from a non-explicit-ghost clustering.
197  // (because in that case it can be erroneous due to the lack of
198  // information about empty areas)
199  if (ca_optimisation_used){
200  if (reclustered.has_area() &&
201  (incljets.size() > 0) &&
202  (! incljets[0].validated_csab()->has_explicit_ghosts())){
203  CompositeJetStructure *css = dynamic_cast<CompositeJetStructure *>(reclustered.structure_non_const_ptr());
204  assert(css);
205  css->discard_area();
206  }
207  }
208 
209  return reclustered;
210 }
211 
212 //----------------------------------------------------------------------
213 // the parts that really do the reclustering
214 //----------------------------------------------------------------------
215 
216 // get the subjets in the simple case of C/A+C/A
217 void Recluster::_recluster_ca(const vector<PseudoJet> & all_pieces,
218  vector<PseudoJet> & subjets,
219  double Rfilt) const{
220  subjets.clear();
221 
222  // each individual piece should have a C/A cluster sequence
223  // associated to it. So a simple loop would do the job
224  for (vector<PseudoJet>::const_iterator piece_it = all_pieces.begin();
225  piece_it!=all_pieces.end(); piece_it++){
226  // just extract the exclusive subjets of 'jet'
227  const ClusterSequence *cs = piece_it->associated_cluster_sequence();
228  vector<PseudoJet> local_subjets;
229 
230  double dcut = Rfilt / cs->jet_def().R();
231  if (dcut>=1.0){
232  // remember that in this case all the pairwise interpiece
233  // distances are supposed to be larger than Rfilt (this was
234  // tested in _check_ca), which means that they can never
235  // recombine with each other.
236  local_subjets.push_back(*piece_it);
237  } else {
238  local_subjets = piece_it->exclusive_subjets(dcut*dcut);
239  }
240 
241  copy(local_subjets.begin(), local_subjets.end(), back_inserter(subjets));
242  }
243 }
244 
245 
246 // perform the reclustering itself for all cases where the "C/A trick"
247 // does not apply
248 void Recluster::_recluster_generic(const PseudoJet & jet,
249  vector<PseudoJet> & incljets,
250  const JetDefinition & new_jet_def,
251  bool do_areas) const{
252  // create a new, internal, ClusterSequence from the jet constituents
253  // get the incljets directly from there
254  //
255  // If the jet has area support then we separate the ghosts from the
256  // "regular" particles so the incljets will also have area
257  // support. Note that we do this regardless of whether rho is zero
258  // or not.
259  //
260  // Note that to be able to separate the ghosts, one needs explicit
261  // ghosts!!
262  // ---------------------------------------------------------------
263  if (do_areas){
264  //vector<PseudoJet> all_constituents = jet.constituents();
265  vector<PseudoJet> regular_constituents, ghosts;
266  SelectorIsPureGhost().sift(jet.constituents(), ghosts, regular_constituents);
267 
268  // figure out the ghost area from the 1st ghost (if none, any value
269  // would probably do as the area will be 0 and subtraction will have
270  // no effect!)
271  double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
272  ClusterSequenceActiveAreaExplicitGhosts * csa
273  = new ClusterSequenceActiveAreaExplicitGhosts(regular_constituents,
274  new_jet_def,
275  ghosts, ghost_area);
276 
277  incljets = csa->inclusive_jets();
278 
279  // allow the cs to be deleted when it's no longer used
280  //
281  // Note that there is at least one constituent in the jet so there
282  // is in principle at least one incljet. But one may have used a
283  // nasty recombiner or jet def that left an empty set of incljets,
284  // so we'd rather play it safe (e.g. GridJetPlugin, with the
285  // constituents outside the range of the grid)
286  if (incljets.size())
287  csa->delete_self_when_unused();
288  else
289  delete csa;
290  } else {
291  ClusterSequence * cs = new ClusterSequence(jet.constituents(), new_jet_def);
292  incljets = cs->inclusive_jets();
293  // allow the cs to be deleted when it's no longer used (again, we
294  // add an extra safety check)
295  if (incljets.size())
296  cs->delete_self_when_unused();
297  else
298  delete cs;
299  }
300 }
301 
302 //----------------------------------------------------------------------
303 // various checks and internal constructs
304 //----------------------------------------------------------------------
305 
306 // fundamental info for CompositeJets
307 //----------------------------------------------------------------------
308 
309 // get the pieces down to the fundamental pieces
310 //
311 // Note that this just checks that there is an associated CS to the
312 // fundamental pieces, not that it is still valid
313 bool Recluster::_get_all_pieces(const PseudoJet &jet, vector<PseudoJet> &all_pieces) const{
314  if (jet.has_associated_cluster_sequence()){
315  all_pieces.push_back(jet);
316  return true;
317  }
318 
319  if (jet.has_pieces()){
320  const vector<PseudoJet> pieces = jet.pieces();
321  for (vector<PseudoJet>::const_iterator it=pieces.begin(); it!=pieces.end(); it++)
322  if (!_get_all_pieces(*it, all_pieces)) return false;
323  return true;
324  }
325 
326  return false;
327 }
328 
329 // construct the re-clustering jet definition using the recombiner
330 // from whatever definition has been used to obtain the original jet
331 //----------------------------------------------------------------------
332 void Recluster::_acquire_recombiner_from_pieces(const vector<PseudoJet> &all_pieces,
333  JetDefinition &new_jet_def) const{
334  // a quick safety check
335  assert(_acquire_recombiner);
336 
337  // check that all the pieces have the same recombiner
338  //
339  // Note that if the jet has an associated cluster sequence that is no
340  // longer valid, an error will be thrown (needed since it could be the
341  // 1st check called after the enumeration of the pieces)
342  const JetDefinition & jd_ref = all_pieces[0].validated_cs()->jet_def();
343  for (unsigned int i=1; i<all_pieces.size(); i++){
344  if (!all_pieces[i].validated_cs()->jet_def().has_same_recombiner(jd_ref)){
345  throw Error("Recluster instance is configured to determine the recombination scheme (or recombiner) from the original jet, but different pieces of the jet were found to have non-equivalent recombiners.");
346  }
347  }
348 
349  // get the recombiner from the original jet_def
350  new_jet_def.set_recombiner(jd_ref);
351 }
352 
353 // area support
354 //----------------------------------------------------------------------
355 
356 // check if the jet (or all its pieces) have explicit ghosts
357 // (assuming the jet has area support).
358 //
359 // Note that if the jet has an associated cluster sequence that is no
360 // longer valid, an error will be thrown (needed since it could be the
361 // 1st check called after the enumeration of the pieces)
362 bool Recluster::_check_explicit_ghosts(const vector<PseudoJet> &all_pieces) const{
363  for (vector<PseudoJet>::const_iterator it=all_pieces.begin(); it!=all_pieces.end(); it++)
364  if (! it->validated_csab()->has_explicit_ghosts()) return false;
365  return true;
366 }
367 
368 // C/A specific tests
369 //----------------------------------------------------------------------
370 
371 // check if one can apply the simplification for C/A incljets
372 //
373 // This includes:
374 // - the incljet definition asks for C/A incljets
375 // - all the pieces share the same CS
376 // - that CS is C/A with the same recombiner as the incljet def
377 // - the re-clustering radius is not larger than any of the pairwise
378 // distance between the pieces
379 //
380 // Note that if the jet has an associated cluster sequence that is no
381 // longer valid, an error will be thrown (needed since it could be the
382 // 1st check called after the enumeration of the pieces)
383 bool Recluster::_check_ca(const vector<PseudoJet> &all_pieces,
384  const JetDefinition &new_jet_def) const{
385  // check that optimisation is enabled
386  if (!_cambridge_optimisation_enabled) return false;
387 
388  // check that we're reclustering with C/A
389  if (new_jet_def.jet_algorithm() != cambridge_algorithm) return false;
390 
391  // check that the 1st of all the pieces (we're sure there is at
392  // least one) is coming from a C/A clustering. Then check that all
393  // the following pieces share the same ClusterSequence
394  const ClusterSequence * cs_ref = all_pieces[0].validated_cs();
395  if (cs_ref->jet_def().jet_algorithm() != cambridge_algorithm) return false;
396  for (unsigned int i=1; i<all_pieces.size(); i++)
397  if (all_pieces[i].validated_cs() != cs_ref) return false;
398 
399  // check that the 1st piece has the same recombiner as the one used
400  // for the incljet clustering
401  // Note that since they share the same CS, checking the 1st one is enough
402  if (!cs_ref->jet_def().has_same_recombiner(new_jet_def)) return false;
403 
404  // we also have to make sure that the reclustering radius is not larger
405  // than any of the inter-piece distances
406  double Rnew2 = new_jet_def.R();
407  Rnew2 *= Rnew2;
408  for (unsigned int i=0; i<all_pieces.size()-1; i++){
409  for (unsigned int j=i+1; j<all_pieces.size(); j++){
410  if (all_pieces[i].squared_distance(all_pieces[j]) < Rnew2) return false;
411  }
412  }
413 
414  return true;
415 }
416 
417 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh