ClusterSequence.cc

Go to the documentation of this file.

//STARTHEADER
// $Id: ClusterSequence.cc 3162 2013-07-17 14:40:02Z salam $
//
// Copyright (c) 2005-2009, Matteo Cacciari, Gavin Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// 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 "fastjet/Error.hh"
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/version.hh" // stores the current version number
#include<iostream>
#include<sstream>
#include<fstream>
#include<cmath>
#include<cstdlib>
#include<cassert>
#include<string>
#include<set>

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

using namespace std;

JetAlgorithm ClusterSequence::_default_jet_algorithm = kt_algorithm;
//


// destructor that does nothing
ClusterSequence::~ClusterSequence () {}

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

  // don't run anything if the event is empty
  if (n_particles() == 0) return;

  // ----- deal with special cases: plugins & e+e- ------
  if (_jet_algorithm == plugin_algorithm) {
    // allows plugin_xyz() functions to modify cluster sequence
    _plugin_activated = true;
    // let the plugin do its work here
    _jet_def.plugin()->run_clustering( (*this) );
    _plugin_activated = false;
    return;
  } else if (_jet_algorithm == ee_kt_algorithm ||
             _jet_algorithm == ee_genkt_algorithm) {
    // ignore requested strategy
    _strategy = N2Plain;
    if (_jet_algorithm == ee_kt_algorithm) {
      // make sure that R is large enough so that "beam" recomb only
      // occurs when a single particle is left
      // Normally, this should be automatically set to 4 from JetDefinition
      assert(_Rparam > 2.0); 
      // this is used to renormalise the dij to get a "standard" form
      // and our convention in e+e- will be different from that
      // in long.inv case; NB: _invR2 name should be changed -> _renorm_dij?
      _invR2 = 1.0;
    } else {
      // as of 2009-01-09, choose R to be an angular distance, in
      // radians.  Since the algorithm uses 2(1-cos(theta)) as its
      // squared angular measure, make sure that the _R2 is defined
      // in a similar way.
      if (_Rparam > pi) {
        // choose a value that ensures that back-to-back particles will
        // always recombine 
        //_R2 = 4.0000000000001;
        _R2 = 2 * ( 3.0 + cos(_Rparam) );
      } else {
        _R2    = 2 * ( 1.0 - cos(_Rparam) );
      }
      _invR2 = 1.0/_R2;
    }
    //_simple_N2_cluster<EEBriefJet>();
    _simple_N2_cluster_EEBriefJet();
    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) && jet_def.jet_algorithm() != antikt_algorithm)
        || N > 35000/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) {
    // BriefJet provides standard long.invariant kt alg.
    //this->_simple_N2_cluster<BriefJet>();
    this->_simple_N2_cluster_BriefJet();
  } 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.fastjet.fr                         \n"; 
  cout << "#                                                                       \n";
  cout << "# Longitudinally invariant Kt, anti-Kt, and inclusive Cambridge/Aachen  \n";
  cout << "# clustering using fast geometric algorithms, with jet areas and optional\n";
  cout << "# external jet-finder plugins. If you use this code towards a scientific \n";
  cout << "# publication please cite EPJC72(2012)1896 [arXiv:1111.6097] and        \n";
  cout << "# optionally PLB641(2006)57 [hep-ph/0512210].                           \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";
  // make sure we really have the output done.
  cout.flush();
}

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

  _Qtot = 0;

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

    // determine the total energy in the event
    _Qtot += _jets[i].E();
  }
  _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 == genkt_algorithm) {
    double kt2 = jet.kt2();
    double p   = jet_def().extra_param();
    if (p <= 0 && kt2 < 1e-300) kt2 = 1e-300; // dodgy safety check
    return pow(kt2, p);
  } 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 algorithm");}
}


//----------------------------------------------------------------------
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 == ee_kt_algorithm
             || _jet_algorithm == antikt_algorithm
             || _jet_algorithm == genkt_algorithm
             || _jet_algorithm == ee_genkt_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("cs::inclusive_jets(...): 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 for algorithms 
  // that does not support it explicitly.
  // Native algorithm that support it are: kt, ee_kt, cambridge, 
  //   genkt and ee_genkt (both with p>=0)
  // For plugins, we check Plugin::exclusive_sequence_meaningful()
  if (( _jet_def.jet_algorithm() != kt_algorithm) &&
      ( _jet_def.jet_algorithm() != cambridge_algorithm) &&
      ( _jet_def.jet_algorithm() != ee_kt_algorithm) &&
      (((_jet_def.jet_algorithm() != genkt_algorithm) && 
        (_jet_def.jet_algorithm() != ee_genkt_algorithm)) || 
       (_jet_def.extra_param() <0)) &&
      ((_jet_def.jet_algorithm() != plugin_algorithm) ||
       (!_jet_def.plugin()->exclusive_sequence_meaningful())) &&
      (_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;
}

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


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


//----------------------------------------------------------------------
std::vector<PseudoJet> ClusterSequence::exclusive_subjets 
   (const PseudoJet & jet, const double & dcut) const {

  set<const history_element*> subhist;

  // get the set of history elements that correspond to subjets at
  // scale dcut
  get_subhist_set(subhist, jet, dcut, 0);

  // now transfer this into a sequence of jets
  vector<PseudoJet> subjets;
  subjets.reserve(subhist.size());
  for (set<const history_element*>::iterator elem = subhist.begin(); 
       elem != subhist.end(); elem++) {
    subjets.push_back(_jets[(*elem)->jetp_index]);
  }
  return subjets;
}

//----------------------------------------------------------------------
int ClusterSequence::n_exclusive_subjets(const PseudoJet & jet, 
                        const double & dcut) const {
  set<const history_element*> subhist;
  // get the set of history elements that correspond to subjets at
  // scale dcut
  get_subhist_set(subhist, jet, dcut, 0);
  return subhist.size();
}

//----------------------------------------------------------------------
std::vector<PseudoJet> ClusterSequence::exclusive_subjets 
   (const PseudoJet & jet, int n) const {

  set<const history_element*> subhist;

  // get the set of history elements that correspond to subjets at
  // scale dcut
  get_subhist_set(subhist, jet, -1.0, n);

  // now transfer this into a sequence of jets
  vector<PseudoJet> subjets;
  subjets.reserve(subhist.size());
  for (set<const history_element*>::iterator elem = subhist.begin(); 
       elem != subhist.end(); elem++) {
    subjets.push_back(_jets[(*elem)->jetp_index]);
  }
  return subjets;
}


//----------------------------------------------------------------------
double ClusterSequence::exclusive_subdmerge(const PseudoJet & jet, int nsub) const {
  set<const history_element*> subhist;

  // get the set of history elements that correspond to subjets at
  // scale dcut
  get_subhist_set(subhist, jet, -1.0, nsub);
  
  set<const history_element*>::iterator highest = subhist.end();
  highest--;
  return (*highest)->dij;
}


//----------------------------------------------------------------------
double ClusterSequence::exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const {

  set<const history_element*> subhist;

  // get the set of history elements that correspond to subjets at
  // scale dcut
  get_subhist_set(subhist, jet, -1.0, nsub);
  
  set<const history_element*>::iterator highest = subhist.end();
  highest--;
  return (*highest)->max_dij_so_far;
}



//----------------------------------------------------------------------
void ClusterSequence::get_subhist_set(set<const history_element*> & subhist,
                                     const  PseudoJet & jet, 
                                     double dcut, int maxjet) const {
  subhist.clear();
  subhist.insert(&(_history[jet.cluster_hist_index()]));

  // establish the set of jets that are relevant
  int njet = 1;
  while (true) {
    // first find out if we need to probe deeper into jet.
    // Get history element closest to end of sequence
    set<const history_element*>::iterator highest = subhist.end();
    assert (highest != subhist.begin()); 
    highest--;
    const history_element* elem = *highest;
    // make sure we haven't got too many jets
    if (njet == maxjet) break;
    // make sure it has parents
    if (elem->parent1 < 0)            break;
    // make sure that we still resolve it at scale dcut
    if (elem->max_dij_so_far <= dcut) break;

    // then do so: replace "highest" with its two parents
    subhist.erase(highest);
    subhist.insert(&(_history[elem->parent1]));
    subhist.insert(&(_history[elem->parent2]));
    njet++;
  }
}

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

//----------------------------------------------------------------------
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()) std::swap(parent1,parent2);
    return true;
  }
}

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


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

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

void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets, 
                                          const std::string & filename,
                                          const std::string & comment ) const {
  std::ofstream ostr(filename.c_str());
  if (comment != "") ostr << "# " << comment << endl;
  print_jets_for_root(jets, ostr);
}


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

//----------------------------------------------------------------------
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
    // Note: we add the initial particle and not simply 'jet' so that
    //       calling add_constituents with a subtracted jet containing
    //       only one particle will work.
    subjet_vector.push_back(_jets[i]);
    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]) 
        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 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);

}


//======================================================================
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