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
double rap() const
returns the rapidity or some large value when the rapidity is infinite
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
bool applies_jet_by_jet() const
returns true if this can be applied jet by jet
mean_pt_strategies
enum providing a variety of tentative strategies for estimating the background (e.g.
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...
int jetp_index
index in _history where the current jet is recombined with another jet to form its child...
virtual bool is_pure_ghost(const PseudoJet &jet) const
true if a jet is made exclusively of ghosts
int parent2
index in _history where first parent of this jet was created (InexistentParent if this jet is an orig...
std::vector< PseudoJet > unclustered_particles() const
return the set of particles that have not been clustered.
double dij
index in the _jets vector where we will find the
unsigned int n_particles() const
returns the number of particles that were provided to the clustering algorithm (helps the user find t...
const std::vector< history_element > & history() const
allow the user to access the raw internal history.
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...
base class corresponding to errors that can be thrown by FastJet
Like ClusterSequence with computation of the active jet area with the addition of explicit ghosts...
virtual double area(const PseudoJet &jet) const
returns the area of a jet
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
a single element in the clustering history
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...
Strategy strategy_used() const
return the enum value of the strategy used to cluster the event
Parameters to configure the computation of jet areas using ghosts.
bool pass(const PseudoJet &jet) const
return true if the jet passes the selection
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
double perp2() const
returns the squared transverse momentum
class that is intended to hold a full definition of the jet clusterer