FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: TrackJetPlugin.cc 2776 2011-11-25 14:59:58Z salam $ 00003 // 00004 // Copyright (c) 2007-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 00005 // 00006 //---------------------------------------------------------------------- 00007 // This file is part of FastJet. It contains code that has been 00008 // obtained from the Rivet project by Leif Lonnblad, Andy Buckley and 00009 // Jon Butterworth. See http://www.hepforge.org/downloads/rivet. 00010 // Rivet is free software released under the terms of the GNU Public 00011 // License(v2). 00012 // Changes from the original file are listed below. 00013 // 00014 // FastJet is free software; you can redistribute it and/or modify 00015 // it under the terms of the GNU General Public License as published by 00016 // the Free Software Foundation; either version 2 of the License, or 00017 // (at your option) any later version. 00018 // 00019 // The algorithms that underlie FastJet have required considerable 00020 // development and are described in hep-ph/0512210. If you use 00021 // FastJet as part of work towards a scientific publication, please 00022 // include a citation to the FastJet paper. 00023 // 00024 // FastJet is distributed in the hope that it will be useful, 00025 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00026 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00027 // GNU General Public License for more details. 00028 // 00029 // You should have received a copy of the GNU General Public License 00030 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00031 //---------------------------------------------------------------------- 00032 //ENDHEADER 00033 00034 // History of changes from the original TrackJet.cc file in Rivet <=1.1.2 00035 // 00036 // 2011-01-28 Gregory Soyez <soyez@fastjet.fr> 00037 // 00038 // * Replaced the use of sort by stable_sort (see BUGS in the top 00039 // FastJet dir) 00040 // 00041 // 00042 // 2009-01-17 Gregory Soyez <soyez@fastjet.fr> 00043 // 00044 // * Aligned the var names with the previous conventions 00045 // 00046 // * Put the plugin in the fastjet::trackjet namespace 00047 // 00048 // 00049 // 2009-01-06 Gregory Soyez <soyez@fastjet.fr> 00050 // 00051 // * Adapted the original code in a FastJet plugin class. 00052 // 00053 // * Allowed for arbitrary recombination schemes (one for the 00054 // recomstruction of the 'jet' --- i.e. summing the particles 00055 // into a jet --- and one for the accumulation of particles in 00056 // a 'track' --- i.e. the dynamics of the clustering) 00057 00058 00059 // fastjet stuff 00060 #include "fastjet/ClusterSequence.hh" 00061 #include "fastjet/TrackJetPlugin.hh" 00062 00063 // other stuff 00064 #include <list> 00065 #include <memory> 00066 #include <cmath> 00067 #include <vector> 00068 #include <sstream> 00069 00070 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00071 00072 using namespace std; 00073 00074 //------------------------------------------------------------------ 00075 // helper class to sort the particles in pt 00076 //------------------------------------------------------------------ 00077 class TrackJetParticlePtr{ 00078 public: 00079 TrackJetParticlePtr(int i_index, double i_perp2) 00080 : index(i_index), perp2(i_perp2){} 00081 00082 int index; 00083 double perp2; 00084 00085 bool operator <(const TrackJetParticlePtr &other) const { 00086 return perp2>other.perp2; 00087 } 00088 }; 00089 00090 //------------------------------------------------------------------ 00091 // implementation of the TrackJet plugin 00092 //------------------------------------------------------------------ 00093 00094 bool TrackJetPlugin::_first_time = true; 00095 00096 string TrackJetPlugin::description () const { 00097 ostringstream desc; 00098 desc << "TrackJet algorithm with R = " << R(); 00099 return desc.str(); 00100 } 00101 00102 void TrackJetPlugin::run_clustering(ClusterSequence & clust_seq) const { 00103 // print a banner if we run this for the first time 00104 _print_banner(clust_seq.fastjet_banner_stream()); 00105 00106 // we first need to sort the particles in pt 00107 vector<TrackJetParticlePtr> particle_list; 00108 00109 const vector<PseudoJet> & jets = clust_seq.jets(); 00110 int index=0; 00111 for (vector<PseudoJet>::const_iterator mom_it = jets.begin(); mom_it != jets.end(); mom_it++){ 00112 particle_list.push_back(TrackJetParticlePtr(index, mom_it->perp2())); 00113 index++; 00114 } 00115 00116 // sort the particles into decreasing pt 00117 stable_sort(particle_list.begin(), particle_list.end()); 00118 00119 00120 // if we're using a recombination scheme different from the E scheme, 00121 // we first need to update the particles' energy so that they 00122 // are massless (and rapidity = pseudorapidity) 00123 vector<PseudoJet> tuned_particles = clust_seq.jets(); 00124 vector<PseudoJet> tuned_tracks = clust_seq.jets(); 00125 for (vector<PseudoJet>::iterator pit = tuned_particles.begin(); 00126 pit != tuned_particles.end(); pit++) 00127 _jet_recombiner.preprocess(*pit); 00128 for (vector<PseudoJet>::iterator pit = tuned_tracks.begin(); 00129 pit != tuned_tracks.end(); pit++) 00130 _track_recombiner.preprocess(*pit); 00131 00132 00133 // we'll just need the particle indices for what follows 00134 list<int> sorted_pt_index; 00135 for (vector<TrackJetParticlePtr>::iterator mom_it = particle_list.begin(); 00136 mom_it != particle_list.end(); mom_it++) 00137 sorted_pt_index.push_back(mom_it->index); 00138 00139 // now start building the jets 00140 while (sorted_pt_index.size()){ 00141 // note that here 'track' refers to the direction we're using to test if a particle belongs to the jet 00142 // 'jet' refers to the momentum of the jet 00143 // the difference between the two is in the recombination scheme used to compute the sum of 4-vectors 00144 int current_jet_index = sorted_pt_index.front(); 00145 PseudoJet current_jet = tuned_particles[current_jet_index]; 00146 PseudoJet current_track = tuned_tracks[current_jet_index]; 00147 00148 // remove the first particle from the available ones 00149 list<int>::iterator index_it = sorted_pt_index.begin(); 00150 sorted_pt_index.erase(index_it); 00151 00152 // start browsing the remaining ones 00153 index_it = sorted_pt_index.begin(); 00154 while (index_it != sorted_pt_index.end()){ 00155 const PseudoJet & current_particle = tuned_particles[*index_it]; 00156 const PseudoJet & current_particle_track = tuned_tracks[*index_it]; 00157 00158 // check if the particle is within a distance R of the jet 00159 double distance2 = current_track.plain_distance(current_particle_track); 00160 if (distance2 <= _radius2){ 00161 // add the particle to the jet 00162 PseudoJet new_track; 00163 PseudoJet new_jet; 00164 _jet_recombiner.recombine(current_jet, current_particle, new_jet); 00165 _track_recombiner.recombine(current_track, current_particle_track, new_track); 00166 00167 int new_jet_index; 00168 clust_seq.plugin_record_ij_recombination(current_jet_index, *index_it, distance2, new_jet, new_jet_index); 00169 00170 current_jet = new_jet; 00171 current_track = new_track; 00172 current_jet_index = new_jet_index; 00173 00174 // particle has been clustered so remove it from the list 00175 sorted_pt_index.erase(index_it); 00176 00177 // and don't forget to start again from the beginning 00178 // as the jet axis may have changed 00179 index_it = sorted_pt_index.begin(); 00180 } else { 00181 index_it++; 00182 } 00183 } 00184 00185 // now we have a final jet, so cluster it with the beam 00186 clust_seq.plugin_record_iB_recombination(current_jet_index, _radius2); 00187 } 00188 00189 } 00190 00191 // print a banner for reference to the 3rd-party code 00192 void TrackJetPlugin::_print_banner(ostream *ostr) const{ 00193 if (! _first_time) return; 00194 _first_time=false; 00195 00196 // make sure the user has not set the banner stream to NULL 00197 if (!ostr) return; 00198 00199 (*ostr) << "#-------------------------------------------------------------------------" << endl; 00200 (*ostr) << "# You are running the TrackJet plugin for FastJet. It is based on " << endl; 00201 (*ostr) << "# the implementation by Andy Buckley and Manuel Bahr that is to be " << endl; 00202 (*ostr) << "# found in Rivet 1.1.2. See http://www.hepforge.org/downloads/rivet. " << endl; 00203 (*ostr) << "#-------------------------------------------------------------------------" << endl; 00204 00205 // make sure we really have the output done. 00206 ostr->flush(); 00207 } 00208 00209 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh