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.

Classes
class	GhostJet
	a class for our internal storage of ghost jets More...
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 &ghost_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.
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 &ghost_spec, const bool &writeout_combinations, bool &continue_running)
void	_run_AA (const GhostedAreaSpec &ghost_spec)
void	_postprocess_AA (const GhostedAreaSpec &ghost_spec)
	run the postprocessing for the active area (and derived classes)
void	_initialise_and_run_AA (const JetDefinition &jet_def, const GhostedAreaSpec &ghost_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 std::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
std::valarray< double >	_average_area
	child classes benefit from having these at their disposal
std::valarray< double >	_average_area2
std::valarray< PseudoJet >	_average_area_4vector
Private Member Functions
void	_extract_tree (std::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, std::valarray< bool > &, const std::valarray< int > &, std::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, std::valarray< bool > &, const std::valarray< int > &, std::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	_ghost_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

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 &	ghost_spec,
		const bool &	writeout_combinations = `false`
	)

constructor based on JetDefinition and GhostedAreaSpec

Definition at line 204 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, ghost_spec, writeout_combinations);

}

Member Function Documentation

void fastjet::ClusterSequenceActiveArea::_extract_tree ( std::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::ClusterSequence::_extract_tree_children	(	int	pos,
		std::valarray< bool > &	,
		const std::valarray< int > &	,
		std::vector< int > &
	)		const `[private]`

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

Reimplemented from fastjet::ClusterSequence.

Definition at line 1019 of file ClusterSequence.cc.

                                        {
  if (!extracted[position]) {
    // that means we may have unidentified parents around, so go and
    // collect them (extracted[position]) will then be made true)
    _extract_tree_parents(position,extracted,lowest_constituent,unique_tree);
  } 
  
  // now look after the children...
  int child = _history[position].child;
  if (child  >= 0) _extract_tree_children(child,extracted,lowest_constituent,unique_tree);
}

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

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

Reimplemented from fastjet::ClusterSequence.

Definition at line 1051 of file ClusterSequence.cc.

                                        {

  if (!extracted[position]) {
    int parent1 = _history[position].parent1;
    int parent2 = _history[position].parent2;
    // where relevant order parents so that we will first treat the
    // one containing the smaller "lowest_constituent"
    if (parent1 >= 0 && parent2 >= 0) {
      if (lowest_constituent[parent1] > lowest_constituent[parent2]) 
        std::swap(parent1, parent2);
    }
    // then actually run through the parents to extract the constituents...
    if (parent1 >= 0 && !extracted[parent1]) 
      _extract_tree_parents(parent1,extracted,lowest_constituent,unique_tree);
    if (parent2 >= 0 && !extracted[parent2]) 
      _extract_tree_parents(parent2,extracted,lowest_constituent,unique_tree);
    // finally declare this position to be accounted for and push it
    // onto our list.
    unique_tree.push_back(position);
    extracted[position] = true;
  }
}

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

Definition at line 77 of file ClusterSequenceActiveArea.cc.

References fastjet::GhostedAreaSpec::ghost_maxrap(), fastjet::JetDefinition::R(), and fastjet::GhostedAreaSpec::repeat().

{

  // store this for future use
  _ghost_spec_repeat = ghost_spec.repeat();

  // make sure placeholders are there & zeroed
  _resize_and_zero_AA();
     
  // for future reference...
  _maxrap_for_area = ghost_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 (ghost_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::_initialise_and_run_AA	(	const JetDefinition &	jet_def,
		const GhostedAreaSpec &	ghost_spec,
		const bool &	writeout_combinations = `false`
	)		`[protected]`

does the initialisation and running specific to the active areas class

global routine for running active area

Definition at line 53 of file ClusterSequenceActiveArea.cc.

                                                    {

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

void fastjet::ClusterSequenceActiveArea::_postprocess_AA ( const GhostedAreaSpec & ghost_spec ) [protected]

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

Definition at line 152 of file ClusterSequenceActiveArea.cc.

References fastjet::GhostedAreaSpec::repeat().

                                                                                   {
  _average_area  /= ghost_spec.repeat();
  _average_area2 /= ghost_spec.repeat();
  if (ghost_spec.repeat() > 1) {
    // the VC compiler complains if one puts everything on a single line.
    // An alternative solution would be to use -1.0 (+single line)
    const double tmp = ghost_spec.repeat()-1;
    _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/tmp);
  } else {
    _average_area2 = 0.0;
  }

  _non_jet_area  /= ghost_spec.repeat();
  _non_jet_area2 /= ghost_spec.repeat();
  _non_jet_area2  = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/
                         ghost_spec.repeat());
  _non_jet_number /= ghost_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/ghost_spec.repeat()) * _average_area_4vector[i];
  }
  //cerr << "Non-jet area = " << _non_jet_area << " +- " << _non_jet_area2<<endl;
}

void fastjet::ClusterSequenceActiveArea::_resize_and_zero_AA ( ) [protected]

Definition at line 67 of file ClusterSequenceActiveArea.cc.

                                                     {
  // 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::_run_AA ( const GhostedAreaSpec & ghost_spec ) [protected]

Definition at line 122 of file ClusterSequenceActiveArea.cc.

References fastjet::ClusterSequenceActiveAreaExplicitGhosts::has_dangerous_particles(), and fastjet::GhostedAreaSpec::repeat().

                                                                           {
  // 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 < ghost_spec.repeat(); irepeat++) {

    ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets, 
                                                      jet_def(), ghost_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::_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 749 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().

        {

  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. See FAQ for possible explanations." << 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());
  }
}

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

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

Definition at line 559 of file ClusterSequenceActiveArea.cc.

                                                                           {

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

  /*
   * WARNING:
   *   _average_area has explicitly been sized initially to 2*jets().size()
   *   which can be bigger than our_areas (of size _history.size()
   *   if there are some unclustered particles. 
   *   So we must take care about boundaries
   */

  for (unsigned int area_index = 0; area_index<our_areas.size(); area_index++){
    _average_area[area_index]  += our_areas[area_index]; 
    _average_area2[area_index] += our_areas[area_index]*our_areas[area_index]; 
  }

  //_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 < our_area_4vectors.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 ) [protected]

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

Definition at line 476 of file ClusterSequenceActiveArea.cc.

References fastjet::ClusterSequence::history_element::dij, fastjet::ClusterSequence::history(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::is_pure_ghost(), fastjet::ClusterSequence::history_element::parent1, fastjet::ClusterSequence::history_element::parent2, and fastjet::ClusterSequence::strategy_used().

                                                                          {
  
  const vector<history_element> & gs_history  = ghosted_seq.history();
  vector<int> gs2self_hist_map(gs_history.size());

  // first transfer info about strategy used (which isn't necessarily
  // always the one that got asked for...)
  _strategy = ghosted_seq.strategy_used();

  // work our way through to first non-trivial combination
  unsigned igs = 0;
  unsigned iself = 0;
  while (igs < gs_history.size() && 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());

  // if there was no clustering in this event (e.g. SISCone passive
  // area with zero input particles, or with a pt cut on stable cones
  // that kills all jets), then don't bother with the rest (which
  // would crash!)
  if (igs == gs_history.size()) return;
  
  // 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());

}

virtual double fastjet::ClusterSequenceActiveArea::area ( const PseudoJet & ) 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().

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

virtual PseudoJet fastjet::ClusterSequenceActiveArea::area_4vector ( const PseudoJet & ) 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().

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

virtual double fastjet::ClusterSequenceActiveArea::area_error ( const PseudoJet & ) 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()];};

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

References fastjet::RangeDefinition::is_in_range().

                                                                                {
  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 /= _ghost_spec_repeat;
  return empty;
}

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 149 of file ClusterSequenceActiveArea.hh.

{return _has_dangerous_particles;}

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

References 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 /= _ghost_spec_repeat;
  return inrange;
}

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.

NB: This call is OBSOLETE; use media_pt_per_unit_area from the

Definition at line 276 of file ClusterSequenceActiveArea.cc.

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

}