29 #include "fastjet/tools/Filter.hh"
30 #include <fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh>
39 FASTJET_BEGIN_NAMESPACE
46 string Filter::description()
const {
48 ostr <<
"Filter with subjet_def = ";
50 ostr <<
"Cambridge/Aachen algorithm with dynamic Rfilt"
51 <<
" (recomb. scheme deduced from jet, or E-scheme if not unique)";
52 }
else if (_Rfilt > 0) {
53 ostr <<
"Cambridge/Aachen algorithm with Rfilt = "
55 <<
" (recomb. scheme deduced from jet, or E-scheme if not unique)";
57 ostr << _subjet_def.description();
59 ostr<<
", selection " << _selector.description();
61 ostr <<
", subtractor: " << _subtractor->description();
62 }
else if (_rho != 0) {
63 ostr <<
", subtracting with rho = " << _rho;
79 vector<PseudoJet> subjets;
85 _set_filtered_elements(jet, subjets, subjet_def, discard_area);
88 vector<PseudoJet> kept, rejected;
93 selector_copy.
sift(subjets, kept, rejected);
96 return _finalise(jet, kept, rejected, subjet_def, discard_area);
101 void Filter::_set_filtered_elements(
const PseudoJet & jet,
102 vector<PseudoJet> & filtered_elements,
104 bool & discard_area)
const {
109 throw Error(
"Filter can only be applied on jets having constituents");
115 vector<PseudoJet> all_pieces;
116 if ((!_get_all_pieces(jet, all_pieces)) || (all_pieces.size()==0))
117 throw Error(
"Attempt to filter a jet that has no associated ClusterSequence or is not a superposition of jets associated with a ClusterSequence");
121 if (_uses_subtraction()) {
123 throw Error(
"Attempt to filter and subtract (non-zero rho or subtractor) without area info for the original jet");
125 if (!_check_explicit_ghosts(all_pieces))
126 throw Error(
"Attempt to filter and subtract (non-zero rho or subtractor) without explicit ghosts");
131 if ((_Rfilt>=0) || (_Rfiltfunc)){
132 double Rfilt = (_Rfiltfunc) ? (*_Rfiltfunc)(jet) : _Rfilt;
133 const JetDefinition::Recombiner * common_recombiner = _get_common_recombiner(all_pieces);
134 if (common_recombiner) {
135 if (
typeid(*common_recombiner) ==
typeid(JetDefinition::DefaultRecombiner)) {
136 RecombinationScheme scheme =
137 static_cast<const JetDefinition::DefaultRecombiner *
>(common_recombiner)->scheme();
138 subjet_def = JetDefinition(cambridge_algorithm, Rfilt, scheme);
140 subjet_def = JetDefinition(cambridge_algorithm, Rfilt, common_recombiner);
143 subjet_def = JetDefinition(cambridge_algorithm, Rfilt);
146 subjet_def = _subjet_def;
158 bool simple_cafilt = _check_ca(all_pieces);
162 discard_area =
false;
165 filtered_elements.clear();
166 _set_filtered_elements_cafilt(jet, filtered_elements, subjet_def.R());
168 discard_area = (!_uses_subtraction()) && (jet.
has_area()) && (!_check_explicit_ghosts(all_pieces));
175 bool do_areas = (jet.
has_area()) && ((_uses_subtraction()) || (_check_explicit_ghosts(all_pieces)));
176 _set_filtered_elements_generic(jet, filtered_elements, subjet_def, do_areas);
187 void Filter::_set_filtered_elements_cafilt(
const PseudoJet & jet,
188 vector<PseudoJet> & filtered_elements,
192 if (jet.has_associated_cluster_sequence()){
194 const ClusterSequence *cs = jet.associated_cluster_sequence();
195 vector<PseudoJet> local_fe;
197 double dcut = Rfilt / cs->jet_def().R();
199 local_fe.push_back(jet);
202 local_fe = jet.exclusive_subjets(dcut);
208 if (_uses_subtraction()){
209 const ClusterSequenceAreaBase * csab = jet.validated_csab();
210 for (
unsigned int i=0;i<local_fe.size();i++) {
212 local_fe[i] = (*_subtractor)(local_fe[i]);
214 local_fe[i] = csab->subtracted_jet(local_fe[i], _rho);
219 copy(local_fe.begin(), local_fe.end(), back_inserter(filtered_elements));
224 const vector<PseudoJet> & pieces = jet.pieces();
225 for (vector<PseudoJet>::const_iterator it = pieces.begin();
226 it!=pieces.end(); it++)
227 _set_filtered_elements_cafilt(*it, filtered_elements, Rfilt);
233 void Filter::_set_filtered_elements_generic(
const PseudoJet & jet,
234 vector<PseudoJet> & filtered_elements,
235 const JetDefinition & subjet_def,
236 bool do_areas)
const{
249 vector<PseudoJet> all_constituents = jet.constituents();
250 vector<PseudoJet> regular_constituents, ghosts;
252 for (vector<PseudoJet>::iterator it = all_constituents.begin();
253 it != all_constituents.end(); it++){
254 if (it->is_pure_ghost())
255 ghosts.push_back(*it);
257 regular_constituents.push_back(*it);
263 double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
264 ClusterSequenceActiveAreaExplicitGhosts * csa
265 =
new ClusterSequenceActiveAreaExplicitGhosts(regular_constituents,
271 if (_uses_subtraction()) {
273 filtered_elements = (*_subtractor)(csa->inclusive_jets());
275 filtered_elements = csa->subtracted_jets(_rho);
278 filtered_elements = csa->inclusive_jets();
282 csa->delete_self_when_unused();
284 ClusterSequence * cs =
new ClusterSequence(jet.constituents(), subjet_def);
285 filtered_elements = cs->inclusive_jets();
287 cs->delete_self_when_unused();
294 PseudoJet Filter::_finalise(
const PseudoJet & ,
295 vector<PseudoJet> & kept,
296 vector<PseudoJet> & rejected,
297 const JetDefinition & subjet_def,
298 const bool discard_area)
const {
300 const JetDefinition::Recombiner &rec = *(subjet_def.recombiner());
303 PseudoJet filtered_jet = join<StructureType>(kept, rec);
304 StructureType *fs = (StructureType*) filtered_jet.structure_non_const_ptr();
306 fs->_rejected = rejected;
310 assert(fs->_area_4vector_ptr);
311 delete fs->_area_4vector_ptr;
312 fs->_area_4vector_ptr=0;
325 bool Filter::_get_all_pieces(
const PseudoJet &jet, vector<PseudoJet> &all_pieces)
const{
326 if (jet.has_associated_cluster_sequence()){
327 all_pieces.push_back(jet);
331 if (jet.has_pieces()){
332 const vector<PseudoJet> pieces = jet.pieces();
333 for (vector<PseudoJet>::const_iterator it=pieces.begin(); it!=pieces.end(); it++)
334 if (!_get_all_pieces(*it, all_pieces))
return false;
346 const JetDefinition::Recombiner* Filter::_get_common_recombiner(
const vector<PseudoJet> &all_pieces)
const{
347 const JetDefinition & jd_ref = all_pieces[0].validated_cs()->jet_def();
348 for (
unsigned int i=1; i<all_pieces.size(); i++)
349 if (!all_pieces[i].validated_cs()->jet_def().has_same_recombiner(jd_ref))
return NULL;
351 return jd_ref.recombiner();
360 bool Filter::_check_explicit_ghosts(
const vector<PseudoJet> &all_pieces)
const{
361 for (vector<PseudoJet>::const_iterator it=all_pieces.begin(); it!=all_pieces.end(); it++)
362 if (! it->validated_csab()->has_explicit_ghosts())
return false;
378 bool Filter::_check_ca(
const vector<PseudoJet> &all_pieces)
const{
384 const ClusterSequence * cs_ref = all_pieces[0].validated_cs();
386 for (
unsigned int i=1; i<all_pieces.size(); i++)
387 if (all_pieces[i].validated_cs() != cs_ref)
return false;
392 if (!cs_ref->jet_def().has_same_recombiner(_subjet_def))
return false;
396 double Rfilt2 = _subjet_def.R();
398 for (
unsigned int i=0; i<all_pieces.size()-1; i++){
399 for (
unsigned int j=i+1; j<all_pieces.size(); j++){
400 if (all_pieces[i].squared_distance(all_pieces[j]) < Rfilt2)
return false;
412 FASTJET_END_NAMESPACE
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
void sift(const std::vector< PseudoJet > &jets, std::vector< PseudoJet > &jets_that_pass, std::vector< PseudoJet > &jets_that_fail) const
sift the input jets into two vectors – those that pass the selector and those that do not ...
bool takes_reference() const
returns true if this can be applied jet by jet
virtual bool has_constituents() const
returns true if the PseudoJet has constituents
base class corresponding to errors that can be thrown by FastJet
const Selector & set_reference(const PseudoJet &reference)
set the reference jet for this Selector
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
virtual bool has_area() const
check if it has a defined area
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
class that is intended to hold a full definition of the jet clusterer