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
	ClusterSequenceActiveArea ()
	default constructor
template<class L>
	ClusterSequenceActiveArea (const std::vector< L > &pseudojets, const JetDefinition &jet_def, const GhostedAreaSpec &area_spec, const bool &writeout_combinations=false)
	constructor based on JetDefinition and GhostedAreaSpec
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).
virtual double	empty_area (const RangeDefinition &range) const
	rewrite the empty area from the parent class, so as to use all info at our disposal return the total area, in the given y-phi range, that consists of ghost jets or unclustered ghosts
virtual double	n_empty_jets (const RangeDefinition &range) const
	return the true number of empty jets (replaces ClusterSequenceAreaBase::n_empty_jets(.
Protected Member Functions
void	_resize_and_zero_AA ()
void	_initialise_AA (const JetDefinition &jet_def, const GhostedAreaSpec &area_spec, const bool &writeout_combinations, bool &continue_running)
void	_run_AA (const GhostedAreaSpec &area_spec)
void	_postprocess_AA (const GhostedAreaSpec &area_spec)
	run the postprocessing for the active area (and derived classes)
void	_initialise_and_run_AA (const JetDefinition &jet_def, const GhostedAreaSpec &area_spec, const bool &writeout_combinations=false)
	global routine for running active area
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.
bool	has_dangerous_particles () const
	returns true if there are any particles whose transverse momentum if so low that there's a risk of the ghosts having modified the clustering sequence
Protected Attributes
valarray< double >	_average_area
	child classes benefit from having these at their disposal
valarray< double >	_average_area2
valarray< PseudoJet >	_average_area_4vector
Private Member Functions
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.
void	_throw_unless_jets_have_same_perp_or_E (const PseudoJet &jet, const PseudoJet &refjet, double tolerance, const ClusterSequenceActiveAreaExplicitGhosts &jets_ghosted_seq) const
	check if two jets have the same momentum to within the tolerance (and if pt's are not the same we're forgiving and look to see if the energy is the same)
Private Attributes
double	_non_jet_area
double	_non_jet_area2
double	_non_jet_number
double	_maxrap_for_area
double	_safe_rap_for_area
bool	_has_dangerous_particles
int	_area_spec_repeat
	since we are playing nasty games with seeds, we should warn the user a few times
std::vector< GhostJet >	_ghost_jets
std::vector< GhostJet >	_unclustered_ghosts
Classes
class	GhostJet
	a class for our internal storage of ghost jets More...

---

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 56 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.

Enumerator:

median
non_ghost_median
pttot_over_areatot
pttot_over_areatot_cut
mean_ratio_cut
play
median_4vector

Definition at line 80 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

fastjet::ClusterSequenceActiveArea::ClusterSequenceActiveArea

(

)

[inline]

default constructor

Definition at line 60 of file ClusterSequenceActiveArea.hh.

{}

template<class L>

fastjet::ClusterSequenceActiveArea::ClusterSequenceActiveArea	(	const std::vector< L > &	pseudojets,
		const JetDefinition &	jet_def,
		const GhostedAreaSpec &	area_spec,
		const bool &	writeout_combinations = false
	)			[inline]

constructor based on JetDefinition and GhostedAreaSpec

Definition at line 199 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

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::ClusterSequenceAreaBase.

Definition at line 69 of file ClusterSequenceActiveArea.hh.

References fastjet::PseudoJet::cluster_hist_index().

Referenced by _transfer_areas(), empty_area(), parabolic_pt_per_unit_area(), and pt_per_unit_area().

                                                    {
                             return _average_area[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::ClusterSequenceAreaBase.

Definition at line 71 of file ClusterSequenceActiveArea.hh.

References fastjet::PseudoJet::cluster_hist_index().

                                                          {
                             return _average_area2[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::ClusterSequenceAreaBase.

Definition at line 74 of file ClusterSequenceActiveArea.hh.

References fastjet::PseudoJet::cluster_hist_index().

Referenced by parabolic_pt_per_unit_area(), and pt_per_unit_area().

                                                               {
                    return _average_area_4vector[jet.cluster_hist_index()];};

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 271 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, fastjet::PseudoJet::perp(), play, pttot_over_areatot, and pttot_over_areatot_cut.

                                                                     {
  
  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;
  }

}

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 383 of file ClusterSequenceActiveArea.cc.

References _safe_rap_for_area, area(), area_4vector(), fastjet::ClusterSequence::inclusive_jets(), and fastjet::PseudoJet::perp().

                                    {
  
  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::empty_area

(

const RangeDefinition &

range

)

const [virtual]

rewrite the empty area from the parent class, so as to use all info at our disposal return the total area, in the given y-phi range, that consists of ghost jets or unclustered ghosts

Reimplemented from fastjet::ClusterSequenceAreaBase.

Reimplemented in fastjet::ClusterSequencePassiveArea.

Definition at line 439 of file ClusterSequenceActiveArea.cc.

References _area_spec_repeat, _ghost_jets, _unclustered_ghosts, area(), and fastjet::RangeDefinition::is_in_range().

Referenced by fastjet::ClusterSequencePassiveArea::empty_area().

                                                                                {
  double empty = 0.0;
  // first deal with ghost jets
  for (unsigned  i = 0; i < _ghost_jets.size(); i++) {
    if (range.is_in_range(_ghost_jets[i])) {
      empty += _ghost_jets[i].area;
    }
  }
  // then deal with unclustered ghosts
  for (unsigned  i = 0; i < _unclustered_ghosts.size(); i++) {
    if (range.is_in_range(_unclustered_ghosts[i])) {
      empty += _unclustered_ghosts[i].area;
    }
  }
  empty /= _area_spec_repeat;
  return empty;
}

double fastjet::ClusterSequenceActiveArea::n_empty_jets

(

const RangeDefinition &

range

)

const [virtual]

return the true number of empty jets (replaces ClusterSequenceAreaBase::n_empty_jets(.

..))

Reimplemented from fastjet::ClusterSequenceAreaBase.

Reimplemented in fastjet::ClusterSequence1GhostPassiveArea.

Definition at line 458 of file ClusterSequenceActiveArea.cc.

References _area_spec_repeat, _ghost_jets, and fastjet::RangeDefinition::is_in_range().

                                                                                  {
  double inrange = 0;
  for (unsigned  i = 0; i < _ghost_jets.size(); i++) {
    if (range.is_in_range(_ghost_jets[i])) inrange++;
  }
  inrange /= _area_spec_repeat;
  return inrange;
}

void fastjet::ClusterSequenceActiveArea::_resize_and_zero_AA

(

)

[protected]

Definition at line 64 of file ClusterSequenceActiveArea.cc.

References _average_area, _average_area2, _average_area_4vector, fastjet::ClusterSequence::_jets, _non_jet_area, _non_jet_area2, and _non_jet_number.

Referenced by _initialise_AA(), and fastjet::ClusterSequencePassiveArea::_initialise_and_run_PA().

                                                     {
  // 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;
}

void fastjet::ClusterSequenceActiveArea::_initialise_AA	(	const JetDefinition &	jet_def,
		const GhostedAreaSpec &	area_spec,
		const bool &	writeout_combinations,
		bool &	continue_running
	)			[protected]

Definition at line 74 of file ClusterSequenceActiveArea.cc.

References _area_spec_repeat, fastjet::ClusterSequence::_decant_options(), fastjet::ClusterSequence::_fill_initial_history(), _has_dangerous_particles, fastjet::ClusterSequence::_initialise_and_run(), _maxrap_for_area, _resize_and_zero_AA(), _safe_rap_for_area, fastjet::GhostedAreaSpec::ghost_maxrap(), fastjet::ClusterSequence::jet_def(), fastjet::JetDefinition::R(), and fastjet::GhostedAreaSpec::repeat().

Referenced by fastjet::ClusterSequence1GhostPassiveArea::_initialise_and_run_1GPA(), and _initialise_and_run_AA().

{

  // store this for future use
  _area_spec_repeat = area_spec.repeat();

  // make sure placeholders are there & zeroed
  _resize_and_zero_AA();
     
  // 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);
    continue_running = false;
    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();

  // by default it does not...
  _has_dangerous_particles = false;
  
  continue_running = true;
}

void fastjet::ClusterSequenceActiveArea::_run_AA

(

const GhostedAreaSpec &

area_spec

)

[protected]

Definition at line 119 of file ClusterSequenceActiveArea.cc.

References _has_dangerous_particles, fastjet::ClusterSequence::_jets, _transfer_areas(), _transfer_ghost_free_history(), fastjet::ClusterSequence::jet_def(), fastjet::GhostedAreaSpec::repeat(), and fastjet::ClusterSequence::unique_history_order().

Referenced by _initialise_and_run_AA().

                                                                          {
  // 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);

    _has_dangerous_particles |= clust_seq.has_dangerous_particles();
    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);
  }
}

void fastjet::ClusterSequenceActiveArea::_postprocess_AA

(

const GhostedAreaSpec &

area_spec

)

[protected]

run the postprocessing for the active area (and derived classes)

Definition at line 149 of file ClusterSequenceActiveArea.cc.

References _average_area, _average_area2, _average_area_4vector, _non_jet_area, _non_jet_area2, _non_jet_number, and fastjet::GhostedAreaSpec::repeat().

Referenced by fastjet::ClusterSequence1GhostPassiveArea::_initialise_and_run_1GPA(), and _initialise_and_run_AA().

                                                                                  {
  _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::_initialise_and_run_AA	(	const JetDefinition &	jet_def,
		const GhostedAreaSpec &	area_spec,
		const bool &	writeout_combinations = false
	)			[protected]

global routine for running active area

Definition at line 50 of file ClusterSequenceActiveArea.cc.

References _initialise_AA(), _postprocess_AA(), _run_AA(), and fastjet::ClusterSequence::jet_def().

Referenced by fastjet::ClusterSequencePassiveArea::_initialise_and_run_PA().

                                                    {

  bool continue_running;
  _initialise_AA(jet_def,  area_spec, writeout_combinations, continue_running);
  if (continue_running) {
    _run_AA(area_spec);
    _postprocess_AA(area_spec);
  }
}

void fastjet::ClusterSequenceActiveArea::_transfer_ghost_free_history

(

const ClusterSequenceActiveAreaExplicitGhosts &

clust_seq

)

[protected]

transfer the history (and jet-momenta) from clust_seq to our own internal structure while removing ghosts

Definition at line 470 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 fastjet::ClusterSequence1GhostPassiveArea::_run_1GPA(), 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
      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();
}

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

transfer areas from the ClusterSequenceActiveAreaExplicitGhosts object into our internal area bookkeeping.

..

Definition at line 546 of file ClusterSequenceActiveArea.cc.

References _average_area, _average_area2, _average_area_4vector, _ghost_jets, 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, _throw_unless_jets_have_same_perp_or_E(), _unclustered_ghosts, fastjet::ClusterSequenceActiveAreaExplicitGhosts::area(), area(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::area_4vector(), fastjet::ClusterSequence::BeamJet, fastjet::ClusterSequence::history(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::is_pure_ghost(), fastjet::ClusterSequence::jets(), fastjet::ClusterSequence::n_particles(), fastjet::JetDefinition::recombiner(), fastjet::ClusterSequence::unclustered_particles(), and fastjet::ClusterSequence::unique_history_order().

Referenced by fastjet::ClusterSequence1GhostPassiveArea::_run_1GPA(), 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-11; // 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)) {
        // record the existence of the pure ghost jet for future use
        _ghost_jets.push_back(GhostJet(jet,area));
        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 its bounds (if we do
        // then we're in big trouble anyway...)
        while (++j < int(_history.size())) {
          hist_index = unique_hist_order[j];
          if (hist_index >= _initial_n) break;}

        // sanity checking -- do not overrun
        if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in diB matching");

        // sanity check -- make sure we are taking about the same 
        // jet in reference and new sequences
        int refjet_index = _history[_history[hist_index].parent1].jetp_index;
        assert(refjet_index >= 0 && refjet_index < int(_jets.size()));
        const PseudoJet & refjet = _jets[refjet_index];

      //cerr << "Inclusive" << endl;
      //cerr << gs_history[parent1].jetp_index << " " << gs_jets.size() << endl;
      //cerr << _history[_history[hist_index].parent1].jetp_index << " " << _jets.size() << endl;

        // 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
        _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
                                               ghosted_seq);

        // 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 < int(_history.size())) {
        hist_index = unique_hist_order[j];
        if (hist_index >= _initial_n) break;}
      
      // sanity checking -- do not overrun
      if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in dij matching");

      // make sure that our reference history entry is also for
      // an exclusive (dij) clustering (otherwise the comparison jet
      // will not exist)
      if (_history[hist_index].parent2 == BeamJet) throw Error("ClusterSequenceActiveArea: could not match clustering sequences (encountered dij matched with diB)");

      //cerr << "Exclusive: hist_index,hist_size: " << hist_index << " " << _history.size()<< endl;
      //cerr << gs_hist.jetp_index << " " << gs_jets.size() << endl;
      //cerr << _history[hist_index].jetp_index << " " << _jets.size() << endl;

      const PseudoJet & jet = gs_jets[gs_hist.jetp_index];
      const PseudoJet & refjet = _jets[_history[hist_index].jetp_index];

      // run sanity check 
      _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
                                             ghosted_seq);

      // 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);

      // GPS TMP debugging (jetclu) -----------------------
      //ext_area = PseudoJet(1e-100,1e-100,1e-100,4e-100);
      //our_area_4vectors[hist_index] = ext_area;
      //cout << "aa " 
      //     << our_area_4vectors[hist_index].px() << " "
      //     << our_area_4vectors[hist_index].py() << " "
      //     << our_area_4vectors[hist_index].pz() << " "
      //     << our_area_4vectors[hist_index].E() << endl;
      //cout << "bb " 
      //     << ext_area.px() << " "
      //     << ext_area.py() << " "
      //     << ext_area.pz() << " "
      //     << ext_area.E() << endl;
      //---------------------------------------------------

      _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);
    }

  }

  // now add unclustered ghosts to the relevant list so that we can
  // calculate empty area later.
  vector<PseudoJet> unclust = ghosted_seq.unclustered_particles();
  for (unsigned iu = 0; iu < unclust.size();  iu++) {
    if (ghosted_seq.is_pure_ghost(unclust[iu])) {
      double area = ghosted_seq.area(unclust[iu]);
      _unclustered_ghosts.push_back(GhostJet(unclust[iu],area));
    }
  }

  _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]);
  }
}

bool fastjet::ClusterSequenceActiveArea::has_dangerous_particles

(

)

const [inline, protected]

returns true if there are any particles whose transverse momentum if so low that there's a risk of the ghosts having modified the clustering sequence

Definition at line 144 of file ClusterSequenceActiveArea.hh.

{return _has_dangerous_particles;}

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::_throw_unless_jets_have_same_perp_or_E	(	const PseudoJet &	jet,
		const PseudoJet &	refjet,
		double	tolerance,
		const ClusterSequenceActiveAreaExplicitGhosts &	jets_ghosted_seq
	)			const [private]

check if two jets have the same momentum to within the tolerance (and if pt's are not the same we're forgiving and look to see if the energy is the same)

Definition at line 726 of file ClusterSequenceActiveArea.cc.

References fastjet::PseudoJet::E(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::has_dangerous_particles(), fastjet::PseudoJet::perp2(), fastjet::PseudoJet::px(), fastjet::PseudoJet::py(), and fastjet::PseudoJet::pz().

Referenced by _transfer_areas().

        {

  if (abs(jet.perp2()-refjet.perp2()) > 
      tolerance*max(jet.perp2(),refjet.perp2())
      && abs(jet.E()-refjet.E()) > tolerance*max(jet.E(),refjet.E())) {
    ostringstream ostr;
    ostr << "Could not match clustering sequence for an inclusive/exclusive jet when reconstructing areas" << endl;
    ostr << "  Ref-Jet: "
         << refjet.px() << " " 
         << refjet.py() << " " 
         << refjet.pz() << " " 
         << refjet.E() << endl;
    ostr << "  New-Jet: "
         << jet.px() << " " 
         << jet.py() << " " 
         << jet.pz() << " " 
         << jet.E() << endl;
    if (jets_ghosted_seq.has_dangerous_particles()) {
      ostr << "  NB: some particles have pt too low wrt ghosts -- this may be the cause" << endl;}
    //ostr << jet.perp() << " " << refjet.perp() << " "
    //     << jet.perp() - refjet.perp() << endl;
    throw Error(ostr.str());
  }
}

---

Member Data Documentation

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

child classes benefit from having these at their disposal

Definition at line 138 of file ClusterSequenceActiveArea.hh.

Referenced by fastjet::ClusterSequencePassiveArea::_initialise_and_run_PA(), _postprocess_AA(), _resize_and_zero_AA(), fastjet::ClusterSequence1GhostPassiveArea::_run_1GPA(), and _transfer_areas().

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

Definition at line 138 of file ClusterSequenceActiveArea.hh.

Referenced by _postprocess_AA(), _resize_and_zero_AA(), fastjet::ClusterSequence1GhostPassiveArea::_run_1GPA(), and _transfer_areas().

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

Definition at line 139 of file ClusterSequenceActiveArea.hh.

Referenced by fastjet::ClusterSequencePassiveArea::_initialise_and_run_PA(), _postprocess_AA(), _resize_and_zero_AA(), and _transfer_areas().

double fastjet::ClusterSequenceActiveArea::_non_jet_area [private]

Definition at line 149 of file ClusterSequenceActiveArea.hh.

Referenced by _postprocess_AA(), _resize_and_zero_AA(), _transfer_areas(), and pt_per_unit_area().

double fastjet::ClusterSequenceActiveArea::_non_jet_area2 [private]

Definition at line 149 of file ClusterSequenceActiveArea.hh.

Referenced by _postprocess_AA(), _resize_and_zero_AA(), and _transfer_areas().

double fastjet::ClusterSequenceActiveArea::_non_jet_number [private]

Definition at line 149 of file ClusterSequenceActiveArea.hh.

Referenced by _postprocess_AA(), _resize_and_zero_AA(), _transfer_areas(), and pt_per_unit_area().

double fastjet::ClusterSequenceActiveArea::_maxrap_for_area [private]

Definition at line 151 of file ClusterSequenceActiveArea.hh.

Referenced by _initialise_AA().

double fastjet::ClusterSequenceActiveArea::_safe_rap_for_area [private]

Definition at line 152 of file ClusterSequenceActiveArea.hh.

Referenced by _initialise_AA(), _transfer_areas(), parabolic_pt_per_unit_area(), and pt_per_unit_area().

bool fastjet::ClusterSequenceActiveArea::_has_dangerous_particles [private]

Definition at line 154 of file ClusterSequenceActiveArea.hh.

Referenced by _initialise_AA(), and _run_AA().

int fastjet::ClusterSequenceActiveArea::_area_spec_repeat [private]

since we are playing nasty games with seeds, we should warn the user a few times

Definition at line 182 of file ClusterSequenceActiveArea.hh.

Referenced by _initialise_AA(), empty_area(), and n_empty_jets().

std::vector<GhostJet> fastjet::ClusterSequenceActiveArea::_ghost_jets [private]

Definition at line 191 of file ClusterSequenceActiveArea.hh.

Referenced by _transfer_areas(), empty_area(), and n_empty_jets().

std::vector<GhostJet> fastjet::ClusterSequenceActiveArea::_unclustered_ghosts [private]

Definition at line 192 of file ClusterSequenceActiveArea.hh.

Referenced by _transfer_areas(), and empty_area().

---

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

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

---