Index: external/fastjet/PseudoJet.cc
===================================================================
--- external/fastjet/PseudoJet.cc	(revision 4564cba9fa0180ca6c293da8cea545291a6d0116)
+++ external/fastjet/PseudoJet.cc	(revision d7d2da362406b9e37fdf238970417c94779b54b0)
@@ -1,6 +1,6 @@
-//FJSTARTHEADER
-// $Id: PseudoJet.cc 4442 2020-05-05 07:50:11Z soyez $
-//
-// Copyright (c) 2005-2020, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//STARTHEADER
+// $Id$
+//
+// 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
 
 
@@ -33,7 +31,5 @@
 #include "fastjet/PseudoJet.hh"
 #include "fastjet/ClusterSequence.hh"
-#ifndef __FJCORE__
 #include "fastjet/ClusterSequenceAreaBase.hh"
-#endif  // __FJCORE__
 #include "fastjet/CompositeJetStructure.hh"
 #include<valarray>
@@ -71,12 +67,4 @@
   _kt2 = this->px()*this->px() + this->py()*this->py();
   _phi = pseudojet_invalid_phi;
-  // strictly speaking, _rap does not need initialising, because
-  // it's never used as long as _phi == pseudojet_invalid_phi
-  // (and gets set when _phi is requested). However ATLAS
-  // 2013-03-28 complained that they sometimes have NaN's in
-  // _rap and this interferes with some of their internal validation. 
-  // So we initialise it; penalty is about 0.3ns per PseudoJet out of about
-  // 10ns total initialisation time (on a intel Core i7 2.7GHz)
-  _rap = pseudojet_invalid_rap;
 }
 
@@ -178,8 +166,4 @@
 // return the product, coeff * jet
 PseudoJet operator* (double coeff, const PseudoJet & jet) {
-  // see the comment in operator*= about ensuring valid rap phi
-  // before a multiplication to handle case of multiplication by
-  // zero, while maintaining rapidity and phi
-  jet._ensure_valid_rap_phi(); 
   //return PseudoJet(coeff*jet.four_mom());
   // the following code is hopefully more efficient
@@ -203,13 +187,5 @@
 //----------------------------------------------------------------------
 /// multiply the jet's momentum by the coefficient
-PseudoJet & PseudoJet::operator*=(double coeff) {
-  // operator*= aims to maintain the rapidity and azimuth
-  // for the PseudoJet; if they have already been evaluated
-  // this is fine, but if they haven't and coeff is sufficiently
-  // small as to cause a zero or underflow result, then a subsequent
-  // invocation of rap or phi will lead to a non-sensical result. 
-  // So, here, we preemptively ensure that rapidity and phi
-  // are correctly cached
-  _ensure_valid_rap_phi(); 
+void PseudoJet::operator*=(double coeff) {
   _px *= coeff;
   _py *= coeff;
@@ -218,12 +194,10 @@
   _kt2*= coeff*coeff;
   // phi and rap are unchanged
-  return *this;
 }
 
 //----------------------------------------------------------------------
 /// divide the jet's momentum by the coefficient
-PseudoJet & PseudoJet::operator/=(double coeff) {
+void PseudoJet::operator/=(double coeff) {
   (*this) *= 1.0/coeff;
-  return *this;
 }
 
@@ -231,5 +205,5 @@
 //----------------------------------------------------------------------
 /// add the other jet's momentum to this jet
-PseudoJet & PseudoJet::operator+=(const PseudoJet & other_jet) {
+void PseudoJet::operator+=(const PseudoJet & other_jet) {
   _px += other_jet._px;
   _py += other_jet._py;
@@ -237,5 +211,4 @@
   _E  += other_jet._E ;
   _finish_init(); // we need to recalculate phi,rap,kt2
-  return *this;
 }
 
@@ -243,5 +216,5 @@
 //----------------------------------------------------------------------
 /// subtract the other jet's momentum from this jet
-PseudoJet & PseudoJet::operator-=(const PseudoJet & other_jet) {
+void PseudoJet::operator-=(const PseudoJet & other_jet) {
   _px -= other_jet._px;
   _py -= other_jet._py;
@@ -249,5 +222,4 @@
   _E  -= other_jet._E ;
   _finish_init(); // we need to recalculate phi,rap,kt2
-  return *this;
 }
 
@@ -279,6 +251,6 @@
 
 //----------------------------------------------------------------------
-/// transform this jet (given in the rest frame of prest) into a jet
-/// in the lab frame 
+/// transform this jet (given in lab) into a jet in the rest
+/// frame of prest
 //
 // NB: code adapted from that in herwig f77 (checked how it worked
@@ -306,6 +278,6 @@
 
 //----------------------------------------------------------------------
-/// transform this jet (given in lab) into a jet in the rest
-/// frame of prest  
+/// transform this jet (given in the rest frame of prest) into a jet
+/// in the lab frame;
 //
 // NB: code adapted from that in herwig f77 (checked how it worked
@@ -413,9 +385,9 @@
 string PseudoJet::description() const{
   // the "default" case of a PJ which does not belong to any cluster sequence
-  if (!_structure)
+  if (!_structure())
     return "standard PseudoJet (with no associated clustering information)";
   
   // for all the other cases, the description comes from the structure
-  return _structure->description();
+  return _structure()->description();
 }
 
@@ -433,5 +405,5 @@
 // ClusterSequence
 bool PseudoJet::has_associated_cluster_sequence() const{
-  return (_structure) && (_structure->has_associated_cluster_sequence());
+  return (_structure()) && (_structure->has_associated_cluster_sequence());
 }
 
@@ -450,5 +422,5 @@
 // ClusterSequence that is still valid
 bool PseudoJet::has_valid_cluster_sequence() const{
-  return (_structure) && (_structure->has_valid_cluster_sequence());
+  return (_structure()) && (_structure->has_valid_cluster_sequence());
 }
 
@@ -471,7 +443,7 @@
 
 //----------------------------------------------------------------------
-// return true if there is some structure associated with this PseudoJet
+// return true if there is some strusture associated with this PseudoJet
 bool PseudoJet::has_structure() const{
-  return bool(_structure);
+  return _structure();
 }
 
@@ -482,6 +454,6 @@
 // return NULL if there is no associated structure
 const PseudoJetStructureBase* PseudoJet::structure_ptr() const {
-  //if (!_structure) return NULL;
-  return _structure.get();
+  if (!_structure()) return NULL;
+  return _structure();
 }
   
@@ -497,6 +469,6 @@
 // underlying structure.
 PseudoJetStructureBase* PseudoJet::structure_non_const_ptr(){
-  //if (!_structure) return NULL;
-  return _structure.get();
+  if (!_structure()) return NULL;
+  return _structure();
 }
   
@@ -507,7 +479,7 @@
 // throw an error if there is no associated structure
 const PseudoJetStructureBase* PseudoJet::validated_structure_ptr() const {
-  if (!_structure) 
+  if (!_structure()) 
     throw Error("Trying to access the structure of a PseudoJet which has no associated structure");
-  return _structure.get();
+  return _structure();
 }
   
@@ -577,5 +549,5 @@
 // returns true if the PseudoJet has constituents
 bool PseudoJet::has_constituents() const{
-  return (_structure) && (_structure->has_constituents());
+  return (_structure()) && (_structure->has_constituents());
 }
 
@@ -590,5 +562,5 @@
 // returns true if the PseudoJet has support for exclusive subjets
 bool PseudoJet::has_exclusive_subjets() const{
-  return (_structure) && (_structure->has_exclusive_subjets());
+  return (_structure()) && (_structure->has_exclusive_subjets());
 }
 
@@ -605,5 +577,5 @@
 // an Error is thrown if this PseudoJet has no currently valid
 // associated ClusterSequence
-std::vector<PseudoJet> PseudoJet::exclusive_subjets (const double dcut) const {
+std::vector<PseudoJet> PseudoJet::exclusive_subjets (const double & dcut) const {
   return validated_structure_ptr()->exclusive_subjets(*this, dcut);
 }
@@ -616,5 +588,5 @@
 // an Error is thrown if this PseudoJet has no currently valid
 // associated ClusterSequence
-int PseudoJet::n_exclusive_subjets(const double dcut) const {
+int PseudoJet::n_exclusive_subjets(const double & dcut) const {
   return validated_structure_ptr()->n_exclusive_subjets(*this, dcut);
 }
@@ -674,5 +646,5 @@
 // ClusterSequence have no pieces and this methos will return false.
 bool PseudoJet::has_pieces() const{
-  return ((_structure) && (_structure->has_pieces(*this)));
+  return ((_structure()) && (_structure->has_pieces(*this)));
 }
 
@@ -696,6 +668,4 @@
 //----------------------------------------------------------------------
 
-#ifndef __FJCORE__
-
 //----------------------------------------------------------------------
 // if possible, return a valid ClusterSequenceAreaBase pointer; otherwise
@@ -745,5 +715,4 @@
 }
 
-#endif  // __FJCORE__
 
 //----------------------------------------------------------------------
@@ -769,4 +738,33 @@
 }
 
+
+
+//----------------------------------------------------------------------
+/// given a vector of values with a one-to-one correspondence with the
+/// vector of objects, sort objects into an order such that the
+/// associated values would be in increasing order
+template<class T> vector<T>  objects_sorted_by_values(
+                       const vector<T> & objects, 
+		       const vector<double> & values) {
+
+  assert(objects.size() == values.size());
+
+  // get a vector of indices
+  vector<int> indices(values.size());
+  for (size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
+  
+  // sort the indices
+  sort_indices(indices, values);
+  
+  // copy the objects 
+  vector<T> objects_sorted(objects.size());
+  
+  // place the objects in the correct order
+  for (size_t i = 0; i < indices.size(); i++) {
+    objects_sorted[i] = objects[indices[i]];
+  }
+
+  return objects_sorted;
+}
 
 //----------------------------------------------------------------------
@@ -836,5 +834,4 @@
 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2){
   vector<PseudoJet> pieces;
-  pieces.reserve(2);
   pieces.push_back(j1);
   pieces.push_back(j2);
@@ -845,5 +842,4 @@
 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3){
   vector<PseudoJet> pieces;
-  pieces.reserve(3);
   pieces.push_back(j1);
   pieces.push_back(j2);
@@ -855,5 +851,4 @@
 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4){
   vector<PseudoJet> pieces;
-  pieces.reserve(4);
   pieces.push_back(j1);
   pieces.push_back(j2);
