FastJet  3.3.4
Filter.cc
1 //FJSTARTHEADER
2 // $Id: Filter.cc 4442 2020-05-05 07:50:11Z soyez $
3 //
4 // Copyright (c) 2005-2020, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
5 //
6 //----------------------------------------------------------------------
7 // This file is part of FastJet.
8 //
9 // FastJet is free software; you can redistribute it and/or modify
10 // it under the terms of the GNU General Public License as published by
11 // the Free Software Foundation; either version 2 of the License, or
12 // (at your option) any later version.
13 //
14 // The algorithms that underlie FastJet have required considerable
15 // development. They are described in the original FastJet paper,
16 // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
17 // FastJet as part of work towards a scientific publication, please
18 // quote the version you use and include a citation to the manual and
19 // optionally also to hep-ph/0512210.
20 //
21 // FastJet is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
28 //----------------------------------------------------------------------
29 //FJENDHEADER
30 
31 #include "fastjet/tools/Filter.hh"
32 #include "fastjet/tools/Recluster.hh"
33 #include "fastjet/tools/Subtractor.hh"
34 #include <fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh>
35 #include <cassert>
36 #include <algorithm>
37 #include <sstream>
38 #include <typeinfo>
39 
40 using namespace std;
41 
42 
43 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
44 
45 //----------------------------------------------------------------------
46 // Filter class implementation
47 //----------------------------------------------------------------------
48 
49 // class description
50 string Filter::description() const {
51  if (!_initialised){
52  return "uninitialised Filter";
53  }
54 
55  ostringstream ostr;
56  ostr << "Filter with subjet_def = ";
57  if (_Rfiltfunc) {
58  ostr << "Cambridge/Aachen algorithm with dynamic Rfilt"
59  << " (recomb. scheme deduced from jet, or E-scheme if not unique)";
60  } else if (_Rfilt > 0) {
61  ostr << "Cambridge/Aachen algorithm with Rfilt = "
62  << _Rfilt
63  << " (recomb. scheme deduced from jet, or E-scheme if not unique)";
64  } else {
65  ostr << _subjet_def.description();
66  }
67  ostr<< ", selection " << _selector.description();
68  if (_subtractor) {
69  ostr << ", subtractor: " << _subtractor->description();
70  } else if (_rho != 0) {
71  ostr << ", subtracting with rho = " << _rho;
72  }
73  return ostr.str();
74 }
75 
76 
77 // core functions
78 //----------------------------------------------------------------------
79 
80 // return a vector of subjets, which are the ones that would be kept
81 // by the filtering
82 PseudoJet Filter::result(const PseudoJet &jet) const {
83  if (!_initialised){
84  //Q: do we throw or do we return an empty PJ?
85  throw Error("uninitialised Filter");
86  }
87 
88  // start by getting the list of subjets (including a list of sanity
89  // checks)
90  // NB: subjets is empty to begin with (see the comment for
91  // _set_filtered_elements_cafilt)
92  vector<PseudoJet> subjets;
93  //JetDefinition subjet_def;
94  bool ca_optimised = _set_filtered_elements(jet, subjets);
95 
96  // apply subtraction if needed:
97  if (_subtractor){
98  subjets = (*_subtractor)(subjets);
99  } else if (_rho!=0){
100  if (subjets.size()>0){
101  //const ClusterSequenceAreaBase *csab = subjets[0].validated_csab();
102  for (unsigned int i=0;i<subjets.size();i++){
103  //subjets[i]=csab->subtracted_jet(subjets[i], _rho);
104  subjets[i]=Subtractor(_rho)(subjets[i]);
105  }
106  }
107  }
108 
109  // now build the vector of kept and rejected subjets
110  vector<PseudoJet> kept, rejected;
111  // Note that in the following line we make a copy of the _selector
112  // to avoid issues with needing a mutable _selector
113  Selector selector_copy = _selector;
114  if (selector_copy.takes_reference()) selector_copy.set_reference(jet);
115  selector_copy.sift(subjets, kept, rejected);
116 
117  // gather the info under the form of a PseudoJet
118  return _finalise(jet, kept, rejected, ca_optimised);
119 }
120 
121 
122 // sets filtered_elements to be all the subjets on which filtering will work
123 //
124 // return true when the subjets have been optained using teh optimised
125 // method for C/A
126 bool Filter::_set_filtered_elements(const PseudoJet & jet,
127  vector<PseudoJet> & filtered_elements) const {
128  // create the recluster instance
129  Recluster recluster;
130  if ((_Rfilt>=0) || (_Rfiltfunc))
131  recluster = Recluster(cambridge_algorithm, (_Rfiltfunc) ? (*_Rfiltfunc)(jet) : _Rfilt, Recluster::keep_all);
132  else
133  recluster = Recluster(_subjet_def, false, Recluster::keep_all);
134 
135  // get the subjets
136  //JetDefinition subjet_def;
137  return recluster.get_new_jets_and_def(jet, filtered_elements);
138 }
139 
140 // gather the information about what is kept and rejected under the
141 // form of a PseudoJet with a special ClusterSequenceInfo
142 PseudoJet Filter::_finalise(const PseudoJet & /*jet*/,
143  vector<PseudoJet> & kept,
144  vector<PseudoJet> & rejected,
145  bool ca_optimisation_used) const {
146  PseudoJet filtered_jet;
147 
148  if (kept.size()+rejected.size()>0){
149  // figure out which recombiner to use
150  const JetDefinition::Recombiner &rec = (kept.size()>0)
151  ? *(kept[0].associated_cs()->jet_def().recombiner())
152  : *(rejected[0].associated_cs()->jet_def().recombiner());
153 
154  // create an appropriate structure and transfer the info to it
155  filtered_jet = join<StructureType>(kept, rec);
156  } else {
157  filtered_jet = join<StructureType>(kept);
158  }
159  StructureType *fs = (StructureType*) filtered_jet.structure_non_const_ptr();
160  fs->_rejected = rejected;
161 
162  // if we've used C/A optimisation, we need to get rid of the area
163  // information if it comes from a non-explicit-ghost clustering.
164  // (because in that case it can be erroneous due the lack of
165  // information about empty areas)
166  if ((ca_optimisation_used) && (kept.size()+rejected.size()>0)){
167  bool has_non_explicit_ghost_area = (kept.size()>0)
168  ? (kept[0].has_area() && (!(kept[0].validated_csab()->has_explicit_ghosts())))
169  : (rejected[0].has_area() && (!(rejected[0].validated_csab()->has_explicit_ghosts())));
170  if (has_non_explicit_ghost_area)
171  fs->discard_area();
172  }
173 
174  return filtered_jet;
175 }
176 
177 
178 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
fastjet::Selector::sift
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
Definition: Selector.cc:153
fastjet::Selector::takes_reference
bool takes_reference() const
returns true if this can be applied jet by jet
Definition: Selector.hh:287
fastjet::Recluster
Definition: Recluster.hh:73
fastjet::Selector::set_reference
const Selector & set_reference(const PseudoJet &reference)
set the reference jet for this Selector
Definition: Selector.hh:292
fastjet::PseudoJet
Definition: PseudoJet.hh:67
fastjet::Selector
Definition: Selector.hh:149
fastjet::Subtractor
Definition: Subtractor.hh:62
fastjet::Recluster::get_new_jets_and_def
bool get_new_jets_and_def(const PseudoJet &input_jet, std::vector< PseudoJet > &output_jets) const
A lower-level method that does the actual work of reclustering the input jet.
Definition: Recluster.cc:108
fastjet::Error
Definition: Error.hh:47