Fork me on GitHub

source: git/external/fastjet/ClusterSequenceActiveArea.cc@ 0e6b14e

ImprovedOutputFile Timing
Last change on this file since 0e6b14e was b7b836a, checked in by Pavel Demin <pavel-demin@…>, 7 years ago

update FastJet library to 3.3.1 and FastJet Contrib library to 1.036

  • Property mode set to 100644
File size: 29.2 KB
Line 
1//FJSTARTHEADER
2// $Id: ClusterSequenceActiveArea.cc 4354 2018-04-22 07:12:37Z salam $
3//
4// Copyright (c) 2005-2018, 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
42FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
43
44
45using 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
53void 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//----------------------------------------------------------------------
67void 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-------------------------------------
77void 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//----------------------------------------------------------------------
122void 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)
152void 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//----------------------------------------------------------------------
276double 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//----------------------------------------------------------------------
445double 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//----------------------------------------------------------------------
469double 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
483void 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//----------------------------------------------------------------------
566void 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)
756void 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
786FASTJET_END_NAMESPACE
787
Note: See TracBrowser for help on using the repository browser.