FastJet 3.0.2
Filter.cc
00001 //STARTHEADER
00002 // $Id: Filter.cc 2695 2011-11-14 22:33:07Z salam $
00003 //
00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00005 //
00006 //----------------------------------------------------------------------
00007 // This file is part of FastJet.
00008 //
00009 //  FastJet is free software; you can redistribute it and/or modify
00010 //  it under the terms of the GNU General Public License as published by
00011 //  the Free Software Foundation; either version 2 of the License, or
00012 //  (at your option) any later version.
00013 //
00014 //  The algorithms that underlie FastJet have required considerable
00015 //  development and are described in hep-ph/0512210. If you use
00016 //  FastJet as part of work towards a scientific publication, please
00017 //  include a citation to the FastJet paper.
00018 //
00019 //  FastJet is distributed in the hope that it will be useful,
00020 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00021 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00022 //  GNU General Public License for more details.
00023 //
00024 //  You should have received a copy of the GNU General Public License
00025 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00026 //----------------------------------------------------------------------
00027 //ENDHEADER
00028 
00029 #include "fastjet/tools/Filter.hh"
00030 #include <fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh>
00031 #include <cassert>
00032 #include <algorithm>
00033 #include <sstream>
00034 #include <typeinfo>
00035 
00036 using namespace std;
00037 
00038 
00039 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00040 
00041 //----------------------------------------------------------------------
00042 // Filter class implementation
00043 //----------------------------------------------------------------------
00044 
00045 // class description
00046 string Filter::description() const {
00047   ostringstream ostr;
00048   ostr << "Filter with subjet_def = ";
00049   if (_Rfiltfunc) {
00050     ostr << "Cambridge/Aachen algorithm with dynamic Rfilt"
00051          << " (recomb. scheme deduced from jet, or E-scheme if not unique)";
00052   } else if (_Rfilt > 0) {
00053     ostr << "Cambridge/Aachen algorithm with Rfilt = " 
00054          << _Rfilt 
00055          << " (recomb. scheme deduced from jet, or E-scheme if not unique)";
00056   } else {
00057     ostr << _subjet_def.description();
00058   }
00059   ostr<< ", selection " << _selector.description();
00060   if (_subtractor) {
00061     ostr << ", subtractor: " << _subtractor->description();
00062   } else if (_rho != 0) {
00063     ostr << ", subtracting with rho = " << _rho;
00064   }
00065   return ostr.str();
00066 }
00067 
00068 
00069 // core functions
00070 //----------------------------------------------------------------------
00071 
00072 // return a vector of subjets, which are the ones that would be kept
00073 // by the filtering
00074 PseudoJet Filter::result(const PseudoJet &jet) const {
00075   // start by getting the list of subjets (including a list of sanity
00076   // checks)
00077   // NB: subjets is empty to begin with (see the comment for
00078   //     _set_filtered_elements_cafilt)
00079   vector<PseudoJet> subjets; 
00080   JetDefinition subjet_def;
00081   bool discard_area;
00082   // NB: on return, subjet_def is set to the jet definition actually
00083   //     used (so that we can make use of its recombination scheme 
00084   //     when joining the jets to be kept).
00085   _set_filtered_elements(jet, subjets, subjet_def, discard_area);
00086 
00087   // now build the vector of kept and rejected subjets
00088   vector<PseudoJet> kept, rejected;
00089   // Note that in the following line we make a copy of the _selector
00090   // to avoid issues with needing a mutable _selector
00091   Selector selector_copy = _selector;
00092   if (selector_copy.takes_reference()) selector_copy.set_reference(jet);
00093   selector_copy.sift(subjets, kept, rejected);
00094 
00095   // gather the info under the form of a PseudoJet
00096   return _finalise(jet, kept, rejected, subjet_def, discard_area);
00097 }
00098 
00099 
00100 // sets filtered_elements to be all the subjets on which filtering will work
00101 void Filter::_set_filtered_elements(const PseudoJet & jet,
00102                                     vector<PseudoJet> & filtered_elements,
00103                                     JetDefinition & subjet_def,
00104                                     bool & discard_area) const {
00105   // sanity checks
00106   //-------------------------------------------------------------------
00107   // make sure that the jet has constituents
00108   if (! jet.has_constituents())
00109     throw Error("Filter can only be applied on jets having constituents");
00110 
00111   // for a whole variety of checks, we shall need the "recursive"
00112   // pieces of the jet (the jet itself or recursing down to its most
00113   // fundamental pieces).
00114   // So we do compute these once and for all
00115   vector<PseudoJet> all_pieces; //.clear();
00116   if ((!_get_all_pieces(jet, all_pieces)) || (all_pieces.size()==0))
00117     throw Error("Attempt to filter a jet that has no associated ClusterSequence or is not a superposition of jets associated with a ClusterSequence");
00118   
00119   // if the filter uses subtraction, make sure we have a CS that supports area and has
00120   // explicit ghosts 
00121   if (_uses_subtraction()) {
00122     if (!jet.has_area())   
00123       throw Error("Attempt to filter and subtract (non-zero rho or subtractor) without area info for the original jet");
00124 
00125     if (!_check_explicit_ghosts(all_pieces))
00126       throw Error("Attempt to filter and subtract (non-zero rho or subtractor) without explicit ghosts");
00127   }
00128 
00129   // if we're dealing with a dynamic determination of the filtering
00130   // radius, do it now
00131   if ((_Rfilt>=0) || (_Rfiltfunc)){
00132     double Rfilt = (_Rfiltfunc) ? (*_Rfiltfunc)(jet) : _Rfilt;
00133     const JetDefinition::Recombiner * common_recombiner = _get_common_recombiner(all_pieces);
00134     if (common_recombiner) {
00135       if (typeid(*common_recombiner) == typeid(JetDefinition::DefaultRecombiner)) {
00136         RecombinationScheme scheme = 
00137           static_cast<const JetDefinition::DefaultRecombiner *>(common_recombiner)->scheme();
00138         subjet_def = JetDefinition(cambridge_algorithm, Rfilt, scheme);
00139       } else {
00140         subjet_def = JetDefinition(cambridge_algorithm, Rfilt, common_recombiner);
00141       }
00142     } else {
00143       subjet_def = JetDefinition(cambridge_algorithm, Rfilt);
00144     }
00145   } else {
00146     subjet_def = _subjet_def;
00147   }
00148 
00149   // get the jet definition to be use and whether we can apply our
00150   // simplified C/A+C/A filter
00151   //
00152   // we apply C/A clustering iff
00153   //  - the request subjet_def is C/A
00154   //  - the jet is either directly coming from C/A or if it is a
00155   //    superposition of C/A jets
00156   //  - the pieces agree with the recombination scheme of subjet_def
00157   //------------------------------------------------------------------
00158   bool simple_cafilt = _check_ca(all_pieces);
00159 
00160   // extract the subjets
00161   //-------------------------------------------------------------------
00162   discard_area = false;
00163   if (simple_cafilt){
00164     // first make sure that 'filtered_elements' is empty
00165     filtered_elements.clear();
00166     _set_filtered_elements_cafilt(jet, filtered_elements, subjet_def.R());
00167     // in the following case, areas can be erroneous and will be discarded
00168     discard_area = (!_uses_subtraction()) && (jet.has_area()) && (!_check_explicit_ghosts(all_pieces));
00169   } else {
00170     // here we'll simply recluster the jets.
00171     //
00172     // To include an area support we need
00173     //  - the jet to have an area
00174     //  - subtraction requested or explicit ghosts
00175     bool do_areas = (jet.has_area()) && ((_uses_subtraction()) || (_check_explicit_ghosts(all_pieces))); 
00176     _set_filtered_elements_generic(jet, filtered_elements, subjet_def, do_areas);
00177   }
00178 
00179   // order the filtered elements in pt
00180   filtered_elements = sorted_by_pt(filtered_elements);
00181 }
00182 
00183 // set the filtered elements in the simple case of C/A+C/A
00184 //
00185 // WATCH OUT: this could be recursively called, so filtered elements
00186 //            of 'jet' are APPENDED to 'filtered_elements'
00187 void Filter::_set_filtered_elements_cafilt(const PseudoJet & jet, 
00188                                            vector<PseudoJet> & filtered_elements, 
00189                                            double Rfilt) const{
00190   // we know that the jet is either a C/A jet or a superposition of
00191   // such pieces
00192   if (jet.has_associated_cluster_sequence()){
00193     // just extract the exclusive subjets of 'jet'
00194     const ClusterSequence *cs = jet.associated_cluster_sequence(); 
00195     vector<PseudoJet> local_fe;
00196 
00197     double dcut = Rfilt / cs->jet_def().R();
00198     if (dcut>=1.0){
00199       local_fe.push_back(jet);
00200     } else {
00201       dcut *= dcut;
00202       local_fe = jet.exclusive_subjets(dcut);
00203     }
00204 
00205     // subtract the jets if needed
00206     // Note that this one would work on pieces!!
00207     //-----------------------------------------------------------------
00208     if (_uses_subtraction()){
00209       const ClusterSequenceAreaBase * csab = jet.validated_csab();
00210       for (unsigned int i=0;i<local_fe.size();i++) {
00211         if (_subtractor) {
00212           local_fe[i] = (*_subtractor)(local_fe[i]);
00213         } else {
00214           local_fe[i] = csab->subtracted_jet(local_fe[i], _rho);
00215         }
00216       }
00217     }
00218 
00219     copy(local_fe.begin(), local_fe.end(), back_inserter(filtered_elements));
00220     return;
00221   }
00222 
00223   // just recurse into the pieces
00224   const vector<PseudoJet> & pieces = jet.pieces();
00225   for (vector<PseudoJet>::const_iterator it = pieces.begin(); 
00226        it!=pieces.end(); it++)
00227     _set_filtered_elements_cafilt(*it, filtered_elements, Rfilt);
00228 }
00229 
00230 
00231 // set the filtered elements in the generic re-clustering case (wo
00232 // subtraction)
00233 void Filter::_set_filtered_elements_generic(const PseudoJet & jet, 
00234                                             vector<PseudoJet> & filtered_elements,
00235                                             const JetDefinition & subjet_def,
00236                                             bool do_areas) const{
00237   // create a new, internal, ClusterSequence from the jet constituents
00238   // get the subjets directly from there
00239   //
00240   // If the jet has area support then we separate the ghosts from the
00241   // "regular" particles so the subjets will also have area
00242   // support. Note that we do this regardless of whether rho is zero
00243   // or not.
00244   //
00245   // Note that to be able to separate the ghosts, one needs explicit
00246   // ghosts!!
00247   // ---------------------------------------------------------------
00248   if (do_areas){
00249     vector<PseudoJet> all_constituents = jet.constituents();
00250     vector<PseudoJet> regular_constituents, ghosts;  
00251 
00252     for (vector<PseudoJet>::iterator it = all_constituents.begin(); 
00253          it != all_constituents.end(); it++){
00254       if (it->is_pure_ghost())
00255         ghosts.push_back(*it);
00256       else
00257         regular_constituents.push_back(*it);
00258     }
00259 
00260     // figure out the ghost area from the 1st ghost (if none, any value
00261     // would probably do as the area will be 0 and subtraction will have
00262     // no effect!)
00263     double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
00264     ClusterSequenceActiveAreaExplicitGhosts * csa
00265       = new ClusterSequenceActiveAreaExplicitGhosts(regular_constituents, 
00266                                                     subjet_def, 
00267                                                     ghosts, ghost_area);
00268 
00269     // get the subjets: we use the subtracted or unsubtracted ones
00270     // depending on rho or _subtractor being non-zero
00271     if (_uses_subtraction()) {
00272       if (_subtractor) {
00273         filtered_elements = (*_subtractor)(csa->inclusive_jets());
00274       } else {
00275         filtered_elements = csa->subtracted_jets(_rho);
00276       }
00277     } else {
00278       filtered_elements = csa->inclusive_jets();
00279     }
00280 
00281     // allow the cs to be deleted when it's no longer used
00282     csa->delete_self_when_unused();
00283   } else {
00284     ClusterSequence * cs = new ClusterSequence(jet.constituents(), subjet_def);
00285     filtered_elements = cs->inclusive_jets();
00286     // allow the cs to be deleted when it's no longer used
00287     cs->delete_self_when_unused();
00288   }
00289 }
00290 
00291 
00292 // gather the information about what is kept and rejected under the
00293 // form of a PseudoJet with a special ClusterSequenceInfo
00294 PseudoJet Filter::_finalise(const PseudoJet & /*jet*/, 
00295                             vector<PseudoJet> & kept, 
00296                             vector<PseudoJet> & rejected,
00297                             const JetDefinition & subjet_def,
00298                             const bool discard_area) const {
00299   // figure out which recombiner to use
00300   const JetDefinition::Recombiner &rec = *(subjet_def.recombiner());
00301 
00302   // create an appropriate structure and transfer the info to it
00303   PseudoJet filtered_jet = join<StructureType>(kept, rec);
00304   StructureType *fs = (StructureType*) filtered_jet.structure_non_const_ptr();
00305 //  fs->_original_jet = jet;
00306   fs->_rejected = rejected;
00307 
00308   if (discard_area){
00309     // safety check: make sure there is an area to discard!!!
00310     assert(fs->_area_4vector_ptr);
00311     delete fs->_area_4vector_ptr;
00312     fs->_area_4vector_ptr=0;
00313   }
00314   
00315   return filtered_jet;
00316 }
00317 
00318 // various checks
00319 //----------------------------------------------------------------------
00320 
00321 // get the pieces down to the fundamental pieces
00322 // 
00323 // Note that this just checks that there is an associated CS to the
00324 // fundamental pieces, not that it is still valid
00325 bool Filter::_get_all_pieces(const PseudoJet &jet, vector<PseudoJet> &all_pieces) const{
00326   if (jet.has_associated_cluster_sequence()){
00327     all_pieces.push_back(jet);
00328     return true;
00329   }
00330 
00331   if (jet.has_pieces()){
00332     const vector<PseudoJet> pieces = jet.pieces();
00333     for (vector<PseudoJet>::const_iterator it=pieces.begin(); it!=pieces.end(); it++)
00334       if (!_get_all_pieces(*it, all_pieces)) return false;
00335     return true;
00336   }
00337 
00338   return false;
00339 }
00340 
00341 // get the common recombiner to all pieces (NULL if none)
00342 //
00343 // Note that if the jet has an associated cluster sequence that is no
00344 // longer valid, an error will be thrown (needed since it could be the
00345 // 1st check called after the enumeration of the pieces)
00346 const JetDefinition::Recombiner* Filter::_get_common_recombiner(const vector<PseudoJet> &all_pieces) const{
00347   const JetDefinition & jd_ref = all_pieces[0].validated_cs()->jet_def();
00348   for (unsigned int i=1; i<all_pieces.size(); i++)
00349     if (!all_pieces[i].validated_cs()->jet_def().has_same_recombiner(jd_ref)) return NULL;
00350 
00351   return jd_ref.recombiner();
00352 }
00353 
00354 // check if the jet (or all its pieces) have explicit ghosts
00355 // (assuming the jet has area support
00356 //
00357 // Note that if the jet has an associated cluster sequence that is no
00358 // longer valid, an error will be thrown (needed since it could be the
00359 // 1st check called after the enumeration of the pieces)
00360 bool Filter::_check_explicit_ghosts(const vector<PseudoJet> &all_pieces) const{
00361   for (vector<PseudoJet>::const_iterator it=all_pieces.begin(); it!=all_pieces.end(); it++)
00362     if (! it->validated_csab()->has_explicit_ghosts()) return false;
00363   return true;
00364 }
00365 
00366 // check if one can apply the simplification for C/A subjets
00367 //
00368 // This includes:
00369 //  - the subjet definition asks for C/A subjets
00370 //  - all the pieces share the same CS
00371 //  - that CS is C/A with the same recombiner as the subjet def
00372 //  - the filtering radius is not larger than any of the pairwise
00373 //    distance between the pieces
00374 //
00375 // Note that if the jet has an associated cluster sequence that is no
00376 // longer valid, an error will be thrown (needed since it could be the
00377 // 1st check called after the enumeration of the pieces)
00378 bool Filter::_check_ca(const vector<PseudoJet> &all_pieces) const{
00379   if (_subjet_def.jet_algorithm() != cambridge_algorithm) return false;
00380 
00381   // check that the 1st of all the pieces (we're sure there is at
00382   // least one) is coming from a C/A clustering. Then check that all
00383   // the following pieces share the same ClusterSequence
00384   const ClusterSequence * cs_ref = all_pieces[0].validated_cs();
00385   if (cs_ref->jet_def().jet_algorithm() != cambridge_algorithm) return false;
00386   for (unsigned int i=1; i<all_pieces.size(); i++)
00387     if (all_pieces[i].validated_cs() != cs_ref) return false;
00388 
00389   // check that the 1st peice has the same recombiner as the one used
00390   // for the subjet clustering
00391   // Note that since they share the same CS, checking the 2st one is enough
00392   if (!cs_ref->jet_def().has_same_recombiner(_subjet_def)) return false;
00393 
00394   // we also have to make sure that the filtering radius is not larger
00395   // than any of the inter-pieces distance
00396   double Rfilt2 = _subjet_def.R();
00397   Rfilt2 *= Rfilt2;
00398   for (unsigned int i=0; i<all_pieces.size()-1; i++){
00399     for (unsigned int j=i+1; j<all_pieces.size(); j++){
00400       if (all_pieces[i].squared_distance(all_pieces[j]) <  Rfilt2) return false;
00401     }
00402   }
00403 
00404   return true;
00405 }
00406 
00407 //----------------------------------------------------------------------
00408 // FilterInterface implementation 
00409 //----------------------------------------------------------------------
00410 
00411 
00412 FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends