29 #include "fastjet/PseudoJet.hh"
30 #include "fastjet/ClusterSequence.hh"
31 #include "fastjet/ClusterSequenceActiveArea.hh"
32 #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
40 FASTJET_BEGIN_NAMESPACE
51 void ClusterSequenceActiveArea::_initialise_and_run_AA (
54 const bool & writeout_combinations) {
56 bool continue_running;
57 _initialise_AA(jet_def_in, ghost_spec, writeout_combinations, continue_running);
58 if (continue_running) {
60 _postprocess_AA(ghost_spec);
65 void ClusterSequenceActiveArea::_resize_and_zero_AA () {
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;
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)
83 _ghost_spec_repeat = ghost_spec.repeat();
86 _resize_and_zero_AA();
89 _maxrap_for_area = ghost_spec.ghost_maxrap();
90 _safe_rap_for_area = _maxrap_for_area - jet_def_in.R();
99 if (ghost_spec.repeat() <= 0) {
100 _initialise_and_run(jet_def_in, writeout_combinations);
101 continue_running =
false;
106 _decant_options(jet_def_in, writeout_combinations);
110 _fill_initial_history();
113 _has_dangerous_particles =
false;
115 continue_running =
true;
120 void ClusterSequenceActiveArea::_run_AA (
const GhostedAreaSpec & ghost_spec) {
122 vector<PseudoJet> input_jets(_jets);
125 vector<int> unique_tree;
128 for (
int irepeat = 0; irepeat < ghost_spec.repeat(); irepeat++) {
130 ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets,
131 jet_def(), ghost_spec);
133 _has_dangerous_particles |= clust_seq.has_dangerous_particles();
137 _transfer_ghost_free_history(clust_seq);
139 unique_tree = unique_history_order();
143 _transfer_areas(unique_tree, clust_seq);
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) {
156 const double tmp = ghost_spec.repeat()-1;
157 _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/tmp);
159 _average_area2 = 0.0;
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();
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];
274 double ClusterSequenceActiveArea::pt_per_unit_area(
277 vector<PseudoJet> incl_jets = inclusive_jets();
278 vector<double> pt_over_areas;
280 for (
unsigned i = 0; i < incl_jets.size(); i++) {
281 if (abs(incl_jets[i].rap()) < _safe_rap_for_area) {
283 if ( strat == median_4vector ) {
284 this_area = area_4vector(incl_jets[i]).perp();
286 this_area = area(incl_jets[i]);
288 pt_over_areas.push_back(incl_jets[i].perp()/this_area);
293 if (pt_over_areas.size() == 0) {
return 0.0;}
298 sort(pt_over_areas.begin(), pt_over_areas.end());
299 double non_ghost_median_ratio = pt_over_areas[pt_over_areas.size()/2];
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);
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);
312 nj_median_ratio = 0.0;
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) {
324 if ( strat == median_4vector ) {
325 this_area = area_4vector(incl_jets[i]).perp();
327 this_area = area(incl_jets[i]);
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++;
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];
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;
352 cout <<
"-----------------------------------"<<endl;
353 for (
unsigned i = 0; i <= pt_over_areas.size()/2 ; i++ ) {
354 cout <<
"Median "<< i <<
" = " << pt_over_areas[i]<<endl;
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;
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;
375 return ratio_sum/ratio_n;
377 return nj_median_ratio;
443 double ClusterSequenceActiveArea::empty_area(
const Selector & selector)
const {
446 throw Error(
"ClusterSequenceActiveArea: empty area can only be computed from selectors applying jet by jet");
451 for (
unsigned i = 0; i < _ghost_jets.size(); i++) {
452 if (selector.
pass(_ghost_jets[i])) {
453 empty += _ghost_jets[i].area;
457 for (
unsigned i = 0; i < _unclustered_ghosts.size(); i++) {
458 if (selector.
pass(_unclustered_ghosts[i])) {
459 empty += _unclustered_ghosts[i].area;
462 empty /= _ghost_spec_repeat;
467 double ClusterSequenceActiveArea::n_empty_jets(
const Selector & selector)
const {
468 _check_selector_good_for_median(selector);
471 for (
unsigned i = 0; i < _ghost_jets.size(); i++) {
472 if (selector.
pass(_ghost_jets[i])) inrange++;
474 inrange /= _ghost_spec_repeat;
481 void ClusterSequenceActiveArea::_transfer_ghost_free_history(
484 const vector<history_element> & gs_history = ghosted_seq.
history();
485 vector<int> gs2self_hist_map(gs_history.size());
494 while (igs < gs_history.size() && gs_history[igs].parent1 == InexistentParent) {
497 gs2self_hist_map[igs] = iself++;
499 gs2self_hist_map[igs] = Invalid;
506 assert(iself == _history.size());
512 if (igs == gs_history.size())
return;
518 gs2self_hist_map[igs] = Invalid;
524 bool parent1_is_ghost = ghosted_seq.
is_pure_ghost(gs_hist_el.parent1);
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];
534 if (!parent1_is_ghost && parent2_is_ghost) {
535 gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent1];
542 gs2self_hist_map[igs] = _history.size();
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;
548 _do_ij_recombination_step(jet_i, jet_j, gs_hist_el.
dij, newjet_k);
551 assert(gs_history[igs].parent2 == BeamJet);
553 gs2self_hist_map[igs] = _history.size();
555 _do_iB_recombination_step(
556 _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index,
559 }
while (++igs < gs_history.size());
564 void ClusterSequenceActiveArea::_transfer_areas(
565 const vector<int> & unique_hist_order,
568 const vector<history_element> & gs_history = ghosted_seq.
history();
569 const vector<PseudoJet> & gs_jets = ghosted_seq.
jets();
572 const double tolerance = 1e-11;
577 valarray<double> our_areas(_history.size());
580 valarray<PseudoJet> our_area_4vectors(_history.size());
581 our_area_4vectors =
PseudoJet(0.0,0.0,0.0,0.0);
583 for (
unsigned i = 0; i < gs_history.size(); i++) {
585 unsigned gs_hist_index = gs_unique_hist_order[i];
586 if (gs_hist_index < ghosted_seq.
n_particles())
continue;
588 int parent1 = gs_hist.parent1;
591 if (parent2 == BeamJet) {
594 gs_jets[gs_history[parent1].jetp_index];
595 double area_local = ghosted_seq.
area(jet);
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;
611 while (++j <
int(_history.size())) {
612 hist_index = unique_hist_order[j];
613 if (hist_index >= _initial_n)
break;}
616 if (j >=
int(_history.size()))
throw Error(
"ClusterSequenceActiveArea: overran reference array in diB matching");
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];
632 _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
636 our_areas[hist_index] = area_local;
637 our_area_4vectors[hist_index] = ext_area;
643 our_areas[_history[hist_index].parent1] = area_local;
644 our_area_4vectors[_history[hist_index].parent1] = ext_area;
652 while (++j <
int(_history.size())) {
653 hist_index = unique_hist_order[j];
654 if (hist_index >= _initial_n)
break;}
657 if (j >=
int(_history.size()))
throw Error(
"ClusterSequenceActiveArea: overran reference array in dij matching");
662 if (_history[hist_index].parent2 == BeamJet)
throw Error(
"ClusterSequenceActiveArea: could not match clustering sequences (encountered dij matched with diB)");
669 const PseudoJet & refjet = _jets[_history[hist_index].jetp_index];
672 _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
677 double area_local = ghosted_seq.
area(jet);
678 our_areas[hist_index] += area_local;
697 _jet_def.recombiner()->plus_equal(our_area_4vectors[hist_index], ext_area);
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);
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);
721 for (
unsigned iu = 0; iu < unclust.size(); iu++) {
723 double area_local = ghosted_seq.
area(unclust[iu]);
724 _unclustered_ghosts.push_back(GhostJet(unclust[iu],area_local));
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];
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]);
754 void ClusterSequenceActiveArea::_throw_unless_jets_have_same_perp_or_E(
763 && abs(jet.E()-refjet.E()) > tolerance*max(jet.E(),refjet.E())) {
765 ostr <<
"Could not match clustering sequence for an inclusive/exclusive jet when reconstructing areas. See FAQ for possible explanations." << endl;
767 << refjet.px() <<
" "
768 << refjet.py() <<
" "
769 << refjet.pz() <<
" "
770 << refjet.E() << endl;
777 ostr <<
" NB: some particles have pt too low wrt ghosts -- this may be the cause" << endl;}
780 throw Error(ostr.str());
784 FASTJET_END_NAMESPACE