fastjet::ClusterSequenceActiveArea Class Reference

Class that behaves essentially like ClusterSequence except that it also provides access to the area of a jet (which will be a random quantity. More...

#include <ClusterSequenceActiveArea.hh>

Inheritance diagram for fastjet::ClusterSequenceActiveArea:

[legend]Collaboration diagram for fastjet::ClusterSequenceActiveArea:

[legend]List of all members.

Public Types
enum	mean_pt_strategies {   median = 0, non_ghost_median, pttot_over_areatot, pttot_over_areatot_cut,   mean_ratio_cut, play, median_4vector }
	enum providing a variety of tentative strategies for estimating the background (e.g. More...
Public Member Functions
template<class L>
	ClusterSequenceActiveArea (const std::vector< L > &pseudojets, const JetDefinition &jet_def, const ActiveAreaSpec &area_spec, const bool &writeout_combinations=false)
	constructor based on JetDefinition and ActiveAreaSpec
virtual double	area (const PseudoJet &jet) const
	return the area associated with the given jet; this base class returns 0.
virtual double	area_error (const PseudoJet &jet) const
	return the error (uncertainty) associated with the determination of the area of this jet; this base class returns 0.
virtual PseudoJet	area_4vector (const PseudoJet &jet) const
	return a PseudoJet whose 4-vector is defined by the following integral
double	pt_per_unit_area (mean_pt_strategies strat=median, double range=2.0) const
	return the transverse momentum per unit area according to one of the above strategies; for some strategies (those with "cut" in their name) the parameter "range" allows one to exclude a subset of the jets for the background estimation, those that have pt/area > median(pt/area)*range.
void	parabolic_pt_per_unit_area (double &a, double &b, double raprange=-1.0, double exclude_above=-1.0, bool use_area_4vector=false) const
	fits a form pt_per_unit_area(y) = a + b*y^2 in the range abs(y)<raprange (for negative raprange, it defaults to _safe_rap_for_area).
Private Member Functions
void	_initialise_and_run_AA (const JetDefinition &jet_def, const ActiveAreaSpec &area_spec, const bool &writeout_combinations=false)
	does the initialisation and running specific to the active areas class
void	_transfer_ghost_free_history (const ClusterSequenceActiveAreaExplicitGhosts &clust_seq)
	transfer the history (and jet-momenta) from clust_seq to our own internal structure while removing ghosts
void	_transfer_areas (const vector< int > &unique_hist_order, const ClusterSequenceActiveAreaExplicitGhosts &)
	transfer areas from the ClusterSequenceActiveAreaExplicitGhosts object into our internal area bookkeeping.
void	_extract_tree (vector< int > &) const
	routine for extracting the tree in an order that will be independent of any degeneracies in the recombination sequence that don't affect the composition of the final jets
void	_extract_tree_children (int pos, valarray< bool > &, const valarray< int > &, vector< int > &) const
	do the part of the extraction associated with pos, working through its children and their parents
void	_extract_tree_parents (int pos, valarray< bool > &, const valarray< int > &, vector< int > &) const
	do the part of the extraction associated with the parents of pos.
Private Attributes
valarray< double >	_average_area
valarray< double >	_average_area2
valarray< PseudoJet >	_average_area_4vector
double	_non_jet_area
double	_non_jet_area2
double	_non_jet_number
double	_maxrap_for_area
double	_safe_rap_for_area

---

Detailed Description

Class that behaves essentially like ClusterSequence except that it also provides access to the area of a jet (which will be a random quantity.

.. Figure out what to do about seeds later...)

Definition at line 49 of file ClusterSequenceActiveArea.hh.

---

Member Enumeration Documentation

enum fastjet::ClusterSequenceActiveArea::mean_pt_strategies

enum providing a variety of tentative strategies for estimating the background (e.g. non-jet) activity in a highly populated event; the one that has been most extensively tested is median. Enumeration values: median  non_ghost_median  pttot_over_areatot  pttot_over_areatot_cut  mean_ratio_cut  play  median_4vector  Definition at line 70 of file ClusterSequenceActiveArea.hh. {median=0, non_ghost_median, pttot_over_areatot, pttot_over_areatot_cut, mean_ratio_cut, play, median_4vector};

---

Constructor & Destructor Documentation

template<class L> fastjet::ClusterSequenceActiveArea::ClusterSequenceActiveArea (  const std::vector< L > &  pseudojets, const JetDefinition &  jet_def, const ActiveAreaSpec &  area_spec, const bool &  writeout_combinations = false )

constructor based on JetDefinition and ActiveAreaSpec Definition at line 137 of file ClusterSequenceActiveArea.hh. { // transfer the initial jets (type L) into our own array _transfer_input_jets(pseudojets); // run the clustering for active areas _initialise_and_run_AA(jet_def, area_spec, writeout_combinations); }

---

Member Function Documentation

void fastjet::ClusterSequenceActiveArea::_extract_tree (  vector< int > &   )  const [private]

routine for extracting the tree in an order that will be independent of any degeneracies in the recombination sequence that don't affect the composition of the final jets

void fastjet::ClusterSequenceActiveArea::_extract_tree_children (  int  pos, valarray< bool > &  , const valarray< int > &  , vector< int > &  )  const [private]

do the part of the extraction associated with pos, working through its children and their parents

void fastjet::ClusterSequenceActiveArea::_extract_tree_parents (  int  pos, valarray< bool > &  , const valarray< int > &  , vector< int > &  )  const [private]

do the part of the extraction associated with the parents of pos.

void fastjet::ClusterSequenceActiveArea::_initialise_and_run_AA (  const JetDefinition &  jet_def, const ActiveAreaSpec &  area_spec, const bool &  writeout_combinations = false )  [private]

does the initialisation and running specific to the active areas class Definition at line 47 of file ClusterSequenceActiveArea.cc. References _average_area, _average_area2, _average_area_4vector, fastjet::ClusterSequence::_decant_options(), fastjet::ClusterSequence::_fill_initial_history(), fastjet::ClusterSequence::_initialise_and_run(), fastjet::ClusterSequence::_jets, _maxrap_for_area, _non_jet_area, _non_jet_area2, _non_jet_number, _safe_rap_for_area, _transfer_areas(), _transfer_ghost_free_history(), fastjet::ActiveAreaSpec::ghost_maxrap(), fastjet::JetDefinition::R(), fastjet::ActiveAreaSpec::repeat(), and fastjet::ClusterSequence::unique_history_order(). { // initialize our local area information _average_area.resize(2*_jets.size()); _average_area = 0.0; _average_area2.resize(2*_jets.size()); _average_area2 = 0.0; _average_area_4vector.resize(2*_jets.size()); _average_area_4vector = PseudoJet(0.0,0.0,0.0,0.0); _non_jet_area = 0.0; _non_jet_area2 = 0.0; _non_jet_number=0.0; // for future reference... _maxrap_for_area = area_spec.ghost_maxrap(); _safe_rap_for_area = _maxrap_for_area - jet_def.R(); // Make sure we'll have at least one repetition -- then we can // deduce the unghosted clustering sequence from one of the ghosted // sequences. If we do not have any repetitions, then get the // unghosted sequence from the plain unghosted clustering. // // NB: all decanting and filling of initial history will then // be carried out by base-class routine if (area_spec.repeat() <= 0) { _initialise_and_run(jet_def, writeout_combinations); return; } // transfer all relevant info into internal variables _decant_options(jet_def, writeout_combinations); // set up the history entries for the initial particles (those // currently in _jets) _fill_initial_history(); // record the input jets as they are currently vector<PseudoJet> input_jets(_jets); // code for testing the unique tree vector<int> unique_tree; // run the clustering multiple times so as to get areas of all the jets for (int irepeat = 0; irepeat < area_spec.repeat(); irepeat++) { ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets, jet_def, area_spec); if (irepeat == 0) { // take the non-ghost part of the history and put into our own // history. _transfer_ghost_free_history(clust_seq); // get the "unique" order that will be used for transferring all areas. unique_tree = unique_history_order(); } // transfer areas from clust_seq into our object _transfer_areas(unique_tree, clust_seq); } _average_area /= area_spec.repeat(); _average_area2 /= area_spec.repeat(); if (area_spec.repeat() > 1) { _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/ (area_spec.repeat()-1)); } else { _average_area2 = 0.0; } _non_jet_area /= area_spec.repeat(); _non_jet_area2 /= area_spec.repeat(); _non_jet_area2 = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/ area_spec.repeat()); _non_jet_number /= area_spec.repeat(); // following bizarre way of writing things is related to // poverty of operations on PseudoJet objects (as well as some confusion // in one or two places) for (unsigned i = 0; i < _average_area_4vector.size(); i++) { _average_area_4vector[i] = (1.0/area_spec.repeat()) * _average_area_4vector[i]; } //cerr << "Non-jet area = " << _non_jet_area << " +- " << _non_jet_area2<<endl; }

void fastjet::ClusterSequenceActiveArea::_transfer_areas (  const vector< int > &  unique_hist_order, const ClusterSequenceActiveAreaExplicitGhosts &  )  [private]

transfer areas from the ClusterSequenceActiveAreaExplicitGhosts object into our internal area bookkeeping. .. Definition at line 386 of file ClusterSequenceActiveArea.cc. References _average_area, _average_area2, _average_area_4vector, fastjet::ClusterSequence::_history, fastjet::ClusterSequence::_initial_n, fastjet::ClusterSequence::_jet_def, fastjet::ClusterSequence::_jets, _non_jet_area, _non_jet_area2, _non_jet_number, _safe_rap_for_area, fastjet::ClusterSequenceActiveAreaExplicitGhosts::area(), area(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::area_4vector(), fastjet::ClusterSequence::BeamJet, fastjet::PseudoJet::E(), fastjet::ClusterSequence::history(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::is_pure_ghost(), fastjet::ClusterSequence::jets(), fastjet::ClusterSequence::n_particles(), fastjet::PseudoJet::perp(), fastjet::PseudoJet::perp2(), fastjet::PseudoJet::px(), fastjet::PseudoJet::py(), fastjet::PseudoJet::pz(), fastjet::JetDefinition::recombiner(), and fastjet::ClusterSequence::unique_history_order(). Referenced by _initialise_and_run_AA(). { const vector<history_element> & gs_history = ghosted_seq.history(); const vector<PseudoJet> & gs_jets = ghosted_seq.jets(); vector<int> gs_unique_hist_order = ghosted_seq.unique_history_order(); const double tolerance = 1e-13; // to decide when two jets are the same int j = -1; int hist_index = -1; valarray<double> our_areas(_history.size()); our_areas = 0.0; valarray<PseudoJet> our_area_4vectors(_history.size()); our_area_4vectors = PseudoJet(0.0,0.0,0.0,0.0); for (unsigned i = 0; i < gs_history.size(); i++) { // only consider composite particles unsigned gs_hist_index = gs_unique_hist_order[i]; if (gs_hist_index < ghosted_seq.n_particles()) continue; const history_element & gs_hist = gs_history[gs_unique_hist_order[i]]; int parent1 = gs_hist.parent1; int parent2 = gs_hist.parent2; if (parent2 == BeamJet) { // need to look at parent to get the actual jet const PseudoJet & jet = gs_jets[gs_history[parent1].jetp_index]; double area = ghosted_seq.area(jet); PseudoJet ext_area = ghosted_seq.area_4vector(jet); if (ghosted_seq.is_pure_ghost(parent1)) { if (abs(jet.rap()) < _safe_rap_for_area) { _non_jet_area += area; _non_jet_area2 += area*area; _non_jet_number += 1; } } else { // get next "combined-particle" index in our own history // making sure we don't go beyond it's bounds (if we do // then we're in big trouble anyway...) while (++j < static_cast<int>(_history.size())) { hist_index = unique_hist_order[j]; if (hist_index >= _initial_n) break;} // sanity check const PseudoJet & refjet = _jets[_history[_history[hist_index].parent1].jetp_index]; //if (jet.perp2() != refjet.perp2()) { //if (abs(jet.perp2()-refjet.perp2()) > // tolerance*max(jet.perp2(),refjet.perp2())) { // If pt disagrees check E; if they both disagree there's a // problem here... NB: a massive particle with zero pt may // have its pt changed when a ghost is added -- this is why we // also require the energy to be wrong before complaining if (abs(jet.perp2()-refjet.perp2()) > tolerance*max(jet.perp2(),refjet.perp2()) && abs(jet.E()-refjet.E()) > tolerance*max(jet.E(),refjet.E())) { cerr << jet.perp() << " " << refjet.perp() << " " << jet.perp() - refjet.perp() << endl; cerr << refjet.px() << " " << refjet.py() << " " << refjet.pz() << " " << refjet.E() << endl; throw Error("Could not match clustering sequence for an inclusive jet when reconstructing areas"); } // set the area at this clustering stage our_areas[hist_index] = area; our_area_4vectors[hist_index] = ext_area; // update the parent as well -- that way its area is the area // immediately before clustering (i.e. resolve an ambiguity in // the Cambridge case and ensure in the kt case that the original // particles get a correct area) our_areas[_history[hist_index].parent1] = area; our_area_4vectors[_history[hist_index].parent1] = ext_area; } } else if (!ghosted_seq.is_pure_ghost(parent1) && !ghosted_seq.is_pure_ghost(parent2)) { // get next "combined-particle" index in our own history while (++j < static_cast<int>(_history.size())) { hist_index = unique_hist_order[j]; if (hist_index >= _initial_n) break;} const PseudoJet & jet = gs_jets[gs_hist.jetp_index]; const PseudoJet & refjet = _jets[_history[hist_index].jetp_index]; // run sanity check if (abs(jet.perp2()-refjet.perp2()) > tolerance*max(jet.perp2(),refjet.perp2())) { cerr << jet.perp() << " " << refjet.perp() << " "<< jet.perp() - refjet.perp() << endl; throw Error("Could not match clustering sequence for an exclusive jet when reconstructing areas"); } // update area and our local index (maybe redundant since later // the descendants will reupdate it?) double area = ghosted_seq.area(jet); our_areas[hist_index] += area; PseudoJet ext_area = ghosted_seq.area_4vector(jet); //our_area_4vectors[hist_index] = our_area_4vectors[hist_index] + ext_area; _jet_def.recombiner()->plus_equal(our_area_4vectors[hist_index], ext_area); // now update areas of parents (so that they becomes areas // immediately before clustering occurred). This is of use // because it allows us to set the areas of the original hard // particles in the kt algorithm; for the Cambridge case it // means a jet's area will be the area just before it clusters // with another hard jet. const PseudoJet & jet1 = gs_jets[gs_history[parent1].jetp_index]; int our_parent1 = _history[hist_index].parent1; our_areas[our_parent1] = ghosted_seq.area(jet1); our_area_4vectors[our_parent1] = ghosted_seq.area_4vector(jet1); const PseudoJet & jet2 = gs_jets[gs_history[parent2].jetp_index]; int our_parent2 = _history[hist_index].parent2; our_areas[our_parent2] = ghosted_seq.area(jet2); our_area_4vectors[our_parent2] = ghosted_seq.area_4vector(jet2); } } _average_area += our_areas; _average_area2 += our_areas*our_areas; //_average_area_4vector += our_area_4vectors; // Use the proper recombination scheme when averaging the area_4vectors // over multiple ghost runs (i.e. the repeat stage); for (unsigned i = 0; i < _average_area_4vector.size(); i++) { _jet_def.recombiner()->plus_equal(_average_area_4vector[i], our_area_4vectors[i]); } }

void fastjet::ClusterSequenceActiveArea::_transfer_ghost_free_history (  const ClusterSequenceActiveAreaExplicitGhosts &  clust_seq  )  [private]

transfer the history (and jet-momenta) from clust_seq to our own internal structure while removing ghosts Definition at line 308 of file ClusterSequenceActiveArea.cc. References fastjet::ClusterSequence::_do_iB_recombination_step(), fastjet::ClusterSequence::_do_ij_recombination_step(), fastjet::ClusterSequence::_history, fastjet::ClusterSequence::_strategy, fastjet::ClusterSequence::BeamJet, fastjet::ClusterSequence::history(), fastjet::ClusterSequence::InexistentParent, fastjet::ClusterSequence::Invalid, fastjet::ClusterSequenceActiveAreaExplicitGhosts::is_pure_ghost(), and fastjet::ClusterSequence::strategy_used(). Referenced by _initialise_and_run_AA(). { const vector<history_element> & gs_history = ghosted_seq.history(); vector<int> gs2self_hist_map(gs_history.size()); // work our way through to first non-trivial combination unsigned igs = 0; unsigned iself = 0; while (gs_history[igs].parent1 == InexistentParent) { // record correspondence if (!ghosted_seq.is_pure_ghost(igs)) { gs2self_hist_map[igs] = iself++; } else { gs2self_hist_map[igs] = Invalid; } igs++; }; // make sure the count of non-ghost initial jets is equal to // what we already have in terms of initial jets assert(iself == _history.size()); // now actually transfer things do { // if we are a pure ghost, then go on to next round if (ghosted_seq.is_pure_ghost(igs)) { gs2self_hist_map[igs] = Invalid; continue; } const history_element & gs_hist_el = gs_history[igs]; bool parent1_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent1); bool parent2_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent2); // if exactly one parent is a ghost then maintain info about the // non-ghost correspondence for this jet, and then go on to next // recombination in the ghosted sequence if (parent1_is_ghost && !parent2_is_ghost && gs_hist_el.parent2 >= 0) { gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent2]; continue; } if (!parent1_is_ghost && parent2_is_ghost) { gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent1]; continue; } // no parents are ghosts... if (gs_hist_el.parent2 >= 0) { // recombination of two non-ghosts gs2self_hist_map[igs] = _history.size(); // record the recombination in our own sequence int newjet_k; // dummy var -- not used //cerr << igs << " " << gs_hist_el.parent1 << " " << gs_hist_el.parent2 << endl; //cerr << gs2self_hist_map[gs_hist_el.parent1] << " " << gs2self_hist_map[gs_hist_el.parent2] << endl; int jet_i = _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index; int jet_j = _history[gs2self_hist_map[gs_hist_el.parent2]].jetp_index; //cerr << "recombining "<< jet_i << " and "<< jet_j << endl; _do_ij_recombination_step(jet_i, jet_j, gs_hist_el.dij, newjet_k); } else { // we have a non-ghost that has become a beam-jet assert(gs_history[igs].parent2 == BeamJet); // record position gs2self_hist_map[igs] = _history.size(); // record the recombination in our own sequence _do_iB_recombination_step( _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index, gs_hist_el.dij); } } while (++igs < gs_history.size()); // finally transfer info about strategy used (which isn't necessarily // always the one that got asked for...) _strategy = ghosted_seq.strategy_used(); }

virtual double fastjet::ClusterSequenceActiveArea::area (  const PseudoJet &  jet  )  const [inline, virtual]

return the area associated with the given jet; this base class returns 0. Reimplemented from fastjet::ClusterSequenceWithArea. Definition at line 59 of file ClusterSequenceActiveArea.hh. References fastjet::PseudoJet::cluster_hist_index(). Referenced by _transfer_areas(), parabolic_pt_per_unit_area(), print_jets(), and pt_per_unit_area(). { return _average_area[jet.cluster_hist_index()];};

virtual PseudoJet fastjet::ClusterSequenceActiveArea::area_4vector (  const PseudoJet &  jet  )  const [inline, virtual]

return a PseudoJet whose 4-vector is defined by the following integral drap d PseudoJet("rap,phi,pt=one") * Theta("rap,phi inside jet boundary") where PseudoJet("rap,phi,pt=one") is a 4-vector with the given rapidity (rap), azimuth (phi) and pt=1, while Theta("rap,phi inside jet boundary") is a function that is 1 when rap,phi define a direction inside the jet boundary and 0 otherwise. This base class returns a null 4-vector. Reimplemented from fastjet::ClusterSequenceWithArea. Definition at line 64 of file ClusterSequenceActiveArea.hh. References fastjet::PseudoJet::cluster_hist_index(). Referenced by parabolic_pt_per_unit_area(), print_jets(), and pt_per_unit_area(). { return _average_area_4vector[jet.cluster_hist_index()];};

virtual double fastjet::ClusterSequenceActiveArea::area_error (  const PseudoJet &  jet  )  const [inline, virtual]

return the error (uncertainty) associated with the determination of the area of this jet; this base class returns 0. Reimplemented from fastjet::ClusterSequenceWithArea. Definition at line 61 of file ClusterSequenceActiveArea.hh. References fastjet::PseudoJet::cluster_hist_index(). { return _average_area2[jet.cluster_hist_index()];};

void fastjet::ClusterSequenceActiveArea::parabolic_pt_per_unit_area (  double &  a, double &  b, double  raprange = -1.0, double  exclude_above = -1.0, bool  use_area_4vector = false )  const

fits a form pt_per_unit_area(y) = a + b*y^2 in the range abs(y)<raprange (for negative raprange, it defaults to _safe_rap_for_area). Definition at line 250 of file ClusterSequenceActiveArea.cc. References _safe_rap_for_area, area(), area_4vector(), and fastjet::ClusterSequence::inclusive_jets(). { double this_raprange; if (raprange <= 0) {this_raprange = _safe_rap_for_area;} else {this_raprange = raprange;} int n=0; int n_excluded = 0; double mean_f=0, mean_x2=0, mean_x4=0, mean_fx2=0; vector<PseudoJet> incl_jets = inclusive_jets(); for (unsigned i = 0; i < incl_jets.size(); i++) { if (abs(incl_jets[i].rap()) < this_raprange) { double this_area; if ( use_area_4vector ) { this_area = area_4vector(incl_jets[i]).perp(); } else { this_area = area(incl_jets[i]); } double f = incl_jets[i].perp()/this_area; if (exclude_above <= 0.0 || f < exclude_above) { double x = incl_jets[i].rap(); double x2 = x*x; mean_f += f; mean_x2 += x2; mean_x4 += x2*x2; mean_fx2 += f*x2; n++; } else { n_excluded++; } } } if (n <= 1) { // meaningful results require at least two jets inside the // area -- mind you if there are empty jets we should be in // any case doing something special... a = 0.0; b = 0.0; } else { mean_f /= n; mean_x2 /= n; mean_x4 /= n; mean_fx2 /= n; b = (mean_f*mean_x2 - mean_fx2)/(mean_x2*mean_x2 - mean_x4); a = mean_f - b*mean_x2; } //cerr << "n_excluded = "<< n_excluded << endl; }

double fastjet::ClusterSequenceActiveArea::pt_per_unit_area (  mean_pt_strategies  strat = median, double  range = 2.0 )  const

return the transverse momentum per unit area according to one of the above strategies; for some strategies (those with "cut" in their name) the parameter "range" allows one to exclude a subset of the jets for the background estimation, those that have pt/area > median(pt/area)*range. Definition at line 138 of file ClusterSequenceActiveArea.cc. References _non_jet_area, _non_jet_number, _safe_rap_for_area, area(), area_4vector(), fastjet::ClusterSequence::inclusive_jets(), mean_ratio_cut, median, median_4vector, non_ghost_median, play, pttot_over_areatot, and pttot_over_areatot_cut. Referenced by print_jets(). { vector<PseudoJet> incl_jets = inclusive_jets(); vector<double> pt_over_areas; for (unsigned i = 0; i < incl_jets.size(); i++) { if (abs(incl_jets[i].rap()) < _safe_rap_for_area) { double this_area; if ( strat == median_4vector ) { this_area = area_4vector(incl_jets[i]).perp(); } else { this_area = area(incl_jets[i]); } pt_over_areas.push_back(incl_jets[i].perp()/this_area); } } // there is nothing inside our region, so answer will always be zero if (pt_over_areas.size() == 0) {return 0.0;} // get median (pt/area) [this is the "old" median definition. It considers // only the "real" jets in calculating the median, i.e. excluding the // only-ghost ones] sort(pt_over_areas.begin(), pt_over_areas.end()); double non_ghost_median_ratio = pt_over_areas[pt_over_areas.size()/2]; // new median definition that takes into account non-jet area (i.e. // jets composed only of ghosts), and for fractional median position // interpolates between the corresponding entries in the pt_over_areas array double nj_median_pos = (pt_over_areas.size()-1 - _non_jet_number)/2.0; double nj_median_ratio; if (nj_median_pos >= 0 && pt_over_areas.size() > 1) { int int_nj_median = int(nj_median_pos); nj_median_ratio = pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos) + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median); } else { nj_median_ratio = 0.0; } // get various forms of mean (pt/area) double pt_sum = 0.0, pt_sum_with_cut = 0.0; double area_sum = _non_jet_area, area_sum_with_cut = _non_jet_area; double ratio_sum = 0.0; double ratio_n = _non_jet_number; for (unsigned i = 0; i < incl_jets.size(); i++) { if (abs(incl_jets[i].rap()) < _safe_rap_for_area) { double this_area; if ( strat == median_4vector ) { this_area = area_4vector(incl_jets[i]).perp(); } else { this_area = area(incl_jets[i]); } pt_sum += incl_jets[i].perp(); area_sum += this_area; double ratio = incl_jets[i].perp()/this_area; if (ratio < range*nj_median_ratio) { pt_sum_with_cut += incl_jets[i].perp(); area_sum_with_cut += this_area; ratio_sum += ratio; ratio_n++; } } } if (strat == play) { double trunc_sum = 0, trunc_sumsqr = 0; vector<double> means(pt_over_areas.size()), sd(pt_over_areas.size()); for (unsigned i = 0; i < pt_over_areas.size() ; i++ ) { double ratio = pt_over_areas[i]; trunc_sum += ratio; trunc_sumsqr += ratio*ratio; means[i] = trunc_sum / (i+1); sd[i] = sqrt(abs(means[i]*means[i] - trunc_sumsqr/(i+1))); cerr << "i, means, sd: " <<i<<", "<< means[i] <<", "<<sd[i]<<", "<< sd[i]/sqrt(i+1.0)<<endl; } cout << "-----------------------------------"<<endl; for (unsigned i = 0; i <= pt_over_areas.size()/2 ; i++ ) { cout << "Median "<< i <<" = " << pt_over_areas[i]<<endl; } cout << "Number of non-jets: "<<_non_jet_number<<endl; cout << "Area of non-jets: "<<_non_jet_area<<endl; cout << "Default median position: " << (pt_over_areas.size()-1)/2.0<<endl; cout << "NJ median position: " << nj_median_pos <<endl; cout << "NJ median value: " << nj_median_ratio <<endl; return 0.0; } switch(strat) { case median: case median_4vector: return nj_median_ratio; case non_ghost_median: return non_ghost_median_ratio; case pttot_over_areatot: return pt_sum / area_sum; case pttot_over_areatot_cut: return pt_sum_with_cut / area_sum_with_cut; case mean_ratio_cut: return ratio_sum/ratio_n; default: return nj_median_ratio; } }

---

Member Data Documentation

valarray<double> fastjet::ClusterSequenceActiveArea::_average_area [private]

Definition at line 100 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_and_run_AA(), and _transfer_areas().

valarray<double> fastjet::ClusterSequenceActiveArea::_average_area2 [private]

Definition at line 100 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_and_run_AA(), and _transfer_areas().

valarray<PseudoJet> fastjet::ClusterSequenceActiveArea::_average_area_4vector [private]

Definition at line 101 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_and_run_AA(), and _transfer_areas().

double fastjet::ClusterSequenceActiveArea::_maxrap_for_area [private]

Definition at line 104 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_and_run_AA().

double fastjet::ClusterSequenceActiveArea::_non_jet_area [private]

Definition at line 102 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_and_run_AA(), _transfer_areas(), and pt_per_unit_area().

double fastjet::ClusterSequenceActiveArea::_non_jet_area2 [private]

Definition at line 102 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_and_run_AA(), and _transfer_areas().

double fastjet::ClusterSequenceActiveArea::_non_jet_number [private]

Definition at line 102 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_and_run_AA(), _transfer_areas(), and pt_per_unit_area().

double fastjet::ClusterSequenceActiveArea::_safe_rap_for_area [private]

Definition at line 105 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_and_run_AA(), _transfer_areas(), parabolic_pt_per_unit_area(), and pt_per_unit_area().

---

The documentation for this class was generated from the following files:

• include/fastjet/ClusterSequenceActiveArea.hh
• src/ClusterSequenceActiveArea.cc

---