FastJet 3.0.0
|
00001 #ifndef __FASTJET_TOOLS_FILTER_HH__ 00002 #define __FASTJET_TOOLS_FILTER_HH__ 00003 00004 //STARTHEADER 00005 // $Id: Filter.hh 2616 2011-09-30 18:03:40Z 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 class Filter : public Transformer{ 00092 public: 00093 /// trivial ctor 00094 /// Note: this is just for derived classes 00095 /// a Filter initialised through this constructor will not work! 00096 Filter() : _Rfiltfunc(0){}; 00097 00098 /// define a filter that decomposes a jet into subjets using a 00099 /// generic JetDefinition and then keeps only a subset of these 00100 /// subjets according to a Selector. Optionally, each subjet may be 00101 /// internally bakground-subtracted prior to selection. 00102 /// 00103 /// \param subjet_def the jet definition applied to obtain the subjets 00104 /// \param selector the Selector applied to compute the kept subjets 00105 /// \param rho if non-zero, backgruond-subtract each subjet befor selection 00106 /// 00107 /// Note: internal subtraction only applies on jets that are 00108 /// obtained with a cluster sequence with area support and explicit 00109 /// ghosts 00110 Filter(JetDefinition subjet_def, Selector selector, double rho = 0.0) : 00111 _subjet_def(subjet_def), _Rfiltfunc(0), _Rfilt(-1), _selector(selector), _rho(rho), _subtractor(0) {} 00112 00113 /// Same as the full constructor (see above) but just specifying the radius 00114 /// By default, Cambridge-Aachen is used 00115 /// If the jet (or all its pieces) is obtained with a non-default 00116 /// recombiner, that one will be used 00117 /// \param Rfilt the filtering radius 00118 Filter(double Rfilt, Selector selector, double rho = 0.0) : 00119 _Rfiltfunc(0), _Rfilt(Rfilt), _selector(selector), _rho(rho), _subtractor(0) { 00120 if (_Rfilt<0) 00121 throw Error("Attempt to create a Filter with a negative filtering radius"); 00122 } 00123 00124 /// Same as the full constructor (see above) but just specifying a 00125 /// filtering radius that will depend on the jet being filtered 00126 /// As for the previous case, Cambridge-Aachen is used 00127 /// If the jet (or all its pieces) is obtained with a non-default 00128 /// recombiner, that one will be used 00129 /// \param Rfilt_func the filtering radius function of a PseudoJet 00130 Filter(FunctionOfPseudoJet<double> *Rfilt_func, Selector selector, double rho = 0.0) : 00131 _Rfiltfunc(Rfilt_func), _Rfilt(-1), _selector(selector), _rho(rho), _subtractor(0) {} 00132 00133 /// default dtor 00134 virtual ~Filter(){}; 00135 00136 /// Set a subtractor that is applied to all individual subjets before 00137 /// deciding which ones to keep. It takes precedence over a non-zero rho. 00138 void set_subtractor(const Transformer * subtractor) {_subtractor = subtractor;} 00139 00140 /// runs the filtering and sets kept and rejected to be the jets of interest 00141 /// (with non-zero rho, they will have been subtracted). 00142 /// 00143 /// \param jet the jet that gets filtered 00144 /// \return the filtered jet 00145 virtual PseudoJet result(const PseudoJet & jet) const; 00146 00147 /// class description 00148 virtual std::string description() const; 00149 00150 // the type of the associated structure 00151 typedef FilterStructure StructureType; 00152 00153 private: 00154 /// sets filtered_elements to be all the subjets on which filtering will work 00155 /// [NB: this routine is work in progress as part of a transition to a Filter 00156 /// that also works on jet collections] 00157 void _set_filtered_elements(const PseudoJet & jet, 00158 std::vector<PseudoJet> & filtered_elements, 00159 bool & discard_area) const; 00160 00161 /// set the filtered elements in the simple case of C/A+C/A 00162 void _set_filtered_elements_cafilt(const PseudoJet & jet, 00163 std::vector<PseudoJet> & filtered_elements, 00164 double Rfilt) const; 00165 00166 /// set the filtered elements in the generic re-clustering case 00167 void _set_filtered_elements_generic(const PseudoJet & jet, 00168 std::vector<PseudoJet> & filtered_elements) const; 00169 00170 /// gather the information about what is kept and rejected under the 00171 /// form of a PseudoJet with a special ClusterSequenceInfo 00172 PseudoJet _finalise(const PseudoJet & jet, 00173 std::vector<PseudoJet> & kept, 00174 std::vector<PseudoJet> & rejected, 00175 const bool discard_area) const; 00176 00177 // a series of checks 00178 //-------------------------------------------------------------------- 00179 /// get the pieces down to the fundamental pieces 00180 bool _get_all_pieces(const PseudoJet &jet, std::vector<PseudoJet> &all_pieces) const; 00181 00182 /// get the common recombiner to all pieces (NULL if none) 00183 const JetDefinition::Recombiner* _get_common_recombiner() const; 00184 00185 /// check if one can apply the simplified trick for C/A subjets 00186 bool _check_ca() const; 00187 00188 /// check if the jet (or all its pieces) have explicit ghosts 00189 /// (assuming the jet has area support 00190 /// 00191 /// Note that if the jet has an associated cluster sequence that is no 00192 /// longer valid, an error will be thrown 00193 bool _check_explicit_ghosts() const; 00194 00195 bool _uses_subtraction() const {return (_subtractor || _rho != 0);} 00196 00197 mutable JetDefinition _subjet_def; 00198 ///< the jet definition to use to extract the subjets 00199 FunctionOfPseudoJet<double> *_Rfiltfunc; 00200 ///< a dynamic filtering radius function of the jet being filtered 00201 double _Rfilt; ///< a constant specifying the subjet radius (with C/A) 00202 mutable Selector _selector; ///< the subjet selection criterium 00203 double _rho; ///< the background density (used for subtraction when possible) 00204 const Transformer * _subtractor; ///< for subtracting bkgd density from subjets 00205 00206 protected: 00207 // internal useful variables 00208 mutable std::vector<PseudoJet> all_pieces; 00209 }; 00210 00211 00212 00213 //---------------------------------------------------------------------- 00214 /// @ingroup tools_generic 00215 /// \class FilterStructure 00216 /// Class to contain structure information for a filtered jet. 00217 class FilterStructure : public CompositeJetStructure { 00218 public: 00219 /// constructor from an original ClusterSequenceInfo 00220 /// We just share the original ClusterSequenceWrapper and initialise 00221 /// the rest 00222 FilterStructure(const std::vector<PseudoJet> & pieces, 00223 const JetDefinition::Recombiner *rec = 0) 00224 : CompositeJetStructure(pieces, rec){} 00225 00226 /// virtual dtor to allow further overloading 00227 virtual ~FilterStructure(){} 00228 00229 /// description 00230 virtual std::string description() const { return "Filtered PseudoJet"; } 00231 00232 //------------------------------------------------------------------ 00233 /// @name The filter-specific information 00234 //------------------------------------------------------------------ 00235 00236 // /// returns the original jet (the first of the original jets 00237 // /// if you filtered a collection of jets) 00238 // const PseudoJet & original() const {return _original_jet;} 00239 00240 /// returns the subjets that were not kept during the filtering procedure 00241 /// (subtracted if the filter requests it, and valid in the original cs) 00242 const std::vector<PseudoJet> & rejected() const {return _rejected;} 00243 00244 friend class Filter; // allow the filter to change the protected/private members 00245 00246 protected: 00247 // PseudoJet _original_jet; ///< the original jet 00248 std::vector<PseudoJet> _rejected; ///< the subjets rejected by the filter 00249 }; 00250 00251 00252 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh 00253 00254 #endif // __FASTJET_TOOLS_FILTER_HH__