FastJet  3.0.6
ClusterSequenceActiveArea.cc
1 //STARTHEADER
2 // $Id: ClusterSequenceActiveArea.cc 2687 2011-11-14 11:17:51Z soyez $
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 
fastjet::ClusterSequenceActiveAreaExplicitGhosts
Definition: ClusterSequenceActiveAreaExplicitGhosts.hh:55
fastjet::Selector::applies_jet_by_jet
bool applies_jet_by_jet() const
returns true if this can be applied jet by jet
Definition: Selector.hh:196
fastjet::JetDefinition
Definition: JetDefinition.hh:173
fastjet::ClusterSequenceActiveAreaExplicitGhosts::area
virtual double area(const PseudoJet &jet) const
returns the area of a jet
Definition: ClusterSequenceActiveAreaExplicitGhosts.cc:64
fastjet::ClusterSequence::strategy_used
Strategy strategy_used() const
return the enum value of the strategy used to cluster the event
Definition: ClusterSequence.hh:281
fastjet::ClusterSequence::history_element::parent2
int parent2
index in _history where first parent of this jet was created (InexistentParent if this jet is an orig...
Definition: ClusterSequence.hh:409
fastjet::ClusterSequence::history_element::jetp_index
int jetp_index
index in _history where the current jet is recombined with another jet to form its child.
Definition: ClusterSequence.hh:420
fastjet::ClusterSequence::history
const std::vector< history_element > & history() const
allow the user to access the raw internal history.
Definition: ClusterSequence.hh:963
fastjet::ClusterSequence::history_element
Definition: ClusterSequence.hh:404
fastjet::ClusterSequenceActiveAreaExplicitGhosts::is_pure_ghost
virtual bool is_pure_ghost(const PseudoJet &jet) const
true if a jet is made exclusively of ghosts
Definition: ClusterSequenceActiveAreaExplicitGhosts.cc:83
fastjet::ClusterSequenceActiveAreaExplicitGhosts::area_4vector
virtual PseudoJet area_4vector(const PseudoJet &jet) const
returns a four vector corresponding to the sum (E-scheme) of the ghost four-vectors composing the jet...
Definition: ClusterSequenceActiveAreaExplicitGhosts.cc:78
fastjet::ClusterSequence::jets
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
Definition: ClusterSequence.hh:959
fastjet::GhostedAreaSpec
Definition: GhostedAreaSpec.hh:62
fastjet::PseudoJet
Definition: PseudoJet.hh:65
fastjet::ClusterSequence::n_particles
unsigned int n_particles() const
returns the number of particles that were provided to the clustering algorithm (helps the user find t...
Definition: ClusterSequence.hh:967
fastjet::ClusterSequence::history_element::dij
double dij
index in the _jets vector where we will find the
Definition: ClusterSequence.hh:427
fastjet::PseudoJet::rap
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:121
fastjet::ClusterSequenceActiveAreaExplicitGhosts::has_dangerous_particles
bool has_dangerous_particles() const
returns true if there are any particles whose transverse momentum if so low that there's a risk of th...
Definition: ClusterSequenceActiveAreaExplicitGhosts.hh:132
fastjet::Selector
Definition: Selector.hh:141
fastjet::ClusterSequence::unique_history_order
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...
Definition: ClusterSequence.cc:1303
fastjet::ClusterSequenceActiveArea::mean_pt_strategies
mean_pt_strategies
enum providing a variety of tentative strategies for estimating the background (e....
Definition: ClusterSequenceActiveArea.hh:85
fastjet::Selector::pass
bool pass(const PseudoJet &jet) const
return true if the jet passes the selection
Definition: Selector.hh:171
fastjet::PseudoJet::perp2
double perp2() const
returns the squared transverse momentum
Definition: PseudoJet.hh:139
fastjet::ClusterSequence::unclustered_particles
std::vector< PseudoJet > unclustered_particles() const
return the set of particles that have not been clustered.
Definition: ClusterSequence.cc:1358
fastjet::Error
Definition: Error.hh:41