FastJet 3.0.2
TrackJetPlugin.cc
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends