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 |
---|
155 | typedef FilterStructure StructureType; |
---|
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. |
---|
220 | class FilterStructure : public CompositeJetStructure { |
---|
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__ |
---|