| 1 | //STARTHEADER
|
|---|
| 2 | // $Id$
|
|---|
| 3 | //
|
|---|
| 4 | // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
|
|---|
| 5 | //
|
|---|
| 6 | //----------------------------------------------------------------------
|
|---|
| 7 | // This file is part of FastJet.
|
|---|
| 8 | //
|
|---|
| 9 | // FastJet is free software; you can redistribute it and/or modify
|
|---|
| 10 | // it under the terms of the GNU General Public License as published by
|
|---|
| 11 | // the Free Software Foundation; either version 2 of the License, or
|
|---|
| 12 | // (at your option) any later version.
|
|---|
| 13 | //
|
|---|
| 14 | // The algorithms that underlie FastJet have required considerable
|
|---|
| 15 | // development and are described in hep-ph/0512210. If you use
|
|---|
| 16 | // FastJet as part of work towards a scientific publication, please
|
|---|
| 17 | // include a citation to the FastJet paper.
|
|---|
| 18 | //
|
|---|
| 19 | // FastJet is distributed in the hope that it will be useful,
|
|---|
| 20 | // but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|---|
| 21 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|---|
| 22 | // GNU General Public License for more details.
|
|---|
| 23 | //
|
|---|
| 24 | // You should have received a copy of the GNU General Public License
|
|---|
| 25 | // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
|
|---|
| 26 | //----------------------------------------------------------------------
|
|---|
| 27 | //ENDHEADER
|
|---|
| 28 |
|
|---|
| 29 | #include "fastjet/PseudoJet.hh"
|
|---|
| 30 | #include "fastjet/ClusterSequence.hh"
|
|---|
| 31 | #include "fastjet/ClusterSequenceActiveArea.hh"
|
|---|
| 32 | #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
|
|---|
| 33 | #include<iostream>
|
|---|
| 34 | #include<vector>
|
|---|
| 35 | #include<sstream>
|
|---|
| 36 | #include<algorithm>
|
|---|
| 37 | #include<cmath>
|
|---|
| 38 | #include<valarray>
|
|---|
| 39 |
|
|---|
| 40 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
|
|---|
| 41 |
|
|---|
| 42 |
|
|---|
| 43 | using namespace std;
|
|---|
| 44 |
|
|---|
| 45 |
|
|---|
| 46 | //int ClusterSequenceActiveArea::_n_seed_warnings = 0;
|
|---|
| 47 | //const int _max_seed_warnings = 10;
|
|---|
| 48 |
|
|---|
| 49 | //----------------------------------------------------------------------
|
|---|
| 50 | /// global routine for running active area
|
|---|
| 51 | void ClusterSequenceActiveArea::_initialise_and_run_AA (
|
|---|
| 52 | const JetDefinition & jet_def_in,
|
|---|
| 53 | const GhostedAreaSpec & ghost_spec,
|
|---|
| 54 | const bool & writeout_combinations) {
|
|---|
| 55 |
|
|---|
| 56 | bool continue_running;
|
|---|
| 57 | _initialise_AA(jet_def_in, ghost_spec, writeout_combinations, continue_running);
|
|---|
| 58 | if (continue_running) {
|
|---|
| 59 | _run_AA(ghost_spec);
|
|---|
| 60 | _postprocess_AA(ghost_spec);
|
|---|
| 61 | }
|
|---|
| 62 | }
|
|---|
| 63 |
|
|---|
| 64 | //----------------------------------------------------------------------
|
|---|
| 65 | void ClusterSequenceActiveArea::_resize_and_zero_AA () {
|
|---|
| 66 | // initialize our local area information
|
|---|
| 67 | _average_area.resize(2*_jets.size()); _average_area = 0.0;
|
|---|
| 68 | _average_area2.resize(2*_jets.size()); _average_area2 = 0.0;
|
|---|
| 69 | _average_area_4vector.resize(2*_jets.size());
|
|---|
| 70 | _average_area_4vector = PseudoJet(0.0,0.0,0.0,0.0);
|
|---|
| 71 | _non_jet_area = 0.0; _non_jet_area2 = 0.0; _non_jet_number=0.0;
|
|---|
| 72 | }
|
|---|
| 73 |
|
|---|
| 74 | //---------------------------------a-------------------------------------
|
|---|
| 75 | void ClusterSequenceActiveArea::_initialise_AA (
|
|---|
| 76 | const JetDefinition & jet_def_in,
|
|---|
| 77 | const GhostedAreaSpec & ghost_spec,
|
|---|
| 78 | const bool & writeout_combinations,
|
|---|
| 79 | bool & continue_running)
|
|---|
| 80 | {
|
|---|
| 81 |
|
|---|
| 82 | // store this for future use
|
|---|
| 83 | _ghost_spec_repeat = ghost_spec.repeat();
|
|---|
| 84 |
|
|---|
| 85 | // make sure placeholders are there & zeroed
|
|---|
| 86 | _resize_and_zero_AA();
|
|---|
| 87 |
|
|---|
| 88 | // for future reference...
|
|---|
| 89 | _maxrap_for_area = ghost_spec.ghost_maxrap();
|
|---|
| 90 | _safe_rap_for_area = _maxrap_for_area - jet_def_in.R();
|
|---|
| 91 |
|
|---|
| 92 | // Make sure we'll have at least one repetition -- then we can
|
|---|
| 93 | // deduce the unghosted clustering sequence from one of the ghosted
|
|---|
| 94 | // sequences. If we do not have any repetitions, then get the
|
|---|
| 95 | // unghosted sequence from the plain unghosted clustering.
|
|---|
| 96 | //
|
|---|
| 97 | // NB: all decanting and filling of initial history will then
|
|---|
| 98 | // be carried out by base-class routine
|
|---|
| 99 | if (ghost_spec.repeat() <= 0) {
|
|---|
| 100 | _initialise_and_run(jet_def_in, writeout_combinations);
|
|---|
| 101 | continue_running = false;
|
|---|
| 102 | return;
|
|---|
| 103 | }
|
|---|
| 104 |
|
|---|
| 105 | // transfer all relevant info into internal variables
|
|---|
| 106 | _decant_options(jet_def_in, writeout_combinations);
|
|---|
| 107 |
|
|---|
| 108 | // set up the history entries for the initial particles (those
|
|---|
| 109 | // currently in _jets)
|
|---|
| 110 | _fill_initial_history();
|
|---|
| 111 |
|
|---|
| 112 | // by default it does not...
|
|---|
| 113 | _has_dangerous_particles = false;
|
|---|
| 114 |
|
|---|
| 115 | continue_running = true;
|
|---|
| 116 | }
|
|---|
| 117 |
|
|---|
| 118 |
|
|---|
| 119 | //----------------------------------------------------------------------
|
|---|
| 120 | void ClusterSequenceActiveArea::_run_AA (const GhostedAreaSpec & ghost_spec) {
|
|---|
| 121 | // record the input jets as they are currently
|
|---|
| 122 | vector<PseudoJet> input_jets(_jets);
|
|---|
| 123 |
|
|---|
| 124 | // code for testing the unique tree
|
|---|
| 125 | vector<int> unique_tree;
|
|---|
| 126 |
|
|---|
| 127 | // run the clustering multiple times so as to get areas of all the jets
|
|---|
| 128 | for (int irepeat = 0; irepeat < ghost_spec.repeat(); irepeat++) {
|
|---|
| 129 |
|
|---|
| 130 | ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets,
|
|---|
| 131 | jet_def(), ghost_spec);
|
|---|
| 132 |
|
|---|
| 133 | _has_dangerous_particles |= clust_seq.has_dangerous_particles();
|
|---|
| 134 | if (irepeat == 0) {
|
|---|
| 135 | // take the non-ghost part of the history and put into our own
|
|---|
| 136 | // history.
|
|---|
| 137 | _transfer_ghost_free_history(clust_seq);
|
|---|
| 138 | // get the "unique" order that will be used for transferring all areas.
|
|---|
| 139 | unique_tree = unique_history_order();
|
|---|
| 140 | }
|
|---|
| 141 |
|
|---|
| 142 | // transfer areas from clust_seq into our object
|
|---|
| 143 | _transfer_areas(unique_tree, clust_seq);
|
|---|
| 144 | }
|
|---|
| 145 | }
|
|---|
| 146 |
|
|---|
| 147 |
|
|---|
| 148 | //----------------------------------------------------------------------
|
|---|
| 149 | /// run the postprocessing for the active area (and derived classes)
|
|---|
| 150 | void ClusterSequenceActiveArea::_postprocess_AA (const GhostedAreaSpec & ghost_spec) {
|
|---|
| 151 | _average_area /= ghost_spec.repeat();
|
|---|
| 152 | _average_area2 /= ghost_spec.repeat();
|
|---|
| 153 | if (ghost_spec.repeat() > 1) {
|
|---|
| 154 | // the VC compiler complains if one puts everything on a single line.
|
|---|
| 155 | // An alternative solution would be to use -1.0 (+single line)
|
|---|
| 156 | const double tmp = ghost_spec.repeat()-1;
|
|---|
| 157 | _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/tmp);
|
|---|
| 158 | } else {
|
|---|
| 159 | _average_area2 = 0.0;
|
|---|
| 160 | }
|
|---|
| 161 |
|
|---|
| 162 | _non_jet_area /= ghost_spec.repeat();
|
|---|
| 163 | _non_jet_area2 /= ghost_spec.repeat();
|
|---|
| 164 | _non_jet_area2 = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/
|
|---|
| 165 | ghost_spec.repeat());
|
|---|
| 166 | _non_jet_number /= ghost_spec.repeat();
|
|---|
| 167 |
|
|---|
| 168 | // following bizarre way of writing things is related to
|
|---|
| 169 | // poverty of operations on PseudoJet objects (as well as some confusion
|
|---|
| 170 | // in one or two places)
|
|---|
| 171 | for (unsigned i = 0; i < _average_area_4vector.size(); i++) {
|
|---|
| 172 | _average_area_4vector[i] = (1.0/ghost_spec.repeat()) * _average_area_4vector[i];
|
|---|
| 173 | }
|
|---|
| 174 | //cerr << "Non-jet area = " << _non_jet_area << " +- " << _non_jet_area2<<endl;
|
|---|
| 175 | }
|
|---|
| 176 |
|
|---|
| 177 |
|
|---|
| 178 | // //----------------------------------------------------------------------
|
|---|
| 179 | // void ClusterSequenceActiveArea::_initialise_and_run_AA (
|
|---|
| 180 | // const JetDefinition & jet_def,
|
|---|
| 181 | // const GhostedAreaSpec & ghost_spec,
|
|---|
| 182 | // const bool & writeout_combinations)
|
|---|
| 183 | // {
|
|---|
| 184 | //
|
|---|
| 185 | // // store this for future use
|
|---|
| 186 | // _ghost_spec_repeat = ghost_spec.repeat();
|
|---|
| 187 | //
|
|---|
| 188 | // // initialize our local area information
|
|---|
| 189 | // _average_area.resize(2*_jets.size()); _average_area = 0.0;
|
|---|
| 190 | // _average_area2.resize(2*_jets.size()); _average_area2 = 0.0;
|
|---|
| 191 | // _average_area_4vector.resize(2*_jets.size());
|
|---|
| 192 | // _average_area_4vector = PseudoJet(0.0,0.0,0.0,0.0);
|
|---|
| 193 | // _non_jet_area = 0.0; _non_jet_area2 = 0.0; _non_jet_number=0.0;
|
|---|
| 194 | //
|
|---|
| 195 | // // for future reference...
|
|---|
| 196 | // _maxrap_for_area = ghost_spec.ghost_maxrap();
|
|---|
| 197 | // _safe_rap_for_area = _maxrap_for_area - jet_def.R();
|
|---|
| 198 | //
|
|---|
| 199 | // // Make sure we'll have at least one repetition -- then we can
|
|---|
| 200 | // // deduce the unghosted clustering sequence from one of the ghosted
|
|---|
| 201 | // // sequences. If we do not have any repetitions, then get the
|
|---|
| 202 | // // unghosted sequence from the plain unghosted clustering.
|
|---|
| 203 | // //
|
|---|
| 204 | // // NB: all decanting and filling of initial history will then
|
|---|
| 205 | // // be carried out by base-class routine
|
|---|
| 206 | // if (ghost_spec.repeat() <= 0) {
|
|---|
| 207 | // _initialise_and_run(jet_def, writeout_combinations);
|
|---|
| 208 | // return;
|
|---|
| 209 | // }
|
|---|
| 210 | //
|
|---|
| 211 | // // transfer all relevant info into internal variables
|
|---|
| 212 | // _decant_options(jet_def, writeout_combinations);
|
|---|
| 213 | //
|
|---|
| 214 | // // set up the history entries for the initial particles (those
|
|---|
| 215 | // // currently in _jets)
|
|---|
| 216 | // _fill_initial_history();
|
|---|
| 217 | //
|
|---|
| 218 | // // record the input jets as they are currently
|
|---|
| 219 | // vector<PseudoJet> input_jets(_jets);
|
|---|
| 220 | //
|
|---|
| 221 | // // code for testing the unique tree
|
|---|
| 222 | // vector<int> unique_tree;
|
|---|
| 223 | //
|
|---|
| 224 | //
|
|---|
| 225 | //
|
|---|
| 226 | //
|
|---|
| 227 | // // run the clustering multiple times so as to get areas of all the jets
|
|---|
| 228 | // for (int irepeat = 0; irepeat < ghost_spec.repeat(); irepeat++) {
|
|---|
| 229 | //
|
|---|
| 230 | // ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets,
|
|---|
| 231 | // jet_def, ghost_spec);
|
|---|
| 232 | //
|
|---|
| 233 | // if (irepeat == 0) {
|
|---|
| 234 | // // take the non-ghost part of the history and put into our own
|
|---|
| 235 | // // history.
|
|---|
| 236 | // _transfer_ghost_free_history(clust_seq);
|
|---|
| 237 | // // get the "unique" order that will be used for transferring all areas.
|
|---|
| 238 | // unique_tree = unique_history_order();
|
|---|
| 239 | // }
|
|---|
| 240 | //
|
|---|
| 241 | // // transfer areas from clust_seq into our object
|
|---|
| 242 | // _transfer_areas(unique_tree, clust_seq);
|
|---|
| 243 | // }
|
|---|
| 244 | //
|
|---|
| 245 | // _average_area /= ghost_spec.repeat();
|
|---|
| 246 | // _average_area2 /= ghost_spec.repeat();
|
|---|
| 247 | // if (ghost_spec.repeat() > 1) {
|
|---|
| 248 | // _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/
|
|---|
| 249 | // (ghost_spec.repeat()-1));
|
|---|
| 250 | // } else {
|
|---|
| 251 | // _average_area2 = 0.0;
|
|---|
| 252 | // }
|
|---|
| 253 | //
|
|---|
| 254 | // _non_jet_area /= ghost_spec.repeat();
|
|---|
| 255 | // _non_jet_area2 /= ghost_spec.repeat();
|
|---|
| 256 | // _non_jet_area2 = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/
|
|---|
| 257 | // ghost_spec.repeat());
|
|---|
| 258 | // _non_jet_number /= ghost_spec.repeat();
|
|---|
| 259 | //
|
|---|
| 260 | // // following bizarre way of writing things is related to
|
|---|
| 261 | // // poverty of operations on PseudoJet objects (as well as some confusion
|
|---|
| 262 | // // in one or two places)
|
|---|
| 263 | // for (unsigned i = 0; i < _average_area_4vector.size(); i++) {
|
|---|
| 264 | // _average_area_4vector[i] = (1.0/ghost_spec.repeat()) * _average_area_4vector[i];
|
|---|
| 265 | // }
|
|---|
| 266 | // //cerr << "Non-jet area = " << _non_jet_area << " +- " << _non_jet_area2<<endl;
|
|---|
| 267 | //
|
|---|
| 268 | //
|
|---|
| 269 | // }
|
|---|
| 270 | //
|
|---|
| 271 |
|
|---|
| 272 |
|
|---|
| 273 | //----------------------------------------------------------------------
|
|---|
| 274 | double ClusterSequenceActiveArea::pt_per_unit_area(
|
|---|
| 275 | mean_pt_strategies strat, double range) const {
|
|---|
| 276 |
|
|---|
| 277 | vector<PseudoJet> incl_jets = inclusive_jets();
|
|---|
| 278 | vector<double> pt_over_areas;
|
|---|
| 279 |
|
|---|
| 280 | for (unsigned i = 0; i < incl_jets.size(); i++) {
|
|---|
| 281 | if (abs(incl_jets[i].rap()) < _safe_rap_for_area) {
|
|---|
| 282 | double this_area;
|
|---|
| 283 | if ( strat == median_4vector ) {
|
|---|
| 284 | this_area = area_4vector(incl_jets[i]).perp();
|
|---|
| 285 | } else {
|
|---|
| 286 | this_area = area(incl_jets[i]);
|
|---|
| 287 | }
|
|---|
| 288 | pt_over_areas.push_back(incl_jets[i].perp()/this_area);
|
|---|
| 289 | }
|
|---|
| 290 | }
|
|---|
| 291 |
|
|---|
| 292 | // there is nothing inside our region, so answer will always be zero
|
|---|
| 293 | if (pt_over_areas.size() == 0) {return 0.0;}
|
|---|
| 294 |
|
|---|
| 295 | // get median (pt/area) [this is the "old" median definition. It considers
|
|---|
| 296 | // only the "real" jets in calculating the median, i.e. excluding the
|
|---|
| 297 | // only-ghost ones]
|
|---|
| 298 | sort(pt_over_areas.begin(), pt_over_areas.end());
|
|---|
| 299 | double non_ghost_median_ratio = pt_over_areas[pt_over_areas.size()/2];
|
|---|
| 300 |
|
|---|
| 301 | // new median definition that takes into account non-jet area (i.e.
|
|---|
| 302 | // jets composed only of ghosts), and for fractional median position
|
|---|
| 303 | // interpolates between the corresponding entries in the pt_over_areas array
|
|---|
| 304 | double nj_median_pos = (pt_over_areas.size()-1 - _non_jet_number)/2.0;
|
|---|
| 305 | double nj_median_ratio;
|
|---|
| 306 | if (nj_median_pos >= 0 && pt_over_areas.size() > 1) {
|
|---|
| 307 | int int_nj_median = int(nj_median_pos);
|
|---|
| 308 | nj_median_ratio =
|
|---|
| 309 | pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos)
|
|---|
| 310 | + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median);
|
|---|
| 311 | } else {
|
|---|
| 312 | nj_median_ratio = 0.0;
|
|---|
| 313 | }
|
|---|
| 314 |
|
|---|
| 315 |
|
|---|
| 316 | // get various forms of mean (pt/area)
|
|---|
| 317 | double pt_sum = 0.0, pt_sum_with_cut = 0.0;
|
|---|
| 318 | double area_sum = _non_jet_area, area_sum_with_cut = _non_jet_area;
|
|---|
| 319 | double ratio_sum = 0.0;
|
|---|
| 320 | double ratio_n = _non_jet_number;
|
|---|
| 321 | for (unsigned i = 0; i < incl_jets.size(); i++) {
|
|---|
| 322 | if (abs(incl_jets[i].rap()) < _safe_rap_for_area) {
|
|---|
| 323 | double this_area;
|
|---|
| 324 | if ( strat == median_4vector ) {
|
|---|
| 325 | this_area = area_4vector(incl_jets[i]).perp();
|
|---|
| 326 | } else {
|
|---|
| 327 | this_area = area(incl_jets[i]);
|
|---|
| 328 | }
|
|---|
| 329 | pt_sum += incl_jets[i].perp();
|
|---|
| 330 | area_sum += this_area;
|
|---|
| 331 | double ratio = incl_jets[i].perp()/this_area;
|
|---|
| 332 | if (ratio < range*nj_median_ratio) {
|
|---|
| 333 | pt_sum_with_cut += incl_jets[i].perp();
|
|---|
| 334 | area_sum_with_cut += this_area;
|
|---|
| 335 | ratio_sum += ratio; ratio_n++;
|
|---|
| 336 | }
|
|---|
| 337 | }
|
|---|
| 338 | }
|
|---|
| 339 |
|
|---|
| 340 | if (strat == play) {
|
|---|
| 341 | double trunc_sum = 0, trunc_sumsqr = 0;
|
|---|
| 342 | vector<double> means(pt_over_areas.size()), sd(pt_over_areas.size());
|
|---|
| 343 | for (unsigned i = 0; i < pt_over_areas.size() ; i++ ) {
|
|---|
| 344 | double ratio = pt_over_areas[i];
|
|---|
| 345 | trunc_sum += ratio;
|
|---|
| 346 | trunc_sumsqr += ratio*ratio;
|
|---|
| 347 | means[i] = trunc_sum / (i+1);
|
|---|
| 348 | sd[i] = sqrt(abs(means[i]*means[i] - trunc_sumsqr/(i+1)));
|
|---|
| 349 | cerr << "i, means, sd: " <<i<<", "<< means[i] <<", "<<sd[i]<<", "<<
|
|---|
| 350 | sd[i]/sqrt(i+1.0)<<endl;
|
|---|
| 351 | }
|
|---|
| 352 | cout << "-----------------------------------"<<endl;
|
|---|
| 353 | for (unsigned i = 0; i <= pt_over_areas.size()/2 ; i++ ) {
|
|---|
| 354 | cout << "Median "<< i <<" = " << pt_over_areas[i]<<endl;
|
|---|
| 355 | }
|
|---|
| 356 | cout << "Number of non-jets: "<<_non_jet_number<<endl;
|
|---|
| 357 | cout << "Area of non-jets: "<<_non_jet_area<<endl;
|
|---|
| 358 | cout << "Default median position: " << (pt_over_areas.size()-1)/2.0<<endl;
|
|---|
| 359 | cout << "NJ median position: " << nj_median_pos <<endl;
|
|---|
| 360 | cout << "NJ median value: " << nj_median_ratio <<endl;
|
|---|
| 361 | return 0.0;
|
|---|
| 362 | }
|
|---|
| 363 |
|
|---|
| 364 | switch(strat) {
|
|---|
| 365 | case median:
|
|---|
| 366 | case median_4vector:
|
|---|
| 367 | return nj_median_ratio;
|
|---|
| 368 | case non_ghost_median:
|
|---|
| 369 | return non_ghost_median_ratio;
|
|---|
| 370 | case pttot_over_areatot:
|
|---|
| 371 | return pt_sum / area_sum;
|
|---|
| 372 | case pttot_over_areatot_cut:
|
|---|
| 373 | return pt_sum_with_cut / area_sum_with_cut;
|
|---|
| 374 | case mean_ratio_cut:
|
|---|
| 375 | return ratio_sum/ratio_n;
|
|---|
| 376 | default:
|
|---|
| 377 | return nj_median_ratio;
|
|---|
| 378 | }
|
|---|
| 379 |
|
|---|
| 380 | }
|
|---|
| 381 |
|
|---|
| 382 |
|
|---|
| 383 | // The following functionality is now provided by the base class
|
|---|
| 384 | // //----------------------------------------------------------------------
|
|---|
| 385 | // // fit a parabola to pt/area as a function of rapidity, using the
|
|---|
| 386 | // // formulae of CCN28-36 (which actually fits f = a+b*x^2)
|
|---|
| 387 | // void ClusterSequenceActiveArea::parabolic_pt_per_unit_area(
|
|---|
| 388 | // double & a, double & b, double raprange, double exclude_above,
|
|---|
| 389 | // bool use_area_4vector) const {
|
|---|
| 390 | //
|
|---|
| 391 | // double this_raprange;
|
|---|
| 392 | // if (raprange <= 0) {this_raprange = _safe_rap_for_area;}
|
|---|
| 393 | // else {this_raprange = raprange;}
|
|---|
| 394 | //
|
|---|
| 395 | // int n=0;
|
|---|
| 396 | // int n_excluded = 0;
|
|---|
| 397 | // double mean_f=0, mean_x2=0, mean_x4=0, mean_fx2=0;
|
|---|
| 398 | //
|
|---|
| 399 | // vector<PseudoJet> incl_jets = inclusive_jets();
|
|---|
| 400 | //
|
|---|
| 401 | // for (unsigned i = 0; i < incl_jets.size(); i++) {
|
|---|
| 402 | // if (abs(incl_jets[i].rap()) < this_raprange) {
|
|---|
| 403 | // double this_area;
|
|---|
| 404 | // if ( use_area_4vector ) {
|
|---|
| 405 | // this_area = area_4vector(incl_jets[i]).perp();
|
|---|
| 406 | // } else {
|
|---|
| 407 | // this_area = area(incl_jets[i]);
|
|---|
| 408 | // }
|
|---|
| 409 | // double f = incl_jets[i].perp()/this_area;
|
|---|
| 410 | // if (exclude_above <= 0.0 || f < exclude_above) {
|
|---|
| 411 | // double x = incl_jets[i].rap(); double x2 = x*x;
|
|---|
| 412 | // mean_f += f;
|
|---|
| 413 | // mean_x2 += x2;
|
|---|
| 414 | // mean_x4 += x2*x2;
|
|---|
| 415 | // mean_fx2 += f*x2;
|
|---|
| 416 | // n++;
|
|---|
| 417 | // } else {
|
|---|
| 418 | // n_excluded++;
|
|---|
| 419 | // }
|
|---|
| 420 | // }
|
|---|
| 421 | // }
|
|---|
| 422 | //
|
|---|
| 423 | // if (n <= 1) {
|
|---|
| 424 | // // meaningful results require at least two jets inside the
|
|---|
| 425 | // // area -- mind you if there are empty jets we should be in
|
|---|
| 426 | // // any case doing something special...
|
|---|
| 427 | // a = 0.0;
|
|---|
| 428 | // b = 0.0;
|
|---|
| 429 | // } else {
|
|---|
| 430 | // mean_f /= n;
|
|---|
| 431 | // mean_x2 /= n;
|
|---|
| 432 | // mean_x4 /= n;
|
|---|
| 433 | // mean_fx2 /= n;
|
|---|
| 434 | //
|
|---|
| 435 | // b = (mean_f*mean_x2 - mean_fx2)/(mean_x2*mean_x2 - mean_x4);
|
|---|
| 436 | // a = mean_f - b*mean_x2;
|
|---|
| 437 | // }
|
|---|
| 438 | // //cerr << "n_excluded = "<< n_excluded << endl;
|
|---|
| 439 | // }
|
|---|
| 440 |
|
|---|
| 441 |
|
|---|
| 442 | //----------------------------------------------------------------------
|
|---|
| 443 | double ClusterSequenceActiveArea::empty_area(const Selector & selector) const {
|
|---|
| 444 | // make sure that the selector applies jet by jet
|
|---|
| 445 | if (! selector.applies_jet_by_jet()){
|
|---|
| 446 | throw Error("ClusterSequenceActiveArea: empty area can only be computed from selectors applying jet by jet");
|
|---|
| 447 | }
|
|---|
| 448 |
|
|---|
| 449 | double empty = 0.0;
|
|---|
| 450 | // first deal with ghost jets
|
|---|
| 451 | for (unsigned i = 0; i < _ghost_jets.size(); i++) {
|
|---|
| 452 | if (selector.pass(_ghost_jets[i])) {
|
|---|
| 453 | empty += _ghost_jets[i].area;
|
|---|
| 454 | }
|
|---|
| 455 | }
|
|---|
| 456 | // then deal with unclustered ghosts
|
|---|
| 457 | for (unsigned i = 0; i < _unclustered_ghosts.size(); i++) {
|
|---|
| 458 | if (selector.pass(_unclustered_ghosts[i])) {
|
|---|
| 459 | empty += _unclustered_ghosts[i].area;
|
|---|
| 460 | }
|
|---|
| 461 | }
|
|---|
| 462 | empty /= _ghost_spec_repeat;
|
|---|
| 463 | return empty;
|
|---|
| 464 | }
|
|---|
| 465 |
|
|---|
| 466 | //----------------------------------------------------------------------
|
|---|
| 467 | double ClusterSequenceActiveArea::n_empty_jets(const Selector & selector) const {
|
|---|
| 468 | _check_selector_good_for_median(selector);
|
|---|
| 469 |
|
|---|
| 470 | double inrange = 0;
|
|---|
| 471 | for (unsigned i = 0; i < _ghost_jets.size(); i++) {
|
|---|
| 472 | if (selector.pass(_ghost_jets[i])) inrange++;
|
|---|
| 473 | }
|
|---|
| 474 | inrange /= _ghost_spec_repeat;
|
|---|
| 475 | return inrange;
|
|---|
| 476 | }
|
|---|
| 477 |
|
|---|
| 478 | //----------------------------------------------------------------------
|
|---|
| 479 | /// transfer the history (and jet-momenta) from clust_seq to our
|
|---|
| 480 | /// own internal structure while removing ghosts
|
|---|
| 481 | void ClusterSequenceActiveArea::_transfer_ghost_free_history(
|
|---|
| 482 | const ClusterSequenceActiveAreaExplicitGhosts & ghosted_seq) {
|
|---|
| 483 |
|
|---|
| 484 | const vector<history_element> & gs_history = ghosted_seq.history();
|
|---|
| 485 | vector<int> gs2self_hist_map(gs_history.size());
|
|---|
| 486 |
|
|---|
| 487 | // first transfer info about strategy used (which isn't necessarily
|
|---|
| 488 | // always the one that got asked for...)
|
|---|
| 489 | _strategy = ghosted_seq.strategy_used();
|
|---|
| 490 |
|
|---|
| 491 | // work our way through to first non-trivial combination
|
|---|
| 492 | unsigned igs = 0;
|
|---|
| 493 | unsigned iself = 0;
|
|---|
| 494 | while (igs < gs_history.size() && gs_history[igs].parent1 == InexistentParent) {
|
|---|
| 495 | // record correspondence
|
|---|
| 496 | if (!ghosted_seq.is_pure_ghost(igs)) {
|
|---|
| 497 | gs2self_hist_map[igs] = iself++;
|
|---|
| 498 | } else {
|
|---|
| 499 | gs2self_hist_map[igs] = Invalid;
|
|---|
| 500 | }
|
|---|
| 501 | igs++;
|
|---|
| 502 | };
|
|---|
| 503 |
|
|---|
| 504 | // make sure the count of non-ghost initial jets is equal to
|
|---|
| 505 | // what we already have in terms of initial jets
|
|---|
| 506 | assert(iself == _history.size());
|
|---|
| 507 |
|
|---|
| 508 | // if there was no clustering in this event (e.g. SISCone passive
|
|---|
| 509 | // area with zero input particles, or with a pt cut on stable cones
|
|---|
| 510 | // that kills all jets), then don't bother with the rest (which
|
|---|
| 511 | // would crash!)
|
|---|
| 512 | if (igs == gs_history.size()) return;
|
|---|
| 513 |
|
|---|
| 514 | // now actually transfer things
|
|---|
| 515 | do {
|
|---|
| 516 | // if we are a pure ghost, then go on to next round
|
|---|
| 517 | if (ghosted_seq.is_pure_ghost(igs)) {
|
|---|
| 518 | gs2self_hist_map[igs] = Invalid;
|
|---|
| 519 | continue;
|
|---|
| 520 | }
|
|---|
| 521 |
|
|---|
| 522 | const history_element & gs_hist_el = gs_history[igs];
|
|---|
| 523 |
|
|---|
| 524 | bool parent1_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent1);
|
|---|
| 525 | bool parent2_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent2);
|
|---|
| 526 |
|
|---|
| 527 | // if exactly one parent is a ghost then maintain info about the
|
|---|
| 528 | // non-ghost correspondence for this jet, and then go on to next
|
|---|
| 529 | // recombination in the ghosted sequence
|
|---|
| 530 | if (parent1_is_ghost && !parent2_is_ghost && gs_hist_el.parent2 >= 0) {
|
|---|
| 531 | gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent2];
|
|---|
| 532 | continue;
|
|---|
| 533 | }
|
|---|
| 534 | if (!parent1_is_ghost && parent2_is_ghost) {
|
|---|
| 535 | gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent1];
|
|---|
| 536 | continue;
|
|---|
| 537 | }
|
|---|
| 538 |
|
|---|
| 539 | // no parents are ghosts...
|
|---|
| 540 | if (gs_hist_el.parent2 >= 0) {
|
|---|
| 541 | // recombination of two non-ghosts
|
|---|
| 542 | gs2self_hist_map[igs] = _history.size();
|
|---|
| 543 | // record the recombination in our own sequence
|
|---|
| 544 | int newjet_k; // dummy var -- not used
|
|---|
| 545 | int jet_i = _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index;
|
|---|
| 546 | int jet_j = _history[gs2self_hist_map[gs_hist_el.parent2]].jetp_index;
|
|---|
| 547 | //cerr << "recombining "<< jet_i << " and "<< jet_j << endl;
|
|---|
| 548 | _do_ij_recombination_step(jet_i, jet_j, gs_hist_el.dij, newjet_k);
|
|---|
| 549 | } else {
|
|---|
| 550 | // we have a non-ghost that has become a beam-jet
|
|---|
| 551 | assert(gs_history[igs].parent2 == BeamJet);
|
|---|
| 552 | // record position
|
|---|
| 553 | gs2self_hist_map[igs] = _history.size();
|
|---|
| 554 | // record the recombination in our own sequence
|
|---|
| 555 | _do_iB_recombination_step(
|
|---|
| 556 | _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index,
|
|---|
| 557 | gs_hist_el.dij);
|
|---|
| 558 | }
|
|---|
| 559 | } while (++igs < gs_history.size());
|
|---|
| 560 |
|
|---|
| 561 | }
|
|---|
| 562 |
|
|---|
| 563 | //----------------------------------------------------------------------
|
|---|
| 564 | void ClusterSequenceActiveArea::_transfer_areas(
|
|---|
| 565 | const vector<int> & unique_hist_order,
|
|---|
| 566 | const ClusterSequenceActiveAreaExplicitGhosts & ghosted_seq ) {
|
|---|
| 567 |
|
|---|
| 568 | const vector<history_element> & gs_history = ghosted_seq.history();
|
|---|
| 569 | const vector<PseudoJet> & gs_jets = ghosted_seq.jets();
|
|---|
| 570 | vector<int> gs_unique_hist_order = ghosted_seq.unique_history_order();
|
|---|
| 571 |
|
|---|
| 572 | const double tolerance = 1e-11; // to decide when two jets are the same
|
|---|
| 573 |
|
|---|
| 574 | int j = -1;
|
|---|
| 575 | int hist_index = -1;
|
|---|
| 576 |
|
|---|
| 577 | valarray<double> our_areas(_history.size());
|
|---|
| 578 | our_areas = 0.0;
|
|---|
| 579 |
|
|---|
| 580 | valarray<PseudoJet> our_area_4vectors(_history.size());
|
|---|
| 581 | our_area_4vectors = PseudoJet(0.0,0.0,0.0,0.0);
|
|---|
| 582 |
|
|---|
| 583 | for (unsigned i = 0; i < gs_history.size(); i++) {
|
|---|
| 584 | // only consider composite particles
|
|---|
| 585 | unsigned gs_hist_index = gs_unique_hist_order[i];
|
|---|
| 586 | if (gs_hist_index < ghosted_seq.n_particles()) continue;
|
|---|
| 587 | const history_element & gs_hist = gs_history[gs_unique_hist_order[i]];
|
|---|
| 588 | int parent1 = gs_hist.parent1;
|
|---|
| 589 | int parent2 = gs_hist.parent2;
|
|---|
| 590 |
|
|---|
| 591 | if (parent2 == BeamJet) {
|
|---|
| 592 | // need to look at parent to get the actual jet
|
|---|
| 593 | const PseudoJet & jet =
|
|---|
| 594 | gs_jets[gs_history[parent1].jetp_index];
|
|---|
| 595 | double area_local = ghosted_seq.area(jet);
|
|---|
| 596 | PseudoJet ext_area = ghosted_seq.area_4vector(jet);
|
|---|
| 597 |
|
|---|
| 598 | if (ghosted_seq.is_pure_ghost(parent1)) {
|
|---|
| 599 | // record the existence of the pure ghost jet for future use
|
|---|
| 600 | _ghost_jets.push_back(GhostJet(jet,area_local));
|
|---|
| 601 | if (abs(jet.rap()) < _safe_rap_for_area) {
|
|---|
| 602 | _non_jet_area += area_local;
|
|---|
| 603 | _non_jet_area2 += area_local*area_local;
|
|---|
| 604 | _non_jet_number += 1;
|
|---|
| 605 | }
|
|---|
| 606 | } else {
|
|---|
| 607 |
|
|---|
| 608 | // get next "combined-particle" index in our own history
|
|---|
| 609 | // making sure we don't go beyond its bounds (if we do
|
|---|
| 610 | // then we're in big trouble anyway...)
|
|---|
| 611 | while (++j < int(_history.size())) {
|
|---|
| 612 | hist_index = unique_hist_order[j];
|
|---|
| 613 | if (hist_index >= _initial_n) break;}
|
|---|
| 614 |
|
|---|
| 615 | // sanity checking -- do not overrun
|
|---|
| 616 | if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in diB matching");
|
|---|
| 617 |
|
|---|
| 618 | // sanity check -- make sure we are taking about the same
|
|---|
| 619 | // jet in reference and new sequences
|
|---|
| 620 | int refjet_index = _history[_history[hist_index].parent1].jetp_index;
|
|---|
| 621 | assert(refjet_index >= 0 && refjet_index < int(_jets.size()));
|
|---|
| 622 | const PseudoJet & refjet = _jets[refjet_index];
|
|---|
| 623 |
|
|---|
| 624 | //cerr << "Inclusive" << endl;
|
|---|
| 625 | //cerr << gs_history[parent1].jetp_index << " " << gs_jets.size() << endl;
|
|---|
| 626 | //cerr << _history[_history[hist_index].parent1].jetp_index << " " << _jets.size() << endl;
|
|---|
| 627 |
|
|---|
| 628 | // If pt disagrees check E; if they both disagree there's a
|
|---|
| 629 | // problem here... NB: a massive particle with zero pt may
|
|---|
| 630 | // have its pt changed when a ghost is added -- this is why we
|
|---|
| 631 | // also require the energy to be wrong before complaining
|
|---|
| 632 | _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
|
|---|
| 633 | ghosted_seq);
|
|---|
| 634 |
|
|---|
| 635 | // set the area at this clustering stage
|
|---|
| 636 | our_areas[hist_index] = area_local;
|
|---|
| 637 | our_area_4vectors[hist_index] = ext_area;
|
|---|
| 638 |
|
|---|
| 639 | // update the parent as well -- that way its area is the area
|
|---|
| 640 | // immediately before clustering (i.e. resolve an ambiguity in
|
|---|
| 641 | // the Cambridge case and ensure in the kt case that the original
|
|---|
| 642 | // particles get a correct area)
|
|---|
| 643 | our_areas[_history[hist_index].parent1] = area_local;
|
|---|
| 644 | our_area_4vectors[_history[hist_index].parent1] = ext_area;
|
|---|
| 645 |
|
|---|
| 646 | }
|
|---|
| 647 | }
|
|---|
| 648 | else if (!ghosted_seq.is_pure_ghost(parent1) &&
|
|---|
| 649 | !ghosted_seq.is_pure_ghost(parent2)) {
|
|---|
| 650 |
|
|---|
| 651 | // get next "combined-particle" index in our own history
|
|---|
| 652 | while (++j < int(_history.size())) {
|
|---|
| 653 | hist_index = unique_hist_order[j];
|
|---|
| 654 | if (hist_index >= _initial_n) break;}
|
|---|
| 655 |
|
|---|
| 656 | // sanity checking -- do not overrun
|
|---|
| 657 | if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in dij matching");
|
|---|
| 658 |
|
|---|
| 659 | // make sure that our reference history entry is also for
|
|---|
| 660 | // an exclusive (dij) clustering (otherwise the comparison jet
|
|---|
| 661 | // will not exist)
|
|---|
| 662 | if (_history[hist_index].parent2 == BeamJet) throw Error("ClusterSequenceActiveArea: could not match clustering sequences (encountered dij matched with diB)");
|
|---|
| 663 |
|
|---|
| 664 | //cerr << "Exclusive: hist_index,hist_size: " << hist_index << " " << _history.size()<< endl;
|
|---|
| 665 | //cerr << gs_hist.jetp_index << " " << gs_jets.size() << endl;
|
|---|
| 666 | //cerr << _history[hist_index].jetp_index << " " << _jets.size() << endl;
|
|---|
| 667 |
|
|---|
| 668 | const PseudoJet & jet = gs_jets[gs_hist.jetp_index];
|
|---|
| 669 | const PseudoJet & refjet = _jets[_history[hist_index].jetp_index];
|
|---|
| 670 |
|
|---|
| 671 | // run sanity check
|
|---|
| 672 | _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
|
|---|
| 673 | ghosted_seq);
|
|---|
| 674 |
|
|---|
| 675 | // update area and our local index (maybe redundant since later
|
|---|
| 676 | // the descendants will reupdate it?)
|
|---|
| 677 | double area_local = ghosted_seq.area(jet);
|
|---|
| 678 | our_areas[hist_index] += area_local;
|
|---|
| 679 |
|
|---|
| 680 | PseudoJet ext_area = ghosted_seq.area_4vector(jet);
|
|---|
| 681 |
|
|---|
| 682 | // GPS TMP debugging (jetclu) -----------------------
|
|---|
| 683 | //ext_area = PseudoJet(1e-100,1e-100,1e-100,4e-100);
|
|---|
| 684 | //our_area_4vectors[hist_index] = ext_area;
|
|---|
| 685 | //cout << "aa "
|
|---|
| 686 | // << our_area_4vectors[hist_index].px() << " "
|
|---|
| 687 | // << our_area_4vectors[hist_index].py() << " "
|
|---|
| 688 | // << our_area_4vectors[hist_index].pz() << " "
|
|---|
| 689 | // << our_area_4vectors[hist_index].E() << endl;
|
|---|
| 690 | //cout << "bb "
|
|---|
| 691 | // << ext_area.px() << " "
|
|---|
| 692 | // << ext_area.py() << " "
|
|---|
| 693 | // << ext_area.pz() << " "
|
|---|
| 694 | // << ext_area.E() << endl;
|
|---|
| 695 | //---------------------------------------------------
|
|---|
| 696 |
|
|---|
| 697 | _jet_def.recombiner()->plus_equal(our_area_4vectors[hist_index], ext_area);
|
|---|
| 698 |
|
|---|
| 699 | // now update areas of parents (so that they becomes areas
|
|---|
| 700 | // immediately before clustering occurred). This is of use
|
|---|
| 701 | // because it allows us to set the areas of the original hard
|
|---|
| 702 | // particles in the kt algorithm; for the Cambridge case it
|
|---|
| 703 | // means a jet's area will be the area just before it clusters
|
|---|
| 704 | // with another hard jet.
|
|---|
| 705 | const PseudoJet & jet1 = gs_jets[gs_history[parent1].jetp_index];
|
|---|
| 706 | int our_parent1 = _history[hist_index].parent1;
|
|---|
| 707 | our_areas[our_parent1] = ghosted_seq.area(jet1);
|
|---|
| 708 | our_area_4vectors[our_parent1] = ghosted_seq.area_4vector(jet1);
|
|---|
| 709 |
|
|---|
| 710 | const PseudoJet & jet2 = gs_jets[gs_history[parent2].jetp_index];
|
|---|
| 711 | int our_parent2 = _history[hist_index].parent2;
|
|---|
| 712 | our_areas[our_parent2] = ghosted_seq.area(jet2);
|
|---|
| 713 | our_area_4vectors[our_parent2] = ghosted_seq.area_4vector(jet2);
|
|---|
| 714 | }
|
|---|
| 715 |
|
|---|
| 716 | }
|
|---|
| 717 |
|
|---|
| 718 | // now add unclustered ghosts to the relevant list so that we can
|
|---|
| 719 | // calculate empty area later.
|
|---|
| 720 | vector<PseudoJet> unclust = ghosted_seq.unclustered_particles();
|
|---|
| 721 | for (unsigned iu = 0; iu < unclust.size(); iu++) {
|
|---|
| 722 | if (ghosted_seq.is_pure_ghost(unclust[iu])) {
|
|---|
| 723 | double area_local = ghosted_seq.area(unclust[iu]);
|
|---|
| 724 | _unclustered_ghosts.push_back(GhostJet(unclust[iu],area_local));
|
|---|
| 725 | }
|
|---|
| 726 | }
|
|---|
| 727 |
|
|---|
| 728 | /*
|
|---|
| 729 | * WARNING:
|
|---|
| 730 | * _average_area has explicitly been sized initially to 2*jets().size()
|
|---|
| 731 | * which can be bigger than our_areas (of size _history.size()
|
|---|
| 732 | * if there are some unclustered particles.
|
|---|
| 733 | * So we must take care about boundaries
|
|---|
| 734 | */
|
|---|
| 735 |
|
|---|
| 736 | for (unsigned int area_index = 0; area_index<our_areas.size(); area_index++){
|
|---|
| 737 | _average_area[area_index] += our_areas[area_index];
|
|---|
| 738 | _average_area2[area_index] += our_areas[area_index]*our_areas[area_index];
|
|---|
| 739 | }
|
|---|
| 740 |
|
|---|
| 741 | //_average_area_4vector += our_area_4vectors;
|
|---|
| 742 | // Use the proper recombination scheme when averaging the area_4vectors
|
|---|
| 743 | // over multiple ghost runs (i.e. the repeat stage);
|
|---|
| 744 | for (unsigned i = 0; i < our_area_4vectors.size(); i++) {
|
|---|
| 745 | _jet_def.recombiner()->plus_equal(_average_area_4vector[i],
|
|---|
| 746 | our_area_4vectors[i]);
|
|---|
| 747 | }
|
|---|
| 748 | }
|
|---|
| 749 |
|
|---|
| 750 |
|
|---|
| 751 | /// check if two jets have the same momentum to within the
|
|---|
| 752 | /// tolerance (and if pt's are not the same we're forgiving and
|
|---|
| 753 | /// look to see if the energy is the same)
|
|---|
| 754 | void ClusterSequenceActiveArea::_throw_unless_jets_have_same_perp_or_E(
|
|---|
| 755 | const PseudoJet & jet,
|
|---|
| 756 | const PseudoJet & refjet,
|
|---|
| 757 | double tolerance,
|
|---|
| 758 | const ClusterSequenceActiveAreaExplicitGhosts & jets_ghosted_seq
|
|---|
| 759 | ) const {
|
|---|
| 760 |
|
|---|
| 761 | if (abs(jet.perp2()-refjet.perp2()) >
|
|---|
| 762 | tolerance*max(jet.perp2(),refjet.perp2())
|
|---|
| 763 | && abs(jet.E()-refjet.E()) > tolerance*max(jet.E(),refjet.E())) {
|
|---|
| 764 | ostringstream ostr;
|
|---|
| 765 | ostr << "Could not match clustering sequence for an inclusive/exclusive jet when reconstructing areas. See FAQ for possible explanations." << endl;
|
|---|
| 766 | ostr << " Ref-Jet: "
|
|---|
| 767 | << refjet.px() << " "
|
|---|
| 768 | << refjet.py() << " "
|
|---|
| 769 | << refjet.pz() << " "
|
|---|
| 770 | << refjet.E() << endl;
|
|---|
| 771 | ostr << " New-Jet: "
|
|---|
| 772 | << jet.px() << " "
|
|---|
| 773 | << jet.py() << " "
|
|---|
| 774 | << jet.pz() << " "
|
|---|
| 775 | << jet.E() << endl;
|
|---|
| 776 | if (jets_ghosted_seq.has_dangerous_particles()) {
|
|---|
| 777 | ostr << " NB: some particles have pt too low wrt ghosts -- this may be the cause" << endl;}
|
|---|
| 778 | //ostr << jet.perp() << " " << refjet.perp() << " "
|
|---|
| 779 | // << jet.perp() - refjet.perp() << endl;
|
|---|
| 780 | throw Error(ostr.str());
|
|---|
| 781 | }
|
|---|
| 782 | }
|
|---|
| 783 |
|
|---|
| 784 | FASTJET_END_NAMESPACE
|
|---|
| 785 |
|
|---|
