FastJet 3.0.2
|
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