FastJet 3.0.2
|
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