FastJet  3.0.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
NestedDefsPlugin.cc
1 //STARTHEADER
2 // $Id: NestedDefsPlugin.cc 2577 2011-09-13 15:11:38Z salam $
3 //
4 // Copyright (c) 2007-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 // TODO
30 // ? Maybe one could provide additional recomb. dists as an "extra".;
31 
32 // fastjet stuff
33 #include "fastjet/ClusterSequence.hh"
34 #include "fastjet/NestedDefsPlugin.hh"
35 
36 // other stuff
37 #include <vector>
38 #include <sstream>
39 
40 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
41 
42 using namespace std;
43 
44 string NestedDefsPlugin::description () const {
45  ostringstream desc;
46 
47  desc << "NestedDefs: successive application of " ;
48  unsigned int i=1;
49  for (list<JetDefinition>::const_iterator it=_defs.begin();it!=_defs.end();it++){
50  desc << "Definition " << i++ << " [" << it->description() << "] - ";
51  }
52 
53  return desc.str();
54 }
55 
56 void NestedDefsPlugin::run_clustering(ClusterSequence & clust_seq) const {
57  vector<PseudoJet> momenta;
58 
59  // build the initial list of particles
60  momenta = clust_seq.jets();
61  unsigned int step_n = momenta.size();
62 
63  // initialise the conversion table, which works as follows
64  // conversion_table[step_cs_jet_index] = main_cs_jet_index
65  vector<unsigned int> conversion_table(2*step_n);
66  vector<unsigned int> new_conversion_table;
67  for (unsigned int i=0;i<step_n;i++)
68  conversion_table[i]=i;
69 
70  // Now the steps go as follows:
71  // for each definition in the list,
72  // - do the clustering,
73  // - copy the history into the main one
74  // - update the list of momenta and the index conversion table
75  list<JetDefinition>::const_iterator def_iterator = _defs.begin();
76  unsigned int def_index=0;
77  bool last_def=false;
78 
79  while (def_iterator!=_defs.end()){
80  last_def = (def_index == (_defs.size()-1));
81 
82  // do the clustering
83  ClusterSequence step_cs(momenta, *def_iterator);
84 
85  // clear the momenta as we shall fill them again
86  momenta.clear();
87  new_conversion_table.clear();
88 
89  // retrieve the history
90  const vector<ClusterSequence::history_element> & step_history = step_cs.history();
91 
92  // copy the history
93  // note that we skip the initial steps which are just the
94  // declaration of the particles.
95  vector<ClusterSequence::history_element>::const_iterator
96  hist_iterator = step_history.begin();
97 
98  for (unsigned int i=step_n;i!=0;i--)
99  hist_iterator++;
100 
101  while (hist_iterator != step_history.end()){
102  // check if it is a recombination with the beam or a simple recombination
103  if (hist_iterator->parent2 == ClusterSequence::BeamJet){
104  // save this jet for future clustering
105  // unless we've reached the last def in which case, record the clustering
106  unsigned int step_jet_index = step_cs.history()[hist_iterator->parent1].jetp_index;
107  if (last_def){
108  clust_seq.plugin_record_iB_recombination(conversion_table[step_jet_index],
109  hist_iterator->dij);
110  } else {
111  momenta.push_back(step_cs.jets()[step_jet_index]);
112  new_conversion_table.push_back(conversion_table[step_jet_index]);
113  }
114  } else {
115  // record combination
116  // note that we set the recombination distance to 0 except for the last alg
117  unsigned int step_jet1_index = step_cs.history()[hist_iterator->parent1].jetp_index;
118  unsigned int step_jet2_index = step_cs.history()[hist_iterator->parent2].jetp_index;
119  PseudoJet newjet = step_cs.jets()[hist_iterator->jetp_index];
120  int jet_k;
121  clust_seq.plugin_record_ij_recombination(conversion_table[step_jet1_index],
122  conversion_table[step_jet2_index],
123  last_def ? hist_iterator->dij : 0.0,
124  newjet, jet_k);
125 
126  // save info in the conversion table for tracking purposes
127  conversion_table[hist_iterator->jetp_index]=jet_k;
128  }
129 
130  // go to the next history element
131  hist_iterator++;
132  }
133 
134  // finalise this step:
135  // - update nr of particles
136  // - update conversion table
137  step_n = momenta.size();
138  for (unsigned int i=0;i<step_n;i++)
139  conversion_table[i] = new_conversion_table[i];
140 
141  // go to the next alg
142  def_index++;
143  def_iterator++;
144  }
145 
146 }
147 
148 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh