FastJet 3.0.2
Filter.hh
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__
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends