FastJet 3.0.2
GridJetPlugin.cc
00001 //STARTHEADER
00002 // $Id: GridJetPlugin.cc 2268 2011-06-20 15:12:26Z salam $
00003 //
00004 // Copyright (c) 2011, Matteo Cacciari, Gavin Salam and Gregory Soyez
00005 //
00006 //----------------------------------------------------------------------
00007 // This file is part of FastJet.
00008 //
00009 //  FastJet is free software; you can redistribute it and/or modify
00010 //  it under the terms of the GNU General Public License as published by
00011 //  the Free Software Foundation; either version 2 of the License, or
00012 //  (at your option) any later version.
00013 //
00014 //  The algorithms that underlie FastJet have required considerable
00015 //  development and are described in hep-ph/0512210. If you use
00016 //  FastJet as part of work towards a scientific publication, please
00017 //  include a citation to the FastJet paper.
00018 //
00019 //  FastJet is distributed in the hope that it will be useful,
00020 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00021 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00022 //  GNU General Public License for more details.
00023 //
00024 //  You should have received a copy of the GNU General Public License
00025 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00026 //----------------------------------------------------------------------
00027 //ENDHEADER
00028 
00029 // fastjet stuff
00030 #include "fastjet/ClusterSequence.hh"
00031 #include "fastjet/GridJetPlugin.hh"
00032 
00033 // other stuff
00034 #include <vector>
00035 #include <sstream>
00036 
00037 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00038 
00039 using namespace std;
00040 
00041 //----------------------------------------------------------------------
00042 GridJetPlugin::GridJetPlugin (double ymax,
00043                               double requested_grid_spacing, 
00044                               const JetDefinition & post_jet_def) :
00045   _ymin(-ymax), _ymax(ymax), 
00046   _requested_grid_spacing(requested_grid_spacing) ,
00047   _post_jet_def(post_jet_def)
00048 {
00049   setup_grid();
00050 }
00051 
00052 void GridJetPlugin::setup_grid() {
00053   // since we've exchanged the arguments of the constructor,
00054   // there's a danger of calls with exchanged ymax,spacing arguments -- 
00055   // the following check should catch most such situations.
00056   assert(_ymax>0 && _ymax - _ymin >= _requested_grid_spacing);
00057 
00058   double ny_double = (_ymax-_ymin) / _requested_grid_spacing;
00059   _ny = int(ny_double+0.49999);
00060   _dy = (_ymax-_ymin) / _ny;
00061   
00062   _nphi = int (twopi / _requested_grid_spacing + 0.5);
00063   _dphi = twopi / _nphi;
00064 
00065   // some sanity checking (could throw a fastjet::Error)
00066   assert(_ny >= 1 && _nphi >= 1);
00067 
00068   _ntotal = _nphi * _ny;
00069 }
00070 
00071 
00072 //----------------------------------------------------------------------
00073 string GridJetPlugin::description () const {
00074   ostringstream desc;
00075   desc << "GridJetPlugin plugin with ymax = " << _ymax << ", dy = " << _dy << ", dphi = " << _dphi << " (requested grid spacing was " << _requested_grid_spacing << ")";
00076   if (_post_jet_def.jet_algorithm() != undefined_jet_algorithm) {
00077     desc << ", followed by " << _post_jet_def.description();
00078   }
00079   return desc.str();
00080 }
00081 
00082 
00083 //----------------------------------------------------------------------
00084 double GridJetPlugin::R() const {return sqrt(_dy*_dphi/pi);}
00085 
00086 
00087 //----------------------------------------------------------------------
00088 int GridJetPlugin::igrid(const PseudoJet & p) const {
00089   // directly taking int does not work for values between -1 and 0
00090   // so use floor instead
00091   // double iy_double = (p.rap() - _ymin) / _dy;
00092   // if (iy_double < 0.0) return -1;
00093   // int iy = int(iy_double);
00094   // if (iy >= _ny) return -1;
00095 
00096   // writing it as below gives a huge speed gain (factor two!). Even
00097   // though answers are identical and the routine here is not the
00098   // speed-critical step. It's not at all clear why.
00099   int iy = int(floor( (p.rap() - _ymin) / _dy ));
00100   if (iy < 0 || iy >= _ny) return -1;
00101 
00102   int iphi = int( p.phi()/_dphi );
00103   assert(iphi >= 0 && iphi <= _nphi);
00104   if (iphi == _nphi) iphi = 0; // just in case of rounding errors
00105 
00106   int igrid_res = iy*_nphi + iphi;
00107   assert (igrid_res >= 0 && igrid_res < _ny*_nphi);
00108   return igrid_res;
00109 }
00110 
00111 
00112 //----------------------------------------------------------------------
00113 void GridJetPlugin::run_clustering(ClusterSequence & cs) const {
00114   
00115   // we will create a grid; 
00116   //  * -1 will indicate there is no jet here currently
00117   //  * a number >= 0 will mean that particle indicated by the index
00118   //    is currently the jet on the grid
00119   vector<int> grid(_ntotal, -1);
00120  
00121   int nparticles = cs.jets().size();
00122   double dij_or_diB = 1.0;
00123 
00124   int ngrid_active = 0;
00125 
00126   // combine particles with whatever is in the grid
00127   for (int i = 0; i < nparticles; i++) {
00128     int igrd = igrid(cs.jets()[i]);
00129     //cout << i << " " << cs.jets()[i].rap() << " " << cs.jets()[i].phi() 
00130     //   << " " << igrd << " " << grid.size() << " " << _ntotal << endl;
00131     if (igrd < 0) continue;
00132     assert(igrd <= _ntotal);
00133     if (grid[igrd] == -1) {
00134       grid[igrd] = i; // jet index of initial particle i is i
00135       ngrid_active++;
00136     } else {
00137       int k;
00138       cs.plugin_record_ij_recombination(grid[igrd], i, dij_or_diB, k);
00139       grid[igrd] = k; // grid takes jet index of new particle
00140       //cout << "  res: " << cs.jets()[k].rap() << " " << cs.jets()[k].phi() << endl;
00141     }
00142   }
00143 
00144   if (_post_jet_def.jet_algorithm() == undefined_jet_algorithm) {
00145     // make the final jets via iB recombinations
00146     for (unsigned igrd = 0; igrd < grid.size(); igrd++) {
00147       if (grid[igrd] != -1) cs.plugin_record_iB_recombination(grid[igrd], 
00148                                                               dij_or_diB);
00149     }
00150   } else {
00151     // otherwise post-cluster the grid elements with a normal jet algorithm
00152     vector<PseudoJet> inputs;
00153     vector<int>       cs_indices;
00154     inputs.reserve(ngrid_active);
00155     cs_indices.reserve(2*ngrid_active);
00156     for (unsigned igrd = 0; igrd < grid.size(); igrd++) {
00157       if (grid[igrd] != -1) {
00158         inputs.push_back(cs.jets()[grid[igrd]]);
00159         cs_indices.push_back(grid[igrd]);
00160       }
00161     }
00162     ClusterSequence post_cs(inputs, _post_jet_def);
00163     const vector<ClusterSequence::history_element> & post_history = post_cs.history();
00164     const vector<PseudoJet>                        & post_jets = post_cs.jets();
00165     for (unsigned ihist = ngrid_active; ihist < post_history.size(); ihist++) {
00166       const ClusterSequence::history_element & hist = post_history[ihist];
00167       int post_ij1 = post_history[hist.parent1].jetp_index;
00168       int ij1 = cs_indices[post_ij1];
00169       if (hist.parent2 >= 0) {
00170         int post_ij2 = post_history[hist.parent2].jetp_index;
00171         int ij2 = cs_indices[post_ij2];
00172         int k;
00173         cs.plugin_record_ij_recombination(ij1, ij2, hist.dij, post_jets[hist.jetp_index], k); 
00174         assert(int(cs_indices.size()) == hist.jetp_index);
00175         cs_indices.push_back(k);
00176       } else {
00177         cs.plugin_record_iB_recombination(ij1, hist.dij);
00178       }
00179     }
00180     
00181   }
00182 }
00183 
00184 FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends