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