FastJet  3.0.2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Filter.hh
1 #ifndef __FASTJET_TOOLS_FILTER_HH__
2 #define __FASTJET_TOOLS_FILTER_HH__
3 
4 //STARTHEADER
5 // $Id: Filter.hh 2694 2011-11-14 22:27:51Z salam $
6 //
7 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
8 //
9 //----------------------------------------------------------------------
10 // This file is part of FastJet.
11 //
12 // FastJet is free software; you can redistribute it and/or modify
13 // it under the terms of the GNU General Public License as published by
14 // the Free Software Foundation; either version 2 of the License, or
15 // (at your option) any later version.
16 //
17 // The algorithms that underlie FastJet have required considerable
18 // development and are described in hep-ph/0512210. If you use
19 // FastJet as part of work towards a scientific publication, please
20 // include a citation to the FastJet paper.
21 //
22 // FastJet is distributed in the hope that it will be useful,
23 // but WITHOUT ANY WARRANTY; without even the implied warranty of
24 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 // GNU General Public License for more details.
26 //
27 // You should have received a copy of the GNU General Public License
28 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
29 //----------------------------------------------------------------------
30 //ENDHEADER
31 
32 #include <fastjet/ClusterSequence.hh>
33 #include <fastjet/Selector.hh>
34 #include <fastjet/CompositeJetStructure.hh> // to derive the FilterStructure from CompositeJetStructure
35 #include <fastjet/tools/Transformer.hh> // to derive Filter from Transformer
36 #include <iostream>
37 #include <string>
38 
39 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
40 
41 // fwd declarations
42 class Filter;
43 class FilterStructure;
44 
45 //----------------------------------------------------------------------
46 /// @ingroup tools_generic
47 /// \class Filter
48 /// Class that helps perform filtering (Butterworth, Davison, Rubin
49 /// and Salam, arXiv:0802.2470) and trimming (Krohn, Thaler and Wang,
50 /// arXiv:0912.1342) on jets, optionally in conjunction with
51 /// subtraction (Cacciari and Salam, arXiv:0707.1378).
52 ///
53 /// For example, to apply filtering that reclusters a jet's
54 /// constituents with the Cambridge/Aachen jet algorithm with R=0.3
55 /// and then selects the 3 hardest subjets, one can use the following
56 /// code:
57 /// \code
58 /// Filter filter(JetDefinition(cambridge_algorithm, 0.3), SelectorNHardest(3));
59 /// PseudoJet filtered_jet = filter(original_jet);
60 /// \endcode
61 ///
62 /// To obtain trimming, involving for example the selection of all
63 /// subjets carrying at least 3% of the original jet's pt, the
64 /// selector would be replaced by SelectorPtFractionMin(0.03).
65 ///
66 /// To additionally perform subtraction on the subjets prior to
67 /// selection, either include a 3rd argument specifying the background
68 /// density rho, or call the set_subtractor(...) member function. If
69 /// subtraction is requested, the original jet must be the result of a
70 /// clustering with active area with explicit ghosts support or a
71 /// merging of such pieces.
72 ///
73 /// The information on the subjets that were kept and rejected can be
74 /// obtained using:
75 /// \code
76 /// vector<PseudoJet> kept_subjets = filtered_jet.pieces();
77 /// vector<PseudoJet> rejected_subjets = filtered_jet.structure_of<Filter>().rejected();
78 /// \endcode
79 ///
80 /// \section impl Implementation Note
81 ///
82 /// If the original jet was defined with the Cambridge/Aachen
83 /// algorithm (or is made of pieces each of which comes from the C/A
84 /// alg) and the filtering definition is C/A, then the filter does not
85 /// rerun the C/A algorithm on the constituents, but instead makes use
86 /// of the existent C/A cluster sequence in the original jet. This
87 /// increases the speed of the filter.
88 ///
89 /// See also \subpage Example11 for a further usage example.
90 ///
91 /// Support for areas, reuse of C/A cluster sequences, etc.,
92 /// considerably complicates the implementation of Filter. For an
93 /// explanation of how a simpler filter might be coded, see the
94 /// "User-defined transformers" appendix of the manual.
95 class Filter : public Transformer{
96 public:
97  /// trivial ctor
98  /// Note: this is just for derived classes
99  /// a Filter initialised through this constructor will not work!
100  Filter() : _Rfiltfunc(0){};
101 
102  /// define a filter that decomposes a jet into subjets using a
103  /// generic JetDefinition and then keeps only a subset of these
104  /// subjets according to a Selector. Optionally, each subjet may be
105  /// internally bakground-subtracted prior to selection.
106  ///
107  /// \param subjet_def the jet definition applied to obtain the subjets
108  /// \param selector the Selector applied to compute the kept subjets
109  /// \param rho if non-zero, backgruond-subtract each subjet befor selection
110  ///
111  /// Note: internal subtraction only applies on jets that are
112  /// obtained with a cluster sequence with area support and explicit
113  /// ghosts
114  Filter(JetDefinition subjet_def, Selector selector, double rho = 0.0) :
115  _subjet_def(subjet_def), _Rfiltfunc(0), _Rfilt(-1), _selector(selector), _rho(rho), _subtractor(0) {}
116 
117  /// Same as the full constructor (see above) but just specifying the radius
118  /// By default, Cambridge-Aachen is used
119  /// If the jet (or all its pieces) is obtained with a non-default
120  /// recombiner, that one will be used
121  /// \param Rfilt the filtering radius
122  Filter(double Rfilt, Selector selector, double rho = 0.0) :
123  _Rfiltfunc(0), _Rfilt(Rfilt), _selector(selector), _rho(rho), _subtractor(0) {
124  if (_Rfilt<0)
125  throw Error("Attempt to create a Filter with a negative filtering radius");
126  }
127 
128  /// Same as the full constructor (see above) but just specifying a
129  /// filtering radius that will depend on the jet being filtered
130  /// As for the previous case, Cambridge-Aachen is used
131  /// If the jet (or all its pieces) is obtained with a non-default
132  /// recombiner, that one will be used
133  /// \param Rfilt_func the filtering radius function of a PseudoJet
134  Filter(FunctionOfPseudoJet<double> *Rfilt_func, Selector selector, double rho = 0.0) :
135  _Rfiltfunc(Rfilt_func), _Rfilt(-1), _selector(selector), _rho(rho), _subtractor(0) {}
136 
137  /// default dtor
138  virtual ~Filter(){};
139 
140  /// Set a subtractor that is applied to all individual subjets before
141  /// deciding which ones to keep. It takes precedence over a non-zero rho.
142  void set_subtractor(const Transformer * subtractor) {_subtractor = subtractor;}
143 
144  /// runs the filtering and sets kept and rejected to be the jets of interest
145  /// (with non-zero rho, they will have been subtracted).
146  ///
147  /// \param jet the jet that gets filtered
148  /// \return the filtered jet
149  virtual PseudoJet result(const PseudoJet & jet) const;
150 
151  /// class description
152  virtual std::string description() const;
153 
154  // the type of the associated structure
156 
157 private:
158  /// Sets filtered_elements to be all the subjets on which filtering will work.
159  /// It also sets the subjet_def to be used in joining things (the bit of
160  /// subjet def that is of interest for later is the recombiner).
161  void _set_filtered_elements(const PseudoJet & jet,
162  std::vector<PseudoJet> & filtered_elements,
163  JetDefinition & subjet_def,
164  bool & discard_area) const;
165 
166  /// set the filtered elements in the simple case of C/A+C/A
167  void _set_filtered_elements_cafilt(const PseudoJet & jet,
168  std::vector<PseudoJet> & filtered_elements,
169  double Rfilt) const;
170 
171  /// set the filtered elements in the generic re-clustering case
172  void _set_filtered_elements_generic(const PseudoJet & jet,
173  std::vector<PseudoJet> & filtered_elements,
174  const JetDefinition & subjet_def,
175  bool do_areas) const;
176 
177  /// gather the information about what is kept and rejected under the
178  /// form of a PseudoJet with a special ClusterSequenceInfo
179  PseudoJet _finalise(const PseudoJet & jet,
180  std::vector<PseudoJet> & kept,
181  std::vector<PseudoJet> & rejected,
182  const JetDefinition & subjet_def,
183  const bool discard_area) const;
184 
185  // a series of checks
186  //--------------------------------------------------------------------
187  /// get the pieces down to the fundamental pieces
188  bool _get_all_pieces(const PseudoJet &jet, std::vector<PseudoJet> &all_pieces) const;
189 
190  /// get the common recombiner to all pieces (NULL if none)
191  const JetDefinition::Recombiner* _get_common_recombiner(const std::vector<PseudoJet> &all_pieces) const;
192 
193  /// check if one can apply the simplified trick for C/A subjets
194  bool _check_ca(const std::vector<PseudoJet> &all_pieces) const;
195 
196  /// check if the jet (or all its pieces) have explicit ghosts
197  /// (assuming the jet has area support
198  ///
199  /// Note that if the jet has an associated cluster sequence that is no
200  /// longer valid, an error will be thrown
201  bool _check_explicit_ghosts(const std::vector<PseudoJet> &all_pieces) const;
202 
203  bool _uses_subtraction() const {return (_subtractor || _rho != 0);}
204 
205  JetDefinition _subjet_def; ///< the jet definition to use to extract the subjets
206  FunctionOfPseudoJet<double> *_Rfiltfunc;
207  ///< a dynamic filtering radius function of the jet being filtered
208  double _Rfilt; ///< a constant specifying the subjet radius (with C/A)
209  Selector _selector; ///< the subjet selection criterium
210  double _rho; ///< the background density (used for subtraction when possible)
211  const Transformer * _subtractor; ///< for subtracting bkgd density from subjets
212 };
213 
214 
215 
216 //----------------------------------------------------------------------
217 /// @ingroup tools_generic
218 /// \class FilterStructure
219 /// Class to contain structure information for a filtered jet.
221 public:
222  /// constructor from an original ClusterSequenceInfo
223  /// We just share the original ClusterSequenceWrapper and initialise
224  /// the rest
225  FilterStructure(const std::vector<PseudoJet> & pieces_in,
226  const JetDefinition::Recombiner *rec = 0)
227  : CompositeJetStructure(pieces_in, rec){}
228 
229  /// virtual dtor to allow further overloading
230  virtual ~FilterStructure(){}
231 
232  /// description
233  virtual std::string description() const { return "Filtered PseudoJet"; }
234 
235  //------------------------------------------------------------------
236  /// @name The filter-specific information
237  //------------------------------------------------------------------
238 
239 // /// returns the original jet (the first of the original jets
240 // /// if you filtered a collection of jets)
241 // const PseudoJet & original() const {return _original_jet;}
242 
243  /// returns the subjets that were not kept during the filtering procedure
244  /// (subtracted if the filter requests it, and valid in the original cs)
245  const std::vector<PseudoJet> & rejected() const {return _rejected;}
246 
247  friend class Filter; // allow the filter to change the protected/private members
248 
249 protected:
250 // PseudoJet _original_jet; ///< the original jet
251  std::vector<PseudoJet> _rejected; ///< the subjets rejected by the filter
252 };
253 
254 
255 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
256 
257 #endif // __FASTJET_TOOLS_FILTER_HH__