FastJet  3.0.6
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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
192  if (jet.has_associated_cluster_sequence()){
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;
264  ClusterSequenceActiveAreaExplicitGhosts * csa
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
282  csa->delete_self_when_unused();
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
287  cs->delete_self_when_unused();
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{
326  if (jet.has_associated_cluster_sequence()){
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
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition: PseudoJet.cc:785
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
bool takes_reference() const
returns true if this can be applied jet by jet
Definition: Selector.hh:268
virtual bool has_constituents() const
returns true if the PseudoJet has constituents
Definition: PseudoJet.cc:560
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:41
const Selector & set_reference(const PseudoJet &reference)
set the reference jet for this Selector
Definition: Selector.hh:273
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition: Selector.hh:141
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
virtual bool has_area() const
check if it has a defined area
Definition: PseudoJet.cc:694
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