source: svn/trunk/Utilities/Fastjet/src/ClusterSequence.cc@ 196

//STARTHEADER
// $Id: ClusterSequence.cc,v 1.1 2008-11-06 14:32:13 ovyn Exp $
//
// Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
//  FastJet is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.
//
//  The algorithms that underlie FastJet have required considerable
//  development and are described in hep-ph/0512210. If you use
//  FastJet as part of work towards a scientific publication, please
//  include a citation to the FastJet paper.
//
//  FastJet is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with FastJet; if not, write to the Free Software
//  Foundation, Inc.:
//      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//----------------------------------------------------------------------
//ENDHEADER

#include "../include/fastjet/Error.hh"
#include "../include/fastjet/PseudoJet.hh"
#include "../include/fastjet/ClusterSequence.hh"
#include "../include/fastjet/version.hh" // stores the current version number
#include<iostream>
#include<sstream>
#include<cmath>
#include<cstdlib>
#include<cassert>
#include<string>

FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh

using namespace std;

//// initialised static member has to go in the .cc code
JetAlgorithm ClusterSequence::_default_jet_algorithm = kt_algorithm;
//


//----------------------------------------------------------------------
void ClusterSequence::_initialise_and_run (
                                  const double & R,
                                  const Strategy & strategy,
                                  const bool & writeout_combinations) {

  JetDefinition jet_def(_default_jet_algorithm, R, strategy);
  _initialise_and_run(jet_def, writeout_combinations);
}


//----------------------------------------------------------------------
void ClusterSequence::_initialise_and_run (
                                  const JetDefinition & jet_def,
                                  const bool & writeout_combinations) {

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

  // run the plugin if that's what's decreed
  if (_jet_algorithm == plugin_algorithm) {
    _plugin_activated = true;
    _jet_def.plugin()->run_clustering( (*this) );
    _plugin_activated = false;
    return;
  }


  // automatically redefine the strategy according to N if that is
  // what the user requested -- transition points (and especially
  // their R-dependence) are based on empirical observations for a
  // R=0.4, 0.7 and 1.0, running on toth (3.4GHz, Pentium IV D [dual
  // core] with 2MB of cache).
  if (_strategy == Best) {
    int N = _jets.size();
    if (N > 6200/pow(_Rparam,2.0) 
        && jet_def.jet_algorithm() == cambridge_algorithm) {
      _strategy = NlnNCam;}
    else
#ifndef DROP_CGAL
    if (N > 16000/pow(_Rparam,1.15)) {
      _strategy = NlnN; }   
    else                    
#endif  // DROP_CGAL
      if (N > 450) {
      _strategy = N2MinHeapTiled;
    }
    else if (N > 55*max(0.5,min(1.0,_Rparam))) {// empirical scaling with R
      _strategy = N2Tiled;
    } else {
      _strategy = N2Plain;
    }
  }


  // run the code containing the selected strategy
  if (_strategy == NlnN || _strategy == NlnN3pi 
      || _strategy == NlnN4pi ) {
    this->_delaunay_cluster();
  } else if (_strategy ==  N3Dumb ) {
    this->_really_dumb_cluster();
  } else if (_strategy == N2Tiled) {
    this->_faster_tiled_N2_cluster();
  } else if (_strategy == N2PoorTiled) {
    this->_tiled_N2_cluster();
  } else if (_strategy == N2Plain) {
    this->_simple_N2_cluster();
  } else if (_strategy == N2MinHeapTiled) {
    this->_minheap_faster_tiled_N2_cluster();
  } else if (_strategy == NlnNCam4pi) {
    this->_CP2DChan_cluster();
  } else if (_strategy == NlnNCam2pi2R) {
    this->_CP2DChan_cluster_2pi2R();
  } else if (_strategy == NlnNCam) {
    this->_CP2DChan_cluster_2piMultD();
  } else {
    ostringstream err;
    err << "Unrecognised value for strategy: "<<_strategy;
    throw Error(err.str());
    //assert(false);
  }
}


// these needs to be defined outside the class definition.
bool ClusterSequence::_first_time = true;
int ClusterSequence::_n_exclusive_warnings = 0;


//----------------------------------------------------------------------
// the version string
string fastjet_version_string() {
  return "FastJet version "+string(fastjet_version);
}


//----------------------------------------------------------------------
// prints a banner on the first call
void ClusterSequence::_print_banner() {

  if (!_first_time) {return;}
  _first_time = false;
  
  
  //Symp. Discr. Alg, p.472 (2002) and  CGAL (http://www.cgal.org);

  cout << "#--------------------------------------------------------------------------\n";
  cout << "#                      FastJet release " << fastjet_version << endl;
  cout << "#            Written by M. Cacciari, G.P. Salam and G. Soyez            \n"; 
  cout << "#              http://www.lpthe.jussieu.fr/~salam/fastjet               \n"; 
  cout << "#                                                                       \n";
  cout << "# Longitudinally invariant Kt, anti-Kt, and inclusive Cambridge/Aachen  \n";
  cout << "# clustering using fast geometric algorithms, with area measures and optional\n";
  cout << "# external jet-finder plugins.                                          \n";
  cout << "# Please cite Phys. Lett. B641 (2006) [hep-ph/0512210] if you use this code.\n";
  cout << "#                                                                       \n";
  cout << "# This package uses T.Chan's closest pair algorithm, Proc.13th ACM-SIAM \n";
  cout << "# Symp. Discr. Alg, p.472 (2002), S.Fortune's Voronoi algorithm and code " ;
#ifndef DROP_CGAL
  cout << endl << "# and CGAL: http://www.cgal.org/";
#endif  // DROP_CGAL
  cout << ".\n";
  cout << "#-------------------------------------------------------------------------\n";
}

//----------------------------------------------------------------------
// transfer all relevant info into internal variables
void ClusterSequence::_decant_options(const JetDefinition & jet_def,
                                      const bool & writeout_combinations) {

  // let the user know what's going on
  _print_banner();

  // make a local copy of the jet definition (for future use?)
  _jet_def = jet_def;
  
  _writeout_combinations = writeout_combinations;
  _jet_algorithm = jet_def.jet_algorithm();
  _Rparam = jet_def.R();  _R2 = _Rparam*_Rparam; _invR2 = 1.0/_R2;
  _strategy = jet_def.strategy();

  // disallow interference from the plugin
  _plugin_activated = false;
  
}


//----------------------------------------------------------------------
// initialise the history in a standard way
void ClusterSequence::_fill_initial_history () {

  if (_jets.size() == 0) {throw Error("Cannot run jet-finder on empty event");}

  // reserve sufficient space for everything
  _jets.reserve(_jets.size()*2);
  _history.reserve(_jets.size()*2);

  for (int i = 0; i < static_cast<int>(_jets.size()) ; i++) {
    history_element element;
    element.parent1 = InexistentParent;
    element.parent2 = InexistentParent;
    element.child   = Invalid;
    element.jetp_index = i;
    element.dij     = 0.0;
    element.max_dij_so_far = 0.0;

    _history.push_back(element);
    
    // do any momentum preprocessing needed by the recombination scheme
    _jet_def.recombiner()->preprocess(_jets[i]);

    // get cross-referencing right from PseudoJets
    _jets[i].set_cluster_hist_index(i);
  }
  _initial_n = _jets.size();
}


//----------------------------------------------------------------------
// Return the component corresponding to the specified index.
// taken from CLHEP
string ClusterSequence::strategy_string ()  const {
  string strategy;
  switch(_strategy) {
  case NlnN:
    strategy = "NlnN"; break;
  case NlnN3pi:
    strategy = "NlnN3pi"; break;
  case NlnN4pi:
    strategy = "NlnN4pi"; break;
  case N2Plain:
    strategy = "N2Plain"; break;
  case N2Tiled:
    strategy = "N2Tiled"; break;
  case N2MinHeapTiled:
    strategy = "N2MinHeapTiled"; break;
  case N2PoorTiled:
    strategy = "N2PoorTiled"; break;
  case N3Dumb:
    strategy = "N3Dumb"; break;
  case NlnNCam4pi:
    strategy = "NlnNCam4pi"; break;
  case NlnNCam2pi2R:
    strategy = "NlnNCam2pi2R"; break;
  case NlnNCam:
    strategy = "NlnNCam"; break; // 2piMultD
  case plugin_strategy:
    strategy = "plugin strategy"; break;
  default:
    strategy = "Unrecognized";
  }
  return strategy;
}  


double ClusterSequence::jet_scale_for_algorithm(
                                  const PseudoJet & jet) const {
  if (_jet_algorithm == kt_algorithm)             {return jet.kt2();}
  else if (_jet_algorithm == cambridge_algorithm) {return 1.0;}
  else if (_jet_algorithm == antikt_algorithm) {
    double kt2=jet.kt2();
    return kt2 > 1e-300 ? 1.0/kt2 : 1e300;
  } else if (_jet_algorithm == cambridge_for_passive_algorithm) {
    double kt2 = jet.kt2();
    double lim = _jet_def.extra_param();
    if (kt2 < lim*lim && kt2 != 0.0) {
      return 1.0/kt2;
    } else {return 1.0;}
  } else {throw Error("Unrecognised jet finder");}
}


//----------------------------------------------------------------------
/// transfer the sequence contained in other_seq into our own;
/// any plugin "extras" contained in the from_seq will be lost
/// from there.
void ClusterSequence::transfer_from_sequence(ClusterSequence & from_seq) {

  // the metadata
  _jet_def                 = from_seq._jet_def                ;
  _writeout_combinations   = from_seq._writeout_combinations  ;
  _initial_n               = from_seq._initial_n              ;
  _Rparam                  = from_seq._Rparam                 ;
  _R2                      = from_seq._R2                     ;
  _invR2                   = from_seq._invR2                  ;
  _strategy                = from_seq._strategy               ;
  _jet_algorithm           = from_seq._jet_algorithm          ;
  _plugin_activated        = from_seq._plugin_activated       ;

  // the data
  _jets     = from_seq._jets;
  _history  = from_seq._history;
  // the following transferse ownership of the extras from the from_seq
  _extras   = from_seq._extras;

}

//----------------------------------------------------------------------
// record an ij recombination and reset the _jets[newjet_k] momentum and
// user index to be those of newjet
void ClusterSequence::plugin_record_ij_recombination(
           int jet_i, int jet_j, double dij, 
           const PseudoJet & newjet, int & newjet_k) {

  plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k);

  // now transfer newjet into place
  int tmp_index = _jets[newjet_k].cluster_hist_index();
  _jets[newjet_k] = newjet;
  _jets[newjet_k].set_cluster_hist_index(tmp_index);
}


//----------------------------------------------------------------------
// return all inclusive jets with pt > ptmin
vector<PseudoJet> ClusterSequence::inclusive_jets (const double & ptmin) const{
  double dcut = ptmin*ptmin;
  int i = _history.size() - 1; // last jet
  vector<PseudoJet> jets;
  if (_jet_algorithm == kt_algorithm) {
    while (i >= 0) {
      // with our specific definition of dij and diB (i.e. R appears only in 
      // dij), then dij==diB is the same as the jet.perp2() and we can exploit
      // this in selecting the jets...
      if (_history[i].max_dij_so_far < dcut) {break;}
      if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) {
        // for beam jets
        int parent1 = _history[i].parent1;
        jets.push_back(_jets[_history[parent1].jetp_index]);}
      i--;
    }
  } else if (_jet_algorithm == cambridge_algorithm) {
    while (i >= 0) {
      // inclusive jets are all at end of clustering sequence in the
      // Cambridge algorithm -- so if we find a non-exclusive jet, then
      // we can exit
      if (_history[i].parent2 != BeamJet) {break;}
      int parent1 = _history[i].parent1;
      const PseudoJet & jet = _jets[_history[parent1].jetp_index];
      if (jet.perp2() >= dcut) {jets.push_back(jet);}
      i--;
    }
  } else if (_jet_algorithm == plugin_algorithm 
             || _jet_algorithm == antikt_algorithm
             || _jet_algorithm == cambridge_for_passive_algorithm) {
    // for inclusive jets with a plugin algorithm, we make no
    // assumptions about anything (relation of dij to momenta,
    // ordering of the dij, etc.)
    while (i >= 0) {
      if (_history[i].parent2 == BeamJet) {
        int parent1 = _history[i].parent1;
        const PseudoJet & jet = _jets[_history[parent1].jetp_index];
        if (jet.perp2() >= dcut) {jets.push_back(jet);}
      }
      i--;
    }
  } else {throw Error("Unrecognized jet algorithm");}
  return jets;
}


//----------------------------------------------------------------------
// return the number of exclusive jets that would have been obtained
// running the algorithm in exclusive mode with the given dcut
int ClusterSequence::n_exclusive_jets (const double & dcut) const {

  // first locate the point where clustering would have stopped (i.e. the
  // first time max_dij_so_far > dcut)
  int i = _history.size() - 1; // last jet
  while (i >= 0) {
    if (_history[i].max_dij_so_far <= dcut) {break;}
    i--;
  }
  int stop_point = i + 1;
  // relation between stop_point, njets assumes one extra jet disappears
  // at each clustering.
  int njets = 2*_initial_n - stop_point;
  return njets;
}

//----------------------------------------------------------------------
// return all exclusive jets that would have been obtained running
// the algorithm in exclusive mode with the given dcut
vector<PseudoJet> ClusterSequence::exclusive_jets (const double & dcut) const {
  int njets = n_exclusive_jets(dcut);
  return exclusive_jets(njets);
}


//----------------------------------------------------------------------
// return the jets obtained by clustering the event to n jets.
vector<PseudoJet> ClusterSequence::exclusive_jets (const int & njets) const {

  // make sure the user does not ask for more than jets than there
  // were particles in the first place.
  assert (njets <= _initial_n);

  // provide a warning when extracting exclusive jets 
  if (_jet_def.jet_algorithm() != kt_algorithm && _n_exclusive_warnings < 5) {
    _n_exclusive_warnings++;
    cerr << "FastJet WARNING: dcut and exclusive jets for jet-finders other than kt should be interpreted with care." << endl;
  }


  // calculate the point where we have to stop the clustering.
  // relation between stop_point, njets assumes one extra jet disappears
  // at each clustering.
  int stop_point = 2*_initial_n - njets;

  // some sanity checking to make sure that e+e- does not give us
  // surprises (should we ever implement e+e-)...
  if (2*_initial_n != static_cast<int>(_history.size())) {
    ostringstream err;
    err << "2*_initial_n != _history.size() -- this endangers internal assumptions!\n";
    throw Error(err.str());
    //assert(false);
  }

  // now go forwards and reconstitute the jets that we have --
  // basically for any history element, see if the parent jets to
  // which it refers were created before the stopping point -- if they
  // were then add them to the list, otherwise they are subsequent
  // recombinations of the jets that we are looking for.
  vector<PseudoJet> jets;
  for (unsigned int i = stop_point; i < _history.size(); i++) {
    int parent1 = _history[i].parent1;
    if (parent1 < stop_point) {
      jets.push_back(_jets[_history[parent1].jetp_index]);
    }
    int parent2 = _history[i].parent2;
    if (parent2 < stop_point && parent2 > 0) {
      jets.push_back(_jets[_history[parent2].jetp_index]);
    }
    
  }

  // sanity check...
  if (static_cast<int>(jets.size()) != njets) {
    ostringstream err;
    err << "ClusterSequence::exclusive_jets: size of returned vector ("
         <<jets.size()<<") does not coincide with requested number of jets ("
         <<njets<<")";
    throw Error(err.str());
  }

  return jets;
}

//----------------------------------------------------------------------
/// return the dmin corresponding to the recombination that went from
/// n+1 to n jets
double ClusterSequence::exclusive_dmerge (const int & njets) const {
  assert(njets >= 0);
  if (njets >= _initial_n) {return 0.0;}
  return _history[2*_initial_n-njets-1].dij;
}


//----------------------------------------------------------------------
/// return the maximum of the dmin encountered during all recombinations 
/// up to the one that led to an n-jet final state; identical to
/// exclusive_dmerge, except in cases where the dmin do not increase
/// monotonically.
double ClusterSequence::exclusive_dmerge_max (const int & njets) const {
  assert(njets >= 0);
  if (njets >= _initial_n) {return 0.0;}
  return _history[2*_initial_n-njets-1].max_dij_so_far;
}


//----------------------------------------------------------------------
// work through the object's history until
bool ClusterSequence::object_in_jet(const PseudoJet & object, 
                                    const PseudoJet & jet) const {

  // make sure the object conceivably belongs to this clustering
  // sequence
  assert(_potentially_valid(object) && _potentially_valid(jet));

  const PseudoJet * this_object = &object;
  const PseudoJet * childp;
  while(true) {
    if (this_object->cluster_hist_index() == jet.cluster_hist_index()) {
      return true;
    } else if (has_child(*this_object, childp)) {this_object = childp;}
    else {return false;}
  }
}

//----------------------------------------------------------------------
/// if the jet has parents in the clustering, it returns true
/// and sets parent1 and parent2 equal to them.
///
/// if it has no parents it returns false and sets parent1 and
/// parent2 to zero
bool ClusterSequence::has_parents(const PseudoJet & jet, PseudoJet & parent1, 
                              PseudoJet & parent2) const {

  const history_element & hist = _history[jet.cluster_hist_index()];

  // make sure we do not run into any unexpected situations --
  // i.e. both parents valid, or neither
  assert ((hist.parent1 >= 0 && hist.parent2 >= 0) || 
          (hist.parent1 < 0 && hist.parent2 < 0));

  if (hist.parent1 < 0) {
    parent1 = PseudoJet(0.0,0.0,0.0,0.0);
    parent2 = parent1;
    return false;
  } else {
    parent1 = _jets[_history[hist.parent1].jetp_index];
    parent2 = _jets[_history[hist.parent2].jetp_index];
    // order the parents in decreasing pt
    if (parent1.perp2() < parent2.perp2()) swap(parent1,parent2);
    return true;
  }
}

//----------------------------------------------------------------------
/// if the jet has a child then return true and give the child jet
/// otherwise return false and set the child to zero
bool ClusterSequence::has_child(const PseudoJet & jet, PseudoJet & child) const {

  //const history_element & hist = _history[jet.cluster_hist_index()];
  //
  //if (hist.child >= 0) {
  //  child = _jets[_history[hist.child].jetp_index];
  //  return true;
  //} else {
  //  child = PseudoJet(0.0,0.0,0.0,0.0);
  //  return false;
  //}
  const PseudoJet * childp;
  bool res = has_child(jet, childp);
  if (res) {
    child = *childp;
    return true;
  } else {
    child = PseudoJet(0.0,0.0,0.0,0.0);
    return false;
  }
}

bool ClusterSequence::has_child(const PseudoJet & jet, const PseudoJet * & childp) const {

  const history_element & hist = _history[jet.cluster_hist_index()];

  // check that this jet has a child and that the child corresponds to
  // a true jet [RETHINK-IF-CHANGE-NUMBERING: what is the right
  // behaviour if the child is the same jet but made inclusive...?]
  if (hist.child >= 0 && _history[hist.child].jetp_index >= 0) {
    childp = &(_jets[_history[hist.child].jetp_index]);
    return true;
  } else {
    childp = NULL;
    return false;
  }
}


//----------------------------------------------------------------------
/// if this jet has a child (and so a partner) return true
/// and give the partner, otherwise return false and set the
/// partner to zero
bool ClusterSequence::has_partner(const PseudoJet & jet, 
                              PseudoJet & partner) const {

  const history_element & hist = _history[jet.cluster_hist_index()];

  // make sure we have a child and that the child does not correspond
  // to a clustering with the beam (or some other invalid quantity)
  if (hist.child >= 0 && _history[hist.child].parent2 >= 0) {
    const history_element & child_hist = _history[hist.child];
    if (child_hist.parent1 == jet.cluster_hist_index()) {
      // partner will be child's parent2 -- for iB clustering
      // parent2 will not be valid
      partner = _jets[_history[child_hist.parent2].jetp_index];
    } else {
      // partner will be child's parent1
      partner = _jets[_history[child_hist.parent1].jetp_index];
    }
    return true;
  } else {
    partner = PseudoJet(0.0,0.0,0.0,0.0);
    return false;
  }
}


//----------------------------------------------------------------------
// return a vector of the particles that make up a jet
vector<PseudoJet> ClusterSequence::constituents (const PseudoJet & jet) const {
  vector<PseudoJet> subjets;
  add_constituents(jet, subjets);
  return subjets;
}

//----------------------------------------------------------------------
/// output the supplied vector of jets in a format that can be read
/// by an appropriate root script; the format is:
/// jet-n jet-px jet-py jet-pz jet-E 
///   particle-n particle-rap particle-phi particle-pt
///   particle-n particle-rap particle-phi particle-pt
///   ...
/// #END
/// ... [i.e. above repeated]
void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets, 
                                          ostream & ostr) const {
  for (unsigned i = 0; i < jets.size(); i++) {
    ostr << i  << " "
         << jets[i].px() << " "
         << jets[i].py() << " "
         << jets[i].pz() << " "
         << jets[i].E() << endl;
    vector<PseudoJet> cst = constituents(jets[i]);
    for (unsigned j = 0; j < cst.size() ; j++) {
      ostr << " " << j << " "
           << cst[j].rap() << " "
           << cst[j].phi() << " "
           << cst[j].perp() << endl;
    }
    ostr << "#END" << endl;
  }
}

// Not yet. Perhaps in a future release
// //----------------------------------------------------------------------
// // print out all inclusive jets with pt > ptmin
// void ClusterSequence::print_jets (const double & ptmin) const{
//     vector<PseudoJet> jets = sorted_by_pt(inclusive_jets(ptmin));
// 
//     for (size_t j = 0; j < jets.size(); j++) {
//        printf("%5u %7.3f %7.3f %9.3f\n",
//        j,jets[j].rap(),jets[j].phi(),jets[j].perp());
//     }
// }

//----------------------------------------------------------------------
/// returns a vector of size n_particles() which indicates, for 
/// each of the initial particles (in the order in which they were
/// supplied), which of the supplied jets it belongs to; if it does
/// not belong to any of the supplied jets, the index is set to -1;
vector<int> ClusterSequence::particle_jet_indices(
                        const vector<PseudoJet> & jets) const {

  vector<int> indices(n_particles());

  // first label all particles as not belonging to any jets
  for (unsigned ipart = 0; ipart < n_particles(); ipart++) 
    indices[ipart] = -1;

  // then for each of the jets relabel its consituents as belonging to
  // that jet
  for (unsigned ijet = 0; ijet < jets.size(); ijet++) {

    vector<PseudoJet> jet_constituents(constituents(jets[ijet]));

    for (unsigned ip = 0; ip < jet_constituents.size(); ip++) {
      // a safe (if slightly redundant) way of getting the particle
      // index (for initial particles it is actually safe to assume
      // ipart=iclust).
      unsigned iclust = jet_constituents[ip].cluster_hist_index();
      unsigned ipart = history()[iclust].jetp_index;
      indices[ipart] = ijet;
    }
  }

  return indices;
}


//----------------------------------------------------------------------
// recursive routine that adds on constituents of jet to the subjet_vector
void ClusterSequence::add_constituents (
           const PseudoJet & jet, vector<PseudoJet> & subjet_vector) const {
  // find out position in cluster history
  int i = jet.cluster_hist_index();
  int parent1 = _history[i].parent1;
  int parent2 = _history[i].parent2;

  if (parent1 == InexistentParent) {
    // It is an original particle (labelled by its parent having value
    // InexistentParent), therefore add it on to the subjet vector
    subjet_vector.push_back(jet);
    return;
  } 

  // add parent 1
  add_constituents(_jets[_history[parent1].jetp_index], subjet_vector);

  // see if parent2 is a real jet; if it is then add its constituents
  if (parent2 != BeamJet) {
    add_constituents(_jets[_history[parent2].jetp_index], subjet_vector);
  }
}



//----------------------------------------------------------------------
// initialise the history in a standard way
void ClusterSequence::_add_step_to_history (
               const int & step_number, const int & parent1, 
               const int & parent2, const int & jetp_index,
               const double & dij) {

  history_element element;
  element.parent1 = parent1;
  element.parent2 = parent2;
  element.jetp_index = jetp_index;
  element.child = Invalid;
  element.dij   = dij;
  element.max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far);
  _history.push_back(element);

  int local_step = _history.size()-1;
  assert(local_step == step_number);

  assert(parent1 >= 0);
  _history[parent1].child = local_step;
  if (parent2 >= 0) {_history[parent2].child = local_step;}

  // get cross-referencing right from PseudoJets
  if (jetp_index != Invalid) {
    assert(jetp_index >= 0);
    //cout << _jets.size() <<" "<<jetp_index<<"\n";
    _jets[jetp_index].set_cluster_hist_index(local_step);
  }

  if (_writeout_combinations) {
    cout << local_step << ": " 
         << parent1 << " with " << parent2
         << "; y = "<< dij<<endl;
  }

}




//======================================================================
// Return an order in which to read the history such that _history[order[i]] 
// will always correspond to the same set of consituent particles if 
// two branching histories are equivalent in terms of the particles
// contained in any given pseudojet.
vector<int> ClusterSequence::unique_history_order() const {

  // first construct an array that will tell us the lowest constituent
  // of a given jet -- this will always be one of the original
  // particles, whose order is well defined and so will help us to
  // follow the tree in a unique manner.
  valarray<int> lowest_constituent(_history.size());
  int hist_n = _history.size();
  lowest_constituent = hist_n; // give it a large number
  for (int i = 0; i < hist_n; i++) {
    // sets things up for the initial partons
    lowest_constituent[i] = min(lowest_constituent[i],i); 
    // propagates them through to the children of this parton
    if (_history[i].child > 0) lowest_constituent[_history[i].child] 
      = min(lowest_constituent[_history[i].child],lowest_constituent[i]);
  }

  // establish an array for what we have and have not extracted so far
  valarray<bool> extracted(_history.size()); extracted = false;
  vector<int> unique_tree;
  unique_tree.reserve(_history.size());

  // now work our way through the tree
  for (unsigned i = 0; i < n_particles(); i++) {
    if (!extracted[i]) {
      unique_tree.push_back(i);
      extracted[i] = true;
      _extract_tree_children(i, extracted, lowest_constituent, unique_tree);
    }
  }

  return unique_tree;
}

//======================================================================
// helper for unique_history_order
void ClusterSequence::_extract_tree_children(
       int position, 
       valarray<bool> & extracted, 
       const valarray<int> & lowest_constituent,
       vector<int> & unique_tree) const {
  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);
}


//======================================================================
// return the list of unclustered particles
vector<PseudoJet> ClusterSequence::unclustered_particles() const {
  vector<PseudoJet> unclustered;
  for (unsigned i = 0; i < n_particles() ; i++) {
    if (_history[i].child == Invalid) 
      unclustered.push_back(_jets[_history[i].jetp_index]);
  }
  return unclustered;
}



//======================================================================
// helper for unique_history_order
void ClusterSequence::_extract_tree_parents(
       int position, 
       valarray<bool> & extracted, 
       const valarray<int> & lowest_constituent,
       vector<int> & unique_tree) const {

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


//======================================================================
/// carries out the bookkeeping associated with the step of recombining
/// jet_i and jet_j (assuming a distance dij) and returns the index
/// of the recombined jet, newjet_k.
void ClusterSequence::_do_ij_recombination_step(
                               const int & jet_i, const int & jet_j, 
                               const double & dij, 
                               int & newjet_k) {

  // create the new jet by recombining the first two
  PseudoJet newjet;
  _jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet);
  _jets.push_back(newjet);
  // original version...
  //_jets.push_back(_jets[jet_i] + _jets[jet_j]);

  // get its index
  newjet_k = _jets.size()-1;

  // get history index
  int newstep_k = _history.size();
  // and provide jet with the info
  _jets[newjet_k].set_cluster_hist_index(newstep_k);

  // finally sort out the history 
  int hist_i = _jets[jet_i].cluster_hist_index();
  int hist_j = _jets[jet_j].cluster_hist_index();

  _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
                       newjet_k, dij);

}


//======================================================================
/// carries out the bookkeeping associated with the step of recombining
/// jet_i with the beam
void ClusterSequence::_do_iB_recombination_step(
                                  const int & jet_i, const double & diB) {
  // get history index
  int newstep_k = _history.size();

  // recombine the jet with the beam
  _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
                       Invalid, diB);

}

FASTJET_END_NAMESPACE