FastJet  3.0.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Pruner.cc
1 //STARTHEADER
2 // $Id: Pruner.cc 2689 2011-11-14 14:51:06Z soyez $
3 //
4 // Copyright (c) 2005-2011, 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 and are described in hep-ph/0512210. If you use
16 // FastJet as part of work towards a scientific publication, please
17 // include a citation to the FastJet paper.
18 //
19 // FastJet is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 // GNU General Public License for more details.
23 //
24 // You should have received a copy of the GNU General Public License
25 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
26 //----------------------------------------------------------------------
27 //ENDHEADER
28 
29 #include "fastjet/tools/Pruner.hh"
30 #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
31 #include "fastjet/Selector.hh"
32 #include <cassert>
33 #include <algorithm>
34 #include <sstream>
35 #include <typeinfo>
36 
37 using namespace std;
38 
39 
40 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
41 
42 
43 //----------------------------------------------------------------------
44 // class Pruner
45 //----------------------------------------------------------------------
46 
47 //----------------------------------------------------------------------
48 // alternative (dynamic) ctor
49 // \param jet_def the jet definition for the internal clustering
50 // \param zcut_dyn dynamic pt-fraction cut in the pruning
51 // \param Rcut_dyn dynamic angular distance cut in the pruning
52 Pruner::Pruner(const JetDefinition &jet_def,
55  : _jet_def(jet_def), _zcut(0), _Rcut_factor(0),
56  _zcut_dyn(zcut_dyn), _Rcut_dyn(Rcut_dyn), _get_recombiner_from_jet(false) {
57  assert(_zcut_dyn != 0 && _Rcut_dyn != 0);
58 }
59 
60 //----------------------------------------------------------------------
61 // action on a single jet
63  // pruning can only be applied to jets that have constituents
64  if (!jet.has_constituents()){
65  throw Error("Pruner: trying to apply the Pruner transformer to a jet that has no constituents");
66  }
67 
68  // if the jet has area support and there are explicit ghosts, we can
69  // transfer that support to the internal re-clustering
70  bool do_areas = jet.has_area() && _check_explicit_ghosts(jet);
71 
72  // build the pruning plugin
73  double Rcut = (_Rcut_dyn) ? (*_Rcut_dyn)(jet) : _Rcut_factor * 2.0*jet.m()/jet.perp();
74  double zcut = (_zcut_dyn) ? (*_zcut_dyn)(jet) : _zcut;
75  PruningPlugin * pruning_plugin;
76  // for some constructors, we get the recombiner from the
77  // input jet -- some acrobatics are needed (see plans for FJ3.1
78  // for a hopefully better solution).
79  if (_get_recombiner_from_jet) {
80  const JetDefinition::Recombiner * common_recombiner =
81  _get_common_recombiner(jet);
82  if (common_recombiner) {
83  JetDefinition jet_def = _jet_def;
84  if (typeid(*common_recombiner) == typeid(JetDefinition::DefaultRecombiner)) {
85  RecombinationScheme scheme =
86  static_cast<const JetDefinition::DefaultRecombiner *>(common_recombiner)->scheme();
87  jet_def.set_recombination_scheme(scheme);
88  } else {
89  jet_def.set_recombiner(common_recombiner);
90  }
91  pruning_plugin = new PruningPlugin(jet_def, zcut, Rcut);
92  } else {
93  // if there wasn't a common recombiner, we just use the default
94  // recombiner that was in _jet_def
95  pruning_plugin = new PruningPlugin(_jet_def, zcut, Rcut);
96  }
97  } else {
98  pruning_plugin = new PruningPlugin(_jet_def, zcut, Rcut);
99  }
100 
101  // now recluster the constituents of the jet with that plugin
102  JetDefinition internal_jet_def(pruning_plugin);
103  // flag the plugin for automatic deletion _before_ we make
104  // copies (so that as long as the copies are also present
105  // it doesn't get deleted).
106  internal_jet_def.delete_plugin_when_unused();
107 
108  ClusterSequence * cs;
109  if (do_areas){
110  vector<PseudoJet> particles, ghosts;
111  SelectorIsPureGhost().sift(jet.constituents(), ghosts, particles);
112  // determine the ghost area from the 1st ghost (if none, any value
113  // will do, as the area will be 0 and subtraction will have
114  // no effect!)
115  double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
116  cs = new ClusterSequenceActiveAreaExplicitGhosts(particles, internal_jet_def,
117  ghosts, ghost_area);
118  } else {
119  cs = new ClusterSequence(jet.constituents(), internal_jet_def);
120  }
121 
122  PseudoJet result_local = SelectorNHardest(1)(cs->inclusive_jets())[0];
123  PrunerStructure * s = new PrunerStructure(result_local);
125 
126  // make sure things remain persistent -- i.e. tell the jet definition
127  // and the cluster sequence that it is their responsibility to clean
128  // up memory once the "result" reaches the end of its life in the user's
129  // code. (The CS deletes itself when the result goes out of scope and
130  // that also triggers deletion of the plugin)
132 
133  return result_local;
134 }
135 
136 // check if the jet has explicit_ghosts (knowing that there is an
137 // area support)
138 bool Pruner::_check_explicit_ghosts(const PseudoJet &jet) const{
139  // if the jet comes from a Clustering check explicit ghosts in that
140  // clustering
142  return jet.validated_csab()->has_explicit_ghosts();
143 
144  // if the jet has pieces, recurse in the pieces
145  if (jet.has_pieces()){
146  vector<PseudoJet> pieces = jet.pieces();
147  for (unsigned int i=0;i<pieces.size(); i++)
148  if (!_check_explicit_ghosts(pieces[i])) return false;
149  // never returned false, so we're OK.
150  return true;
151  }
152 
153  // return false for any other (unknown) structure
154  return false;
155 }
156 
157 // see if there is a common recombiner among the pieces; if there
158 // is return a pointer to it; otherwise, return NULL.
159 //
160 // NB: this way of doing things is not ideal, because quite some work
161 // is needed to get a correct handling of the final recombiner
162 // (e.g. default v. non-default). In future add
163 // set_recombiner(jet_def) to JetDefinition, maybe also add
164 // an invalid_scheme to the default recombiner and then
165 // do all the work below directly with a JetDefinition directly
166 // together with JD::has_same_recombiner(...)
167 const JetDefinition::Recombiner * Pruner::_get_common_recombiner(const PseudoJet &jet) const{
168  if (jet.has_associated_cluster_sequence())
169  return jet.validated_cs()->jet_def().recombiner();
170 
171  // if the jet has pieces, recurse in the pieces
172  if (jet.has_pieces()){
173  vector<PseudoJet> pieces = jet.pieces();
174  if (pieces.size() == 0) return 0;
175  const JetDefinition::Recombiner * reco = _get_common_recombiner(pieces[0]);
176  for (unsigned int i=1;i<pieces.size(); i++)
177  if (_get_common_recombiner(pieces[i]) != reco) return 0;
178  // never returned false, so we're OK.
179  return reco;
180  }
181 
182  // return false for any other (unknown) structure
183  return 0;
184 }
185 
186 
187 // transformer description
188 std::string Pruner::description() const{
189  ostringstream oss;
190  oss << "Pruner with jet_definition = (" << _jet_def.description() << ")";
191  if (_zcut_dyn) {
192  oss << ", dynamic zcut (" << _zcut_dyn->description() << ")"
193  << ", dynamic Rcut (" << _Rcut_dyn->description() << ")";
194  } else {
195  oss << ", zcut = " << _zcut
196  << ", Rcut_factor = " << _Rcut_factor;
197  }
198  return oss.str();
199 }
200 
201 
202 
203 //----------------------------------------------------------------------
204 // class PrunerStructure
205 //----------------------------------------------------------------------
206 
207 // return the other jets that may have been found along with the
208 // result of the pruning
209 // The resulting vector is sorted in pt
210 vector<PseudoJet> PrunerStructure::extra_jets() const{
211  return sorted_by_pt((!SelectorNHardest(1))(validated_cs()->inclusive_jets()));;
212 }
213 
214 
215 //----------------------------------------------------------------------
216 // class PruningRecombiner
217 //----------------------------------------------------------------------
218 
219 // decide whether to recombine things or not
220 void PruningRecombiner::recombine(const PseudoJet &pa,
221  const PseudoJet &pb,
222  PseudoJet &pab) const{
223  PseudoJet p;
224  _recombiner->recombine(pa, pb, p);
225 
226  // if the 2 particles are close enough, do the recombination
227  if (pa.squared_distance(pb)<=_Rcut2){
228  pab=p; return;
229  }
230 
231  double pt2a = pa.perp2();
232  double pt2b = pb.perp2();
233 
234  // check which is the softest
235  if (pt2a < pt2b){
236  if (pt2a<_zcut2*p.perp2()){
237  pab = pb; _rejected.push_back(pa.cluster_hist_index());
238  } else {
239  pab = p;
240  }
241  } else {
242  if (pt2b<_zcut2*p.perp2()) {
243  pab = pa; _rejected.push_back(pb.cluster_hist_index());
244  } else {
245  pab = p;
246  }
247  }
248 }
249 
250 // description
251 string PruningRecombiner::description() const{
252  ostringstream oss;
253  oss << "Pruning recombiner with zcut = " << sqrt(_zcut2)
254  << ", Rcut = " << sqrt(_Rcut2)
255  << ", and underlying recombiner = " << _recombiner->description();
256  return oss.str();
257 }
258 
259 
260 
261 
262 //----------------------------------------------------------------------
263 // class PruningPlugin
264 //----------------------------------------------------------------------
265 // the actual clustering work for the plugin
266 void PruningPlugin::run_clustering(ClusterSequence &input_cs) const{
267  // declare a pruning recombiner
268  PruningRecombiner pruning_recombiner(_zcut, _Rcut, _jet_def.recombiner());
269  JetDefinition jet_def = _jet_def;
270  jet_def.set_recombiner(&pruning_recombiner);
271 
272  // cluster the particles using that recombiner
273  ClusterSequence internal_cs(input_cs.jets(), jet_def);
274  const vector<ClusterSequence::history_element> & internal_hist = internal_cs.history();
275 
276  // transfer the list of "childless" elements into a bool vector
277  vector<bool> kept(internal_hist.size(), true);
278  const vector<unsigned int> &pr_rej = pruning_recombiner.rejected();
279  for (unsigned int i=0;i<pr_rej.size(); i++) kept[pr_rej[i]]=false;
280 
281  // browse the history, keeping only the elements that have not been
282  // vetoed.
283  //
284  // In the process we build a map for the history indices
285  vector<unsigned int> internal2input(internal_hist.size());
286  for (unsigned int i=0; i<input_cs.jets().size(); i++)
287  internal2input[i] = i;
288 
289  for (unsigned int i=input_cs.jets().size(); i<internal_hist.size(); i++){
290  const ClusterSequence::history_element &he = internal_hist[i];
291 
292  // deal with recombinations with the beam
293  if (he.parent2 == ClusterSequence::BeamJet){
294  int internal_jetp_index = internal_hist[he.parent1].jetp_index;
295  int internal_hist_index = internal_cs.jets()[internal_jetp_index].cluster_hist_index();
296 
297  int input_jetp_index = input_cs.history()[internal2input[internal_hist_index]].jetp_index;
298 
299  // cout << "Beam recomb for internal " << internal_hist_index
300  // << " (input jet index=" << input_jetp_index << endl;
301 
302  input_cs.plugin_record_iB_recombination(input_jetp_index, he.dij);
303  continue;
304  }
305 
306  // now, deal with two-body recombinations
307  if (!kept[he.parent1]){ // 1 is rejected, we keep only 2
308  internal2input[i]=internal2input[he.parent2];
309  // cout << "rejecting internal " << he.parent1
310  // << ", mapping internal " << i
311  // << " to internal " << he.parent2
312  // << " i.e. " << internal2input[i] << endl;
313  } else if (!kept[he.parent2]){ // 2 is rejected, we keep only 1
314  internal2input[i]=internal2input[he.parent1];
315  // cout << "rejecting internal " << he.parent2
316  // << ", mapping internal " << i
317  // << " to internal " << he.parent1
318  // << " i.e. " << internal2input[i] << endl;
319  } else { // do the recombination
320  int new_index;
321  input_cs.plugin_record_ij_recombination(input_cs.history()[internal2input[he.parent1]].jetp_index,
322  input_cs.history()[internal2input[he.parent2]].jetp_index,
323  he.dij, internal_cs.jets()[he.jetp_index], new_index);
324  internal2input[i]=input_cs.jets()[new_index].cluster_hist_index();
325  // cout << "merging " << internal2input[he.parent1] << " (int: " << he.parent1 << ")"
326  // << " and " << internal2input[he.parent2] << " (int: " << he.parent2 << ")"
327  // << " into " << internal2input[i] << " (int: " << i << ")" << endl;
328  }
329  }
330 }
331 
332 // returns the plugin description
333 string PruningPlugin::description() const{
334  ostringstream oss;
335  oss << "Pruning plugin with jet_definition = (" << _jet_def.description()
336  <<"), zcut = " << _zcut
337  << ", Rcut = " << _Rcut;
338  return oss.str();
339 }
340 
341 
342 FASTJET_END_NAMESPACE