FastJet 3.0.2
|
00001 #ifndef __FASTJET_TOOLS_FILTER_HH__ 00002 #define __FASTJET_TOOLS_FILTER_HH__ 00003 00004 //STARTHEADER 00005 // $Id: Filter.hh 2694 2011-11-14 22:27:51Z salam $ 00006 // 00007 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 00008 // 00009 //---------------------------------------------------------------------- 00010 // This file is part of FastJet. 00011 // 00012 // FastJet is free software; you can redistribute it and/or modify 00013 // it under the terms of the GNU General Public License as published by 00014 // the Free Software Foundation; either version 2 of the License, or 00015 // (at your option) any later version. 00016 // 00017 // The algorithms that underlie FastJet have required considerable 00018 // development and are described in hep-ph/0512210. If you use 00019 // FastJet as part of work towards a scientific publication, please 00020 // include a citation to the FastJet paper. 00021 // 00022 // FastJet is distributed in the hope that it will be useful, 00023 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00024 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00025 // GNU General Public License for more details. 00026 // 00027 // You should have received a copy of the GNU General Public License 00028 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00029 //---------------------------------------------------------------------- 00030 //ENDHEADER 00031 00032 #include <fastjet/ClusterSequence.hh> 00033 #include <fastjet/Selector.hh> 00034 #include <fastjet/CompositeJetStructure.hh> // to derive the FilterStructure from CompositeJetStructure 00035 #include <fastjet/tools/Transformer.hh> // to derive Filter from Transformer 00036 #include <iostream> 00037 #include <string> 00038 00039 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00040 00041 // fwd declarations 00042 class Filter; 00043 class FilterStructure; 00044 00045 //---------------------------------------------------------------------- 00046 /// @ingroup tools_generic 00047 /// \class Filter 00048 /// Class that helps perform filtering (Butterworth, Davison, Rubin 00049 /// and Salam, arXiv:0802.2470) and trimming (Krohn, Thaler and Wang, 00050 /// arXiv:0912.1342) on jets, optionally in conjunction with 00051 /// subtraction (Cacciari and Salam, arXiv:0707.1378). 00052 /// 00053 /// For example, to apply filtering that reclusters a jet's 00054 /// constituents with the Cambridge/Aachen jet algorithm with R=0.3 00055 /// and then selects the 3 hardest subjets, one can use the following 00056 /// code: 00057 /// \code 00058 /// Filter filter(JetDefinition(cambridge_algorithm, 0.3), SelectorNHardest(3)); 00059 /// PseudoJet filtered_jet = filter(original_jet); 00060 /// \endcode 00061 /// 00062 /// To obtain trimming, involving for example the selection of all 00063 /// subjets carrying at least 3% of the original jet's pt, the 00064 /// selector would be replaced by SelectorPtFractionMin(0.03). 00065 /// 00066 /// To additionally perform subtraction on the subjets prior to 00067 /// selection, either include a 3rd argument specifying the background 00068 /// density rho, or call the set_subtractor(...) member function. If 00069 /// subtraction is requested, the original jet must be the result of a 00070 /// clustering with active area with explicit ghosts support or a 00071 /// merging of such pieces. 00072 /// 00073 /// The information on the subjets that were kept and rejected can be 00074 /// obtained using: 00075 /// \code 00076 /// vector<PseudoJet> kept_subjets = filtered_jet.pieces(); 00077 /// vector<PseudoJet> rejected_subjets = filtered_jet.structure_of<Filter>().rejected(); 00078 /// \endcode 00079 /// 00080 /// \section impl Implementation Note 00081 /// 00082 /// If the original jet was defined with the Cambridge/Aachen 00083 /// algorithm (or is made of pieces each of which comes from the C/A 00084 /// alg) and the filtering definition is C/A, then the filter does not 00085 /// rerun the C/A algorithm on the constituents, but instead makes use 00086 /// of the existent C/A cluster sequence in the original jet. This 00087 /// increases the speed of the filter. 00088 /// 00089 /// See also \subpage Example11 for a further usage example. 00090 /// 00091 /// Support for areas, reuse of C/A cluster sequences, etc., 00092 /// considerably complicates the implementation of Filter. For an 00093 /// explanation of how a simpler filter might be coded, see the 00094 /// "User-defined transformers" appendix of the manual. 00095 class Filter : public Transformer{ 00096 public: 00097 /// trivial ctor 00098 /// Note: this is just for derived classes 00099 /// a Filter initialised through this constructor will not work! 00100 Filter() : _Rfiltfunc(0){}; 00101 00102 /// define a filter that decomposes a jet into subjets using a 00103 /// generic JetDefinition and then keeps only a subset of these 00104 /// subjets according to a Selector. Optionally, each subjet may be 00105 /// internally bakground-subtracted prior to selection. 00106 /// 00107 /// \param subjet_def the jet definition applied to obtain the subjets 00108 /// \param selector the Selector applied to compute the kept subjets 00109 /// \param rho if non-zero, backgruond-subtract each subjet befor selection 00110 /// 00111 /// Note: internal subtraction only applies on jets that are 00112 /// obtained with a cluster sequence with area support and explicit 00113 /// ghosts 00114 Filter(JetDefinition subjet_def, Selector selector, double rho = 0.0) : 00115 _subjet_def(subjet_def), _Rfiltfunc(0), _Rfilt(-1), _selector(selector), _rho(rho), _subtractor(0) {} 00116 00117 /// Same as the full constructor (see above) but just specifying the radius 00118 /// By default, Cambridge-Aachen is used 00119 /// If the jet (or all its pieces) is obtained with a non-default 00120 /// recombiner, that one will be used 00121 /// \param Rfilt the filtering radius 00122 Filter(double Rfilt, Selector selector, double rho = 0.0) : 00123 _Rfiltfunc(0), _Rfilt(Rfilt), _selector(selector), _rho(rho), _subtractor(0) { 00124 if (_Rfilt<0) 00125 throw Error("Attempt to create a Filter with a negative filtering radius"); 00126 } 00127 00128 /// Same as the full constructor (see above) but just specifying a 00129 /// filtering radius that will depend on the jet being filtered 00130 /// As for the previous case, Cambridge-Aachen is used 00131 /// If the jet (or all its pieces) is obtained with a non-default 00132 /// recombiner, that one will be used 00133 /// \param Rfilt_func the filtering radius function of a PseudoJet 00134 Filter(FunctionOfPseudoJet<double> *Rfilt_func, Selector selector, double rho = 0.0) : 00135 _Rfiltfunc(Rfilt_func), _Rfilt(-1), _selector(selector), _rho(rho), _subtractor(0) {} 00136 00137 /// default dtor 00138 virtual ~Filter(){}; 00139 00140 /// Set a subtractor that is applied to all individual subjets before 00141 /// deciding which ones to keep. It takes precedence over a non-zero rho. 00142 void set_subtractor(const Transformer * subtractor) {_subtractor = subtractor;} 00143 00144 /// runs the filtering and sets kept and rejected to be the jets of interest 00145 /// (with non-zero rho, they will have been subtracted). 00146 /// 00147 /// \param jet the jet that gets filtered 00148 /// \return the filtered jet 00149 virtual PseudoJet result(const PseudoJet & jet) const; 00150 00151 /// class description 00152 virtual std::string description() const; 00153 00154 // the type of the associated structure 00155 typedef FilterStructure StructureType; 00156 00157 private: 00158 /// Sets filtered_elements to be all the subjets on which filtering will work. 00159 /// It also sets the subjet_def to be used in joining things (the bit of 00160 /// subjet def that is of interest for later is the recombiner). 00161 void _set_filtered_elements(const PseudoJet & jet, 00162 std::vector<PseudoJet> & filtered_elements, 00163 JetDefinition & subjet_def, 00164 bool & discard_area) const; 00165 00166 /// set the filtered elements in the simple case of C/A+C/A 00167 void _set_filtered_elements_cafilt(const PseudoJet & jet, 00168 std::vector<PseudoJet> & filtered_elements, 00169 double Rfilt) const; 00170 00171 /// set the filtered elements in the generic re-clustering case 00172 void _set_filtered_elements_generic(const PseudoJet & jet, 00173 std::vector<PseudoJet> & filtered_elements, 00174 const JetDefinition & subjet_def, 00175 bool do_areas) const; 00176 00177 /// gather the information about what is kept and rejected under the 00178 /// form of a PseudoJet with a special ClusterSequenceInfo 00179 PseudoJet _finalise(const PseudoJet & jet, 00180 std::vector<PseudoJet> & kept, 00181 std::vector<PseudoJet> & rejected, 00182 const JetDefinition & subjet_def, 00183 const bool discard_area) const; 00184 00185 // a series of checks 00186 //-------------------------------------------------------------------- 00187 /// get the pieces down to the fundamental pieces 00188 bool _get_all_pieces(const PseudoJet &jet, std::vector<PseudoJet> &all_pieces) const; 00189 00190 /// get the common recombiner to all pieces (NULL if none) 00191 const JetDefinition::Recombiner* _get_common_recombiner(const std::vector<PseudoJet> &all_pieces) const; 00192 00193 /// check if one can apply the simplified trick for C/A subjets 00194 bool _check_ca(const std::vector<PseudoJet> &all_pieces) const; 00195 00196 /// check if the jet (or all its pieces) have explicit ghosts 00197 /// (assuming the jet has area support 00198 /// 00199 /// Note that if the jet has an associated cluster sequence that is no 00200 /// longer valid, an error will be thrown 00201 bool _check_explicit_ghosts(const std::vector<PseudoJet> &all_pieces) const; 00202 00203 bool _uses_subtraction() const {return (_subtractor || _rho != 0);} 00204 00205 JetDefinition _subjet_def; ///< the jet definition to use to extract the subjets 00206 FunctionOfPseudoJet<double> *_Rfiltfunc; 00207 ///< a dynamic filtering radius function of the jet being filtered 00208 double _Rfilt; ///< a constant specifying the subjet radius (with C/A) 00209 Selector _selector; ///< the subjet selection criterium 00210 double _rho; ///< the background density (used for subtraction when possible) 00211 const Transformer * _subtractor; ///< for subtracting bkgd density from subjets 00212 }; 00213 00214 00215 00216 //---------------------------------------------------------------------- 00217 /// @ingroup tools_generic 00218 /// \class FilterStructure 00219 /// Class to contain structure information for a filtered jet. 00220 class FilterStructure : public CompositeJetStructure { 00221 public: 00222 /// constructor from an original ClusterSequenceInfo 00223 /// We just share the original ClusterSequenceWrapper and initialise 00224 /// the rest 00225 FilterStructure(const std::vector<PseudoJet> & pieces_in, 00226 const JetDefinition::Recombiner *rec = 0) 00227 : CompositeJetStructure(pieces_in, rec){} 00228 00229 /// virtual dtor to allow further overloading 00230 virtual ~FilterStructure(){} 00231 00232 /// description 00233 virtual std::string description() const { return "Filtered PseudoJet"; } 00234 00235 //------------------------------------------------------------------ 00236 /// @name The filter-specific information 00237 //------------------------------------------------------------------ 00238 00239 // /// returns the original jet (the first of the original jets 00240 // /// if you filtered a collection of jets) 00241 // const PseudoJet & original() const {return _original_jet;} 00242 00243 /// returns the subjets that were not kept during the filtering procedure 00244 /// (subtracted if the filter requests it, and valid in the original cs) 00245 const std::vector<PseudoJet> & rejected() const {return _rejected;} 00246 00247 friend class Filter; // allow the filter to change the protected/private members 00248 00249 protected: 00250 // PseudoJet _original_jet; ///< the original jet 00251 std::vector<PseudoJet> _rejected; ///< the subjets rejected by the filter 00252 }; 00253 00254 00255 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh 00256 00257 #endif // __FASTJET_TOOLS_FILTER_HH__