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