| 1 | // $Id$
|
|---|
| 2 | //
|
|---|
| 3 | // Copyright (c) 2014-, Matteo Cacciari, Gavin. P. Salam and Gregory Soyez
|
|---|
| 4 | //
|
|---|
| 5 | //----------------------------------------------------------------------
|
|---|
| 6 | // This file is part of FastJet contrib.
|
|---|
| 7 | //
|
|---|
| 8 | // It is free software; you can redistribute it and/or modify it under
|
|---|
| 9 | // the terms of the GNU General Public License as published by the
|
|---|
| 10 | // Free Software Foundation; either version 2 of the License, or (at
|
|---|
| 11 | // your option) any later version.
|
|---|
| 12 | //
|
|---|
| 13 | // It is distributed in the hope that it will be useful, but WITHOUT
|
|---|
| 14 | // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
|
|---|
| 15 | // or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
|
|---|
| 16 | // License for more details.
|
|---|
| 17 | //
|
|---|
| 18 | // You should have received a copy of the GNU General Public License
|
|---|
| 19 | // along with this code. If not, see <http://www.gnu.org/licenses/>.
|
|---|
| 20 | //----------------------------------------------------------------------
|
|---|
| 21 |
|
|---|
| 22 | #include "SoftKiller.hh"
|
|---|
| 23 | #include <sstream>
|
|---|
| 24 |
|
|---|
| 25 | using namespace std;
|
|---|
| 26 |
|
|---|
| 27 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
|
|---|
| 28 |
|
|---|
| 29 | namespace contrib{
|
|---|
| 30 |
|
|---|
| 31 | /// ctor with simple initialisation
|
|---|
| 32 | /// \param rapmax the maximal absolute rapidity extent of the grid
|
|---|
| 33 | /// \param cell_size the grid spacing (equivalently, cell size)
|
|---|
| 34 | /// \param sifter when provided, the soft killer is applied
|
|---|
| 35 | /// only to particles that pass the sifter (the
|
|---|
| 36 | /// others are kept untouched)
|
|---|
| 37 | SoftKiller::SoftKiller(double rapmax, double cell_size,
|
|---|
| 38 | Selector sifter) :
|
|---|
| 39 | #ifdef FJCONTRIB_SOFTKILLER_USEFJGRID
|
|---|
| 40 | RectangularGrid(rapmax, cell_size), _sifter(sifter) {}
|
|---|
| 41 | #else // not FJCONTRIB_SOFTKILLER_USEFJGRID
|
|---|
| 42 | _ymax(rapmax), _ymin(-rapmax),
|
|---|
| 43 | _requested_drap(cell_size), _requested_dphi(cell_size),
|
|---|
| 44 | _sifter(sifter) {
|
|---|
| 45 | _setup_grid();
|
|---|
| 46 | }
|
|---|
| 47 | #endif
|
|---|
| 48 |
|
|---|
| 49 | /// ctor with more control over initialisation
|
|---|
| 50 | /// \param rapmin the minimum rapidity extent of the grid
|
|---|
| 51 | /// \param rapmax the maximum rapidity extent of the grid
|
|---|
| 52 | /// \param drap the grid spacing in rapidity
|
|---|
| 53 | /// \param dphi the grid spacing in azimuth
|
|---|
| 54 | /// \param sifter when provided, the soft killer is applied
|
|---|
| 55 | /// onyl to particles that pass the sifter (the
|
|---|
| 56 | /// others are kept untouched)
|
|---|
| 57 | SoftKiller::SoftKiller(double rapmin, double rapmax, double drap, double dphi,
|
|---|
| 58 | Selector sifter) :
|
|---|
| 59 | #ifdef FJCONTRIB_SOFTKILLER_USEFJGRID
|
|---|
| 60 | RectangularGrid(rapmin, rapmax, drap, dphi), _sifter(sifter) {}
|
|---|
| 61 | #else
|
|---|
| 62 | _ymax(rapmax), _ymin(rapmin),
|
|---|
| 63 | _requested_drap(drap), _requested_dphi(dphi),
|
|---|
| 64 | _sifter(sifter) {
|
|---|
| 65 | _setup_grid();
|
|---|
| 66 | }
|
|---|
| 67 | #endif
|
|---|
| 68 |
|
|---|
| 69 | #ifdef FJCONTRIB_SOFTKILLER_USEFJGRID
|
|---|
| 70 | SoftKiller::SoftKiller(const RectangularGrid & grid, Selector sifter) :
|
|---|
| 71 | RectangularGrid(grid), _sifter(sifter) {}
|
|---|
| 72 | #endif
|
|---|
| 73 |
|
|---|
| 74 | /// dummy ctor (will give an unusable SoftKiller)
|
|---|
| 75 | SoftKiller::SoftKiller()
|
|---|
| 76 | #ifdef FJCONTRIB_SOFTKILLER_USEFJGRID
|
|---|
| 77 | {}
|
|---|
| 78 | #else
|
|---|
| 79 | : _ymax(-1.0), _ymin(1.0), _requested_drap(-1.0), _requested_dphi(-1.0) {
|
|---|
| 80 | _ntotal = 0;
|
|---|
| 81 | }
|
|---|
| 82 | #endif // FJCONTRIB_SOFTKILLER_USEFJGRID
|
|---|
| 83 |
|
|---|
| 84 |
|
|---|
| 85 | //------------------------------------------------------------------------
|
|---|
| 86 | // description of the soft killer
|
|---|
| 87 | std::string SoftKiller::description() const{
|
|---|
| 88 | ostringstream oss;
|
|---|
| 89 |
|
|---|
| 90 |
|
|---|
| 91 | #ifdef FJCONTRIB_SOFTKILLER_USEFJGRID
|
|---|
| 92 | oss << "SoftKiller with " << RectangularGrid::description();
|
|---|
| 93 | #else
|
|---|
| 94 | if (_requested_drap < 0 || _requested_dphi < 0)
|
|---|
| 95 | return "Uninitialised SoftKiller";
|
|---|
| 96 |
|
|---|
| 97 | oss << "SoftKiller with rapidity extent " << _ymin << " < rap < " << _ymax
|
|---|
| 98 | << ", cell size drap x dphi = " << _dy << " x " << _dphi;
|
|---|
| 99 | #endif
|
|---|
| 100 | if (_sifter.worker()) {
|
|---|
| 101 | oss << " and applied to particles passing the selection ("
|
|---|
| 102 | << _sifter.description() << ")";
|
|---|
| 103 | }
|
|---|
| 104 | return oss.str();
|
|---|
| 105 | }
|
|---|
| 106 |
|
|---|
| 107 | //------------------------------------------------------------------------
|
|---|
| 108 | // similarly to Transformers in FastJet, introduce a 'result'
|
|---|
| 109 | // method equivalent to the () operator, i.e. returns the event
|
|---|
| 110 | // after the soft killer has been applied
|
|---|
| 111 | //vector<PseudoJet> SoftKiller::result(const vector<PseudoJet> & event) const {
|
|---|
| 112 | void SoftKiller::apply(const vector<PseudoJet> & event,
|
|---|
| 113 | vector<PseudoJet> & reduced_event,
|
|---|
| 114 | double & pt_threshold) const {
|
|---|
| 115 | // a safety check: we impose at least 2 cells (otherwise, this is
|
|---|
| 116 | // equivalent to asking an empty event and there are more
|
|---|
| 117 | // efficient ways to do that)
|
|---|
| 118 | if (n_tiles()<2){
|
|---|
| 119 | throw Error("SoftKiller not properly initialised.");
|
|---|
| 120 | }
|
|---|
| 121 |
|
|---|
| 122 | // we're not set up to handle the case where the event and reduced
|
|---|
| 123 | // event are the same vector; so crash in that case.
|
|---|
| 124 | assert(&event != &reduced_event);
|
|---|
| 125 | // currently we can only handle the case where all tiles have equal
|
|---|
| 126 | // area; that is the case for Rectangular tilings, but in the future
|
|---|
| 127 | // one might imagine having non-rectangular tilings.
|
|---|
| 128 | assert(all_tiles_equal_area());
|
|---|
| 129 |
|
|---|
| 130 | //profiling: CPUTimer t;
|
|---|
| 131 | //profiling: t.start();
|
|---|
| 132 |
|
|---|
| 133 | // init an array to hold the max pt in each grid cell
|
|---|
| 134 | //
|
|---|
| 135 | // This is better (stack) but only C99 (fails with pedantic):
|
|---|
| 136 | // double max_pt2[n_tiles()];
|
|---|
| 137 | // See
|
|---|
| 138 | // http://binglongx.wordpress.com/2011/05/08/create-variable-length-array-on-the-stack-in-c/
|
|---|
| 139 | // for a possible workaround
|
|---|
| 140 | vector<double> max_pt2(n_tiles(), 0.0);
|
|---|
| 141 | //double max_pt2[n_tiles()];
|
|---|
| 142 | //memset(max_pt2, 0, n_tiles()*sizeof(double));
|
|---|
| 143 | //double *max_pt2 = new double[n_tiles()]; // if this is reinstated, reinstate also the delete, below
|
|---|
| 144 | //for (int i = 0; i < n_tiles(); i++) {max_pt2[i] = 0;}
|
|---|
| 145 |
|
|---|
| 146 | //profiling: t.stop();
|
|---|
| 147 | //profiling: cout << " copy ptrs: " << t.total() << endl;
|
|---|
| 148 | //profiling: t.start();
|
|---|
| 149 |
|
|---|
| 150 | // leave away particles that are sifted
|
|---|
| 151 | // vector<PseudoJet> reduced_event, saved_particles;
|
|---|
| 152 | // _sifter.sift(event, reduced_event, saved_particles);
|
|---|
| 153 |
|
|---|
| 154 | vector<const PseudoJet *> event_ptrs(event.size());
|
|---|
| 155 | for (unsigned i = 0; i < event.size(); i++) {
|
|---|
| 156 | event_ptrs[i] = & event[i];
|
|---|
| 157 | }
|
|---|
| 158 | // only run the sifter if it serves a purpose
|
|---|
| 159 | if (_sifter.worker()) _sifter.nullify_non_selected(event_ptrs);
|
|---|
| 160 |
|
|---|
| 161 | //profiling: t.stop();
|
|---|
| 162 | //profiling: cout << " sifter : " << t.total() << endl;
|
|---|
| 163 | //profiling: t.start();
|
|---|
| 164 |
|
|---|
| 165 | //vector<PseudoJet> reduced_event = event;
|
|---|
| 166 |
|
|---|
| 167 | // browse the particles and figure which is the min pt in each cell
|
|---|
| 168 | for(unsigned int i=0;i<event.size(); i++){
|
|---|
| 169 | if (event_ptrs[i] == 0) continue;
|
|---|
| 170 | int idx = tile_index(event[i]);
|
|---|
| 171 | if (idx<0) continue;
|
|---|
| 172 | max_pt2[idx] = max(max_pt2[idx], event[i].pt2());
|
|---|
| 173 | }
|
|---|
| 174 |
|
|---|
| 175 | //profiling: t.stop();
|
|---|
| 176 | //profiling: cout << " browsing : " << t.total() << endl;
|
|---|
| 177 | //profiling: t.start();
|
|---|
| 178 |
|
|---|
| 179 |
|
|---|
| 180 | // if there are some "bad" tiles, then we need to exclude them from
|
|---|
| 181 | // the calculation of the median. We'll do this by condensing the
|
|---|
| 182 | // max_pt2 vector down to just the values for the tiles that are
|
|---|
| 183 | // good.
|
|---|
| 184 | //
|
|---|
| 185 | // tested answers look right in "issue" 2014-08-08-testing-rect-grid
|
|---|
| 186 | if (n_good_tiles() != n_tiles()) {
|
|---|
| 187 | int newn = 0;
|
|---|
| 188 | for (unsigned i = 0; i < max_pt2.size(); i++) {
|
|---|
| 189 | if (tile_is_good(i)) {
|
|---|
| 190 | // clang gets confused with the SharedPtr swap if we don't
|
|---|
| 191 | // have std:: here
|
|---|
| 192 | std::swap(max_pt2[i],max_pt2[newn]);
|
|---|
| 193 | newn++;
|
|---|
| 194 | }
|
|---|
| 195 | }
|
|---|
| 196 | max_pt2.resize(newn);
|
|---|
| 197 | }
|
|---|
| 198 |
|
|---|
| 199 | // sort the list of max values (sort works with arrays)
|
|---|
| 200 | //sort(max_pt2, max_pt2+n_tiles()); // use this one if we have the allocated C-style array above
|
|---|
| 201 | sort(max_pt2.begin(), max_pt2.end());
|
|---|
| 202 |
|
|---|
| 203 | // Note that if we ask that half of the event is empty that means
|
|---|
| 204 | // that we need at least half of the cells empty... so we need to
|
|---|
| 205 | // round up
|
|---|
| 206 | //
|
|---|
| 207 | // For a potentially faster median search, see
|
|---|
| 208 | //
|
|---|
| 209 | // http://ndevilla.free.fr/median/median/index.html?utm_source=feedburner&utm_medium=twitter&utm_campaign=Feed%3A+hnycombinator+%28HN+-+hnycombinator%29
|
|---|
| 210 | int int_median_pos = max_pt2.size()/2;
|
|---|
| 211 |
|
|---|
| 212 | //--------------------------------------------------------------
|
|---|
| 213 | // alternative to what is below:
|
|---|
| 214 | //
|
|---|
| 215 | // apply the pt cut manually (always on pt^2): first get a vector of
|
|---|
| 216 | // indices that we'll want to keep and then when we know the size of
|
|---|
| 217 | // the resulting vector we actually do the transfer of the
|
|---|
| 218 | // PseudoJets.
|
|---|
| 219 | //
|
|---|
| 220 | // This is almost a factor of two faster than doing a
|
|---|
| 221 | // push_back on the PseudoJets themselves (which is logical, since a
|
|---|
| 222 | // push_back probably averages out as doing the copy twice).
|
|---|
| 223 | // (2014-08-11, also tried it with pointers, which seemed
|
|---|
| 224 | // marginally slower).
|
|---|
| 225 | double pt2cut = (1+1e-12)*max_pt2[int_median_pos];
|
|---|
| 226 |
|
|---|
| 227 | vector<int> indices;
|
|---|
| 228 | for(unsigned int i=0;i<event.size(); i++){
|
|---|
| 229 | if ((event_ptrs[i] == 0) || (event[i].pt2() >= pt2cut))
|
|---|
| 230 | indices.push_back(i);
|
|---|
| 231 | }
|
|---|
| 232 | reduced_event.reserve(indices.size());
|
|---|
| 233 | for(unsigned int i=0;i<indices.size(); i++){
|
|---|
| 234 | reduced_event.push_back(event[indices[i]]);
|
|---|
| 235 | }
|
|---|
| 236 |
|
|---|
| 237 | // //vector<PseudoJet> reduced_event;
|
|---|
| 238 | // reduced_event.clear();
|
|---|
| 239 | // for(unsigned int i=0;i<event.size(); i++){
|
|---|
| 240 | // if ((event_ptrs[i] == 0) ||
|
|---|
| 241 | // (event[i].pt2() >= pt2cut))
|
|---|
| 242 | // reduced_event.push_back(event[i]);
|
|---|
| 243 | // }
|
|---|
| 244 |
|
|---|
| 245 | // free memory: reinstate this is max_pt2 becomes an allocated variable again
|
|---|
| 246 | //delete max_pt2;
|
|---|
| 247 |
|
|---|
| 248 | //return reduced_event;
|
|---|
| 249 | pt_threshold = sqrt(pt2cut);
|
|---|
| 250 |
|
|---|
| 251 | // end of alternative
|
|---|
| 252 | //--------------------------------------------------------------
|
|---|
| 253 |
|
|---|
| 254 | // double median_maxpt = sqrt(max_pt2[int_median_pos]);
|
|---|
| 255 | //
|
|---|
| 256 | // //profiling: t.stop();
|
|---|
| 257 | // //profiling: cout << " median : " << t.total() << endl;
|
|---|
| 258 | // //profiling: t.start();
|
|---|
| 259 | //
|
|---|
| 260 | // // apply a cut on pt using a selector
|
|---|
| 261 | // //
|
|---|
| 262 | // // Watch out that the Selector checks pt >= ptcut, and, since it
|
|---|
| 263 | // // uses pt2 to perform the comparison, may lead to rounding
|
|---|
| 264 | // // errors. By multiplying the ptcut by a small amount, we make
|
|---|
| 265 | // // sure the particle with pt=ptcut is killed.
|
|---|
| 266 | // //
|
|---|
| 267 | // // We're actually going to set to null the ones that will be kept
|
|---|
| 268 | // SelectorPtMax((1+1e-12)*median_maxpt).nullify_non_selected(event_ptrs);
|
|---|
| 269 | //
|
|---|
| 270 | // // basic information
|
|---|
| 271 | // _ptcut = median_maxpt; // a good first start
|
|---|
| 272 | //
|
|---|
| 273 | // //profiling: t.stop();
|
|---|
| 274 | // //profiling: cout << " cutting : " << t.total() << endl;
|
|---|
| 275 | // //profiling: t.start();
|
|---|
| 276 | //
|
|---|
| 277 | // // then put back in the saved particles
|
|---|
| 278 | // //
|
|---|
| 279 | // // Note that here we may want to use the killer independently on
|
|---|
| 280 | // // the 2 sets of particles but then the best woulb be to move most
|
|---|
| 281 | // // of the above in a separate method and call it twice (left for
|
|---|
| 282 | // // future experimentation)
|
|---|
| 283 | // //copy(saved_particles.begin(), saved_particles.end(), back_inserter(reduced_event));
|
|---|
| 284 | //
|
|---|
| 285 | // vector<PseudoJet> reduced_event;
|
|---|
| 286 | // // no sizeable speed gain: reduced_event.reserve(event.size());
|
|---|
| 287 | // for (unsigned int i=0;i<event.size();i++)
|
|---|
| 288 | // if (event_ptrs[i] == 0) reduced_event.push_back(event[i]);
|
|---|
| 289 | //
|
|---|
| 290 | // // free memory
|
|---|
| 291 | // delete max_pt2;
|
|---|
| 292 | //
|
|---|
| 293 | // //profiling: t.stop();
|
|---|
| 294 | // //profiling: cout << " finishing: " << t.total() << endl;
|
|---|
| 295 | //
|
|---|
| 296 | // return reduced_event;
|
|---|
| 297 | }
|
|---|
| 298 |
|
|---|
| 299 | #ifndef FJCONTRIB_SOFTKILLER_USEFJGRID
|
|---|
| 300 |
|
|---|
| 301 | // configure the grid
|
|---|
| 302 | void SoftKiller::_setup_grid() {
|
|---|
| 303 | // this grid-definition code is becoming repetitive -- it should
|
|---|
| 304 | // probably be moved somewhere central...
|
|---|
| 305 | double ny_double = (_ymax-_ymin) / _requested_drap;
|
|---|
| 306 | _ny = max(int(ny_double+0.5),1);
|
|---|
| 307 | _dy = (_ymax-_ymin) / _ny;
|
|---|
| 308 | _inverse_dy = _ny/(_ymax-_ymin);
|
|---|
| 309 |
|
|---|
| 310 | _nphi = int (twopi / _requested_dphi + 0.5);
|
|---|
| 311 | _dphi = twopi / _nphi;
|
|---|
| 312 | _inverse_dphi = _nphi/twopi;
|
|---|
| 313 |
|
|---|
| 314 | // some sanity checking (could throw a fastjet::Error)
|
|---|
| 315 | assert(_ny >= 1 && _nphi >= 1);
|
|---|
| 316 |
|
|---|
| 317 | _ntotal = _nphi * _ny;
|
|---|
| 318 | //_max_pt.resize(_ntotal);
|
|---|
| 319 | _cell_area = _dy * _dphi;
|
|---|
| 320 | }
|
|---|
| 321 |
|
|---|
| 322 | // retrieve the grid cell index for a given PseudoJet
|
|---|
| 323 | inline int SoftKiller::tile_index(const PseudoJet & p) const {
|
|---|
| 324 | // writing it as below gives a huge speed gain (factor two!). Even
|
|---|
| 325 | // though answers are identical and the routine here is not the
|
|---|
| 326 | // speed-critical step. It's not at all clear why.
|
|---|
| 327 | int iy = int(floor( (p.rap() - _ymin) * _inverse_dy ));
|
|---|
| 328 | if (iy < 0 || iy >= _ny) return -1;
|
|---|
| 329 |
|
|---|
| 330 | int iphi = int( p.phi() * _inverse_dphi );
|
|---|
| 331 | //assert(iphi >= 0 && iphi <= _nphi);
|
|---|
| 332 | if (iphi == _nphi) iphi = 0; // just in case of rounding errors
|
|---|
| 333 |
|
|---|
| 334 | //int igrid_res = iy*_nphi + iphi;
|
|---|
| 335 | //assert (igrid_res >= 0 && igrid_res < _ny*_nphi);
|
|---|
| 336 | return iy*_nphi + iphi; //igrid_res;
|
|---|
| 337 | }
|
|---|
| 338 |
|
|---|
| 339 | #endif // not FJCONTRIB_SOFTKILLER_USEFJGRID
|
|---|
| 340 |
|
|---|
| 341 |
|
|---|
| 342 |
|
|---|
| 343 | } // namespace contrib
|
|---|
| 344 |
|
|---|
| 345 | FASTJET_END_NAMESPACE
|
|---|
