Index: external/fastjet/ClusterSequence.cc
===================================================================
--- external/fastjet/ClusterSequence.cc	(revision 4564cba9fa0180ca6c293da8cea545291a6d0116)
+++ external/fastjet/ClusterSequence.cc	(revision d7d2da362406b9e37fdf238970417c94779b54b0)
@@ -1,6 +1,6 @@
-//FJSTARTHEADER
-// $Id: ClusterSequence.cc 4442 2020-05-05 07:50:11Z soyez $
+//STARTHEADER
+// $Id$
 //
-// Copyright (c) 2005-2020, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,9 +13,7 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development. They are described in the original FastJet paper,
-//  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
+//  development and are described in hep-ph/0512210. If you use
 //  FastJet as part of work towards a scientific publication, please
-//  quote the version you use and include a citation to the manual and
-//  optionally also to hep-ph/0512210.
+//  include a citation to the FastJet paper.
 //
 //  FastJet is distributed in the hope that it will be useful,
@@ -27,5 +25,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//FJENDHEADER
+//ENDHEADER
 
 #include "fastjet/Error.hh"
@@ -34,10 +32,4 @@
 #include "fastjet/ClusterSequenceStructure.hh"
 #include "fastjet/version.hh" // stores the current version number
-#include "fastjet/internal/LazyTiling9Alt.hh"
-#include "fastjet/internal/LazyTiling9.hh"
-#include "fastjet/internal/LazyTiling25.hh"
-#ifndef __FJCORE__
-#include "fastjet/internal/LazyTiling9SeparateGhosts.hh"
-#endif  // __FJCORE__
 #include<iostream>
 #include<sstream>
@@ -152,6 +144,6 @@
   // set the pointer in the wrapper to this object to NULL to say that
   // we're going out of scope
-  if (_structure_shared_ptr){
-    ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr.get()); 
+  if (_structure_shared_ptr()){
+    ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr()); 
     // normally the csi is purely internal so it really should not be
     // NULL i.e assert should be OK
@@ -195,5 +187,5 @@
 //DEP //----------------------------------------------------------------------
 //DEP void ClusterSequence::_initialise_and_run (
-//DEP 				  const double R,
+//DEP 				  const double & R,
 //DEP 				  const Strategy & strategy,
 //DEP 				  const bool & writeout_combinations) {
@@ -254,10 +246,10 @@
       // 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) );
+	// 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) );
+	_R2    = 2 * ( 1.0 - cos(_Rparam) );
       }
       _invR2 = 1.0/_R2;
@@ -281,10 +273,4 @@
   //             ./fastjet_timing_plugins -kt -nhardest 30 -repeat 50000 -strategy -3 -R 0.5 -nev 1  <  ../../data/Pythia-PtMin1000-LHC-1000ev.dat
   if (_strategy == Best) {
-    _strategy = _best_strategy();
-#ifdef DROP_CGAL
-    // fall back strategy for large N when CGAL is missing
-    if (_strategy == NlnN) _strategy = N2MHTLazy25;
-#endif  // DROP_CGAL
-  } else if (_strategy == BestFJ30) {
     int N = _jets.size();
     //if (N <= 55*max(0.5,min(1.0,_Rparam))) {// old empirical scaling with R
@@ -338,5 +324,5 @@
   // run the code containing the selected strategy
   // 
-  // We order the strategies starting from the ones used by the Best
+  // We order the strategies stqrting from the ones used by the Best
   // strategy in the order of increasing N, then the remaining ones
   // again in the order of increasing N.
@@ -348,40 +334,4 @@
   } else if (_strategy == N2MinHeapTiled) {
     this->_minheap_faster_tiled_N2_cluster();
-  } else if (_strategy == N2MHTLazy9Alt) {
-    // attempt to use an external tiling routine -- it manipulates
-    // the CS history via the plugin mechanism
-    _plugin_activated = true;
-    LazyTiling9Alt tiling(*this);
-    tiling.run();
-    _plugin_activated = false;
-
-  } else if (_strategy == N2MHTLazy25) {
-    // attempt to use an external tiling routine -- it manipulates
-    // the CS history via the plugin mechanism
-    _plugin_activated = true;
-    LazyTiling25 tiling(*this);
-    tiling.run();
-    _plugin_activated = false;
-
-  } else if (_strategy == N2MHTLazy9) {
-    // attempt to use an external tiling routine -- it manipulates
-    // the CS history via the plugin mechanism
-    _plugin_activated = true;
-    LazyTiling9 tiling(*this);
-    tiling.run();
-    _plugin_activated = false;
-
-  } else if (_strategy == N2MHTLazy9AntiKtSeparateGhosts) {
-#ifndef __FJCORE__
-    // attempt to use an external tiling routine -- it manipulates
-    // the CS history via the plugin mechanism
-    _plugin_activated = true;
-    LazyTiling9SeparateGhosts tiling(*this);
-    tiling.run();
-    _plugin_activated = false;
-#else 
-    throw Error("N2MHTLazy9AntiKtSeparateGhosts strategy not supported with FJCORE");
-#endif  // __FJCORE__
-
   } else if (_strategy == NlnN) {
     this->_delaunay_cluster();
@@ -409,5 +359,5 @@
 // these needs to be defined outside the class definition.
 bool ClusterSequence::_first_time = true;
-LimitedWarning ClusterSequence::_exclusive_warnings;
+int ClusterSequence::_n_exclusive_warnings = 0;
 
 
@@ -435,9 +385,9 @@
   (*ostr) << "#     A software package for jet finding and analysis at colliders      \n";
   (*ostr) << "#                           http://fastjet.fr                           \n"; 
-  (*ostr) << "#	                                                                      \n";
+  (*ostr) << "#                                                                       \n";
   (*ostr) << "# Please cite EPJC72(2012)1896 [arXiv:1111.6097] if you use this package\n";
   (*ostr) << "# for scientific work and optionally PLB641(2006)57 [hep-ph/0512210].   \n";
-  (*ostr) << "#                                                                       \n";
-  (*ostr) << "# FastJet is provided without warranty under the GNU GPL v2 or higher.  \n";
+  (*ostr) << "#								      	   \n";
+  (*ostr) << "# FastJet is provided without warranty under the terms of the GNU GPLv2.\n";
   (*ostr) << "# It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code";
 #ifndef DROP_CGAL
@@ -523,4 +473,6 @@
 
 //----------------------------------------------------------------------
+// Return the component corresponding to the specified index.
+// taken from CLHEP
 string ClusterSequence::strategy_string (Strategy strategy_in)  const {
   string strategy;
@@ -540,12 +492,4 @@
   case N2PoorTiled:
     strategy = "N2PoorTiled"; break;
-  case N2MHTLazy9:
-    strategy = "N2MHTLazy9"; break;
-  case N2MHTLazy9Alt:
-    strategy = "N2MHTLazy9Alt"; break;
-  case N2MHTLazy25:
-    strategy = "N2MHTLazy25"; break;
-  case N2MHTLazy9AntiKtSeparateGhosts:
-    strategy = "N2MHTLazy9AntiKtSeparateGhosts"; break;
   case N3Dumb:
     strategy = "N3Dumb"; break;
@@ -584,151 +528,4 @@
     } else {return 1.0;}
   } else {throw Error("Unrecognised jet algorithm");}
-}
-
-//----------------------------------------------------------------------
-// returns a suggestion for the best strategy to use on event
-// multiplicity, algorithm, R, etc.
-//
-// Some of the work to establish the best strategy is collected in
-// issue-tracker/2014-07-auto-strategy-selection;
-// transition_fit_v2.fit indicates the results of the fits that we're
-// using here. (Automatically generated by transition_fit_v2.gp).
-//
-// The transition to NlnN is always present, and it is the the
-// caller's responsibility to drop back down to N2MHTLazy25 if NlnN
-// isn't available.
-//
-// This routine should be called only if the jet alg is one of kt,
-// antikt, cam or genkt.
-Strategy ClusterSequence::_best_strategy() const {
-  int N = _jets.size();
-  // define bounded R, always above 0.1, because we don't trust any
-  // of our parametrizations below R = 0.1
-  double bounded_R = max(_Rparam, 0.1);
-
-  // the very first test thing is a quick hard-coded test to decide
-  // if we immediately opt for N2Plain
-  if (N <= 30 || N <= 39.0/(bounded_R + 0.6)) {
-    return N2Plain;
-  } 
-  
-  // Define objects that describe our various boundaries. A prefix N_
-  // indicates that boundary is for N, while L_ means it's for log(N).
-  //
-  // Hopefully having them static will ensure minimal overhead
-  // in creating them; collecting them in one place should
-  // help with updates?
-  //
-  const static _Parabola N_Tiled_to_MHT_lowR             (-45.4947,54.3528,44.6283);
-  const static _Parabola L_MHT_to_MHTLazy9_lowR          (0.677807,-1.05006,10.6994);
-  const static _Parabola L_MHTLazy9_to_MHTLazy25_akt_lowR(0.169967,-0.512589,12.1572);
-  const static _Parabola L_MHTLazy9_to_MHTLazy25_kt_lowR (0.16237,-0.484612,12.3373);
-  const static _Parabola L_MHTLazy9_to_MHTLazy25_cam_lowR = L_MHTLazy9_to_MHTLazy25_kt_lowR;
-  const static _Parabola L_MHTLazy25_to_NlnN_akt_lowR    (0.0472051,-0.22043,15.9196);
-  const static _Parabola L_MHTLazy25_to_NlnN_kt_lowR     (0.118609,-0.326811,14.8287);
-  const static _Parabola L_MHTLazy25_to_NlnN_cam_lowR    (0.10119,-0.295748,14.3924);
-
-  const static _Line     L_Tiled_to_MHTLazy9_medR         (-1.31304,7.29621);
-  const static _Parabola L_MHTLazy9_to_MHTLazy25_akt_medR = L_MHTLazy9_to_MHTLazy25_akt_lowR;
-  const static _Parabola L_MHTLazy9_to_MHTLazy25_kt_medR  = L_MHTLazy9_to_MHTLazy25_kt_lowR;
-  const static _Parabola L_MHTLazy9_to_MHTLazy25_cam_medR = L_MHTLazy9_to_MHTLazy25_cam_lowR;
-  const static _Parabola L_MHTLazy25_to_NlnN_akt_medR     = L_MHTLazy25_to_NlnN_akt_lowR;
-  const static _Parabola L_MHTLazy25_to_NlnN_kt_medR      = L_MHTLazy25_to_NlnN_kt_lowR;
-  const static _Parabola L_MHTLazy25_to_NlnN_cam_medR     = L_MHTLazy25_to_NlnN_cam_lowR;
-
-  const static double    N_Plain_to_MHTLazy9_largeR         = 75;
-  const static double    N_MHTLazy9_to_MHTLazy25_akt_largeR = 700;
-  const static double    N_MHTLazy9_to_MHTLazy25_kt_largeR  = 1000;
-  const static double    N_MHTLazy9_to_MHTLazy25_cam_largeR = 1000;
-  const static double    N_MHTLazy25_to_NlnN_akt_largeR     = 100000;
-  const static double    N_MHTLazy25_to_NlnN_kt_largeR      = 40000;
-  const static double    N_MHTLazy25_to_NlnN_cam_largeR     = 15000;
-
-  // We have timing studies only for kt, cam and antikt; for other
-  // algorithms we set the local jet_algorithm variable to the one of
-  // kt,cam,antikt that we think will be closest in behaviour to the
-  // other alg.
-  JetAlgorithm jet_algorithm;
-  if (_jet_algorithm == genkt_algorithm) {
-    // for genkt, then we set the local jet_algorithm variable (used
-    // only for strategy choice) to be either kt or antikt, depending on
-    // the p value.
-    double p   = jet_def().extra_param();
-    if (p < 0.0) jet_algorithm = antikt_algorithm;
-    else         jet_algorithm =     kt_algorithm;
-  } else if (_jet_algorithm == cambridge_for_passive_algorithm) {
-    // we assume (but haven't tested) that using the kt-alg timing
-    // transitions should be adequate for cambridge_for_passive_algorithm
-    jet_algorithm = kt_algorithm;
-  } else {
-    jet_algorithm = _jet_algorithm;
-  }
-
-  if (bounded_R < 0.65) {
-    // low R case
-    if          (N    < N_Tiled_to_MHT_lowR(bounded_R))              return N2Tiled;
-    double logN = log(double(N));
-    if          (logN < L_MHT_to_MHTLazy9_lowR(bounded_R))           return N2MinHeapTiled;
-    else {
-      if (jet_algorithm == antikt_algorithm){
-        if      (logN < L_MHTLazy9_to_MHTLazy25_akt_lowR(bounded_R)) return N2MHTLazy9;
-        else if (logN < L_MHTLazy25_to_NlnN_akt_lowR(bounded_R))     return N2MHTLazy25;
-        else                                                         return NlnN;
-      } else if (jet_algorithm == kt_algorithm){
-        if      (logN < L_MHTLazy9_to_MHTLazy25_kt_lowR(bounded_R))  return N2MHTLazy9;
-        else if (logN < L_MHTLazy25_to_NlnN_kt_lowR(bounded_R))      return N2MHTLazy25;
-        else                                                         return NlnN;
-      } else if (jet_algorithm == cambridge_algorithm)  {
-        if      (logN < L_MHTLazy9_to_MHTLazy25_cam_lowR(bounded_R)) return N2MHTLazy9;
-        else if (logN < L_MHTLazy25_to_NlnN_cam_lowR(bounded_R))     return N2MHTLazy25;
-        else                                                         return NlnNCam;
-      }
-    }
-  } else if (bounded_R < 0.5*pi) {
-    // medium R case
-    double logN = log(double(N));
-    if      (logN < L_Tiled_to_MHTLazy9_medR(bounded_R))             return N2Tiled;
-    else {
-      if (jet_algorithm == antikt_algorithm){
-        if      (logN < L_MHTLazy9_to_MHTLazy25_akt_medR(bounded_R)) return N2MHTLazy9;
-        else if (logN < L_MHTLazy25_to_NlnN_akt_medR(bounded_R))     return N2MHTLazy25;
-        else                                                         return NlnN;
-      } else if (jet_algorithm == kt_algorithm){
-        if      (logN < L_MHTLazy9_to_MHTLazy25_kt_medR(bounded_R))  return N2MHTLazy9;
-        else if (logN < L_MHTLazy25_to_NlnN_kt_medR(bounded_R))      return N2MHTLazy25;
-        else                                                         return NlnN;
-      } else if (jet_algorithm == cambridge_algorithm)  {
-        if      (logN < L_MHTLazy9_to_MHTLazy25_cam_medR(bounded_R)) return N2MHTLazy9;
-        else if (logN < L_MHTLazy25_to_NlnN_cam_medR(bounded_R))     return N2MHTLazy25;
-        else                                                         return NlnNCam;
-      }
-    }
-  } else {
-    // large R case (R > pi/2)
-    if      (N    < N_Plain_to_MHTLazy9_largeR)                      return N2Plain;
-    else {
-      if (jet_algorithm == antikt_algorithm){
-        if      (N < N_MHTLazy9_to_MHTLazy25_akt_largeR)             return N2MHTLazy9;
-        else if (N < N_MHTLazy25_to_NlnN_akt_largeR)                 return N2MHTLazy25;
-        else                                                         return NlnN;
-      } else if (jet_algorithm == kt_algorithm){
-        if      (N < N_MHTLazy9_to_MHTLazy25_kt_largeR)              return N2MHTLazy9;
-        else if (N < N_MHTLazy25_to_NlnN_kt_largeR)                  return N2MHTLazy25;
-        else                                                         return NlnN;
-      } else if (jet_algorithm == cambridge_algorithm)  {
-        if      (N < N_MHTLazy9_to_MHTLazy25_cam_largeR)             return N2MHTLazy9;
-        else if (N < N_MHTLazy25_to_NlnN_cam_largeR)                 return N2MHTLazy25;
-        else                                                         return NlnNCam;
-      }
-    }
-  }
-  
-  //bool code_should_never_reach_here = false;
-  //assert(code_should_never_reach_here);
-
-  assert(0 && "Code should never reach here");
-
-  return N2MHTLazy9;
-
 }
 
