FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
Filter.cc
1 //FJSTARTHEADER
2 // $Id: Filter.cc 3760 2014-12-19 10:05:10Z soyez $
3 //
4 // Copyright (c) 2005-2014, 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/Filter.hh"
32 #include "fastjet/tools/Recluster.hh"
33 #include <fastjet/ClusterSequenceActiveAreaExplicitGhosts.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 // Filter class implementation
46 //----------------------------------------------------------------------
47 
48 // class description
49 string Filter::description() const {
50  if (!_initialised){
51  return "uninitialised Filter";
52  }
53 
54  ostringstream ostr;
55  ostr << "Filter with subjet_def = ";
56  if (_Rfiltfunc) {
57  ostr << "Cambridge/Aachen algorithm with dynamic Rfilt"
58  << " (recomb. scheme deduced from jet, or E-scheme if not unique)";
59  } else if (_Rfilt > 0) {
60  ostr << "Cambridge/Aachen algorithm with Rfilt = "
61  << _Rfilt
62  << " (recomb. scheme deduced from jet, or E-scheme if not unique)";
63  } else {
64  ostr << _subjet_def.description();
65  }
66  ostr<< ", selection " << _selector.description();
67  if (_subtractor) {
68  ostr << ", subtractor: " << _subtractor->description();
69  } else if (_rho != 0) {
70  ostr << ", subtracting with rho = " << _rho;
71  }
72  return ostr.str();
73 }
74 
75 
76 // core functions
77 //----------------------------------------------------------------------
78 
79 // return a vector of subjets, which are the ones that would be kept
80 // by the filtering
81 PseudoJet Filter::result(const PseudoJet &jet) const {
82  if (!_initialised){
83  //Q: do we throw or do we return an empty PJ?
84  throw Error("uninitialised Filter");
85  }
86 
87  // start by getting the list of subjets (including a list of sanity
88  // checks)
89  // NB: subjets is empty to begin with (see the comment for
90  // _set_filtered_elements_cafilt)
91  vector<PseudoJet> subjets;
92  //JetDefinition subjet_def;
93  bool ca_optimised = _set_filtered_elements(jet, subjets);
94 
95  // apply subtraction if needed:
96  if (_subtractor){
97  subjets = (*_subtractor)(subjets);
98  } else if (_rho!=0){
99  if (subjets.size()>0){
100  const ClusterSequenceAreaBase *csab = subjets[0].validated_csab();
101  for (unsigned int i=0;i<subjets.size();i++){
102  subjets[i]=csab->subtracted_jet(subjets[i], _rho);
103  }
104  }
105  }
106 
107  // now build the vector of kept and rejected subjets
108  vector<PseudoJet> kept, rejected;
109  // Note that in the following line we make a copy of the _selector
110  // to avoid issues with needing a mutable _selector
111  Selector selector_copy = _selector;
112  if (selector_copy.takes_reference()) selector_copy.set_reference(jet);
113  selector_copy.sift(subjets, kept, rejected);
114 
115  // gather the info under the form of a PseudoJet
116  return _finalise(jet, kept, rejected, ca_optimised);
117 }
118 
119 
120 // sets filtered_elements to be all the subjets on which filtering will work
121 //
122 // return true when the subjets have been optained using teh optimised
123 // method for C/A
124 bool Filter::_set_filtered_elements(const PseudoJet & jet,
125  vector<PseudoJet> & filtered_elements) const {
126  // create the recluster instance
127  Recluster recluster;
128  if ((_Rfilt>=0) || (_Rfiltfunc))
129  recluster = Recluster(cambridge_algorithm, (_Rfiltfunc) ? (*_Rfiltfunc)(jet) : _Rfilt, Recluster::keep_all);
130  else
131  recluster = Recluster(_subjet_def, false, Recluster::keep_all);
132 
133  // get the subjets
134  //JetDefinition subjet_def;
135  return recluster.get_new_jets_and_def(jet, filtered_elements);
136 }
137 
138 // gather the information about what is kept and rejected under the
139 // form of a PseudoJet with a special ClusterSequenceInfo
140 PseudoJet Filter::_finalise(const PseudoJet & /*jet*/,
141  vector<PseudoJet> & kept,
142  vector<PseudoJet> & rejected,
143  bool ca_optimisation_used) const {
144  PseudoJet filtered_jet;
145 
146  if (kept.size()+rejected.size()>0){
147  // figure out which recombiner to use
148  const JetDefinition::Recombiner &rec = (kept.size()>0)
149  ? *(kept[0].associated_cs()->jet_def().recombiner())
150  : *(rejected[0].associated_cs()->jet_def().recombiner());
151 
152  // create an appropriate structure and transfer the info to it
153  filtered_jet = join<StructureType>(kept, rec);
154  } else {
155  filtered_jet = join<StructureType>(kept);
156  }
157  StructureType *fs = (StructureType*) filtered_jet.structure_non_const_ptr();
158  fs->_rejected = rejected;
159 
160  // if we've used C/A optimisation, we need to get rid of the area
161  // information if it comes from a non-explicit-ghost clustering.
162  // (because in that case it can be erroneous due the lack of
163  // information about empty areas)
164  if ((ca_optimisation_used) && (kept.size()+rejected.size()>0)){
165  bool has_non_explicit_ghost_area = (kept.size()>0)
166  ? (kept[0].has_area() && (!(kept[0].validated_csab()->has_explicit_ghosts())))
167  : (rejected[0].has_area() && (!(rejected[0].validated_csab()->has_explicit_ghosts())));
168  if (has_non_explicit_ghost_area)
169  fs->discard_area();
170  }
171 
172  return filtered_jet;
173 }
174 
175 
176 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
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
bool takes_reference() const
returns true if this can be applied jet by jet
Definition: Selector.hh:287
bool get_new_jets_and_def(const PseudoJet &input_jet, std::vector< PseudoJet > &output_jets) const
A lower-level method that does the actual work of reclustering the input jet.
Definition: Recluster.cc:108
base class that sets interface for extensions of ClusterSequence that provide information about the a...
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:47
const Selector & set_reference(const PseudoJet &reference)
set the reference jet for this Selector
Definition: Selector.hh:292
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition: Selector.hh:149
Recluster a jet's constituents with a new jet definition.
Definition: Recluster.hh:73
PseudoJet subtracted_jet(const PseudoJet &jet, const double rho) const
return a subtracted jet, using area_4vector, given rho
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67