FastJet  3.3.3
TrackJetPlugin.cc
1 //FJSTARTHEADER
2 // $Id: TrackJetPlugin.cc 4420 2019-11-29 09:28:20Z soyez $
3 //
4 // Copyright (c) 2007-2019, 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 // History of changes from the original TrackJet.cc file in Rivet <=1.1.2
32 //
33 // 2011-01-28 Gregory Soyez <soyez@fastjet.fr>
34 //
35 // * Replaced the use of sort by stable_sort (see BUGS in the top
36 // FastJet dir)
37 //
38 //
39 // 2009-01-17 Gregory Soyez <soyez@fastjet.fr>
40 //
41 // * Aligned the var names with the previous conventions
42 //
43 // * Put the plugin in the fastjet::trackjet namespace
44 //
45 //
46 // 2009-01-06 Gregory Soyez <soyez@fastjet.fr>
47 //
48 // * Adapted the original code in a FastJet plugin class.
49 //
50 // * Allowed for arbitrary recombination schemes (one for the
51 // recomstruction of the 'jet' --- i.e. summing the particles
52 // into a jet --- and one for the accumulation of particles in
53 // a 'track' --- i.e. the dynamics of the clustering)
54 
55 
56 // fastjet stuff
57 #include "fastjet/ClusterSequence.hh"
58 #include "fastjet/TrackJetPlugin.hh"
59 
60 // other stuff
61 #include <list>
62 #include <memory>
63 #include <cmath>
64 #include <vector>
65 #include <sstream>
66 
67 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
68 
69 using namespace std;
70 
71 //------------------------------------------------------------------
72 // helper class to sort the particles in pt
73 //------------------------------------------------------------------
74 class TrackJetParticlePtr{
75 public:
76  TrackJetParticlePtr(int i_index, double i_perp2)
77  : index(i_index), perp2(i_perp2){}
78 
79  int index;
80  double perp2;
81 
82  bool operator <(const TrackJetParticlePtr &other) const {
83  return perp2>other.perp2;
84  }
85 };
86 
87 //------------------------------------------------------------------
88 // implementation of the TrackJet plugin
89 //------------------------------------------------------------------
90 
91 bool TrackJetPlugin::_first_time = true;
92 
93 string TrackJetPlugin::description () const {
94  ostringstream desc;
95  desc << "TrackJet algorithm with R = " << R();
96  return desc.str();
97 }
98 
99 void TrackJetPlugin::run_clustering(ClusterSequence & clust_seq) const {
100  // print a banner if we run this for the first time
101  _print_banner(clust_seq.fastjet_banner_stream());
102 
103  // we first need to sort the particles in pt
104  vector<TrackJetParticlePtr> particle_list;
105 
106  const vector<PseudoJet> & jets = clust_seq.jets();
107  int index=0;
108  for (vector<PseudoJet>::const_iterator mom_it = jets.begin(); mom_it != jets.end(); mom_it++){
109  particle_list.push_back(TrackJetParticlePtr(index, mom_it->perp2()));
110  index++;
111  }
112 
113  // sort the particles into decreasing pt
114  stable_sort(particle_list.begin(), particle_list.end());
115 
116 
117  // if we're using a recombination scheme different from the E scheme,
118  // we first need to update the particles' energy so that they
119  // are massless (and rapidity = pseudorapidity)
120  vector<PseudoJet> tuned_particles = clust_seq.jets();
121  vector<PseudoJet> tuned_tracks = clust_seq.jets();
122  for (vector<PseudoJet>::iterator pit = tuned_particles.begin();
123  pit != tuned_particles.end(); pit++)
124  _jet_recombiner.preprocess(*pit);
125  for (vector<PseudoJet>::iterator pit = tuned_tracks.begin();
126  pit != tuned_tracks.end(); pit++)
127  _track_recombiner.preprocess(*pit);
128 
129 
130  // we'll just need the particle indices for what follows
131  list<int> sorted_pt_index;
132  for (vector<TrackJetParticlePtr>::iterator mom_it = particle_list.begin();
133  mom_it != particle_list.end(); mom_it++)
134  sorted_pt_index.push_back(mom_it->index);
135 
136  // now start building the jets
137  while (sorted_pt_index.size()){
138  // note that here 'track' refers to the direction we're using to test if a particle belongs to the jet
139  // 'jet' refers to the momentum of the jet
140  // the difference between the two is in the recombination scheme used to compute the sum of 4-vectors
141  int current_jet_index = sorted_pt_index.front();
142  PseudoJet current_jet = tuned_particles[current_jet_index];
143  PseudoJet current_track = tuned_tracks[current_jet_index];
144 
145  // remove the first particle from the available ones
146  list<int>::iterator index_it = sorted_pt_index.begin();
147  sorted_pt_index.erase(index_it);
148 
149  // start browsing the remaining ones
150  index_it = sorted_pt_index.begin();
151  while (index_it != sorted_pt_index.end()){
152  const PseudoJet & current_particle = tuned_particles[*index_it];
153  const PseudoJet & current_particle_track = tuned_tracks[*index_it];
154 
155  // check if the particle is within a distance R of the jet
156  double distance2 = current_track.plain_distance(current_particle_track);
157  if (distance2 <= _radius2){
158  // add the particle to the jet
159  PseudoJet new_track;
160  PseudoJet new_jet;
161  _jet_recombiner.recombine(current_jet, current_particle, new_jet);
162  _track_recombiner.recombine(current_track, current_particle_track, new_track);
163 
164  int new_jet_index;
165  clust_seq.plugin_record_ij_recombination(current_jet_index, *index_it, distance2, new_jet, new_jet_index);
166 
167  current_jet = new_jet;
168  current_track = new_track;
169  current_jet_index = new_jet_index;
170 
171  // particle has been clustered so remove it from the list
172  sorted_pt_index.erase(index_it);
173 
174  // and don't forget to start again from the beginning
175  // as the jet axis may have changed
176  index_it = sorted_pt_index.begin();
177  } else {
178  index_it++;
179  }
180  }
181 
182  // now we have a final jet, so cluster it with the beam
183  clust_seq.plugin_record_iB_recombination(current_jet_index, _radius2);
184  }
185 
186 }
187 
188 // print a banner for reference to the 3rd-party code
189 void TrackJetPlugin::_print_banner(ostream *ostr) const{
190  if (! _first_time) return;
191  _first_time=false;
192 
193  // make sure the user has not set the banner stream to NULL
194  if (!ostr) return;
195 
196  (*ostr) << "#-------------------------------------------------------------------------" << endl;
197  (*ostr) << "# You are running the TrackJet plugin for FastJet. It is based on " << endl;
198  (*ostr) << "# the implementation by Andy Buckley and Manuel Bahr that is to be " << endl;
199  (*ostr) << "# found in Rivet 1.1.2. See http://www.hepforge.org/downloads/rivet. " << endl;
200  (*ostr) << "#-------------------------------------------------------------------------" << endl;
201 
202  // make sure we really have the output done.
203  ostr->flush();
204 }
205 
206 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
deals with clustering
double plain_distance(const PseudoJet &other) const
returns squared cylinder (rap-phi) distance between this jet and another
Definition: PseudoJet.cc:395
void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, int &newjet_k)
record the fact that there has been a recombination between jets()[jet_i] and jets()[jet_k], with the specified dij, and return the index (newjet_k) allocated to the new jet, whose momentum is assumed to be the 4-vector sum of that of jet_i and jet_j
void plugin_record_iB_recombination(int jet_i, double diB)
record the fact that there has been a recombination between jets()[jet_i] and the beam...
static std::ostream * fastjet_banner_stream()
returns a pointer to the stream to be used to print banners (cout by default).
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...