@@ -778,13 +575,4 @@
 
 
-ClusterSequence & ClusterSequence::operator=(const ClusterSequence & cs) {
-  // self assignment is trivial
-  if (&cs != this) {
-    _deletes_self_when_unused = false;
-    transfer_from_sequence(cs);
-  }
-  return *this;
-}
-
 //----------------------------------------------------------------------
 // transfer the sequence contained in other_seq into our own;
@@ -827,5 +615,5 @@
 
   // clean up existing structure
-  if (_structure_shared_ptr) {
+  if (_structure_shared_ptr()) {
     // If there are jets associated with an old version of the CS and
     // a new one, keeping track of when to delete the CS becomes more
@@ -835,5 +623,5 @@
     // anything that is currently associated with the cluster sequence
     // should be told that its cluster sequence no longer exists
-    ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr.get()); 
+    ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr()); 
     assert(csi != NULL);
     csi->set_associated_cs(NULL);
@@ -874,5 +662,5 @@
 //----------------------------------------------------------------------
 // return all inclusive jets with pt > ptmin
-vector<PseudoJet> ClusterSequence::inclusive_jets (const double ptmin) const{
+vector<PseudoJet> ClusterSequence::inclusive_jets (const double & ptmin) const{
   double dcut = ptmin*ptmin;
   int i = _history.size() - 1; // last jet
@@ -926,5 +714,5 @@
 // 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 {
+int ClusterSequence::n_exclusive_jets (const double & dcut) const {
 
   // first locate the point where clustering would have stopped (i.e. the
@@ -945,5 +733,5 @@
 // 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 {
+vector<PseudoJet> ClusterSequence::exclusive_jets (const double & dcut) const {
   int njets = n_exclusive_jets(dcut);
   return exclusive_jets(njets);
@@ -954,5 +742,5 @@
 // return the jets obtained by clustering the event to n jets.
 // Throw an error if there are fewer than n particles.
-vector<PseudoJet> ClusterSequence::exclusive_jets (const int njets) const {
+vector<PseudoJet> ClusterSequence::exclusive_jets (const int & njets) const {
 
   // make sure the user does not ask for more than jets than there
@@ -971,9 +759,9 @@
 // return the jets obtained by clustering the event to n jets.
 // If there are fewer than n particles, simply return all particles
-vector<PseudoJet> ClusterSequence::exclusive_jets_up_to (const int njets) const {
+vector<PseudoJet> ClusterSequence::exclusive_jets_up_to (const int & njets) const {
 
   // 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/Aachen, 
+  // 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()
@@ -985,6 +773,8 @@
        (_jet_def.extra_param() <0)) &&
       ((_jet_def.jet_algorithm() != plugin_algorithm) ||
-       (!_jet_def.plugin()->exclusive_sequence_meaningful()))) {
-    _exclusive_warnings.warn("dcut and exclusive jets for jet-finders other than kt, C/A or genkt with p>=0 should be interpreted with care.");
+       (!_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;
   }
 
@@ -1039,5 +829,5 @@
 /// return the dmin corresponding to the recombination that went from
 /// n+1 to n jets
-double ClusterSequence::exclusive_dmerge (const int njets) const {
+double ClusterSequence::exclusive_dmerge (const int & njets) const {
   assert(njets >= 0);
   if (njets >= _initial_n) {return 0.0;}
@@ -1051,5 +841,5 @@
 /// exclusive_dmerge, except in cases where the dmin do not increase
 /// monotonically.
-double ClusterSequence::exclusive_dmerge_max (const int njets) const {
+double ClusterSequence::exclusive_dmerge_max (const int & njets) const {
   assert(njets >= 0);
   if (njets >= _initial_n) {return 0.0;}
@@ -1063,5 +853,5 @@
 /// the algorithm with the given dcut.
 std::vector<PseudoJet> ClusterSequence::exclusive_subjets 
-   (const PseudoJet & jet, const double dcut) const {
+   (const PseudoJet & jet, const double & dcut) const {
 
   set<const history_element*> subhist;
@@ -1086,5 +876,5 @@
 /// exclusive_subjets.size()
 int ClusterSequence::n_exclusive_subjets(const PseudoJet & jet, 
-                        const double dcut) const {
+                        const double & dcut) const {
   set<const history_element*> subhist;
   // get the set of history elements that correspond to subjets at
@@ -1391,5 +1181,5 @@
 // //----------------------------------------------------------------------
 // // print out all inclusive jets with pt > ptmin
-// void ClusterSequence::print_jets (const double ptmin) const{
+// void ClusterSequence::print_jets (const double & ptmin) const{
 //     vector<PseudoJet> jets = sorted_by_pt(inclusive_jets(ptmin));
 // 
@@ -1467,8 +1257,7 @@
 // initialise the history in a standard way
 void ClusterSequence::_add_step_to_history (
-               //NO_LONGER_USED: const int step_number,
-               const int parent1, 
-	       const int parent2, const int jetp_index,
-	       const double dij) {
+	       const int & step_number, const int & parent1, 
+	       const int & parent2, const int & jetp_index,
+	       const double & dij) {
 
   history_element element;
@@ -1482,26 +1271,9 @@
 
   int local_step = _history.size()-1;
-  //#ifndef __NO_ASSERTS__
-  //assert(local_step == step_number);
-  //#endif
-
-  // sanity check: make sure the particles have not already been recombined
-  //
-  // Note that good practice would make this an assert (since this is
-  // a serious internal issue). However, we decided to throw an
-  // InternalError so that the end user can decide to catch it and
-  // retry the clustering with a different strategy.
+  assert(local_step == step_number);
 
   assert(parent1 >= 0);
-  if (_history[parent1].child != Invalid){
-    throw InternalError("trying to recomine an object that has previsously been recombined");
-  }
   _history[parent1].child = local_step;
-  if (parent2 >= 0) {
-    if (_history[parent2].child != Invalid){
-      throw InternalError("trying to recomine an object that has previsously been recombined");
-    }
-    _history[parent2].child = local_step;
-  }
+  if (parent2 >= 0) {_history[parent2].child = local_step;}
 
   // get cross-referencing right from PseudoJets
@@ -1657,6 +1429,6 @@
 /// of the recombined jet, newjet_k.
 void ClusterSequence::_do_ij_recombination_step(
-                               const int jet_i, const int jet_j, 
-			       const double dij, 
+                               const int & jet_i, const int & jet_j, 
+			       const double & dij, 
 			       int & newjet_k) {
 
@@ -1684,10 +1456,6 @@
   int hist_j = _jets[jet_j].cluster_hist_index();
 
-  _add_step_to_history(min(hist_i, hist_j), max(hist_i,hist_j),
+  _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
 		       newjet_k, dij);
-
-  //  _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
-  //		       newjet_k, dij);
-
 
 }
@@ -1698,14 +1466,11 @@
 /// jet_i with the beam
 void ClusterSequence::_do_iB_recombination_step(
-				  const int jet_i, const double diB) {
+				  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(_jets[jet_i].cluster_hist_index(),BeamJet,
+  _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
 		       Invalid, diB);
-
-  // // get history index
-  // int newstep_k = _history.size();
-  // 
-  // _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
-  //         	       Invalid, diB);
 
 }
