FastJet  3.0.2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
TrackJetPlugin.cc
1 //STARTHEADER
2 // $Id: TrackJetPlugin.cc 2776 2011-11-25 14:59:58Z salam $
3 //
4 // Copyright (c) 2007-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
5 //
6 //----------------------------------------------------------------------
7 // This file is part of FastJet. It contains code that has been
8 // obtained from the Rivet project by Leif Lonnblad, Andy Buckley and
9 // Jon Butterworth. See http://www.hepforge.org/downloads/rivet.
10 // Rivet is free software released under the terms of the GNU Public
11 // License(v2).
12 // Changes from the original file are listed below.
13 //
14 // FastJet is free software; you can redistribute it and/or modify
15 // it under the terms of the GNU General Public License as published by
16 // the Free Software Foundation; either version 2 of the License, or
17 // (at your option) any later version.
18 //
19 // The algorithms that underlie FastJet have required considerable
20 // development and are described in hep-ph/0512210. If you use
21 // FastJet as part of work towards a scientific publication, please
22 // include a citation to the FastJet paper.
23 //
24 // FastJet is distributed in the hope that it will be useful,
25 // but WITHOUT ANY WARRANTY; without even the implied warranty of
26 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27 // GNU General Public License for more details.
28 //
29 // You should have received a copy of the GNU General Public License
30 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
31 //----------------------------------------------------------------------
32 //ENDHEADER
33 
34 // History of changes from the original TrackJet.cc file in Rivet <=1.1.2
35 //
36 // 2011-01-28 Gregory Soyez <soyez@fastjet.fr>
37 //
38 // * Replaced the use of sort by stable_sort (see BUGS in the top
39 // FastJet dir)
40 //
41 //
42 // 2009-01-17 Gregory Soyez <soyez@fastjet.fr>
43 //
44 // * Aligned the var names with the previous conventions
45 //
46 // * Put the plugin in the fastjet::trackjet namespace
47 //
48 //
49 // 2009-01-06 Gregory Soyez <soyez@fastjet.fr>
50 //
51 // * Adapted the original code in a FastJet plugin class.
52 //
53 // * Allowed for arbitrary recombination schemes (one for the
54 // recomstruction of the 'jet' --- i.e. summing the particles
55 // into a jet --- and one for the accumulation of particles in
56 // a 'track' --- i.e. the dynamics of the clustering)
57 
58 
59 // fastjet stuff
60 #include "fastjet/ClusterSequence.hh"
61 #include "fastjet/TrackJetPlugin.hh"
62 
63 // other stuff
64 #include <list>
65 #include <memory>
66 #include <cmath>
67 #include <vector>
68 #include <sstream>
69 
70 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
71 
72 using namespace std;
73 
74 //------------------------------------------------------------------
75 // helper class to sort the particles in pt
76 //------------------------------------------------------------------
77 class TrackJetParticlePtr{
78 public:
79  TrackJetParticlePtr(int i_index, double i_perp2)
80  : index(i_index), perp2(i_perp2){}
81 
82  int index;
83  double perp2;
84 
85  bool operator <(const TrackJetParticlePtr &other) const {
86  return perp2>other.perp2;
87  }
88 };
89 
90 //------------------------------------------------------------------
91 // implementation of the TrackJet plugin
92 //------------------------------------------------------------------
93 
94 bool TrackJetPlugin::_first_time = true;
95 
96 string TrackJetPlugin::description () const {
97  ostringstream desc;
98  desc << "TrackJet algorithm with R = " << R();
99  return desc.str();
100 }
101 
102 void TrackJetPlugin::run_clustering(ClusterSequence & clust_seq) const {
103  // print a banner if we run this for the first time
104  _print_banner(clust_seq.fastjet_banner_stream());
105 
106  // we first need to sort the particles in pt
107  vector<TrackJetParticlePtr> particle_list;
108 
109  const vector<PseudoJet> & jets = clust_seq.jets();
110  int index=0;
111  for (vector<PseudoJet>::const_iterator mom_it = jets.begin(); mom_it != jets.end(); mom_it++){
112  particle_list.push_back(TrackJetParticlePtr(index, mom_it->perp2()));
113  index++;
114  }
115 
116  // sort the particles into decreasing pt
117  stable_sort(particle_list.begin(), particle_list.end());
118 
119 
120  // if we're using a recombination scheme different from the E scheme,
121  // we first need to update the particles' energy so that they
122  // are massless (and rapidity = pseudorapidity)
123  vector<PseudoJet> tuned_particles = clust_seq.jets();
124  vector<PseudoJet> tuned_tracks = clust_seq.jets();
125  for (vector<PseudoJet>::iterator pit = tuned_particles.begin();
126  pit != tuned_particles.end(); pit++)
127  _jet_recombiner.preprocess(*pit);
128  for (vector<PseudoJet>::iterator pit = tuned_tracks.begin();
129  pit != tuned_tracks.end(); pit++)
130  _track_recombiner.preprocess(*pit);
131 
132 
133  // we'll just need the particle indices for what follows
134  list<int> sorted_pt_index;
135  for (vector<TrackJetParticlePtr>::iterator mom_it = particle_list.begin();
136  mom_it != particle_list.end(); mom_it++)
137  sorted_pt_index.push_back(mom_it->index);
138 
139  // now start building the jets
140  while (sorted_pt_index.size()){
141  // note that here 'track' refers to the direction we're using to test if a particle belongs to the jet
142  // 'jet' refers to the momentum of the jet
143  // the difference between the two is in the recombination scheme used to compute the sum of 4-vectors
144  int current_jet_index = sorted_pt_index.front();
145  PseudoJet current_jet = tuned_particles[current_jet_index];
146  PseudoJet current_track = tuned_tracks[current_jet_index];
147 
148  // remove the first particle from the available ones
149  list<int>::iterator index_it = sorted_pt_index.begin();
150  sorted_pt_index.erase(index_it);
151 
152  // start browsing the remaining ones
153  index_it = sorted_pt_index.begin();
154  while (index_it != sorted_pt_index.end()){
155  const PseudoJet & current_particle = tuned_particles[*index_it];
156  const PseudoJet & current_particle_track = tuned_tracks[*index_it];
157 
158  // check if the particle is within a distance R of the jet
159  double distance2 = current_track.plain_distance(current_particle_track);
160  if (distance2 <= _radius2){
161  // add the particle to the jet
162  PseudoJet new_track;
163  PseudoJet new_jet;
164  _jet_recombiner.recombine(current_jet, current_particle, new_jet);
165  _track_recombiner.recombine(current_track, current_particle_track, new_track);
166 
167  int new_jet_index;
168  clust_seq.plugin_record_ij_recombination(current_jet_index, *index_it, distance2, new_jet, new_jet_index);
169 
170  current_jet = new_jet;
171  current_track = new_track;
172  current_jet_index = new_jet_index;
173 
174  // particle has been clustered so remove it from the list
175  sorted_pt_index.erase(index_it);
176 
177  // and don't forget to start again from the beginning
178  // as the jet axis may have changed
179  index_it = sorted_pt_index.begin();
180  } else {
181  index_it++;
182  }
183  }
184 
185  // now we have a final jet, so cluster it with the beam
186  clust_seq.plugin_record_iB_recombination(current_jet_index, _radius2);
187  }
188 
189 }
190 
191 // print a banner for reference to the 3rd-party code
192 void TrackJetPlugin::_print_banner(ostream *ostr) const{
193  if (! _first_time) return;
194  _first_time=false;
195 
196  // make sure the user has not set the banner stream to NULL
197  if (!ostr) return;
198 
199  (*ostr) << "#-------------------------------------------------------------------------" << endl;
200  (*ostr) << "# You are running the TrackJet plugin for FastJet. It is based on " << endl;
201  (*ostr) << "# the implementation by Andy Buckley and Manuel Bahr that is to be " << endl;
202  (*ostr) << "# found in Rivet 1.1.2. See http://www.hepforge.org/downloads/rivet. " << endl;
203  (*ostr) << "#-------------------------------------------------------------------------" << endl;
204 
205  // make sure we really have the output done.
206  ostr->flush();
207 }
208 
209 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh