1 | //STARTHEADER |
---|
2 | // $Id: TrackJetPlugin.cc 859 2012-11-28 01:49:23Z pavel $ |
---|
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 |
---|