FastJet  3.0.6
Filter.cc
1 //STARTHEADER
2 // $Id: Filter.cc 2695 2011-11-14 22:33:07Z salam $
3 //
4 // Copyright (c) 2005-2011, 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 and are described in hep-ph/0512210. If you use
16 // FastJet as part of work towards a scientific publication, please
17 // include a citation to the FastJet paper.
18 //
19 // FastJet is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 // GNU General Public License for more details.
23 //
24 // You should have received a copy of the GNU General Public License
25 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
26 //----------------------------------------------------------------------
27 //ENDHEADER
28 
29 #include "fastjet/tools/Filter.hh"
30 #include <fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh>
31 #include <cassert>
32 #include <algorithm>
33 #include <sstream>
34 #include <typeinfo>
35 
36 using namespace std;
37 
38 
39 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
40 
41 //----------------------------------------------------------------------
42 // Filter class implementation
43 //----------------------------------------------------------------------
44 
45 // class description
46 string Filter::description() const {
47  ostringstream ostr;
48  ostr << "Filter with subjet_def = ";
49  if (_Rfiltfunc) {
50  ostr << "Cambridge/Aachen algorithm with dynamic Rfilt"
51  << " (recomb. scheme deduced from jet, or E-scheme if not unique)";
52  } else if (_Rfilt > 0) {
53  ostr << "Cambridge/Aachen algorithm with Rfilt = "
54  << _Rfilt
55  << " (recomb. scheme deduced from jet, or E-scheme if not unique)";
56  } else {
57  ostr << _subjet_def.description();
58  }
59  ostr<< ", selection " << _selector.description();
60  if (_subtractor) {
61  ostr << ", subtractor: " << _subtractor->description();
62  } else if (_rho != 0) {
63  ostr << ", subtracting with rho = " << _rho;
64  }
65  return ostr.str();
66 }
67 
68 
69 // core functions
70 //----------------------------------------------------------------------
71 
72 // return a vector of subjets, which are the ones that would be kept
73 // by the filtering
74 PseudoJet Filter::result(const PseudoJet &jet) const {
75  // start by getting the list of subjets (including a list of sanity
76  // checks)
77  // NB: subjets is empty to begin with (see the comment for
78  // _set_filtered_elements_cafilt)
79  vector<PseudoJet> subjets;
80  JetDefinition subjet_def;
81  bool discard_area;
82  // NB: on return, subjet_def is set to the jet definition actually
83  // used (so that we can make use of its recombination scheme
84  // when joining the jets to be kept).
85  _set_filtered_elements(jet, subjets, subjet_def, discard_area);
86 
87  // now build the vector of kept and rejected subjets
88  vector<PseudoJet> kept, rejected;
89  // Note that in the following line we make a copy of the _selector
90  // to avoid issues with needing a mutable _selector
91  Selector selector_copy = _selector;
92  if (selector_copy.takes_reference()) selector_copy.set_reference(jet);
93  selector_copy.sift(subjets, kept, rejected);
94 
95  // gather the info under the form of a PseudoJet
96  return _finalise(jet, kept, rejected, subjet_def, discard_area);
97 }
98 
99 
100 // sets filtered_elements to be all the subjets on which filtering will work
101 void Filter::_set_filtered_elements(const PseudoJet & jet,
102  vector<PseudoJet> & filtered_elements,
103  JetDefinition & subjet_def,
104  bool & discard_area) const {
105  // sanity checks
106  //-------------------------------------------------------------------
107  // make sure that the jet has constituents
108  if (! jet.has_constituents())
109  throw Error("Filter can only be applied on jets having constituents");
110 
111  // for a whole variety of checks, we shall need the "recursive"
112  // pieces of the jet (the jet itself or recursing down to its most
113  // fundamental pieces).
114  // So we do compute these once and for all
115  vector<PseudoJet> all_pieces; //.clear();
116  if ((!_get_all_pieces(jet, all_pieces)) || (all_pieces.size()==0))
117  throw Error("Attempt to filter a jet that has no associated ClusterSequence or is not a superposition of jets associated with a ClusterSequence");
118 
119  // if the filter uses subtraction, make sure we have a CS that supports area and has
120  // explicit ghosts
121  if (_uses_subtraction()) {
122  if (!jet.has_area())
123  throw Error("Attempt to filter and subtract (non-zero rho or subtractor) without area info for the original jet");
124 
125  if (!_check_explicit_ghosts(all_pieces))
126  throw Error("Attempt to filter and subtract (non-zero rho or subtractor) without explicit ghosts");
127  }
128 
129  // if we're dealing with a dynamic determination of the filtering
130  // radius, do it now
131  if ((_Rfilt>=0) || (_Rfiltfunc)){
132  double Rfilt = (_Rfiltfunc) ? (*_Rfiltfunc)(jet) : _Rfilt;
133  const JetDefinition::Recombiner * common_recombiner = _get_common_recombiner(all_pieces);
134  if (common_recombiner) {
135  if (typeid(*common_recombiner) == typeid(JetDefinition::DefaultRecombiner)) {
136  RecombinationScheme scheme =
137  static_cast<const JetDefinition::DefaultRecombiner *>(common_recombiner)->scheme();
138  subjet_def = JetDefinition(cambridge_algorithm, Rfilt, scheme);
139  } else {
140  subjet_def = JetDefinition(cambridge_algorithm, Rfilt, common_recombiner);
141  }
142  } else {
143  subjet_def = JetDefinition(cambridge_algorithm, Rfilt);
144  }
145  } else {
146  subjet_def = _subjet_def;
147  }
148 
149  // get the jet definition to be use and whether we can apply our
150  // simplified C/A+C/A filter
151  //
152  // we apply C/A clustering iff
153  // - the request subjet_def is C/A
154  // - the jet is either directly coming from C/A or if it is a
155  // superposition of C/A jets
156  // - the pieces agree with the recombination scheme of subjet_def
157  //------------------------------------------------------------------
158  bool simple_cafilt = _check_ca(all_pieces);
159 
160  // extract the subjets
161  //-------------------------------------------------------------------
162  discard_area = false;
163  if (simple_cafilt){
164  // first make sure that 'filtered_elements' is empty
165  filtered_elements.clear();
166  _set_filtered_elements_cafilt(jet, filtered_elements, subjet_def.R());
167  // in the following case, areas can be erroneous and will be discarded
168  discard_area = (!_uses_subtraction()) && (jet.has_area()) && (!_check_explicit_ghosts(all_pieces));
169  } else {
170  // here we'll simply recluster the jets.
171  //
172  // To include an area support we need
173  // - the jet to have an area
174  // - subtraction requested or explicit ghosts
175  bool do_areas = (jet.has_area()) && ((_uses_subtraction()) || (_check_explicit_ghosts(all_pieces)));
176  _set_filtered_elements_generic(jet, filtered_elements, subjet_def, do_areas);
177  }
178 
179  // order the filtered elements in pt
180  filtered_elements = sorted_by_pt(filtered_elements);
181 }
182 
183 // set the filtered elements in the simple case of C/A+C/A
184 //
185 // WATCH OUT: this could be recursively called, so filtered elements
186 // of 'jet' are APPENDED to 'filtered_elements'
187 void Filter::_set_filtered_elements_cafilt(const PseudoJet & jet,
188  vector<PseudoJet> & filtered_elements,
189  double Rfilt) const{
190  // we know that the jet is either a C/A jet or a superposition of
191  // such pieces
193  // just extract the exclusive subjets of 'jet'
194  const ClusterSequence *cs = jet.associated_cluster_sequence();
195  vector<PseudoJet> local_fe;
196 
197  double dcut = Rfilt / cs->jet_def().R();
198  if (dcut>=1.0){
199  local_fe.push_back(jet);
200  } else {
201  dcut *= dcut;
202  local_fe = jet.exclusive_subjets(dcut);
203  }
204 
205  // subtract the jets if needed
206  // Note that this one would work on pieces!!
207  //-----------------------------------------------------------------
208  if (_uses_subtraction()){
209  const ClusterSequenceAreaBase * csab = jet.validated_csab();
210  for (unsigned int i=0;i<local_fe.size();i++) {
211  if (_subtractor) {
212  local_fe[i] = (*_subtractor)(local_fe[i]);
213  } else {
214  local_fe[i] = csab->subtracted_jet(local_fe[i], _rho);
215  }
216  }
217  }
218 
219  copy(local_fe.begin(), local_fe.end(), back_inserter(filtered_elements));
220  return;
221  }
222 
223  // just recurse into the pieces
224  const vector<PseudoJet> & pieces = jet.pieces();
225  for (vector<PseudoJet>::const_iterator it = pieces.begin();
226  it!=pieces.end(); it++)
227  _set_filtered_elements_cafilt(*it, filtered_elements, Rfilt);
228 }
229 
230 
231 // set the filtered elements in the generic re-clustering case (wo
232 // subtraction)
233 void Filter::_set_filtered_elements_generic(const PseudoJet & jet,
234  vector<PseudoJet> & filtered_elements,
235  const JetDefinition & subjet_def,
236  bool do_areas) const{
237  // create a new, internal, ClusterSequence from the jet constituents
238  // get the subjets directly from there
239  //
240  // If the jet has area support then we separate the ghosts from the
241  // "regular" particles so the subjets will also have area
242  // support. Note that we do this regardless of whether rho is zero
243  // or not.
244  //
245  // Note that to be able to separate the ghosts, one needs explicit
246  // ghosts!!
247  // ---------------------------------------------------------------
248  if (do_areas){
249  vector<PseudoJet> all_constituents = jet.constituents();
250  vector<PseudoJet> regular_constituents, ghosts;
251 
252  for (vector<PseudoJet>::iterator it = all_constituents.begin();
253  it != all_constituents.end(); it++){
254  if (it->is_pure_ghost())
255  ghosts.push_back(*it);
256  else
257  regular_constituents.push_back(*it);
258  }
259 
260  // figure out the ghost area from the 1st ghost (if none, any value
261  // would probably do as the area will be 0 and subtraction will have
262  // no effect!)
263  double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
265  = new ClusterSequenceActiveAreaExplicitGhosts(regular_constituents,
266  subjet_def,
267  ghosts, ghost_area);
268 
269  // get the subjets: we use the subtracted or unsubtracted ones
270  // depending on rho or _subtractor being non-zero
271  if (_uses_subtraction()) {
272  if (_subtractor) {
273  filtered_elements = (*_subtractor)(csa->inclusive_jets());
274  } else {
275  filtered_elements = csa->subtracted_jets(_rho);
276  }
277  } else {
278  filtered_elements = csa->inclusive_jets();
279  }
280 
281  // allow the cs to be deleted when it's no longer used
283  } else {
284  ClusterSequence * cs = new ClusterSequence(jet.constituents(), subjet_def);
285  filtered_elements = cs->inclusive_jets();
286  // allow the cs to be deleted when it's no longer used
288  }
289 }
290 
291 
292 // gather the information about what is kept and rejected under the
293 // form of a PseudoJet with a special ClusterSequenceInfo
294 PseudoJet Filter::_finalise(const PseudoJet & /*jet*/,
295  vector<PseudoJet> & kept,
296  vector<PseudoJet> & rejected,
297  const JetDefinition & subjet_def,
298  const bool discard_area) const {
299  // figure out which recombiner to use
300  const JetDefinition::Recombiner &rec = *(subjet_def.recombiner());
301 
302  // create an appropriate structure and transfer the info to it
303  PseudoJet filtered_jet = join<StructureType>(kept, rec);
304  StructureType *fs = (StructureType*) filtered_jet.structure_non_const_ptr();
305 // fs->_original_jet = jet;
306  fs->_rejected = rejected;
307 
308  if (discard_area){
309  // safety check: make sure there is an area to discard!!!
310  assert(fs->_area_4vector_ptr);
311  delete fs->_area_4vector_ptr;
312  fs->_area_4vector_ptr=0;
313  }
314 
315  return filtered_jet;
316 }
317 
318 // various checks
319 //----------------------------------------------------------------------
320 
321 // get the pieces down to the fundamental pieces
322 //
323 // Note that this just checks that there is an associated CS to the
324 // fundamental pieces, not that it is still valid
325 bool Filter::_get_all_pieces(const PseudoJet &jet, vector<PseudoJet> &all_pieces) const{
327  all_pieces.push_back(jet);
328  return true;
329  }
330 
331  if (jet.has_pieces()){
332  const vector<PseudoJet> pieces = jet.pieces();
333  for (vector<PseudoJet>::const_iterator it=pieces.begin(); it!=pieces.end(); it++)
334  if (!_get_all_pieces(*it, all_pieces)) return false;
335  return true;
336  }
337 
338  return false;
339 }
340 
341 // get the common recombiner to all pieces (NULL if none)
342 //
343 // Note that if the jet has an associated cluster sequence that is no
344 // longer valid, an error will be thrown (needed since it could be the
345 // 1st check called after the enumeration of the pieces)
346 const JetDefinition::Recombiner* Filter::_get_common_recombiner(const vector<PseudoJet> &all_pieces) const{
347  const JetDefinition & jd_ref = all_pieces[0].validated_cs()->jet_def();
348  for (unsigned int i=1; i<all_pieces.size(); i++)
349  if (!all_pieces[i].validated_cs()->jet_def().has_same_recombiner(jd_ref)) return NULL;
350 
351  return jd_ref.recombiner();
352 }
353 
354 // check if the jet (or all its pieces) have explicit ghosts
355 // (assuming the jet has area support
356 //
357 // Note that if the jet has an associated cluster sequence that is no
358 // longer valid, an error will be thrown (needed since it could be the
359 // 1st check called after the enumeration of the pieces)
360 bool Filter::_check_explicit_ghosts(const vector<PseudoJet> &all_pieces) const{
361  for (vector<PseudoJet>::const_iterator it=all_pieces.begin(); it!=all_pieces.end(); it++)
362  if (! it->validated_csab()->has_explicit_ghosts()) return false;
363  return true;
364 }
365 
366 // check if one can apply the simplification for C/A subjets
367 //
368 // This includes:
369 // - the subjet definition asks for C/A subjets
370 // - all the pieces share the same CS
371 // - that CS is C/A with the same recombiner as the subjet def
372 // - the filtering radius is not larger than any of the pairwise
373 // distance between the pieces
374 //
375 // Note that if the jet has an associated cluster sequence that is no
376 // longer valid, an error will be thrown (needed since it could be the
377 // 1st check called after the enumeration of the pieces)
378 bool Filter::_check_ca(const vector<PseudoJet> &all_pieces) const{
379  if (_subjet_def.jet_algorithm() != cambridge_algorithm) return false;
380 
381  // check that the 1st of all the pieces (we're sure there is at
382  // least one) is coming from a C/A clustering. Then check that all
383  // the following pieces share the same ClusterSequence
384  const ClusterSequence * cs_ref = all_pieces[0].validated_cs();
385  if (cs_ref->jet_def().jet_algorithm() != cambridge_algorithm) return false;
386  for (unsigned int i=1; i<all_pieces.size(); i++)
387  if (all_pieces[i].validated_cs() != cs_ref) return false;
388 
389  // check that the 1st peice has the same recombiner as the one used
390  // for the subjet clustering
391  // Note that since they share the same CS, checking the 2st one is enough
392  if (!cs_ref->jet_def().has_same_recombiner(_subjet_def)) return false;
393 
394  // we also have to make sure that the filtering radius is not larger
395  // than any of the inter-pieces distance
396  double Rfilt2 = _subjet_def.R();
397  Rfilt2 *= Rfilt2;
398  for (unsigned int i=0; i<all_pieces.size()-1; i++){
399  for (unsigned int j=i+1; j<all_pieces.size(); j++){
400  if (all_pieces[i].squared_distance(all_pieces[j]) < Rfilt2) return false;
401  }
402  }
403 
404  return true;
405 }
406 
407 //----------------------------------------------------------------------
408 // FilterInterface implementation
409 //----------------------------------------------------------------------
410 
411 
412 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
const JetDefinition & jet_def() const
return a reference to the jet definition
std::vector< PseudoJet > subtracted_jets(const double rho, const double ptmin=0.0) const
return a vector of all subtracted jets, using area_4vector, given rho.
const ClusterSequenceAreaBase * validated_csab() const
shorthand for validated_cluster_sequence_area_base()
Definition: PseudoJet.cc:685
JetAlgorithm jet_algorithm() const
return information about the definition...
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition: PseudoJet.cc:785
virtual bool has_area() const
check if it has a defined area
Definition: PseudoJet.cc:694
deals with clustering
bool has_associated_cluster_sequence() const
returns true if this PseudoJet has an associated ClusterSequence.
Definition: PseudoJet.cc:416
std::vector< PseudoJet > exclusive_subjets(const double &dcut) const
return a vector of all subjets of the current jet (in the sense of the exclusive algorithm) that woul...
Definition: PseudoJet.cc:589
PseudoJetStructureBase * structure_non_const_ptr()
return a non-const pointer to the structure (of type PseudoJetStructureBase*) associated with this Ps...
Definition: PseudoJet.cc:480
bool takes_reference() const
returns true if this can be applied jet by jet
Definition: Selector.hh:268
virtual bool has_pieces() const
returns true if a jet has pieces
Definition: PseudoJet.cc:657
A class that will provide the recombination scheme facilities and/or allow a user to extend these fac...
virtual std::vector< PseudoJet > constituents() const
retrieve the constituents.
Definition: PseudoJet.cc:566
const Recombiner * recombiner() const
return a pointer to the currently defined recombiner.
base class that sets interface for extensions of ClusterSequence that provide information about the a...
virtual std::vector< PseudoJet > pieces() const
retrieve the pieces that make up the jet.
Definition: PseudoJet.cc:666
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:41
Class to contain structure information for a filtered jet.
Definition: Filter.hh:220
Like ClusterSequence with computation of the active jet area with the addition of explicit ghosts...
RecombinationScheme
the various recombination schemes
virtual bool has_constituents() const
returns true if the PseudoJet has constituents
Definition: PseudoJet.cc:560
const ClusterSequence * associated_cluster_sequence() const
get a (const) pointer to the parent ClusterSequence (NULL if inexistent)
Definition: PseudoJet.cc:423
const Selector & set_reference(const PseudoJet &reference)
set the reference jet for this Selector
Definition: Selector.hh:273
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:101
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition: Selector.hh:141
PseudoJet subtracted_jet(const PseudoJet &jet, const double rho) const
return a subtracted jet, using area_4vector, given rho
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
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.
bool has_same_recombiner(const JetDefinition &other_jd) const
returns true if the current jet definitions shares the same recombiner as teh one passed as an argume...
An abstract base class that will provide the recombination scheme facilities and/or allow a user to e...
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:65
class that is intended to hold a full definition of the jet clusterer