[d7d2da3] | 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__
|
---|