Fork me on GitHub

source: git/external/fastjet/ClusterSequenceActiveArea.cc@ ded7435

ImprovedOutputFile Timing dual_readout llp
Last change on this file since ded7435 was d7d2da3, checked in by pavel <pavel@…>, 12 years ago

move branches/ModularDelphes to trunk

  • Property mode set to 100644
File size: 29.0 KB
Line 
1//STARTHEADER
2// $Id$
3//
4// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
5//
6//----------------------------------------------------------------------
7// This file is part of FastJet.
8//
9// FastJet is free software; you can redistribute it and/or modify
10// it under the terms of the GNU General Public License as published by
11// the Free Software Foundation; either version 2 of the License, or
12// (at your option) any later version.
13//
14// The algorithms that underlie FastJet have required considerable
15// development and are described in hep-ph/0512210. If you use
16// FastJet as part of work towards a scientific publication, please
17// include a citation to the FastJet paper.
18//
19// FastJet is distributed in the hope that it will be useful,
20// but WITHOUT ANY WARRANTY; without even the implied warranty of
21// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22// GNU General Public License for more details.
23//
24// You should have received a copy of the GNU General Public License
25// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
26//----------------------------------------------------------------------
27//ENDHEADER
28
29#include "fastjet/PseudoJet.hh"
30#include "fastjet/ClusterSequence.hh"
31#include "fastjet/ClusterSequenceActiveArea.hh"
32#include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
33#include<iostream>
34#include<vector>
35#include<sstream>
36#include<algorithm>
37#include<cmath>
38#include<valarray>
39
40FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
41
42
43using 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
51void 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//----------------------------------------------------------------------
65void 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-------------------------------------
75void 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//----------------------------------------------------------------------
120void 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)
150void 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//----------------------------------------------------------------------
274double 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//----------------------------------------------------------------------
443double 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//----------------------------------------------------------------------
467double 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
481void 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//----------------------------------------------------------------------
564void 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)
754void 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
784FASTJET_END_NAMESPACE
785
Note: See TracBrowser for help on using the repository browser.