FastJet  3.3.3
ClusterSequenceActiveArea.cc
1 //FJSTARTHEADER
2 // $Id: ClusterSequenceActiveArea.cc 4420 2019-11-29 09:28:20Z soyez $
3 //
4 // Copyright (c) 2005-2019, 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. They are described in the original FastJet paper,
16 // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
17 // FastJet as part of work towards a scientific publication, please
18 // quote the version you use and include a citation to the manual and
19 // optionally also to hep-ph/0512210.
20 //
21 // FastJet is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
28 //----------------------------------------------------------------------
29 //FJENDHEADER
30 
31 #include "fastjet/PseudoJet.hh"
32 #include "fastjet/ClusterSequence.hh"
33 #include "fastjet/ClusterSequenceActiveArea.hh"
34 #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_in,
55  const GhostedAreaSpec & ghost_spec,
56  const bool & writeout_combinations) {
57 
58  bool continue_running;
59  _initialise_AA(jet_def_in, ghost_spec, writeout_combinations, continue_running);
60  if (continue_running) {
61  _run_AA(ghost_spec);
62  _postprocess_AA(ghost_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_in,
79  const GhostedAreaSpec & ghost_spec,
80  const bool & writeout_combinations,
81  bool & continue_running)
82 {
83 
84  // store this for future use
85  _ghost_spec_repeat = ghost_spec.repeat();
86 
87  // make sure placeholders are there & zeroed
88  _resize_and_zero_AA();
89 
90  // for future reference...
91  _maxrap_for_area = ghost_spec.ghost_maxrap();
92  _safe_rap_for_area = _maxrap_for_area - jet_def_in.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 (ghost_spec.repeat() <= 0) {
102  _initialise_and_run(jet_def_in, writeout_combinations);
103  continue_running = false;
104  return;
105  }
106 
107  // transfer all relevant info into internal variables
108  _decant_options(jet_def_in, 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 & ghost_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 < ghost_spec.repeat(); irepeat++) {
131 
132  ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets,
133  jet_def(), ghost_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 & ghost_spec) {
153  _average_area /= ghost_spec.repeat();
154  _average_area2 /= ghost_spec.repeat();
155  if (ghost_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 = ghost_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 /= ghost_spec.repeat();
165  _non_jet_area2 /= ghost_spec.repeat();
166  _non_jet_area2 = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/
167  ghost_spec.repeat());
168  _non_jet_number /= ghost_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/ghost_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 & ghost_spec,
184 // const bool & writeout_combinations)
185 // {
186 //
187 // // store this for future use
188 // _ghost_spec_repeat = ghost_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 = ghost_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 (ghost_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 < ghost_spec.repeat(); irepeat++) {
231 //
232 // ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets,
233 // jet_def, ghost_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 /= ghost_spec.repeat();
248 // _average_area2 /= ghost_spec.repeat();
249 // if (ghost_spec.repeat() > 1) {
250 // _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/
251 // (ghost_spec.repeat()-1));
252 // } else {
253 // _average_area2 = 0.0;
254 // }
255 //
256 // _non_jet_area /= ghost_spec.repeat();
257 // _non_jet_area2 /= ghost_spec.repeat();
258 // _non_jet_area2 = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/
259 // ghost_spec.repeat());
260 // _non_jet_number /= ghost_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/ghost_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 Selector & selector) const {
446  // make sure that the selector applies jet by jet
447  if (! selector.applies_jet_by_jet()){
448  throw Error("ClusterSequenceActiveArea: empty area can only be computed from selectors applying jet by jet");
449  }
450 
451  double empty = 0.0;
452  // first deal with ghost jets
453  for (unsigned i = 0; i < _ghost_jets.size(); i++) {
454  if (selector.pass(_ghost_jets[i])) {
455  empty += _ghost_jets[i].area;
456  }
457  }
458  // then deal with unclustered ghosts
459  for (unsigned i = 0; i < _unclustered_ghosts.size(); i++) {
460  if (selector.pass(_unclustered_ghosts[i])) {
461  empty += _unclustered_ghosts[i].area;
462  }
463  }
464  empty /= _ghost_spec_repeat;
465  return empty;
466 }
467 
468 //----------------------------------------------------------------------
469 double ClusterSequenceActiveArea::n_empty_jets(const Selector & selector) const {
470  _check_selector_good_for_median(selector);
471 
472  double inrange = 0;
473  for (unsigned i = 0; i < _ghost_jets.size(); i++) {
474  if (selector.pass(_ghost_jets[i])) inrange++;
475  }
476  inrange /= _ghost_spec_repeat;
477  return inrange;
478 }
479 
480 //----------------------------------------------------------------------
481 /// transfer the history (and jet-momenta) from clust_seq to our
482 /// own internal structure while removing ghosts
483 void ClusterSequenceActiveArea::_transfer_ghost_free_history(
484  const ClusterSequenceActiveAreaExplicitGhosts & ghosted_seq) {
485 
486  const vector<history_element> & gs_history = ghosted_seq.history();
487  vector<int> gs2self_hist_map(gs_history.size());
488 
489  // first transfer info about strategy used (which isn't necessarily
490  // always the one that got asked for...)
491  _strategy = ghosted_seq.strategy_used();
492 
493  // work our way through to first non-trivial combination
494  unsigned igs = 0;
495  unsigned iself = 0;
496  while (igs < gs_history.size() && gs_history[igs].parent1 == InexistentParent) {
497  // record correspondence
498  if (!ghosted_seq.is_pure_ghost(igs)) {
499  gs2self_hist_map[igs] = iself++;
500  } else {
501  gs2self_hist_map[igs] = Invalid;
502  }
503  igs++;
504  };
505 
506  // make sure the count of non-ghost initial jets is equal to
507  // what we already have in terms of initial jets
508  assert(iself == _history.size());
509 
510  // if there was no clustering in this event (e.g. SISCone passive
511  // area with zero input particles, or with a pt cut on stable cones
512  // that kills all jets), then don't bother with the rest (which
513  // would crash!)
514  if (igs == gs_history.size()) return;
515 
516  // now actually transfer things
517  do {
518  // if we are a pure ghost, then go on to next round
519  if (ghosted_seq.is_pure_ghost(igs)) {
520  gs2self_hist_map[igs] = Invalid;
521  continue;
522  }
523 
524  const history_element & gs_hist_el = gs_history[igs];
525 
526  bool parent1_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent1);
527  bool parent2_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent2);
528 
529  // if exactly one parent is a ghost then maintain info about the
530  // non-ghost correspondence for this jet, and then go on to next
531  // recombination in the ghosted sequence
532  if (parent1_is_ghost && !parent2_is_ghost && gs_hist_el.parent2 >= 0) {
533  gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent2];
534  continue;
535  }
536  if (!parent1_is_ghost && parent2_is_ghost) {
537  gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent1];
538  continue;
539  }
540 
541  // no parents are ghosts...
542  if (gs_hist_el.parent2 >= 0) {
543  // recombination of two non-ghosts
544  gs2self_hist_map[igs] = _history.size();
545  // record the recombination in our own sequence
546  int newjet_k; // dummy var -- not used
547  int jet_i = _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index;
548  int jet_j = _history[gs2self_hist_map[gs_hist_el.parent2]].jetp_index;
549  //cerr << "recombining "<< jet_i << " and "<< jet_j << endl;
550  _do_ij_recombination_step(jet_i, jet_j, gs_hist_el.dij, newjet_k);
551  } else {
552  // we have a non-ghost that has become a beam-jet
553  assert(gs_history[igs].parent2 == BeamJet);
554  // record position
555  gs2self_hist_map[igs] = _history.size();
556  // record the recombination in our own sequence
557  _do_iB_recombination_step(
558  _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index,
559  gs_hist_el.dij);
560  }
561  } while (++igs < gs_history.size());
562 
563 }
564 
565 //----------------------------------------------------------------------
566 void ClusterSequenceActiveArea::_transfer_areas(
567  const vector<int> & unique_hist_order,
568  const ClusterSequenceActiveAreaExplicitGhosts & ghosted_seq ) {
569 
570  const vector<history_element> & gs_history = ghosted_seq.history();
571  const vector<PseudoJet> & gs_jets = ghosted_seq.jets();
572  vector<int> gs_unique_hist_order = ghosted_seq.unique_history_order();
573 
574  const double tolerance = 1e-11; // to decide when two jets are the same
575 
576  int j = -1;
577  int hist_index = -1;
578 
579  valarray<double> our_areas(_history.size());
580  our_areas = 0.0;
581 
582  valarray<PseudoJet> our_area_4vectors(_history.size());
583  our_area_4vectors = PseudoJet(0.0,0.0,0.0,0.0);
584 
585  for (unsigned i = 0; i < gs_history.size(); i++) {
586  // only consider composite particles
587  unsigned gs_hist_index = gs_unique_hist_order[i];
588  if (gs_hist_index < ghosted_seq.n_particles()) continue;
589  const history_element & gs_hist = gs_history[gs_unique_hist_order[i]];
590  int parent1 = gs_hist.parent1;
591  int parent2 = gs_hist.parent2;
592 
593  if (parent2 == BeamJet) {
594  // need to look at parent to get the actual jet
595  const PseudoJet & jet =
596  gs_jets[gs_history[parent1].jetp_index];
597  double area_local = ghosted_seq.area(jet);
598  PseudoJet ext_area = ghosted_seq.area_4vector(jet);
599 
600  if (ghosted_seq.is_pure_ghost(parent1)) {
601  // record the existence of the pure ghost jet for future use
602  _ghost_jets.push_back(GhostJet(jet,area_local));
603  if (abs(jet.rap()) < _safe_rap_for_area) {
604  _non_jet_area += area_local;
605  _non_jet_area2 += area_local*area_local;
606  _non_jet_number += 1;
607  }
608  } else {
609 
610  // get next "combined-particle" index in our own history
611  // making sure we don't go beyond its bounds (if we do
612  // then we're in big trouble anyway...)
613  while (++j < int(_history.size())) {
614  hist_index = unique_hist_order[j];
615  if (hist_index >= _initial_n) break;}
616 
617  // sanity checking -- do not overrun
618  if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in diB matching");
619 
620  // sanity check -- make sure we are taking about the same
621  // jet in reference and new sequences
622  int refjet_index = _history[_history[hist_index].parent1].jetp_index;
623  assert(refjet_index >= 0 && refjet_index < int(_jets.size()));
624  const PseudoJet & refjet = _jets[refjet_index];
625 
626  //cerr << "Inclusive" << endl;
627  //cerr << gs_history[parent1].jetp_index << " " << gs_jets.size() << endl;
628  //cerr << _history[_history[hist_index].parent1].jetp_index << " " << _jets.size() << endl;
629 
630  // If pt disagrees check E; if they both disagree there's a
631  // problem here... NB: a massive particle with zero pt may
632  // have its pt changed when a ghost is added -- this is why we
633  // also require the energy to be wrong before complaining
634  _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
635  ghosted_seq);
636 
637  // set the area at this clustering stage
638  our_areas[hist_index] = area_local;
639  our_area_4vectors[hist_index] = ext_area;
640 
641  // update the parent as well -- that way its area is the area
642  // immediately before clustering (i.e. resolve an ambiguity in
643  // the Cambridge case and ensure in the kt case that the original
644  // particles get a correct area)
645  our_areas[_history[hist_index].parent1] = area_local;
646  our_area_4vectors[_history[hist_index].parent1] = ext_area;
647 
648  }
649  }
650  else if (!ghosted_seq.is_pure_ghost(parent1) &&
651  !ghosted_seq.is_pure_ghost(parent2)) {
652 
653  // get next "combined-particle" index in our own history
654  while (++j < int(_history.size())) {
655  hist_index = unique_hist_order[j];
656  if (hist_index >= _initial_n) break;}
657 
658  // sanity checking -- do not overrun
659  if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in dij matching");
660 
661  // make sure that our reference history entry is also for
662  // an exclusive (dij) clustering (otherwise the comparison jet
663  // will not exist)
664  if (_history[hist_index].parent2 == BeamJet) throw Error("ClusterSequenceActiveArea: could not match clustering sequences (encountered dij matched with diB)");
665 
666  //cerr << "Exclusive: hist_index,hist_size: " << hist_index << " " << _history.size()<< endl;
667  //cerr << gs_hist.jetp_index << " " << gs_jets.size() << endl;
668  //cerr << _history[hist_index].jetp_index << " " << _jets.size() << endl;
669 
670  const PseudoJet & jet = gs_jets[gs_hist.jetp_index];
671  const PseudoJet & refjet = _jets[_history[hist_index].jetp_index];
672 
673  // run sanity check
674  _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
675  ghosted_seq);
676 
677  // update area and our local index (maybe redundant since later
678  // the descendants will reupdate it?)
679  double area_local = ghosted_seq.area(jet);
680  our_areas[hist_index] += area_local;
681 
682  PseudoJet ext_area = ghosted_seq.area_4vector(jet);
683 
684  // GPS TMP debugging (jetclu) -----------------------
685  //ext_area = PseudoJet(1e-100,1e-100,1e-100,4e-100);
686  //our_area_4vectors[hist_index] = ext_area;
687  //cout << "aa "
688  // << our_area_4vectors[hist_index].px() << " "
689  // << our_area_4vectors[hist_index].py() << " "
690  // << our_area_4vectors[hist_index].pz() << " "
691  // << our_area_4vectors[hist_index].E() << endl;
692  //cout << "bb "
693  // << ext_area.px() << " "
694  // << ext_area.py() << " "
695  // << ext_area.pz() << " "
696  // << ext_area.E() << endl;
697  //---------------------------------------------------
698 
699  _jet_def.recombiner()->plus_equal(our_area_4vectors[hist_index], ext_area);
700 
701  // now update areas of parents (so that they becomes areas
702  // immediately before clustering occurred). This is of use
703  // because it allows us to set the areas of the original hard
704  // particles in the kt algorithm; for the Cambridge case it
705  // means a jet's area will be the area just before it clusters
706  // with another hard jet.
707  const PseudoJet & jet1 = gs_jets[gs_history[parent1].jetp_index];
708  int our_parent1 = _history[hist_index].parent1;
709  our_areas[our_parent1] = ghosted_seq.area(jet1);
710  our_area_4vectors[our_parent1] = ghosted_seq.area_4vector(jet1);
711 
712  const PseudoJet & jet2 = gs_jets[gs_history[parent2].jetp_index];
713  int our_parent2 = _history[hist_index].parent2;
714  our_areas[our_parent2] = ghosted_seq.area(jet2);
715  our_area_4vectors[our_parent2] = ghosted_seq.area_4vector(jet2);
716  }
717 
718  }
719 
720  // now add unclustered ghosts to the relevant list so that we can
721  // calculate empty area later.
722  vector<PseudoJet> unclust = ghosted_seq.unclustered_particles();
723  for (unsigned iu = 0; iu < unclust.size(); iu++) {
724  if (ghosted_seq.is_pure_ghost(unclust[iu])) {
725  double area_local = ghosted_seq.area(unclust[iu]);
726  _unclustered_ghosts.push_back(GhostJet(unclust[iu],area_local));
727  }
728  }
729 
730  /*
731  * WARNING:
732  * _average_area has explicitly been sized initially to 2*jets().size()
733  * which can be bigger than our_areas (of size _history.size()
734  * if there are some unclustered particles.
735  * So we must take care about boundaries
736  */
737 
738  for (unsigned int area_index = 0; area_index<our_areas.size(); area_index++){
739  _average_area[area_index] += our_areas[area_index];
740  _average_area2[area_index] += our_areas[area_index]*our_areas[area_index];
741  }
742 
743  //_average_area_4vector += our_area_4vectors;
744  // Use the proper recombination scheme when averaging the area_4vectors
745  // over multiple ghost runs (i.e. the repeat stage);
746  for (unsigned i = 0; i < our_area_4vectors.size(); i++) {
747  _jet_def.recombiner()->plus_equal(_average_area_4vector[i],
748  our_area_4vectors[i]);
749  }
750 }
751 
752 
753 /// check if two jets have the same momentum to within the
754 /// tolerance (and if pt's are not the same we're forgiving and
755 /// look to see if the energy is the same)
756 void ClusterSequenceActiveArea::_throw_unless_jets_have_same_perp_or_E(
757  const PseudoJet & jet,
758  const PseudoJet & refjet,
759  double tolerance,
760  const ClusterSequenceActiveAreaExplicitGhosts & jets_ghosted_seq
761 ) const {
762 
763  if (abs(jet.perp2()-refjet.perp2()) >
764  tolerance*max(jet.perp2(),refjet.perp2())
765  && abs(jet.E()-refjet.E()) > tolerance*max(jet.E(),refjet.E())) {
766  ostringstream ostr;
767  ostr << "Could not match clustering sequence for an inclusive/exclusive jet when reconstructing areas. See FAQ for possible explanations." << endl;
768  ostr << " Ref-Jet: "
769  << refjet.px() << " "
770  << refjet.py() << " "
771  << refjet.pz() << " "
772  << refjet.E() << endl;
773  ostr << " New-Jet: "
774  << jet.px() << " "
775  << jet.py() << " "
776  << jet.pz() << " "
777  << jet.E() << endl;
778  if (jets_ghosted_seq.has_dangerous_particles()) {
779  ostr << " NB: some particles have pt too low wrt ghosts -- this may be the cause" << endl;}
780  //ostr << jet.perp() << " " << refjet.perp() << " "
781  // << jet.perp() - refjet.perp() << endl;
782  throw Error(ostr.str());
783  }
784 }
785 
786 FASTJET_END_NAMESPACE
787 
virtual bool is_pure_ghost(const PseudoJet &jet) const override
true if a jet is made exclusively of ghosts
std::vector< int > unique_history_order() const
routine that returns an order in which to read the history such that clusterings that lead to identic...
mean_pt_strategies
enum providing a variety of tentative strategies for estimating the background (e.g.
const std::vector< history_element > & history() const
allow the user to access the raw internal history.
int jetp_index
index in _history where the current jet is recombined with another jet to form its child...
virtual double area(const PseudoJet &jet) const override
returns the area of a jet
int parent2
index in _history where first parent of this jet was created (InexistentParent if this jet is an orig...
bool applies_jet_by_jet() const
returns true if this can be applied jet by jet
Definition: Selector.hh:215
double dij
index in the _jets vector where we will find the
bool pass(const PseudoJet &jet) const
return true if the jet passes the selection
Definition: Selector.hh:178
Strategy strategy_used() const
return the enum value of the strategy used to cluster the event
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:47
Like ClusterSequence with computation of the active jet area with the addition of explicit ghosts...
double perp2() const
returns the squared transverse momentum
Definition: PseudoJet.hh:143
bool has_dangerous_particles() const
returns true if there are any particles whose transverse momentum if so low that there&#39;s a risk of th...
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:125
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition: Selector.hh:149
a single element in the clustering history
virtual PseudoJet area_4vector(const PseudoJet &jet) const override
returns a four vector corresponding to the sum (E-scheme) of the ghost four-vectors composing the jet...
Parameters to configure the computation of jet areas using ghosts.
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
class that is intended to hold a full definition of the jet clusterer
unsigned int n_particles() const
returns the number of particles that were provided to the clustering algorithm (helps the user find t...
std::vector< PseudoJet > unclustered_particles() const
return the set of particles that have not been clustered.