No OneTemporary
Actions

Size

97 KB

Subscribers

None

View Options

	diff --git a/include/Rivet/Makefile.am b/include/Rivet/Makefile.am
	--- a/include/Rivet/Makefile.am
	+++ b/include/Rivet/Makefile.am
	@@ -1,187 +1,188 @@
	## Internal headers - not to be installed
	nobase_dist_noinst_HEADERS =
	## Public headers - to be installed
	nobase_pkginclude_HEADERS =


	## Rivet interface
	nobase_pkginclude_HEADERS += \
	Rivet.hh \
	Run.hh \
	Event.hh \
	ParticleBase.hh \
	Particle.fhh Particle.hh \
	Jet.fhh Jet.hh \
	Projection.fhh Projection.hh \
	ProjectionApplier.hh \
	ProjectionHandler.hh \
	Analysis.hh \
	AnalysisHandler.hh \
	AnalysisInfo.hh \
	AnalysisBuilder.hh \
	AnalysisLoader.hh


	## Build config stuff
	nobase_pkginclude_HEADERS += \
	Config/RivetCommon.hh \
	Config/RivetConfig.hh \
	Config/BuildOptions.hh


	## Projections
	nobase_pkginclude_HEADERS += \
	Projections/AliceCommon.hh \
	Projections/AxesDefinition.hh \
	Projections/Beam.hh \
	Projections/BeamThrust.hh \
	Projections/CentralEtHCM.hh \
	Projections/CentralityProjection.hh \
	Projections/ChargedFinalState.hh \
	Projections/ChargedLeptons.hh \
	Projections/ConstLossyFinalState.hh \
	Projections/DirectFinalState.hh \
	Projections/DISFinalState.hh \
	Projections/DISKinematics.hh \
	Projections/DISLepton.hh \
	Projections/DressedLeptons.hh \
	Projections/FastJets.hh \
	Projections/PxConePlugin.hh \
	+ Projections/pxcone.h \
	Projections/FinalPartons.hh \
	Projections/FinalState.hh \
	Projections/FoxWolframMoments.hh \
	Projections/FParameter.hh \
	Projections/GeneratedPercentileProjection.hh \
	Projections/HadronicFinalState.hh \
	Projections/HeavyHadrons.hh \
	Projections/Hemispheres.hh \
	Projections/IdentifiedFinalState.hh \
	Projections/ImpactParameterProjection.hh \
	Projections/IndirectFinalState.hh \
	Projections/InitialQuarks.hh \
	Projections/InvMassFinalState.hh \
	Projections/JetAlg.hh \
	Projections/JetShape.hh \
	Projections/LeadingParticlesFinalState.hh \
	Projections/LossyFinalState.hh \
	Projections/MergedFinalState.hh \
	Projections/MissingMomentum.hh \
	Projections/MixedFinalState.hh \
	Projections/NeutralFinalState.hh \
	Projections/NonHadronicFinalState.hh \
	Projections/NonPromptFinalState.hh \
	Projections/ParisiTensor.hh \
	Projections/ParticleFinder.hh \
	Projections/PartonicTops.hh \
	Projections/PercentileProjection.hh \
	Projections/PrimaryHadrons.hh \
	Projections/PrimaryParticles.hh \
	Projections/PromptFinalState.hh \
	Projections/SingleValueProjection.hh \
	Projections/SmearedParticles.hh \
	Projections/SmearedJets.hh \
	Projections/SmearedMET.hh \
	Projections/Sphericity.hh \
	Projections/Spherocity.hh \
	Projections/TauFinder.hh \
	Projections/Thrust.hh \
	Projections/TriggerCDFRun0Run1.hh \
	Projections/TriggerCDFRun2.hh \
	Projections/TriggerProjection.hh \
	Projections/TriggerUA5.hh \
	Projections/UnstableFinalState.hh \
	Projections/UserCentEstimate.hh \
	Projections/VetoedFinalState.hh \
	Projections/VisibleFinalState.hh \
	Projections/WFinder.hh \
	Projections/ZFinder.hh
	## Meta-projection convenience headers
	nobase_pkginclude_HEADERS += \
	Projections/FinalStates.hh \
	Projections/Smearing.hh


	## Analysis base class headers
	# TODO: Move to Rivet/AnalysisTools header dir?
	nobase_pkginclude_HEADERS += \
	Analyses/MC_Cent_pPb.hh \
	Analyses/MC_ParticleAnalysis.hh \
	Analyses/MC_JetAnalysis.hh \
	Analyses/MC_JetSplittings.hh


	## Tools
	nobase_pkginclude_HEADERS += \
	Tools/AliceCommon.hh \
	Tools/AtlasCommon.hh \
	Tools/BeamConstraint.hh \
	Tools/BinnedHistogram.hh \
	Tools/CentralityBinner.hh \
	Tools/Cmp.fhh \
	Tools/Cmp.hh \
	Tools/Correlators.hh \
	Tools/Cutflow.hh \
	Tools/Cuts.fhh \
	Tools/Cuts.hh \
	Tools/Exceptions.hh \
	Tools/JetUtils.hh \
	Tools/Logging.hh \
	Tools/Random.hh \
	Tools/ParticleBaseUtils.hh \
	Tools/ParticleIdUtils.hh \
	Tools/ParticleUtils.hh \
	Tools/ParticleName.hh \
	Tools/Percentile.hh \
	Tools/PrettyPrint.hh \
	Tools/RivetPaths.hh \
	Tools/RivetSTL.hh \
	Tools/RivetFastJet.hh \
	Tools/RivetHepMC.hh \
	Tools/RivetYODA.hh \
	Tools/RivetMT2.hh \
	Tools/SmearingFunctions.hh \
	Tools/MomentumSmearingFunctions.hh \
	Tools/ParticleSmearingFunctions.hh \
	Tools/JetSmearingFunctions.hh \
	Tools/TypeTraits.hh \
	Tools/Utils.hh \
	Tools/Nsubjettiness/AxesDefinition.hh \
	Tools/Nsubjettiness/ExtraRecombiners.hh \
	Tools/Nsubjettiness/MeasureDefinition.hh \
	Tools/Nsubjettiness/Njettiness.hh \
	Tools/Nsubjettiness/NjettinessPlugin.hh \
	Tools/Nsubjettiness/Nsubjettiness.hh \
	Tools/Nsubjettiness/TauComponents.hh \
	Tools/Nsubjettiness/XConePlugin.hh
	nobase_dist_noinst_HEADERS += \
	Tools/osdir.hh


	## Maths
	nobase_pkginclude_HEADERS += \
	Math/Matrices.hh \
	Math/Vector3.hh \
	Math/VectorN.hh \
	Math/MatrixN.hh \
	Math/MatrixDiag.hh \
	Math/MathHeader.hh \
	Math/Vectors.hh \
	Math/LorentzTrans.hh \
	Math/Matrix3.hh \
	Math/MathUtils.hh \
	Math/Vector4.hh \
	Math/Math.hh \
	Math/Units.hh \
	Math/Constants.hh \
	Math/eigen/util.h \
	Math/eigen/regressioninternal.h \
	Math/eigen/regression.h \
	Math/eigen/vector.h \
	Math/eigen/ludecompositionbase.h \
	Math/eigen/ludecomposition.h \
	Math/eigen/linearsolver.h \
	Math/eigen/linearsolverbase.h \
	Math/eigen/matrix.h \
	Math/eigen/vectorbase.h \
	Math/eigen/projective.h \
	Math/eigen/matrixbase.h
	diff --git a/include/Rivet/Projections/PxConePlugin.hh b/include/Rivet/Projections/PxConePlugin.hh
	--- a/include/Rivet/Projections/PxConePlugin.hh
	+++ b/include/Rivet/Projections/PxConePlugin.hh
	@@ -1,165 +1,144 @@
	//FJSTARTHEADER
	-// $Id: PxConePlugin.hh,v 1.1 2019/01/04 09:47:31 leif Exp $
	+// $Id: PxConePlugin.hh 3433 2014-07-23 08:17:03Z salam $
	//
	// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. 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. They are described in the original FastJet paper,
	// hep-ph/0512210 and in the manual, arXiv:1111.6097. 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.
	//
	// 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, see <http://www.gnu.org/licenses/>.
	//----------------------------------------------------------------------
	//FJENDHEADER

	#ifndef __PXCONEPLUGIN_HH__
	#define __PXCONEPLUGIN_HH__

	#include "fastjet/JetDefinition.hh"

	// questionable whether this should be in fastjet namespace or not...

	-// FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
	+//FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
	namespace Rivet {
	-
	//----------------------------------------------------------------------
	//
	/// @ingroup plugins
	/// \class PxConePlugin
	/// Implementation of the PxCone algorithm (plugin for fastjet v2.1 upwards)
	///
	/// PxConePlugin is a plugin for fastjet (v2.1 upwards) that provides
	/// an interface to the fortran pxcone iterative cone algorithm with
	/// midpoint seeds.
	///
	/// Pxcone was written by Luis del Pozo and Michael H. Seymour. It is
	/// not a "supported" program, so if you encounter problems, you are
	/// on your own...
	///
	/// Note that pxcone sometimes encounters non-stable iterations; in
	/// such cases it returns an error -- the plugin propagates this by
	/// throwing a fastjet::Error exception; if the user wishes to have
	/// robust code, they should catch this exception.
	///
	/// Pxcone has a hard-coded limit (by default 4000) on the maximum
	/// number of particles and protojets; if the number of particles or
	/// protojets exceeds this, again a fastjet::Error exception will be
	/// thrown.
	///
	/// The functionality of pxcone is described at
	/// http://www.hep.man.ac.uk/u/wplano/ConeJet.ps
	//
	//----------------------------------------------------------------------
	class PxConePlugin : public fastjet::JetDefinition::Plugin {
	public:

	/// constructor for the PxConePlugin, whose arguments have the
	/// following meaning:
	///
	/// - the cone_radius is as usual in cone algorithms
	///
	/// - stables cones (protojets) below min_jet_energy are discarded
	/// before calling the splitting procedure to resolve overlaps
	/// (called epslon in pxcone).
	///
	/// - when two protojets overlap, if
	/// (overlapping_Et)/(Et_of_softer_protojet) < overlap_threshold
	/// the overlapping energy is split between the two protojets;
	/// otherwise the less energetic protojet is discarded. Called
	/// ovlim in pxcone.
	///
	/// - pxcone carries out p-scheme recombination, and the resulting
	/// jets are massless; setting E_scheme_jets = true (default
	/// false) doesn't change the jet composition, but the final
	/// momentum sum for the jets is carried out by direct
	/// four-vector addition instead of p-scheme recombination.
	///
	PxConePlugin (double cone_radius_in ,
	double min_jet_energy_in = 5.0 ,
	double overlap_threshold_in = 0.5,
	bool E_scheme_jets_in = false) :
	_cone_radius (cone_radius_in ),
	_min_jet_energy (min_jet_energy_in ),
	_overlap_threshold (overlap_threshold_in),
	_E_scheme_jets (E_scheme_jets_in ) {}


	// some functions to return info about parameters ----------------

	/// the cone radius
	double cone_radius () const {return _cone_radius ;}

	/// minimum jet energy (protojets below this are thrown own before
	/// merging/splitting) -- called epslon in pxcone
	double min_jet_energy () const {return _min_jet_energy ;}

	/// Maximum fraction of overlap energy in a jet -- called ovlim in pxcone.
	double overlap_threshold () const {return _overlap_threshold ;}

	/// if true then the final jets are returned as the E-scheme recombination
	/// of the particle momenta (by default, pxcone returns massless jets with
	/// a mean phi,eta type of recombination); regardless of what is
	/// returned, the internal pxcone jet-finding procedure is
	/// unaffected.
	bool E_scheme_jets() const {return _E_scheme_jets ;}


	// the things that are required by base class
	virtual std::string description () const;
	virtual void run_clustering(fastjet::ClusterSequence &) const;
	/// the plugin mechanism's standard way of accessing the jet radius
	virtual double R() const {return cone_radius();}

	private:

	double _cone_radius ;
	double _min_jet_energy ;
	double _overlap_threshold ;

	bool _E_scheme_jets;

	static bool _first_time;

	/// print a banner for reference to the 3rd-party code
	void _print_banner(std::ostream *ostr) const;
	};

	-// FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
	+//FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
	}
	-
	-//extern "C" {
	- void pxcone (
	- const int * mode , // 1=>e+e-, 2=>hadron-hadron
	- const int * ntrak , // Number of particles
	- const int * itkdm , // First dimension of PTRAK array:
	- const double * ptrak , // Array of particle 4-momenta (Px,Py,Pz,E)
	- const double * coner , // Cone size (half angle) in radians
	- const double * epslon , // Minimum Jet energy (GeV)
	- const double * ovlim , // Maximum fraction of overlap energy in a jet
	- const int * mxjet , // Maximum possible number of jets
	- int * njet , // Number of jets found
	- double * pjet, // 5-vectors of jets
	- int * ipass, // Particle k belongs to jet number IPASS(k)-1
	- // IPASS = -1 if not assosciated to a jet
	- int * ijmul, // Jet i contains IJMUL[i] particles
	- int * ierr // = 0 if all is OK ; = -1 otherwise
	- );
	-//}
	-
	#endif // __PXCONEPLUGIN_HH__
	diff --git a/include/Rivet/Projections/pxcone.h b/include/Rivet/Projections/pxcone.h
	new file mode 100644
	--- /dev/null
	+++ b/include/Rivet/Projections/pxcone.h
	@@ -0,0 +1,36 @@
	+
	+// actual physical parameters:
	+//
	+// coner
	+// epsilon
	+// ovlim
	+
	+extern "C" {
	+#ifdef WIN32
	+ void _stdcall PXCONE
	+#else
	+ void pxcone_
	+#endif
	+ (
	+ const int & mode , // 1=>e+e-, 2=>hadron-hadron
	+ const int & ntrak , // Number of particles
	+ const int & itkdm , // First dimension of PTRAK array:
	+ const double * ptrak , // Array of particle 4-momenta (Px,Py,Pz,E)
	+ const double & coner , // Cone size (half angle) in radians
	+ const double & epslon , // Minimum Jet energy (GeV)
	+ const double & ovlim , // Maximum fraction of overlap energy in a jet
	+ const int & mxjet , // Maximum possible number of jets
	+ int & njet , // Number of jets found
	+ double * pjet, // 5-vectors of jets
	+ int * ipass, // Particle k belongs to jet number IPASS(k)-1
	+ // IPASS = -1 if not assosciated to a jet
	+ int * ijmul, // Jet i contains IJMUL[i] particles
	+ int & ierr // = 0 if all is OK ; = -1 otherwise
	+ );
	+}
	+
	+#ifdef WIN32
	+#define pxcone PXCONE
	+#else
	+#define pxcone pxcone_
	+#endif
	diff --git a/src/Projections/Makefile.am b/src/Projections/Makefile.am
	--- a/src/Projections/Makefile.am
	+++ b/src/Projections/Makefile.am
	@@ -1,47 +1,48 @@
	noinst_LTLIBRARIES = libRivetProjections.la
	libRivetProjections_la_SOURCES = \
	Beam.cc \
	BeamThrust.cc \
	ChargedFinalState.cc \
	ChargedLeptons.cc \
	CentralEtHCM.cc \
	DISFinalState.cc \
	DISKinematics.cc \
	DISLepton.cc \
	DressedLeptons.cc \
	FastJets.cc \
	PxConePlugin.cc \
	+ pxcone.cc \
	FinalPartons.cc \
	FinalState.cc \
	FParameter.cc \
	HadronicFinalState.cc \
	HeavyHadrons.cc \
	Hemispheres.cc \
	IdentifiedFinalState.cc \
	InitialQuarks.cc \
	InvMassFinalState.cc \
	JetAlg.cc \
	JetShape.cc \
	LeadingParticlesFinalState.cc \
	MergedFinalState.cc \
	MissingMomentum.cc \
	NeutralFinalState.cc \
	NonHadronicFinalState.cc \
	NonPromptFinalState.cc \
	ParisiTensor.cc \
	ParticleFinder.cc \
	PrimaryHadrons.cc \
	PrimaryParticles.cc \
	PromptFinalState.cc \
	Sphericity.cc \
	Spherocity.cc \
	TauFinder.cc \
	Thrust.cc \
	TriggerCDFRun0Run1.cc \
	TriggerCDFRun2.cc \
	TriggerUA5.cc \
	UnstableFinalState.cc \
	VetoedFinalState.cc \
	VisibleFinalState.cc \
	WFinder.cc \
	ZFinder.cc
	diff --git a/src/Projections/PxConePlugin.cc b/src/Projections/PxConePlugin.cc
	--- a/src/Projections/PxConePlugin.cc
	+++ b/src/Projections/PxConePlugin.cc
	@@ -1,1387 +1,190 @@
	//FJSTARTHEADER
	-// $Id: PxConePlugin.cc,v 1.1 2019/01/04 09:47:31 leif Exp $
	+// $Id: PxConePlugin.cc 3433 2014-07-23 08:17:03Z salam $
	//
	// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. 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. They are described in the original FastJet paper,
	// hep-ph/0512210 and in the manual, arXiv:1111.6097. 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.
	//
	// 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, see <http://www.gnu.org/licenses/>.
	//----------------------------------------------------------------------
	//FJENDHEADER

	#include "Rivet/Projections/PxConePlugin.hh"

	#include "fastjet/ClusterSequence.hh"
	#include <sstream>

	+// pxcone stuff
	+#include "Rivet/Projections/pxcone.h"
	+
	+
	// FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
	namespace Rivet {

	using namespace std;
	-using namespace fastjet;

	bool PxConePlugin::_first_time = true;

	string PxConePlugin::description () const {
	ostringstream desc;

	desc << "PxCone jet algorithm with "
	<< "cone_radius = " << cone_radius () << ", "
	<< "min_jet_energy = " << min_jet_energy () << ", "
	<< "overlap_threshold = " << overlap_threshold () << ", "
	<< "E_scheme_jets = " << E_scheme_jets ()
	<< " (NB: non-standard version of PxCone, containing small bug fixes by Gavin Salam)";

	return desc.str();
	}


	-void PxConePlugin::run_clustering(ClusterSequence & clust_seq) const {
	+void PxConePlugin::run_clustering(fastjet::ClusterSequence & clust_seq) const {
	// print a banner if we run this for the first time
	//_print_banner(clust_seq.fastjet_banner_stream());

	// only have hh mode
	int mode = 2;

	int ntrak = clust_seq.jets().size(), itkdm = 4;
	double ptrak = new double[ntrak4+1];
	for (int i = 0; i < ntrak; i++) {
	ptrak[4*i+0] = clust_seq.jets()[i].px();
	ptrak[4*i+1] = clust_seq.jets()[i].py();
	ptrak[4*i+2] = clust_seq.jets()[i].pz();
	ptrak[4*i+3] = clust_seq.jets()[i].E();
	}

	// max number of allowed jets
	int mxjet = ntrak;
	int njet;
	double pjet = new double[mxjet5+1];
	int *ipass = new int[ntrak+1];
	int *ijmul = new int[mxjet+1];
	int ierr;
	- double coner = cone_radius();
	- double epslon = min_jet_energy();
	- double ovlim = overlap_threshold();

	// run pxcone
	pxcone(
	- &mode , // 1=>e+e-, 2=>hadron-hadron
	- &ntrak , // Number of particles
	- &itkdm , // First dimension of PTRAK array:
	- ptrak , // Array of particle 4-momenta (Px,Py,Pz,E)
	- &coner , // Cone size (half angle) in radians
	- &epslon , // Minimum Jet energy (GeV)
	- &ovlim , // Maximum fraction of overlap energy in a jet
	- &mxjet , // Maximum possible number of jets
	- &njet , // Number of jets found
	+ mode , // 1=>e+e-, 2=>hadron-hadron
	+ ntrak , // Number of particles
	+ itkdm , // First dimension of PTRAK array:
	+ ptrak , // Array of particle 4-momenta (Px,Py,Pz,E)
	+ cone_radius() , // Cone size (half angle) in radians
	+ min_jet_energy() , // Minimum Jet energy (GeV)
	+ overlap_threshold() , // Maximum fraction of overlap energy in a jet
	+ mxjet , // Maximum possible number of jets
	+ njet , // Number of jets found
	pjet , // 5-vectors of jets
	- ipass, // Particle k belongs to jet number IPASS(k)-1
	- // IPASS = -1 if not assosciated to a jet
	- ijmul, // Jet i contains IJMUL[i] particles
	- &ierr // = 0 if all is OK ; = -1 otherwise
	+ ipass, // Particle k belongs to jet number IPASS(k)-1
	+ // IPASS = -1 if not assosciated to a jet
	+ ijmul, // Jet i contains IJMUL[i] particles
	+ ierr // = 0 if all is OK ; = -1 otherwise
	);

	- if (ierr != 0) throw Error("An error occurred while running PXCONE");
	+ if (ierr != 0) throw fastjet::Error("An error occurred while running PXCONE");

	// now transfer information back
	valarray<int> last_index_created(njet);

	vector<vector<int> > jet_particle_content(njet);

	// get a list of particles in each jet
	for (int itrak = 0; itrak < ntrak; itrak++) {
	int jet_i = ipass[itrak] - 1;
	if (jet_i >= 0) jet_particle_content[jet_i].push_back(itrak);
	}

	// now transfer the jets back into our own structure -- we will
	// mimic the cone code with a sequential recombination sequence in
	// which the jets are built up by adding one particle at a time
	for(int ipxjet = njet-1; ipxjet >= 0; ipxjet--) {
	const vector<int> & jet_trak_list = jet_particle_content[ipxjet];
	int jet_k = jet_trak_list[0];

	for (unsigned ilist = 1; ilist < jet_trak_list.size(); ilist++) {
	int jet_i = jet_k;
	// retrieve our misappropriated index for the jet
	int jet_j = jet_trak_list[ilist];
	// do a fake recombination step with dij=0
	double dij = 0.0;
	//clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k);
	if (ilist != jet_trak_list.size()-1 \|\| E_scheme_jets()) {
	// our E-scheme recombination in cases where it doesn't matter
	clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k);
	} else {
	// put in pxcone's momentum for the last recombination so that the
	// final inclusive jet corresponds exactly to PXCONE's
	clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij,
	- PseudoJet(pjet[5ipxjet+0],pjet[5ipxjet+1],
	+ fastjet::PseudoJet(pjet[5ipxjet+0],pjet[5ipxjet+1],
	pjet[5ipxjet+2],pjet[5ipxjet+3]),
	jet_k);
	}
	}

	// NB: put a sensible looking d_iB just to be nice...
	double d_iB = clust_seq.jets()[jet_k].perp2();
	clust_seq.plugin_record_iB_recombination(jet_k, d_iB);
	}


	//// following code is for testing only
	//cout << endl;
	//for (int ijet = 0; ijet < njet; ijet++) {
	// PseudoJet jet(pjet[ijet][0],pjet[ijet][1],pjet[ijet][2],pjet[ijet][3]);
	// cout << jet.perp() << " " << jet.rap() << endl;
	//}
	//cout << "-----------------------------------------------------\n";
	//vector<PseudoJet> ourjets(clust_seq.inclusive_jets());
	//for (vector<PseudoJet>::const_iterator ourjet = ourjets.begin();
	// ourjet != ourjets.end(); ourjet++) {
	// cout << ourjet->perp() << " " << ourjet->rap() << endl;
	//}
	////cout << endl;

	delete[] ptrak;
	delete[] ipass;
	delete[] ijmul;
	delete[] pjet;
	}

	// print a banner for reference to the 3rd-party code
	void PxConePlugin::_print_banner(ostream *ostr) const{
	if (! _first_time) return;
	_first_time=false;

	// make sure the user has not set the banner stream to NULL
	if (!ostr) return;

	(*ostr) << "#-------------------------------------------------------------------------" << endl;
	(*ostr) << "# You are running the PxCone plugin for FastJet " << endl;
	(*ostr) << "# Original code by the Luis Del Pozo, David Ward and Michael H. Seymour " << endl;
	(*ostr) << "# If you use this plugin, please cite " << endl;
	(*ostr) << "# M. H. Seymour and C. Tevlin, JHEP 0611 (2006) 052 [hep-ph/0609100]. " << endl;
	(*ostr) << "# in addition to the usual FastJet reference. " << endl;
	(*ostr) << "#-------------------------------------------------------------------------" << endl;

	// make sure we really have the output done.
	ostr->flush();
	}

	// FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
	}
	-
	-/* pxcone.f -- translated by f2c and hacked by Leif Lönnblad to avoid
	- linking with libf2c.
	-*/
	-
	-#include <cmath>
	-#include <cstdio>
	-#include <cstdlib>
	-
	-/* Table of constant values, which are actually non const to be able
	- to be used as fortran arguments. */
	-static int MAXV = 20000;
	-static int VDIM = 3;
	-
	-// The standard fortran SIGN function for doubles.
	-double d_sign(const double a, const double b) {
	- return b < 0.0? -fabs(a): fabs(*a);
	-}
	-
	-// The standard fortran MOD function for doubles.
	-double d_mod(double a, double p) {
	- return a - int(a/p)*p;
	-}
	-
	-// The main PXCONE function.
	-void pxcone(const int mode, const int ntrak, const int *itkdm,
	- const double ptrak, const double coner, const double epslon, const double
	- const ovlim, const int mxjet, int njet, double pjet, int
	- ipass, int ijmul, int ierr)
	-{
	- /* Initialized data */
	-
	- static int ncall = 0;
	- static int nprint = 0;
	- static double rold = 0.;
	- static double epsold = 0.;
	- static double ovold = 0.;
	-
	- /* System generated locals */
	- int ptrak_dim1, ptrak_offset, i__1, i__2;
	- double d__1, d__2, d__3;
	-
	- /* Local variables */
	- static double cosr, rsep, ppsq, ptsq, cos2r;
	- static int i__, j, n;
	- static double vseed[3];
	- static int iterr;
	- extern /* Subroutine / int pxord_(const double , int , const int ,
	- int , double );
	- static int n1, n2;
	- static double pj[20000] /* was [4][5000] /, pp[20000] / was [4][
	- 5000] */;
	- static int mu;
	- static double pu[15000] /* was [3][5000] */, cosval;
	- extern void pxaddv(int, double , double , double *);
	- static int jetlis[25000000] /* was [5000][5000] */;
	- extern double pxmdpi(double);
	- extern /* Subroutine / int pxsear_(const int , double , const int ,
	- double , double , double , int , int *,
	- double , int , int ), pxolap_(const int , int *,
	- const int , int , double , double , const double *);
	- static int unstbl;
	- extern bool pxuvec(int, double , double );
	- extern int pxzeri(int, int *);
	- extern int pxnorv_(int , double , double , int );
	- extern void pxzerv(int, double *);
	- static double vec1[3], vec2[3];
	-
	-/* .********************************************************* */
	-/* . ------ */
	-/* . PXCONE */
	-/* . ------ */
	-/* . */
	-/* . Code downloaded from the following web page */
	-/* . */
	-/* . http://aliceinfo.cern.ch/alicvs/viewvc/JETAN/pxcone.F?view=markup&pathrev=v4-05-04 */
	-/* . */
	-/* . on 17/10/2006 by G. Salam. Permission subsequently granted by Michael */
	-/* . H. Seymour (on behalf of the PxCone authors) for this code to be */
	-/* . distributed together with FastJet under the terms of the GNU Public */
	-/* . License v2 (see the file COPYING in the main FastJet directory). */
	-/* . */
	-/* .********** Pre Release Version 26.2.93 */
	-/* . */
	-/* . Driver for the Cone Jet finding algorithm of L.A. del Pozo. */
	-/* . Based on algorithm from D.E. Soper. */
	-/* . Finds jets inside cone of half angle CONER with energy > EPSLON. */
	-/* . Jets which receive more than a fraction OVLIM of their energy from */
	-/* . overlaps with other jets are excluded. */
	-/* . Output jets are ordered in energy. */
	-/* . If MODE.EQ.2 momenta are stored as (eta,phi,<empty>,pt) */
	-/* . Usage : */
	-/* . */
	-/* . INTEGER ITKDM,MXTRK */
	-/* . PARAMETER (ITKDM=4.or.more,MXTRK=1.or.more) */
	-/* . INTEGER MXJET, MXTRAK, MXPROT */
	-/* . PARAMETER (MXJET=10,MXTRAK=500,MXPROT=500) */
	-/* . INTEGER IPASS (MXTRAK),IJMUL (MXJET) */
	-/* . INTEGER NTRAK,NJET,IERR,MODE */
	-/* . DOUBLE PRECISION PTRAK (ITKDM,MXTRK),PJET (5,MXJET) */
	-/* . DOUBLE PRECISION CONER, EPSLON, OVLIM */
	-/* . NTRAK = 1.to.MXTRAK */
	-/* . CONER = ... */
	-/* . EPSLON = ... */
	-/* . OVLIM = ... */
	-/* . CALL PXCONE (MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,MXJET, */
	-/* . + NJET,PJET,IPASS,IJMUL,IERR) */
	-/* . */
	-/* . INPUT : MODE 1=>e+e-, 2=>hadron-hadron */
	-/* . INPUT : NTRAK Number of particles */
	-/* . INPUT : ITKDM First dimension of PTRAK array */
	-/* . INPUT : PTRAK Array of particle 4-momenta (Px,Py,Pz,E) */
	-/* . INPUT : CONER Cone size (half angle) in radians */
	-/* . INPUT : EPSLON Minimum Jet energy (GeV) */
	-/* . INPUT : OVLIM Maximum fraction of overlap energy in a jet */
	-/* . INPUT : MXJET Maximum possible number of jets */
	-/* . OUTPUT : NJET Number of jets found */
	-/* . OUTPUT : PJET 5-vectors of jets */
	-/* . OUTPUT : IPASS(k) Particle k belongs to jet number IPASS(k) */
	-/* . IPASS = -1 if not assosciated to a jet */
	-/* . OUTPUT : IJMUL(i) Jet i contains IJMUL(i) particles */
	-/* . OUTPUT : IERR = 0 if all is OK ; = -1 otherwise */
	-/* . */
	-/* . CALLS : PXSEAR, PXSAME, PXNEW, PXTRY, PXORD, PXUVEC, PXOLAP */
	-/* . CALLED : User */
	-/* . */
	-/* . AUTHOR : L.A. del Pozo */
	-/* . CREATED : 26-Feb-93 */
	-/* . LAST MOD : 2-Mar-93 */
	-/* . */
	-/* . Modification Log. */
	-/* . 25-Feb-07: G P Salam - fix bugs concerning 2pi periodicity in eta phi mode */
	-/* . - added commented code to get consistent behaviour */
	-/* . regardless of particle order (replaces n-way */
	-/* . midpoints with 2-way midpoints however...) */
	-/* . 2-Jan-97: M Wobisch - fix bug concerning COS2R in eta phi mode */
	-/* . 4-Apr-93: M H Seymour - Change 2d arrays to 1d in PXTRY & PXNEW */
	-/* . 2-Apr-93: M H Seymour - Major changes to add boost-invariant mode */
	-/* . 1-Apr-93: M H Seymour - Increase all array sizes */
	-/* . 30-Mar-93: M H Seymour - Change all REAL variables to DOUBLE PRECISION */
	-/* . 30-Mar-93: M H Seymour - Change OVLIM into an input parameter */
	-/* . 2-Mar-93: L A del Pozo - Fix Bugs in PXOLAP */
	-/* . 1-Mar-93: L A del Pozo - Remove Cern library routine calls */
	-/* . 1-Mar-93: L A del Pozo - Add Print out of welcome and R and Epsilon */
	-/* . */
	-/* .********************************************************* */
	-/* +SEQ,DECLARE. */
	-/* ** External Arrays */
	-/* ** Internal Arrays */
	-/* ** Used in the routine. */
	-/* MWobisch */
	-/* MWobisch */
	- /* Parameter adjustments */
	- --ipass;
	- ptrak_dim1 = *itkdm;
	- ptrak_offset = 1 + ptrak_dim1 * 1;
	- ptrak -= ptrak_offset;
	- --ijmul;
	- pjet -= 6;
	-
	- /* Function Body */
	-/* MWobisch */
	-/* *************************************** */
	- rsep = 2.;
	-/* *************************************** */
	-/* MWobisch */
	- *ierr = 0;
	-
	-/* ** INITIALIZE */
	- if (ncall <= 0) {
	- rold = (float)0.;
	- epsold = (float)0.;
	- ovold = (float)0.;
	- }
	- ++ncall;
	-
	-/* ** Print welcome and Jetfinder parameters */
	- if ((coner != rold \|\| epslon != epsold \|\| *ovlim != ovold) && nprint <=
	- 10) {
	- printf("%s\n", " ********* PXCONE: Cone Jet-finder *********");
	- printf("%s\n", " Written by Luis Del Pozo of OPAL");
	- printf("%s\n", " Modified for eta-phi by Mike Seymour");
	- printf("%s\n", " Includes bug fixes by Wobisch, Salam");
	- printf("%s%5.2f%s\n", " Cone Size R = ",*coner," Radians");
	- printf("%s%5.2f%s\n", " Min Jet energy Epsilon = ",*epslon," GeV");
	- printf("%s%5.2f\n", " Overlap fraction parameter = ",*ovlim);
	- printf("%s\n", " PXCONE is not a supported product and is");
	- printf("%s\n", " is provided for comparative purposes only");
	- printf("%s\n", " ***********************************************");
	-/* WRITE (6,) /
	-/* WRITE (6,) ' ******** PXCONE: Cone Jet-finder ********' /
	-/* WRITE (6,) ' Written by Luis Del Pozo of OPAL' /
	-/* WRITE (6,) ' Modified for eta-phi by Mike Seymour' /
	-/* WRITE (6,) ' Includes bug fixes by Wobisch, Salam' /
	-/* WRITE(6,1000)' Cone Size R = ',CONER,' Radians' */
	-/* WRITE(6,1001)' Min Jet energy Epsilon = ',EPSLON,' GeV' */
	-/* WRITE(6,1002)' Overlap fraction parameter = ',OVLIM */
	-/* WRITE (6,) ' PXCONE is not a supported product and is' /
	-/* WRITE (6,) ' is provided for comparative purposes only' /
	-/* WRITE (6,) ' *********************************************' /
	-/* MWobisch */
	- if (rsep < (float)1.999) {
	- printf("%s\n", " ******************************************");
	- printf("%s\n", " ******************************************");
	- printf("%s\n", " M Wobisch: private change !!!!!!!!!!!! ");
	- printf("%s%5.2f\n", " Rsep is set to ",rsep);
	- printf("%s\n", " this is ONLY meaningful in a NLO calculation");
	- printf("%s\n", " ------------------------ ");
	- printf("%s\n", " please check what you're doing!!");
	- printf("%s\n", " ******************************************");
	- printf("%s\n", " ******************************************");
	- printf("%s\n", " ******************************************");
	- printf("%s\n", "");
	- printf("%s\n", "");
	- printf("%s\n", "");
	-/* WRITE(6,) ' ' /
	-/* WRITE (6,) ' ****************************************' /
	-/* WRITE (6,) ' ****************************************' /
	-/* WRITE(6,) ' M Wobisch: private change !!!!!!!!!!!! ' /
	-/* WRITE(6,) ' Rsep is set to ',RSEP /
	-/* WRITE(6,) ' this is ONLY meaningful in a NLO calculation' /
	-/* WRITE(6,) ' ------------------------ ' /
	-/* WRITE(6,) ' please check what you''re doing!!' /
	-/* WRITE(6,) ' or ask: Markus.Wobisch@desy.de --' /
	-/* WRITE (6,) ' ****************************************' /
	-/* WRITE (6,) ' ****************************************' /
	-/* WRITE (6,) ' ****************************************' /
	-/* WRITE(6,) ' ' /
	-/* WRITE(6,) ' ' /
	- }
	-/* MWobisch */
	-/* WRITE (6,) /
	-/* 1000 FORMAT(A18,F5.2,A10) */
	-/* 1001 FORMAT(A29,F5.2,A5) */
	-/* 1002 FORMAT(A33,F5.2) */
	- ++nprint;
	- rold = *coner;
	- epsold = *epslon;
	- ovold = *ovlim;
	- }
	-
	-/* ** Copy calling array PTRAK to internal array PP(4,NTRAK) */
	-
	- if (*ntrak > 5000) {
	-/* WRITE (6,) ' PXCONE: Ntrak too large: ',NTRAK /
	- printf("%s%d\n", " PXCONE: Ntrak too large: ", *ntrak);
	- *ierr = -1;
	- return;
	- }
	- if (*mode != 2) {
	- i__1 = *ntrak;
	- for (i__ = 1; i__ <= i__1; ++i__) {
	- for (j = 1; j <= 4; ++j) {
	- pp[j + (i__ << 2) - 5] = ptrak[j + i__ * ptrak_dim1];
	-/* L101: */
	- }
	-/* L100: */
	- }
	- } else {
	-/* ** Converting to eta,phi,pt if necessary */
	- i__1 = *ntrak;
	- for (i__ = 1; i__ <= i__1; ++i__) {
	-/* Computing 2nd power */
	- d__1 = ptrak[i__ * ptrak_dim1 + 1];
	-/* Computing 2nd power */
	- d__2 = ptrak[i__ * ptrak_dim1 + 2];
	- ptsq = d__1 * d__1 + d__2 * d__2;
	-/* Computing 2nd power */
	- d__3 = ptrak[i__ * ptrak_dim1 + 3];
	-/* Computing 2nd power */
	- d__2 = sqrt(ptsq + d__3 * d__3) + (d__1 = ptrak[i__ * ptrak_dim1
	- + 3], abs(d__1));
	- ppsq = d__2 * d__2;
	- if (ptsq <= ppsq * (float)4.25e-18) {
	- pp[(i__ << 2) - 4] = 20.;
	- } else {
	- pp[(i__ << 2) - 4] = log(ppsq / ptsq) * (float).5;
	- }
	- pp[(i__ << 2) - 4] = d_sign(&pp[(i__ << 2) - 4], &ptrak[i__ *
	- ptrak_dim1 + 3]);
	- if (ptsq == 0.) {
	- pp[(i__ << 2) - 3] = 0.;
	- } else {
	- pp[(i__ << 2) - 3] = atan2(ptrak[i__ * ptrak_dim1 + 2], ptrak[
	- i__ * ptrak_dim1 + 1]);
	- }
	- pp[(i__ << 2) - 2] = 0.;
	- pp[(i__ << 2) - 1] = sqrt(ptsq);
	- pu[i__ * 3 - 3] = pp[(i__ << 2) - 4];
	- pu[i__ * 3 - 2] = pp[(i__ << 2) - 3];
	- pu[i__ * 3 - 1] = pp[(i__ << 2) - 2];
	-/* L104: */
	- }
	- }
	-
	-/* ** Zero output variables */
	-
	- *njet = 0;
	- i__1 = *ntrak;
	- for (i__ = 1; i__ <= i__1; ++i__) {
	- for (j = 1; j <= 5000; ++j) {
	- jetlis[j + i__ * 5000 - 5001] = false;
	-/* L103: */
	- }
	-/* L102: */
	- }
	- pxzerv(MAXV, pj);
	- pxzeri(*mxjet, &ijmul[1]);
	-
	- if (*mode != 2) {
	- cosr = cos(*coner);
	- cos2r = cos(*coner);
	- } else {
	-/* Purely for convenience, work in terms of 1-R2 */
	-/* Computing 2nd power */
	- d__1 = *coner;
	- cosr = 1 - d__1 * d__1;
	-/* MW -- select Rsep: 1-(RsepCONER)2 /
	-/* Computing 2nd power */
	- d__1 = rsep * *coner;
	- cos2r = 1 - d__1 * d__1;
	-/* ORIGINAL COS2R = 1-(2CONER)2 /
	- }
	- unstbl = false;
	- if (*mode != 2) {
	- if ( !pxuvec(*ntrak, pp, pu) ) {
	- *ierr =1;
	- return;
	- }
	- }
	-/* ** Look for jets using particle diretions as seed axes */
	-
	- i__1 = *ntrak;
	- for (n = 1; n <= i__1; ++n) {
	- for (mu = 1; mu <= 3; ++mu) {
	- vseed[mu - 1] = pu[mu + n * 3 - 4];
	-/* L120: */
	- }
	- pxsear_(mode, &cosr, ntrak, pu, pp, vseed, njet, jetlis, pj, &unstbl,
	- ierr);
	- if (*ierr != 0) {
	- return;
	- }
	-/* L110: */
	- }
	-/* MW - for Rsep=1 goto 145 */
	-/* GOTO 145 */
	-/* ** Now look between all pairs of jets as seed axes. */
	-/* NJTORG = NJET ! GPS -- to get consistent behaviour (2-way midpnts) */
	-/* DO 140 N1 = 1,NJTORG-1 ! GPS -- to get consistent behaviour (2-way midpnts) */
	- i__1 = *njet - 1;
	- for (n1 = 1; n1 <= i__1; ++n1) {
	- vec1[0] = pj[(n1 << 2) - 4];
	- vec1[1] = pj[(n1 << 2) - 3];
	- vec1[2] = pj[(n1 << 2) - 2];
	- if (*mode != 2) {
	- pxnorv_(&VDIM, vec1, vec1, &iterr);
	- }
	-/* DO 150 N2 = N1+1,NJTORG ! GPS -- to get consistent behaviour */
	- i__2 = *njet;
	- for (n2 = n1 + 1; n2 <= i__2; ++n2) {
	- vec2[0] = pj[(n2 << 2) - 4];
	- vec2[1] = pj[(n2 << 2) - 3];
	- vec2[2] = pj[(n2 << 2) - 2];
	- if (*mode != 2) {
	- pxnorv_(&VDIM, vec2, vec2, &iterr);
	- }
	- pxaddv(VDIM, vec1, vec2, vseed);
	- if (*mode != 2) {
	- pxnorv_(&VDIM, vseed, vseed, &iterr);
	- } else {
	- vseed[0] /= 2;
	-/* VSEED(2)=VSEED(2)/2 */
	-/* GPS 25/02/07 */
	- d__2 = vec2[1] - vec1[1];
	- d__1 = vec1[1] + pxmdpi(d__2) * .5;
	- vseed[1] = pxmdpi(d__1);
	- }
	-/* ---ONLY BOTHER IF THEY ARE BETWEEN 1 AND 2 CONE RADII APART */
	- if (*mode != 2) {
	- cosval = vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] *
	- vec2[2];
	- } else {
	- if (abs(vec1[0]) >= 20. \|\| abs(vec2[0]) >= 20.) {
	- cosval = -1e3;
	- } else {
	-/* Computing 2nd power */
	- d__1 = vec1[0] - vec2[0];
	- d__3 = vec1[1] - vec2[1];
	-/* Computing 2nd power */
	- d__2 = pxmdpi(d__3);
	- cosval = 1 - (d__1 * d__1 + d__2 * d__2);
	- }
	- }
	- if (cosval <= cosr && cosval >= cos2r) {
	- pxsear_(mode, &cosr, ntrak, pu, pp, vseed, njet, jetlis, pj, &
	- unstbl, ierr);
	- }
	-/* CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET, */
	-/* + JETLIS,PJ,UNSTBL,IERR) */
	- if (*ierr != 0) {
	- return;
	- }
	-/* L150: */
	- }
	-/* L140: */
	- }
	- if (unstbl) {
	- *ierr = -1;
	-/* WRITE (6,) ' PXCONE: Too many iterations to find a proto-jet' /
	- printf(" PXCONE: Too many iterations to find a proto-jet\n");
	- return;
	- }
	-/* L145: */
	-/* ** Now put the jet list into order by jet energy, eliminating jets */
	-/* ** with energy less than EPSLON. */
	- pxord_(epslon, njet, ntrak, jetlis, pj);
	-
	-/* ** Take care of jet overlaps */
	- pxolap_(mode, njet, ntrak, jetlis, pj, pp, ovlim);
	-
	-/* ** Order jets again as some have been eliminated, or lost energy. */
	- pxord_(epslon, njet, ntrak, jetlis, pj);
	-
	-/* ** All done!, Copy output into output arrays */
	- if (njet > mxjet) {
	-/* WRITE (6,) ' PXCONE: Found more than MXJET jets' /
	- printf(" PXCONE: Found more than MXJET jets\n");
	- *ierr = -1;
	- goto L99;
	- }
	- if (*mode != 2) {
	- i__1 = *njet;
	- for (i__ = 1; i__ <= i__1; ++i__) {
	- for (j = 1; j <= 4; ++j) {
	- pjet[j + i__ * 5] = pj[j + (i__ << 2) - 5];
	-/* L310: */
	- }
	-/* L300: */
	- }
	- } else {
	- i__1 = *njet;
	- for (i__ = 1; i__ <= i__1; ++i__) {
	- pjet[i__ * 5 + 1] = pj[(i__ << 2) - 1] * cos(pj[(i__ << 2) - 3]);
	- pjet[i__ * 5 + 2] = pj[(i__ << 2) - 1] * sin(pj[(i__ << 2) - 3]);
	- pjet[i__ * 5 + 3] = pj[(i__ << 2) - 1] * sinh(pj[(i__ << 2) - 4]);
	- pjet[i__ * 5 + 4] = pj[(i__ << 2) - 1] * cosh(pj[(i__ << 2) - 4]);
	-/* L315: */
	- }
	- }
	- i__1 = *ntrak;
	- for (i__ = 1; i__ <= i__1; ++i__) {
	- ipass[i__] = -1;
	- i__2 = *njet;
	- for (j = 1; j <= i__2; ++j) {
	- if (jetlis[j + i__ * 5000 - 5001]) {
	- ++ijmul[j];
	- ipass[i__] = j;
	- }
	-/* L330: */
	- }
	-/* L320: */
	- }
	-L99:
	- return;
	-} /* pxcone_ */
	-
	-/* CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper */
	-/* -- Author : */
	-/* ----------------------------------------------------------------------- */
	-/* Subroutine / int pxnorv_(int n, double a, double b,
	- int *iterr)
	-{
	- /* System generated locals */
	- int i__1;
	- double d__1;
	-
	- /* Builtin functions */
	- double sqrt(double);
	-
	- /* Local variables */
	- static double c__;
	- static int i__;
	-
	- /* Parameter adjustments */
	- --b;
	- --a;
	-
	- /* Function Body */
	- c__ = 0.;
	- i__1 = *n;
	- for (i__ = 1; i__ <= i__1; ++i__) {
	-/* Computing 2nd power */
	- d__1 = a[i__];
	- c__ += d__1 * d__1;
	-/* L10: */
	- }
	- if (c__ <= 0.) {
	- return 0;
	- }
	- c__ = 1 / sqrt(c__);
	- i__1 = *n;
	- for (i__ = 1; i__ <= i__1; ++i__) {
	- b[i__] = a[i__] * c__;
	-/* L20: */
	- }
	- return 0;
	-} /* pxnorv_ */
	-
	-/* CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper */
	-/* CMZ : 1.06/00 15/03/94 12.17.46 by P. Schleper */
	-/* -- Author : */
	-
	-/* +DECK,PXOLAP. */
	-/* Subroutine / int pxolap_(const int mode, int njet, const int ntrak,
	- int jetlis, double pj, double pp, const double ovlim)
	-{
	- /* Initialized data */
	-
	- static int ijmin = 0;
	-
	- /* System generated locals */
	- int i__1, i__2, i__3;
	- double d__1, d__2, d__3;
	-
	- /* Local variables */
	- static int ijet[5000];
	- static double thet, cost;
	- static int i__, j, n;
	- static double eover, thmin;
	- extern void pxang3(double , double , double &, double &);
	- static int nj, mu;
	- static int ovelap;
	- extern double pxmdpi(double);
	- static double vec1[3], vec2[3];
	-
	-
	-/* ** Looks for particles assigned to more than 1 jet, and reassigns them */
	-/* ** If more than a fraction OVLIM of a jet's energy is contained in */
	-/* ** higher energy jets, that jet is neglected. */
	-/* ** Particles assigned to the jet closest in angle (a la CDF, Snowmass). */
	-/* +SEQ,DECLARE. */
	- /* Parameter adjustments */
	- pp -= 5;
	- pj -= 5;
	- jetlis -= 5001;
	-
	- /* Function Body */
	-
	- if (*njet <= 1) {
	- return 0;
	- }
	-/* ** Look for jets with large overlaps with higher energy jets. */
	- i__1 = *njet;
	- for (i__ = 2; i__ <= i__1; ++i__) {
	-/* ** Find overlap energy between jets I and all higher energy jets. */
	- eover = (float)0.;
	- i__2 = *ntrak;
	- for (n = 1; n <= i__2; ++n) {
	- ovelap = false;
	- i__3 = i__ - 1;
	- for (j = 1; j <= i__3; ++j) {
	- if (jetlis[i__ + n * 5000] && jetlis[j + n * 5000]) {
	- ovelap = true;
	- }
	-/* L120: */
	- }
	- if (ovelap) {
	- eover += pp[(n << 2) + 4];
	- }
	-/* L110: */
	- }
	-/* ** Is the fraction of energy shared larger than OVLIM? */
	- if (eover > ovlim pj[(i__ << 2) + 4]) {
	-/* ** De-assign all particles from Jet I */
	- i__2 = *ntrak;
	- for (n = 1; n <= i__2; ++n) {
	- jetlis[i__ + n * 5000] = false;
	-/* L130: */
	- }
	- }
	-/* L100: */
	- }
	-/* ** Now there are no big overlaps, assign every particle in */
	-/* ** more than 1 jet to the closet jet. */
	-/* ** Any particles now in more than 1 jet are assigned to the CLOSET */
	-/* ** jet (in angle). */
	- i__1 = *ntrak;
	- for (i__ = 1; i__ <= i__1; ++i__) {
	- nj = 0;
	- i__2 = *njet;
	- for (j = 1; j <= i__2; ++j) {
	- if (jetlis[j + i__ * 5000]) {
	- ++nj;
	- ijet[nj - 1] = j;
	- }
	-/* L150: */
	- }
	- if (nj > 1) {
	-/* ** Particle in > 1 jet - calc angles... */
	- vec1[0] = pp[(i__ << 2) + 1];
	- vec1[1] = pp[(i__ << 2) + 2];
	- vec1[2] = pp[(i__ << 2) + 3];
	- thmin = (float)0.;
	- i__2 = nj;
	- for (j = 1; j <= i__2; ++j) {
	- vec2[0] = pj[(ijet[j - 1] << 2) + 1];
	- vec2[1] = pj[(ijet[j - 1] << 2) + 2];
	- vec2[2] = pj[(ijet[j - 1] << 2) + 3];
	- if (*mode != 2) {
	- pxang3(vec1, vec2, cost, thet);
	- } else {
	-/* Computing 2nd power */
	- d__1 = vec1[0] - vec2[0];
	- d__3 = vec1[1] - vec2[1];
	-/* Computing 2nd power */
	- d__2 = pxmdpi(d__3);
	- thet = d__1 * d__1 + d__2 * d__2;
	- }
	- if (j == 1) {
	- thmin = thet;
	- ijmin = ijet[j - 1];
	- } else if (thet < thmin) {
	- thmin = thet;
	- ijmin = ijet[j - 1];
	- }
	-/* L160: */
	- }
	-/* ** Assign track to IJMIN */
	- i__2 = *njet;
	- for (j = 1; j <= i__2; ++j) {
	- jetlis[j + i__ * 5000] = false;
	-/* L170: */
	- }
	- jetlis[ijmin + i__ * 5000] = true;
	- }
	-/* L140: */
	- }
	-/* ** Recompute PJ */
	- i__1 = *njet;
	- for (i__ = 1; i__ <= i__1; ++i__) {
	- for (mu = 1; mu <= 4; ++mu) {
	- pj[mu + (i__ << 2)] = (float)0.;
	-/* L210: */
	- }
	- i__2 = *ntrak;
	- for (n = 1; n <= i__2; ++n) {
	- if (jetlis[i__ + n * 5000]) {
	- if (*mode != 2) {
	- for (mu = 1; mu <= 4; ++mu) {
	- pj[mu + (i__ << 2)] += pp[mu + (n << 2)];
	-/* L230: */
	- }
	- } else {
	- pj[(i__ << 2) + 1] += pp[(n << 2) + 4] / (pp[(n << 2) + 4]
	- + pj[(i__ << 2) + 4]) * (pp[(n << 2) + 1] - pj[(
	- i__ << 2) + 1]);
	-/* GPS 25/02/07 */
	- d__2 = pp[(n << 2) + 2] - pj[(i__ << 2) + 2];
	- d__1 = pj[(i__ << 2) + 2] + pp[(n << 2) + 4] / (pp[(n <<
	- 2) + 4] + pj[(i__ << 2) + 4]) * pxmdpi(d__2);
	- pj[(i__ << 2) + 2] = pxmdpi(d__1);
	-/* PJ(2,I)=PJ(2,I) */
	-/* + + PP(4,N)/(PP(4,N)+PJ(4,I))PXMDPI(PP(2,N)-PJ(2,I)) /
	- pj[(i__ << 2) + 4] += pp[(n << 2) + 4];
	- }
	- }
	-/* L220: */
	- }
	-/* L200: */
	- }
	- return 0;
	-} /* pxolap_ */
	-
	-/* CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper */
	-/* CMZ : 1.06/00 14/03/94 15.37.45 by P. Schleper */
	-/* -- Author : */
	-
	-/* +DECK,PXORD. */
	-/* Subroutine / int pxord_(const double epslon, int njet, const int ntrak,
	- int jetlis, double pj)
	-{
	- /* System generated locals */
	- int i__1, i__2;
	-
	- /* Local variables */
	- static int i__, j, index[5000];
	- static double elist[5000], ptemp[20000] /* was [4][5000] */;
	- static int logtmp[25000000] /* was [5000][5000] */;
	- extern /* Subroutine / int pxsorv_(int , double , int ,
	- char);
	-
	-
	-/* ** Routine to put jets into order and eliminate tose less than EPSLON */
	-/* +SEQ,DECLARE. */
	-/* ** Puts jets in order of energy: 1 = highest energy etc. */
	-/* ** Then Eliminate jets with energy below EPSLON */
	-
	-/* ** Copy input arrays. */
	- /* Parameter adjustments */
	- pj -= 5;
	- jetlis -= 5001;
	-
	- /* Function Body */
	- i__1 = *njet;
	- for (i__ = 1; i__ <= i__1; ++i__) {
	- for (j = 1; j <= 4; ++j) {
	- ptemp[j + (i__ << 2) - 5] = pj[j + (i__ << 2)];
	-/* L110: */
	- }
	- i__2 = *ntrak;
	- for (j = 1; j <= i__2; ++j) {
	- logtmp[i__ + j * 5000 - 5001] = jetlis[i__ + j * 5000];
	-/* L120: */
	- }
	-/* L100: */
	- }
	- i__1 = *njet;
	- for (i__ = 1; i__ <= i__1; ++i__) {
	- elist[i__ - 1] = pj[(i__ << 2) + 4];
	-/* L150: */
	- }
	-/* ** Sort the energies... */
	- pxsorv_(njet, elist, index, 'I');
	-/* ** Fill PJ and JETLIS according to sort ( sort is in ascending order!!) */
	- i__1 = *njet;
	- for (i__ = 1; i__ <= i__1; ++i__) {
	- for (j = 1; j <= 4; ++j) {
	- pj[j + (i__ << 2)] = ptemp[j + (index[*njet + 1 - i__ - 1] << 2)
	- - 5];
	-/* L210: */
	- }
	- i__2 = *ntrak;
	- for (j = 1; j <= i__2; ++j) {
	- jetlis[i__ + j * 5000] = logtmp[index[njet + 1 - i__ - 1] + j
	- 5000 - 5001];
	-/* L220: */
	- }
	-/* L200: */
	- }
	-/* * Jets are now in order */
	-/* ** Now eliminate jets with less than Epsilon energy */
	- i__1 = *njet;
	- for (i__ = 1; i__ <= i__1; ++i__) {
	- if (pj[(i__ << 2) + 4] < *epslon) {
	- --(*njet);
	- pj[(i__ << 2) + 4] = (float)0.;
	- }
	-/* L300: */
	- }
	- return 0;
	-} /* pxord_ */
	-
	-/* ******************************************************************* */
	-/* CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper */
	-/* CMZ : 1.06/00 14/03/94 15.37.44 by P. Schleper */
	-/* -- Author : */
	-/* +DECK,PXSEAR. */
	-/* Subroutine / int pxsear_(const int mode, double cosr, const int ntrak,
	- double pu, double pp, double vseed, int njet,
	- int jetlis, double pj, int unstbl, int ierr)
	-{
	- /* System generated locals */
	- int i__1;
	-
	- /* Local variables */
	- static int iter;
	- static double pnew[4];
	- static int n;
	- static double naxis[3], oaxis[3];
	- extern int pxnew(int , int , int, int);
	- extern /* Subroutine / int pxtry_(const int , double , const int ,
	- double , double , double , double ,
	- double , int , int *);
	- static int ok;
	- static int mu;
	- static int oldlis[5000];
	- extern int pxsame_(int , int , int);
	- static int newlis[5000];
	-
	-
	-/* +SEQ,DECLARE. */
	-/* ** Using VSEED as a trial axis , look for a stable jet. */
	-/* ** Check stable jets against those already found and add to PJ. */
	-/* ** Will try up to MXITER iterations to get a stable set of particles */
	-/* ** in the cone. */
	-
	- /* Parameter adjustments */
	- pj -= 5;
	- jetlis -= 5001;
	- --vseed;
	- pp -= 5;
	- pu -= 4;
	-
	- /* Function Body */
	- for (mu = 1; mu <= 3; ++mu) {
	- oaxis[mu - 1] = vseed[mu];
	-/* L100: */
	- }
	- i__1 = *ntrak;
	- for (n = 1; n <= i__1; ++n) {
	- oldlis[n - 1] = false;
	-/* L110: */
	- }
	- for (iter = 1; iter <= 30; ++iter) {
	- pxtry_(mode, cosr, ntrak, &pu[4], &pp[5], oaxis, naxis, pnew, newlis,
	- &ok);
	-/* ** Return immediately if there were no particles in the cone. */
	- if (! ok) {
	- return 0;
	- }
	- if (pxsame_(newlis, oldlis, *ntrak)) {
	-/* ** We have a stable jet. */
	- if (pxnew(newlis, &jetlis[5001], ntrak, njet)) {
	-/* ** And the jet is a new one. So add it to our arrays. */
	-/* ** Check arrays are big anough... */
	- if (*njet == 5000) {
	-/* WRITE (6,) ' PXCONE: Found more than MXPROT proto-jets' /
	- printf(" PXCONE: Found more than MXPROT proto-jets\n");
	- *ierr = -1;
	- return 0;
	- }
	- ++(*njet);
	- i__1 = *ntrak;
	- for (n = 1; n <= i__1; ++n) {
	- jetlis[njet + n 5000] = newlis[n - 1];
	-/* L130: */
	- }
	- for (mu = 1; mu <= 4; ++mu) {
	- pj[mu + (*njet << 2)] = pnew[mu - 1];
	-/* L140: */
	- }
	- }
	- return 0;
	- }
	-/* ** The jet was not stable, so we iterate again */
	- i__1 = *ntrak;
	- for (n = 1; n <= i__1; ++n) {
	- oldlis[n - 1] = newlis[n - 1];
	-/* L150: */
	- }
	- for (mu = 1; mu <= 3; ++mu) {
	- oaxis[mu - 1] = naxis[mu - 1];
	-/* L160: */
	- }
	-/* L120: */
	- }
	- *unstbl = true;
	- return 0;
	-} /* pxsear_ */
	-
	-/* CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper */
	-/* -- Author : */
	-/* ----------------------------------------------------------------------- */
	-/* Subroutine / int pxsorv_(int n, double a, int k, char opt)
	-{
	- /* System generated locals */
	- int i__1;
	-
	- /* Builtin functions */
	- // /* Subroutine / int s_stop(char , ftnlen);
	-
	- /* Local variables */
	- static double b[5000];
	- static int i__, j, il[5000], ir[5000];
	-
	-/* Sort A(N) into ascending order */
	-/* OPT = 'I' : return index array K only */
	-/* OTHERWISE : return sorted A and index array K */
	-/* ----------------------------------------------------------------------- */
	-
	-/* INT N,I,J,K(N),IL(NMAX),IR(NMAX) */
	-/* LUND */
	-
	-/* DOUBLE PRECISION A(N),B(NMAX) */
	-/* LUND */
	- /* Parameter adjustments */
	- --k;
	- --a;
	-
	- /* Function Body */
	- if (*n > 5000) {
	- // WRITE s_stop("Sorry, not enough room in Mike's PXSORV", (ftnlen)39);
	- printf("Sorry, not enough room in Mike's PXSORV\n");
	- abort();
	- }
	- il[0] = 0;
	- ir[0] = 0;
	- i__1 = *n;
	- for (i__ = 2; i__ <= i__1; ++i__) {
	- il[i__ - 1] = 0;
	- ir[i__ - 1] = 0;
	- j = 1;
	-L2:
	- if (a[i__] > a[j]) {
	- goto L5;
	- }
	-/* L3: */
	- if (il[j - 1] == 0) {
	- goto L4;
	- }
	- j = il[j - 1];
	- goto L2;
	-L4:
	- ir[i__ - 1] = -j;
	- il[j - 1] = i__;
	- goto L10;
	-L5:
	- if (ir[j - 1] <= 0) {
	- goto L6;
	- }
	- j = ir[j - 1];
	- goto L2;
	-L6:
	- ir[i__ - 1] = ir[j - 1];
	- ir[j - 1] = i__;
	-L10:
	- ;
	- }
	- i__ = 1;
	- j = 1;
	- goto L8;
	-L20:
	- j = il[j - 1];
	-L8:
	- if (il[j - 1] > 0) {
	- goto L20;
	- }
	-L9:
	- k[i__] = j;
	- b[i__ - 1] = a[j];
	- ++i__;
	- if ((i__1 = ir[j - 1]) < 0) {
	- goto L12;
	- } else if (i__1 == 0) {
	- goto L30;
	- } else {
	- goto L13;
	- }
	-L13:
	- j = ir[j - 1];
	- goto L8;
	-L12:
	- j = -ir[j - 1];
	- goto L9;
	-L30:
	- if ( opt == 'I') {
	- return 0;
	- }
	- i__1 = *n;
	- for (i__ = 1; i__ <= i__1; ++i__) {
	-/* L31: */
	- a[i__] = b[i__ - 1];
	- }
	-/* L999: */
	- return 0;
	-} /* pxsorv_ */
	-
	-/* ******************************************************************** */
	-/* CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper */
	-/* CMZ : 1.06/00 14/03/94 15.37.44 by P. Schleper */
	-/* -- Author : */
	-
	-/* +DECK,PXTRY. */
	-/* Subroutine / int pxtry_(const int mode, double cosr, const int ntrak,
	- double pu, double pp, double oaxis, double naxis,
	- double pnew, int newlis, int *ok)
	-{
	- /* System generated locals */
	- int i__1;
	- double d__1, d__2, d__3;
	-
	- /* Builtin functions */
	- double sqrt(double);
	-
	- /* Local variables */
	- static double norm;
	- static int n, mu;
	- static double cosval;
	- extern double pxmdpi(double);
	- static double normsq;
	- static int npp, npu;
	-
	-
	-/* +SEQ,DECLARE. */
	-/* ** Note that although PU and PP are assumed to be 2d arrays, they */
	-/* ** are used as 1d in this routine for efficiency */
	-/* ** Finds all particles in cone of size COSR about OAXIS direction. */
	-/* ** Calculates 4-momentum sum of all particles in cone (PNEW) , and */
	-/* ** returns this as new jet axis NAXIS (Both unit Vectors) */
	-
	- /* Parameter adjustments */
	- --newlis;
	- --pnew;
	- --naxis;
	- --oaxis;
	- --pp;
	- --pu;
	-
	- /* Function Body */
	- *ok = false;
	- for (mu = 1; mu <= 4; ++mu) {
	- pnew[mu] = (float)0.;
	-/* L100: */
	- }
	- npu = -3;
	- npp = -4;
	- i__1 = *ntrak;
	- for (n = 1; n <= i__1; ++n) {
	- npu += 3;
	- npp += 4;
	- if (*mode != 2) {
	- cosval = (float)0.;
	- for (mu = 1; mu <= 3; ++mu) {
	- cosval += oaxis[mu] * pu[mu + npu];
	-/* L120: */
	- }
	- } else {
	- if ((d__1 = pu[npu + 1], abs(d__1)) >= 20. \|\| abs(oaxis[1]) >=
	- 20.) {
	- cosval = -1e3;
	- } else {
	-/* Computing 2nd power */
	- d__1 = oaxis[1] - pu[npu + 1];
	- d__3 = oaxis[2] - pu[npu + 2];
	-/* Computing 2nd power */
	- d__2 = pxmdpi(d__3);
	- cosval = 1 - (d__1 * d__1 + d__2 * d__2);
	- }
	- }
	- if (cosval >= *cosr) {
	- newlis[n] = true;
	- *ok = true;
	- if (*mode != 2) {
	- for (mu = 1; mu <= 4; ++mu) {
	- pnew[mu] += pp[mu + npp];
	-/* L130: */
	- }
	- } else {
	- pnew[1] += pp[npp + 4] / (pp[npp + 4] + pnew[4]) * (pp[npp +
	- 1] - pnew[1]);
	-/* PNEW(2)=PNEW(2) */
	-/* + + PP(4+NPP)/(PP(4+NPP)+PNEW(4)) */
	-/* + PXMDPI(PP(2+NPP)-PNEW(2)) /
	-/* GPS 25/02/07 */
	- d__2 = pp[npp + 2] - pnew[2];
	- d__1 = pnew[2] + pp[npp + 4] / (pp[npp + 4] + pnew[4]) *
	- pxmdpi(d__2);
	- pnew[2] = pxmdpi(d__1);
	- pnew[4] += pp[npp + 4];
	- }
	- } else {
	- newlis[n] = false;
	- }
	-/* L110: */
	- }
	-/* ** If there are particles in the cone, calc new jet axis */
	- if (*ok) {
	- if (*mode != 2) {
	- normsq = (float)0.;
	- for (mu = 1; mu <= 3; ++mu) {
	-/* Computing 2nd power */
	- d__1 = pnew[mu];
	- normsq += d__1 * d__1;
	-/* L140: */
	- }
	- norm = sqrt(normsq);
	- } else {
	- norm = 1.;
	- }
	- for (mu = 1; mu <= 3; ++mu) {
	- naxis[mu] = pnew[mu] / norm;
	-/* L150: */
	- }
	- }
	- return 0;
	-} /* pxtry_ */
	-
	-/* ******************************************************************** */
	-/* CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper */
	-/* CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper */
	-/* -- Author : */
	-/* +DECK,PXUVEC. */
	-
	-/* ** Routine to calculate unit vectors PU of all particles PP return
	- false on error*/
	-bool pxuvec(int ntrak, double pp, double pu) {
	-
	- /* Parameter adjustments */
	- pu -= 4;
	- pp -= 5;
	-
	- for (int n = 1; n <= ntrak; ++n) {
	- double mag = 0.0;
	- for ( int mu = 1; mu <= 3; ++mu)
	- mag += pp[mu + (n << 2)]*pp[mu + (n << 2)];
	- mag = sqrt(mag);
	- if (mag == 0.0 ) {
	- printf(" PXCONE: An input particle has zero mod(p)\n");
	- return false;
	- }
	- for (int mu = 1; mu <= 3; ++mu)
	- pu[mu + n * 3] = pp[mu + (n << 2)] / mag;
	- }
	- return true;
	-}
	-
	-
	-/* Set integer vector to zero */
	-void pxzeri(int n, int *a){
	- for (int i = 0; i < n; ++i) a[i] = 0;
	-}
	-
	-
	-/* This is a set of routines written by Mike Seymour to provide the */
	-/* services presumably normally provided by standard OPAL routines */
	-
	-
	-/* Set vector a to zero */
	-void pxzerv(int n, double *a) {
	- for (int i = 0; i < n; ++i) a[i] = 0.;
	-}
	-
	-/* add vectors c = a + b */
	-void pxaddv(int n, double a, double b, double *c) {
	- for (int i = 0; i < n; ++i) c[i] = a[i] + b[i];
	-}
	-
	-/* calculate angle between two vectors */
	-void pxang3(double a, double b, double &cost, double &thet) {
	-
	- cost = 1.0;
	- thet = 0.0;
	- double c = (a[0]a[0] + a[1]a[1] + a[2]a[2])
	- (b[0]b[0] + b[1]b[1] + b[2]*b[2]);
	- if (c <= 0.) return;
	-
	- c = 1/sqrt(c);
	- cost = (a[0]b[0] + a[1]b[1] + a[2]b[2])c;
	- thet = acos(cost);
	-
	-}
	-
	-
	-
	-/* ** Note that although JETLIS is assumed to be a 2d array, it */
	-/* ** it is used as 1d in this routine for efficiency */
	-/* ** Checks to see if TSTLIS entries correspond to a jet already found */
	-/* ** and entered in JETLIS */
	-int pxnew(int tstlis, int jetlis, int ntrak, int njet) {
	-
	- int match;
	- for (int i = 0; i < njet; ++i) {
	- match = true;
	- int in = i - 5000;
	- for (int n = 0; n < ntrak; ++n) {
	- in += 5000;
	- if (tstlis[n] != jetlis[in]) {
	- match = false;
	- break;
	- }
	- }
	- if (match) return false;
	- }
	- return true;
	-}
	-
	-/* ** Returns T if the first N elements of LIST1 are the same as the */
	-/* ** first N elements of LIST2. */
	-int pxsame_(int list1, int list2, int n) {
	- for (int i = 0; i < n; ++i)
	- if (list1[i] != list2[i])
	- return false;
	- return true;
	-}
	-
	-
	-/* ---RETURNS PHI, MOVED ONTO THE RANGE [-PI,PI) */
	-double pxmdpi(double phi) {
	-
	- if (phi <= M_PI) {
	- if (phi > -M_PI)
	- return abs(phi) < 1e-15? 0.0: phi;
	- else if (phi > -3*M_PI)
	- phi += 2*M_PI;
	- else
	- phi = -d_mod(M_PI - phi, 2*M_PI) + M_PI;
	- } else if (phi <= 3.0*M_PI) {
	- phi -= 2*M_PI;
	- } else {
	- phi = d_mod(phi + M_PI, 2*M_PI) - M_PI;
	- }
	-
	- return abs(phi) < 1e-15? 0.0: phi;
	-
	-}
	-
	diff --git a/src/Projections/pxcone.cc b/src/Projections/pxcone.cc
	new file mode 100644
	--- /dev/null
	+++ b/src/Projections/pxcone.cc
	@@ -0,0 +1,1405 @@
	+/* pxcone.f -- translated by f2c and hacked by Leif Lönnblad to avoid
	+ linking with libf2c.
	+*/
	+
	+#include <cmath>
	+#include <cstdio>
	+#include <cstdlib>
	+//#include "Rivet/Projections/pxcone.h"
	+using namespace std;
	+
	+#ifdef __cplusplus
	+extern "C" {
	+#endif
	+#include "f2c.h"
	+
	+/* Table of constant values, which are actually non const to be able
	+ to be used as fortran arguments. */
	+static int MAXV = 20000;
	+static int VDIM = 3;
	+static double twopi = 6.283185307;
	+
	+// The standard fortran SIGN function for doubles.
	+double d_sign(double a, double b) {
	+ return b < 0.0? -fabs(a): fabs(*a);
	+}
	+
	+// The standard fortran MOD function for doubles.
	+double d_mod(double * a, double * p) {
	+ return (a) - int((a)/(p))(*p);
	+}
	+
	+// The main PXCONE function.
	+void pxcone_(int mode, int ntrak, int *itkdm,
	+ double ptrak, double coner, double epslon, double
	+ ovlim, int mxjet, int njet, double pjet, int
	+ ipass, int ijmul, int ierr)
	+{
	+ /* Initialized data */
	+
	+ static int ncall = 0;
	+ static int nprint = 0;
	+ static double rold = 0.;
	+ static double epsold = 0.;
	+ static double ovold = 0.;
	+
	+ /* System generated locals */
	+ int ptrak_dim1, ptrak_offset, i__1, i__2;
	+ double d__1, d__2, d__3;
	+
	+ /* Local variables */
	+ static double cosr, rsep, ppsq, ptsq, cos2r;
	+ static int i__, j, n;
	+ static double vseed[3];
	+ static int iterr;
	+ extern /* Subroutine / int pxord_(double , int , int ,
	+ int , double );
	+ static int n1, n2;
	+ static double pj[20000] /* was [4][5000] /, pp[20000] / was [4][
	+ 5000] */;
	+ static int mu;
	+ static double pu[15000] /* was [3][5000] */, cosval;
	+ extern /* Subroutine / int pxaddv_(int , double , double ,
	+ double , int );
	+ static int jetlis[25000000] /* was [5000][5000] */;
	+ extern double pxmdpi_(double *);
	+ extern /* Subroutine / int pxsear_(int , double , int ,
	+ double , double , double , int , int *,
	+ double , int , int ), pxolap_(int , int *,
	+ int , int , double , double , double *);
	+ static int unstbl;
	+ extern /* Subroutine / int pxuvec_(int , double , double ,
	+ int ), pxzeri_(int , int ), pxnorv_(int ,
	+ double , double , int ), pxzerv_(int ,
	+ double *);
	+ static double vec1[3], vec2[3];
	+
	+/* .********************************************************* */
	+/* . ------ */
	+/* . PXCONE */
	+/* . ------ */
	+/* . */
	+/* . Code downloaded from the following web page */
	+/* . */
	+/* . http://aliceinfo.cern.ch/alicvs/viewvc/JETAN/pxcone.F?view=markup&pathrev=v4-05-04 */
	+/* . */
	+/* . on 17/10/2006 by G. Salam. Permission subsequently granted by Michael */
	+/* . H. Seymour (on behalf of the PxCone authors) for this code to be */
	+/* . distributed together with FastJet under the terms of the GNU Public */
	+/* . License v2 (see the file COPYING in the main FastJet directory). */
	+/* . */
	+/* .********** Pre Release Version 26.2.93 */
	+/* . */
	+/* . Driver for the Cone Jet finding algorithm of L.A. del Pozo. */
	+/* . Based on algorithm from D.E. Soper. */
	+/* . Finds jets inside cone of half angle CONER with energy > EPSLON. */
	+/* . Jets which receive more than a fraction OVLIM of their energy from */
	+/* . overlaps with other jets are excluded. */
	+/* . Output jets are ordered in energy. */
	+/* . If MODE.EQ.2 momenta are stored as (eta,phi,<empty>,pt) */
	+/* . Usage : */
	+/* . */
	+/* . INTEGER ITKDM,MXTRK */
	+/* . PARAMETER (ITKDM=4.or.more,MXTRK=1.or.more) */
	+/* . INTEGER MXJET, MXTRAK, MXPROT */
	+/* . PARAMETER (MXJET=10,MXTRAK=500,MXPROT=500) */
	+/* . INTEGER IPASS (MXTRAK),IJMUL (MXJET) */
	+/* . INTEGER NTRAK,NJET,IERR,MODE */
	+/* . DOUBLE PRECISION PTRAK (ITKDM,MXTRK),PJET (5,MXJET) */
	+/* . DOUBLE PRECISION CONER, EPSLON, OVLIM */
	+/* . NTRAK = 1.to.MXTRAK */
	+/* . CONER = ... */
	+/* . EPSLON = ... */
	+/* . OVLIM = ... */
	+/* . CALL PXCONE (MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,MXJET, */
	+/* . + NJET,PJET,IPASS,IJMUL,IERR) */
	+/* . */
	+/* . INPUT : MODE 1=>e+e-, 2=>hadron-hadron */
	+/* . INPUT : NTRAK Number of particles */
	+/* . INPUT : ITKDM First dimension of PTRAK array */
	+/* . INPUT : PTRAK Array of particle 4-momenta (Px,Py,Pz,E) */
	+/* . INPUT : CONER Cone size (half angle) in radians */
	+/* . INPUT : EPSLON Minimum Jet energy (GeV) */
	+/* . INPUT : OVLIM Maximum fraction of overlap energy in a jet */
	+/* . INPUT : MXJET Maximum possible number of jets */
	+/* . OUTPUT : NJET Number of jets found */
	+/* . OUTPUT : PJET 5-vectors of jets */
	+/* . OUTPUT : IPASS(k) Particle k belongs to jet number IPASS(k) */
	+/* . IPASS = -1 if not assosciated to a jet */
	+/* . OUTPUT : IJMUL(i) Jet i contains IJMUL(i) particles */
	+/* . OUTPUT : IERR = 0 if all is OK ; = -1 otherwise */
	+/* . */
	+/* . CALLS : PXSEAR, PXSAME, PXNEW, PXTRY, PXORD, PXUVEC, PXOLAP */
	+/* . CALLED : User */
	+/* . */
	+/* . AUTHOR : L.A. del Pozo */
	+/* . CREATED : 26-Feb-93 */
	+/* . LAST MOD : 2-Mar-93 */
	+/* . */
	+/* . Modification Log. */
	+/* . 25-Feb-07: G P Salam - fix bugs concerning 2pi periodicity in eta phi mode */
	+/* . - added commented code to get consistent behaviour */
	+/* . regardless of particle order (replaces n-way */
	+/* . midpoints with 2-way midpoints however...) */
	+/* . 2-Jan-97: M Wobisch - fix bug concerning COS2R in eta phi mode */
	+/* . 4-Apr-93: M H Seymour - Change 2d arrays to 1d in PXTRY & PXNEW */
	+/* . 2-Apr-93: M H Seymour - Major changes to add boost-invariant mode */
	+/* . 1-Apr-93: M H Seymour - Increase all array sizes */
	+/* . 30-Mar-93: M H Seymour - Change all REAL variables to DOUBLE PRECISION */
	+/* . 30-Mar-93: M H Seymour - Change OVLIM into an input parameter */
	+/* . 2-Mar-93: L A del Pozo - Fix Bugs in PXOLAP */
	+/* . 1-Mar-93: L A del Pozo - Remove Cern library routine calls */
	+/* . 1-Mar-93: L A del Pozo - Add Print out of welcome and R and Epsilon */
	+/* . */
	+/* .********************************************************* */
	+/* +SEQ,DECLARE. */
	+/* ** External Arrays */
	+/* ** Internal Arrays */
	+/* ** Used in the routine. */
	+/* MWobisch */
	+/* MWobisch */
	+ /* Parameter adjustments */
	+ --ipass;
	+ ptrak_dim1 = *itkdm;
	+ ptrak_offset = 1 + ptrak_dim1 * 1;
	+ ptrak -= ptrak_offset;
	+ --ijmul;
	+ pjet -= 6;
	+
	+ /* Function Body */
	+/* MWobisch */
	+/* *************************************** */
	+ rsep = 2.;
	+/* *************************************** */
	+/* MWobisch */
	+ *ierr = 0;
	+
	+/* ** INITIALIZE */
	+ if (ncall <= 0) {
	+ rold = (float)0.;
	+ epsold = (float)0.;
	+ ovold = (float)0.;
	+ }
	+ ++ncall;
	+
	+/* ** Print welcome and Jetfinder parameters */
	+ if ((coner != rold \|\| epslon != epsold \|\| *ovlim != ovold) && nprint <=
	+ 10) {
	+ printf("%s\n", " ********* PXCONE: Cone Jet-finder *********");
	+ printf("%s\n", " Written by Luis Del Pozo of OPAL");
	+ printf("%s\n", " Modified for eta-phi by Mike Seymour");
	+ printf("%s\n", " Includes bug fixes by Wobisch, Salam");
	+ printf("%s\n", " Translated to c(++) by Leif Lonnblad");
	+ printf("%s%5.2f%s\n", " Cone Size R = ",*coner," Radians");
	+ printf("%s%5.2f%s\n", " Min Jet energy Epsilon = ",*epslon," GeV");
	+ printf("%s%5.2f\n", " Overlap fraction parameter = ",*ovlim);
	+ printf("%s\n", " PXCONE is not a supported product and is");
	+ printf("%s\n", " is provided for comparative purposes only");
	+ printf("%s\n", " ***********************************************");
	+/* WRITE (6,) /
	+/* WRITE (6,) ' ******** PXCONE: Cone Jet-finder ********' /
	+/* WRITE (6,) ' Written by Luis Del Pozo of OPAL' /
	+/* WRITE (6,) ' Modified for eta-phi by Mike Seymour' /
	+/* WRITE (6,) ' Includes bug fixes by Wobisch, Salam' /
	+/* WRITE(6,1000)' Cone Size R = ',CONER,' Radians' */
	+/* WRITE(6,1001)' Min Jet energy Epsilon = ',EPSLON,' GeV' */
	+/* WRITE(6,1002)' Overlap fraction parameter = ',OVLIM */
	+/* WRITE (6,) ' PXCONE is not a supported product and is' /
	+/* WRITE (6,) ' is provided for comparative purposes only' /
	+/* WRITE (6,) ' *********************************************' /
	+/* MWobisch */
	+ if (rsep < (float)1.999) {
	+ printf("%s\n", " ******************************************");
	+ printf("%s\n", " ******************************************");
	+ printf("%s\n", " M Wobisch: private change !!!!!!!!!!!! ");
	+ printf("%s%5.2f\n", " Rsep is set to ",rsep);
	+ printf("%s\n", " this is ONLY meaningful in a NLO calculation");
	+ printf("%s\n", " ------------------------ ");
	+ printf("%s\n", " please check what you're doing!!");
	+ printf("%s\n", " ******************************************");
	+ printf("%s\n", " ******************************************");
	+ printf("%s\n", " ******************************************");
	+ printf("%s\n", "");
	+ printf("%s\n", "");
	+ printf("%s\n", "");
	+/* WRITE(6,) ' ' /
	+/* WRITE (6,) ' ****************************************' /
	+/* WRITE (6,) ' ****************************************' /
	+/* WRITE(6,) ' M Wobisch: private change !!!!!!!!!!!! ' /
	+/* WRITE(6,) ' Rsep is set to ',RSEP /
	+/* WRITE(6,) ' this is ONLY meaningful in a NLO calculation' /
	+/* WRITE(6,) ' ------------------------ ' /
	+/* WRITE(6,) ' please check what you''re doing!!' /
	+/* WRITE(6,) ' or ask: Markus.Wobisch@desy.de --' /
	+/* WRITE (6,) ' ****************************************' /
	+/* WRITE (6,) ' ****************************************' /
	+/* WRITE (6,) ' ****************************************' /
	+/* WRITE(6,) ' ' /
	+/* WRITE(6,) ' ' /
	+ }
	+/* MWobisch */
	+/* WRITE (6,) /
	+/* 1000 FORMAT(A18,F5.2,A10) */
	+/* 1001 FORMAT(A29,F5.2,A5) */
	+/* 1002 FORMAT(A33,F5.2) */
	+ ++nprint;
	+ rold = *coner;
	+ epsold = *epslon;
	+ ovold = *ovlim;
	+ }
	+
	+/* ** Copy calling array PTRAK to internal array PP(4,NTRAK) */
	+
	+ if (*ntrak > 5000) {
	+/* WRITE (6,) ' PXCONE: Ntrak too large: ',NTRAK /
	+ printf("%s%d\n", " PXCONE: Ntrak too large: ", *ntrak);
	+ *ierr = -1;
	+ return;
	+ }
	+ if (*mode != 2) {
	+ i__1 = *ntrak;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ for (j = 1; j <= 4; ++j) {
	+ pp[j + (i__ << 2) - 5] = ptrak[j + i__ * ptrak_dim1];
	+/* L101: */
	+ }
	+/* L100: */
	+ }
	+ } else {
	+/* ** Converting to eta,phi,pt if necessary */
	+ i__1 = *ntrak;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+/* Computing 2nd power */
	+ d__1 = ptrak[i__ * ptrak_dim1 + 1];
	+/* Computing 2nd power */
	+ d__2 = ptrak[i__ * ptrak_dim1 + 2];
	+ ptsq = d__1 * d__1 + d__2 * d__2;
	+/* Computing 2nd power */
	+ d__3 = ptrak[i__ * ptrak_dim1 + 3];
	+/* Computing 2nd power */
	+ d__2 = sqrt(ptsq + d__3 * d__3) + (d__1 = ptrak[i__ * ptrak_dim1
	+ + 3], abs(d__1));
	+ ppsq = d__2 * d__2;
	+ if (ptsq <= ppsq * (float)4.25e-18) {
	+ pp[(i__ << 2) - 4] = 20.;
	+ } else {
	+ pp[(i__ << 2) - 4] = log(ppsq / ptsq) * (float).5;
	+ }
	+ pp[(i__ << 2) - 4] = d_sign(&pp[(i__ << 2) - 4], &ptrak[i__ *
	+ ptrak_dim1 + 3]);
	+ if (ptsq == 0.) {
	+ pp[(i__ << 2) - 3] = 0.;
	+ } else {
	+ pp[(i__ << 2) - 3] = atan2(ptrak[i__ * ptrak_dim1 + 2], ptrak[
	+ i__ * ptrak_dim1 + 1]);
	+ }
	+ pp[(i__ << 2) - 2] = 0.;
	+ pp[(i__ << 2) - 1] = sqrt(ptsq);
	+ pu[i__ * 3 - 3] = pp[(i__ << 2) - 4];
	+ pu[i__ * 3 - 2] = pp[(i__ << 2) - 3];
	+ pu[i__ * 3 - 1] = pp[(i__ << 2) - 2];
	+/* L104: */
	+ }
	+ }
	+
	+/* ** Zero output variables */
	+
	+ *njet = 0;
	+ i__1 = *ntrak;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ for (j = 1; j <= 5000; ++j) {
	+ jetlis[j + i__ * 5000 - 5001] = false;
	+/* L103: */
	+ }
	+/* L102: */
	+ }
	+ pxzerv_(&MAXV, pj);
	+ pxzeri_(mxjet, &ijmul[1]);
	+
	+ if (*mode != 2) {
	+ cosr = cos(*coner);
	+ cos2r = cos(*coner);
	+ } else {
	+/* Purely for convenience, work in terms of 1-R2 */
	+/* Computing 2nd power */
	+ d__1 = *coner;
	+ cosr = 1 - d__1 * d__1;
	+/* MW -- select Rsep: 1-(RsepCONER)2 /
	+/* Computing 2nd power */
	+ d__1 = rsep * *coner;
	+ cos2r = 1 - d__1 * d__1;
	+/* ORIGINAL COS2R = 1-(2CONER)2 /
	+ }
	+ unstbl = false;
	+ if (*mode != 2) {
	+ pxuvec_(ntrak, pp, pu, ierr);
	+ if (*ierr != 0) {
	+ return;
	+ }
	+ }
	+/* ** Look for jets using particle diretions as seed axes */
	+
	+ i__1 = *ntrak;
	+ for (n = 1; n <= i__1; ++n) {
	+ for (mu = 1; mu <= 3; ++mu) {
	+ vseed[mu - 1] = pu[mu + n * 3 - 4];
	+/* L120: */
	+ }
	+ pxsear_(mode, &cosr, ntrak, pu, pp, vseed, njet, jetlis, pj, &unstbl,
	+ ierr);
	+ if (*ierr != 0) {
	+ return;
	+ }
	+/* L110: */
	+ }
	+/* MW - for Rsep=1 goto 145 */
	+/* GOTO 145 */
	+/* ** Now look between all pairs of jets as seed axes. */
	+/* NJTORG = NJET ! GPS -- to get consistent behaviour (2-way midpnts) */
	+/* DO 140 N1 = 1,NJTORG-1 ! GPS -- to get consistent behaviour (2-way midpnts) */
	+ i__1 = *njet - 1;
	+ for (n1 = 1; n1 <= i__1; ++n1) {
	+ vec1[0] = pj[(n1 << 2) - 4];
	+ vec1[1] = pj[(n1 << 2) - 3];
	+ vec1[2] = pj[(n1 << 2) - 2];
	+ if (*mode != 2) {
	+ pxnorv_(&VDIM, vec1, vec1, &iterr);
	+ }
	+/* DO 150 N2 = N1+1,NJTORG ! GPS -- to get consistent behaviour */
	+ i__2 = *njet;
	+ for (n2 = n1 + 1; n2 <= i__2; ++n2) {
	+ vec2[0] = pj[(n2 << 2) - 4];
	+ vec2[1] = pj[(n2 << 2) - 3];
	+ vec2[2] = pj[(n2 << 2) - 2];
	+ if (*mode != 2) {
	+ pxnorv_(&VDIM, vec2, vec2, &iterr);
	+ }
	+ pxaddv_(&VDIM, vec1, vec2, vseed, &iterr);
	+ if (*mode != 2) {
	+ pxnorv_(&VDIM, vseed, vseed, &iterr);
	+ } else {
	+ vseed[0] /= 2;
	+/* VSEED(2)=VSEED(2)/2 */
	+/* GPS 25/02/07 */
	+ d__2 = vec2[1] - vec1[1];
	+ d__1 = vec1[1] + pxmdpi_(&d__2) * .5;
	+ vseed[1] = pxmdpi_(&d__1);
	+ }
	+/* ---ONLY BOTHER IF THEY ARE BETWEEN 1 AND 2 CONE RADII APART */
	+ if (*mode != 2) {
	+ cosval = vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] *
	+ vec2[2];
	+ } else {
	+ if (abs(vec1[0]) >= 20. \|\| abs(vec2[0]) >= 20.) {
	+ cosval = -1e3;
	+ } else {
	+/* Computing 2nd power */
	+ d__1 = vec1[0] - vec2[0];
	+ d__3 = vec1[1] - vec2[1];
	+/* Computing 2nd power */
	+ d__2 = pxmdpi_(&d__3);
	+ cosval = 1 - (d__1 * d__1 + d__2 * d__2);
	+ }
	+ }
	+ if (cosval <= cosr && cosval >= cos2r) {
	+ pxsear_(mode, &cosr, ntrak, pu, pp, vseed, njet, jetlis, pj, &
	+ unstbl, ierr);
	+ }
	+/* CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET, */
	+/* + JETLIS,PJ,UNSTBL,IERR) */
	+ if (*ierr != 0) {
	+ return;
	+ }
	+/* L150: */
	+ }
	+/* L140: */
	+ }
	+ if (unstbl) {
	+ *ierr = -1;
	+/* WRITE (6,) ' PXCONE: Too many iterations to find a proto-jet' /
	+ printf(" PXCONE: Too many iterations to find a proto-jet\n");
	+ return;
	+ }
	+/* L145: */
	+/* ** Now put the jet list into order by jet energy, eliminating jets */
	+/* ** with energy less than EPSLON. */
	+ pxord_(epslon, njet, ntrak, jetlis, pj);
	+
	+/* ** Take care of jet overlaps */
	+ pxolap_(mode, njet, ntrak, jetlis, pj, pp, ovlim);
	+
	+/* ** Order jets again as some have been eliminated, or lost energy. */
	+ pxord_(epslon, njet, ntrak, jetlis, pj);
	+
	+/* ** All done!, Copy output into output arrays */
	+ if (njet > mxjet) {
	+/* WRITE (6,) ' PXCONE: Found more than MXJET jets' /
	+ printf(" PXCONE: Found more than MXJET jets\n");
	+ *ierr = -1;
	+ goto L99;
	+ }
	+ if (*mode != 2) {
	+ i__1 = *njet;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ for (j = 1; j <= 4; ++j) {
	+ pjet[j + i__ * 5] = pj[j + (i__ << 2) - 5];
	+/* L310: */
	+ }
	+/* L300: */
	+ }
	+ } else {
	+ i__1 = *njet;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ pjet[i__ * 5 + 1] = pj[(i__ << 2) - 1] * cos(pj[(i__ << 2) - 3]);
	+ pjet[i__ * 5 + 2] = pj[(i__ << 2) - 1] * sin(pj[(i__ << 2) - 3]);
	+ pjet[i__ * 5 + 3] = pj[(i__ << 2) - 1] * sinh(pj[(i__ << 2) - 4]);
	+ pjet[i__ * 5 + 4] = pj[(i__ << 2) - 1] * cosh(pj[(i__ << 2) - 4]);
	+/* L315: */
	+ }
	+ }
	+ i__1 = *ntrak;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ ipass[i__] = -1;
	+ i__2 = *njet;
	+ for (j = 1; j <= i__2; ++j) {
	+ if (jetlis[j + i__ * 5000 - 5001]) {
	+ ++ijmul[j];
	+ ipass[i__] = j;
	+ }
	+/* L330: */
	+ }
	+/* L320: */
	+ }
	+L99:
	+ return;
	+} /* pxcone_ */
	+
	+/* CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper */
	+/* -- Author : */
	+/* ----------------------------------------------------------------------- */
	+/* Subroutine / int pxnorv_(int n, double a, double b,
	+ int *iterr)
	+{
	+ /* System generated locals */
	+ int i__1;
	+ double d__1;
	+
	+ /* Builtin functions */
	+ double sqrt(double);
	+
	+ /* Local variables */
	+ static double c__;
	+ static int i__;
	+
	+ /* Parameter adjustments */
	+ --b;
	+ --a;
	+
	+ /* Function Body */
	+ c__ = 0.;
	+ i__1 = *n;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+/* Computing 2nd power */
	+ d__1 = a[i__];
	+ c__ += d__1 * d__1;
	+/* L10: */
	+ }
	+ if (c__ <= 0.) {
	+ return 0;
	+ }
	+ c__ = 1 / sqrt(c__);
	+ i__1 = *n;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ b[i__] = a[i__] * c__;
	+/* L20: */
	+ }
	+ return 0;
	+} /* pxnorv_ */
	+
	+/* CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper */
	+/* CMZ : 1.06/00 15/03/94 12.17.46 by P. Schleper */
	+/* -- Author : */
	+
	+/* +DECK,PXOLAP. */
	+/* Subroutine / int pxolap_(int mode, int njet, int ntrak,
	+ int jetlis, double pj, double pp, double ovlim)
	+{
	+ /* Initialized data */
	+
	+ static int ijmin = 0;
	+
	+ /* System generated locals */
	+ int i__1, i__2, i__3;
	+ double d__1, d__2, d__3;
	+
	+ /* Local variables */
	+ static int ijet[5000];
	+ static double thet, cost;
	+ static int i__, j, n;
	+ static double eover, thmin;
	+ static int iterr;
	+ extern /* Subroutine / int pxang3_(double , double *,
	+ double , double , int *);
	+ static int nj, mu;
	+ static int ovelap;
	+ extern double pxmdpi_(double *);
	+ static double vec1[3], vec2[3];
	+
	+
	+/* ** Looks for particles assigned to more than 1 jet, and reassigns them */
	+/* ** If more than a fraction OVLIM of a jet's energy is contained in */
	+/* ** higher energy jets, that jet is neglected. */
	+/* ** Particles assigned to the jet closest in angle (a la CDF, Snowmass). */
	+/* +SEQ,DECLARE. */
	+ /* Parameter adjustments */
	+ pp -= 5;
	+ pj -= 5;
	+ jetlis -= 5001;
	+
	+ /* Function Body */
	+
	+ if (*njet <= 1) {
	+ return 0;
	+ }
	+/* ** Look for jets with large overlaps with higher energy jets. */
	+ i__1 = *njet;
	+ for (i__ = 2; i__ <= i__1; ++i__) {
	+/* ** Find overlap energy between jets I and all higher energy jets. */
	+ eover = (float)0.;
	+ i__2 = *ntrak;
	+ for (n = 1; n <= i__2; ++n) {
	+ ovelap = false;
	+ i__3 = i__ - 1;
	+ for (j = 1; j <= i__3; ++j) {
	+ if (jetlis[i__ + n * 5000] && jetlis[j + n * 5000]) {
	+ ovelap = true;
	+ }
	+/* L120: */
	+ }
	+ if (ovelap) {
	+ eover += pp[(n << 2) + 4];
	+ }
	+/* L110: */
	+ }
	+/* ** Is the fraction of energy shared larger than OVLIM? */
	+ if (eover > ovlim pj[(i__ << 2) + 4]) {
	+/* ** De-assign all particles from Jet I */
	+ i__2 = *ntrak;
	+ for (n = 1; n <= i__2; ++n) {
	+ jetlis[i__ + n * 5000] = false;
	+/* L130: */
	+ }
	+ }
	+/* L100: */
	+ }
	+/* ** Now there are no big overlaps, assign every particle in */
	+/* ** more than 1 jet to the closet jet. */
	+/* ** Any particles now in more than 1 jet are assigned to the CLOSET */
	+/* ** jet (in angle). */
	+ i__1 = *ntrak;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ nj = 0;
	+ i__2 = *njet;
	+ for (j = 1; j <= i__2; ++j) {
	+ if (jetlis[j + i__ * 5000]) {
	+ ++nj;
	+ ijet[nj - 1] = j;
	+ }
	+/* L150: */
	+ }
	+ if (nj > 1) {
	+/* ** Particle in > 1 jet - calc angles... */
	+ vec1[0] = pp[(i__ << 2) + 1];
	+ vec1[1] = pp[(i__ << 2) + 2];
	+ vec1[2] = pp[(i__ << 2) + 3];
	+ thmin = (float)0.;
	+ i__2 = nj;
	+ for (j = 1; j <= i__2; ++j) {
	+ vec2[0] = pj[(ijet[j - 1] << 2) + 1];
	+ vec2[1] = pj[(ijet[j - 1] << 2) + 2];
	+ vec2[2] = pj[(ijet[j - 1] << 2) + 3];
	+ if (*mode != 2) {
	+ pxang3_(vec1, vec2, &cost, &thet, &iterr);
	+ } else {
	+/* Computing 2nd power */
	+ d__1 = vec1[0] - vec2[0];
	+ d__3 = vec1[1] - vec2[1];
	+/* Computing 2nd power */
	+ d__2 = pxmdpi_(&d__3);
	+ thet = d__1 * d__1 + d__2 * d__2;
	+ }
	+ if (j == 1) {
	+ thmin = thet;
	+ ijmin = ijet[j - 1];
	+ } else if (thet < thmin) {
	+ thmin = thet;
	+ ijmin = ijet[j - 1];
	+ }
	+/* L160: */
	+ }
	+/* ** Assign track to IJMIN */
	+ i__2 = *njet;
	+ for (j = 1; j <= i__2; ++j) {
	+ jetlis[j + i__ * 5000] = false;
	+/* L170: */
	+ }
	+ jetlis[ijmin + i__ * 5000] = true;
	+ }
	+/* L140: */
	+ }
	+/* ** Recompute PJ */
	+ i__1 = *njet;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ for (mu = 1; mu <= 4; ++mu) {
	+ pj[mu + (i__ << 2)] = (float)0.;
	+/* L210: */
	+ }
	+ i__2 = *ntrak;
	+ for (n = 1; n <= i__2; ++n) {
	+ if (jetlis[i__ + n * 5000]) {
	+ if (*mode != 2) {
	+ for (mu = 1; mu <= 4; ++mu) {
	+ pj[mu + (i__ << 2)] += pp[mu + (n << 2)];
	+/* L230: */
	+ }
	+ } else {
	+ pj[(i__ << 2) + 1] += pp[(n << 2) + 4] / (pp[(n << 2) + 4]
	+ + pj[(i__ << 2) + 4]) * (pp[(n << 2) + 1] - pj[(
	+ i__ << 2) + 1]);
	+/* GPS 25/02/07 */
	+ d__2 = pp[(n << 2) + 2] - pj[(i__ << 2) + 2];
	+ d__1 = pj[(i__ << 2) + 2] + pp[(n << 2) + 4] / (pp[(n <<
	+ 2) + 4] + pj[(i__ << 2) + 4]) * pxmdpi_(&d__2);
	+ pj[(i__ << 2) + 2] = pxmdpi_(&d__1);
	+/* PJ(2,I)=PJ(2,I) */
	+/* + + PP(4,N)/(PP(4,N)+PJ(4,I))PXMDPI(PP(2,N)-PJ(2,I)) /
	+ pj[(i__ << 2) + 4] += pp[(n << 2) + 4];
	+ }
	+ }
	+/* L220: */
	+ }
	+/* L200: */
	+ }
	+ return 0;
	+} /* pxolap_ */
	+
	+/* CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper */
	+/* CMZ : 1.06/00 14/03/94 15.37.45 by P. Schleper */
	+/* -- Author : */
	+
	+/* +DECK,PXORD. */
	+/* Subroutine / int pxord_(double epslon, int njet, int ntrak,
	+ int jetlis, double pj)
	+{
	+ /* System generated locals */
	+ int i__1, i__2;
	+
	+ /* Local variables */
	+ static int i__, j, index[5000];
	+ static double elist[5000], ptemp[20000] /* was [4][5000] */;
	+ static int logtmp[25000000] /* was [5000][5000] */;
	+ extern /* Subroutine / int pxsorv_(int , double , int ,
	+ char);
	+
	+
	+/* ** Routine to put jets into order and eliminate tose less than EPSLON */
	+/* +SEQ,DECLARE. */
	+/* ** Puts jets in order of energy: 1 = highest energy etc. */
	+/* ** Then Eliminate jets with energy below EPSLON */
	+
	+/* ** Copy input arrays. */
	+ /* Parameter adjustments */
	+ pj -= 5;
	+ jetlis -= 5001;
	+
	+ /* Function Body */
	+ i__1 = *njet;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ for (j = 1; j <= 4; ++j) {
	+ ptemp[j + (i__ << 2) - 5] = pj[j + (i__ << 2)];
	+/* L110: */
	+ }
	+ i__2 = *ntrak;
	+ for (j = 1; j <= i__2; ++j) {
	+ logtmp[i__ + j * 5000 - 5001] = jetlis[i__ + j * 5000];
	+/* L120: */
	+ }
	+/* L100: */
	+ }
	+ i__1 = *njet;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ elist[i__ - 1] = pj[(i__ << 2) + 4];
	+/* L150: */
	+ }
	+/* ** Sort the energies... */
	+ pxsorv_(njet, elist, index, 'I');
	+/* ** Fill PJ and JETLIS according to sort ( sort is in ascending order!!) */
	+ i__1 = *njet;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ for (j = 1; j <= 4; ++j) {
	+ pj[j + (i__ << 2)] = ptemp[j + (index[*njet + 1 - i__ - 1] << 2)
	+ - 5];
	+/* L210: */
	+ }
	+ i__2 = *ntrak;
	+ for (j = 1; j <= i__2; ++j) {
	+ jetlis[i__ + j * 5000] = logtmp[index[njet + 1 - i__ - 1] + j
	+ 5000 - 5001];
	+/* L220: */
	+ }
	+/* L200: */
	+ }
	+/* * Jets are now in order */
	+/* ** Now eliminate jets with less than Epsilon energy */
	+ i__1 = *njet;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ if (pj[(i__ << 2) + 4] < *epslon) {
	+ --(*njet);
	+ pj[(i__ << 2) + 4] = (float)0.;
	+ }
	+/* L300: */
	+ }
	+ return 0;
	+} /* pxord_ */
	+
	+/* ******************************************************************* */
	+/* CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper */
	+/* CMZ : 1.06/00 14/03/94 15.37.44 by P. Schleper */
	+/* -- Author : */
	+/* +DECK,PXSEAR. */
	+/* Subroutine / int pxsear_(int mode, double cosr, int ntrak,
	+ double pu, double pp, double vseed, int njet,
	+ int jetlis, double pj, int unstbl, int ierr)
	+{
	+ /* System generated locals */
	+ int i__1;
	+
	+ /* Local variables */
	+ static int iter;
	+ static double pnew[4];
	+ static int n;
	+ static double naxis[3], oaxis[3];
	+ extern int pxnew_(int , int , int , int );
	+ extern /* Subroutine / int pxtry_(int , double , int ,
	+ double , double , double , double ,
	+ double , int , int *);
	+ static int ok;
	+ static int mu;
	+ static int oldlis[5000];
	+ extern int pxsame_(int , int , int *);
	+ static int newlis[5000];
	+
	+
	+/* +SEQ,DECLARE. */
	+/* ** Using VSEED as a trial axis , look for a stable jet. */
	+/* ** Check stable jets against those already found and add to PJ. */
	+/* ** Will try up to MXITER iterations to get a stable set of particles */
	+/* ** in the cone. */
	+
	+ /* Parameter adjustments */
	+ pj -= 5;
	+ jetlis -= 5001;
	+ --vseed;
	+ pp -= 5;
	+ pu -= 4;
	+
	+ /* Function Body */
	+ for (mu = 1; mu <= 3; ++mu) {
	+ oaxis[mu - 1] = vseed[mu];
	+/* L100: */
	+ }
	+ i__1 = *ntrak;
	+ for (n = 1; n <= i__1; ++n) {
	+ oldlis[n - 1] = false;
	+/* L110: */
	+ }
	+ for (iter = 1; iter <= 30; ++iter) {
	+ pxtry_(mode, cosr, ntrak, &pu[4], &pp[5], oaxis, naxis, pnew, newlis,
	+ &ok);
	+/* ** Return immediately if there were no particles in the cone. */
	+ if (! ok) {
	+ return 0;
	+ }
	+ if (pxsame_(newlis, oldlis, ntrak)) {
	+/* ** We have a stable jet. */
	+ if (pxnew_(newlis, &jetlis[5001], ntrak, njet)) {
	+/* ** And the jet is a new one. So add it to our arrays. */
	+/* ** Check arrays are big anough... */
	+ if (*njet == 5000) {
	+/* WRITE (6,) ' PXCONE: Found more than MXPROT proto-jets' /
	+ printf(" PXCONE: Found more than MXPROT proto-jets\n");
	+ *ierr = -1;
	+ return 0;
	+ }
	+ ++(*njet);
	+ i__1 = *ntrak;
	+ for (n = 1; n <= i__1; ++n) {
	+ jetlis[njet + n 5000] = newlis[n - 1];
	+/* L130: */
	+ }
	+ for (mu = 1; mu <= 4; ++mu) {
	+ pj[mu + (*njet << 2)] = pnew[mu - 1];
	+/* L140: */
	+ }
	+ }
	+ return 0;
	+ }
	+/* ** The jet was not stable, so we iterate again */
	+ i__1 = *ntrak;
	+ for (n = 1; n <= i__1; ++n) {
	+ oldlis[n - 1] = newlis[n - 1];
	+/* L150: */
	+ }
	+ for (mu = 1; mu <= 3; ++mu) {
	+ oaxis[mu - 1] = naxis[mu - 1];
	+/* L160: */
	+ }
	+/* L120: */
	+ }
	+ *unstbl = true;
	+ return 0;
	+} /* pxsear_ */
	+
	+/* CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper */
	+/* -- Author : */
	+/* ----------------------------------------------------------------------- */
	+/* Subroutine / int pxsorv_(int n, double a, int k, char opt)
	+{
	+ /* System generated locals */
	+ int i__1;
	+
	+ /* Builtin functions */
	+ // /* Subroutine / int s_stop(char , ftnlen);
	+
	+ /* Local variables */
	+ static double b[5000];
	+ static int i__, j, il[5000], ir[5000];
	+
	+/* Sort A(N) into ascending order */
	+/* OPT = 'I' : return index array K only */
	+/* OTHERWISE : return sorted A and index array K */
	+/* ----------------------------------------------------------------------- */
	+
	+/* INT N,I,J,K(N),IL(NMAX),IR(NMAX) */
	+/* LUND */
	+
	+/* DOUBLE PRECISION A(N),B(NMAX) */
	+/* LUND */
	+ /* Parameter adjustments */
	+ --k;
	+ --a;
	+
	+ /* Function Body */
	+ if (*n > 5000) {
	+ // WRITE s_stop("Sorry, not enough room in Mike's PXSORV", (ftnlen)39);
	+ printf("Sorry, not enough room in Mike's PXSORV\n");
	+ abort();
	+ }
	+ il[0] = 0;
	+ ir[0] = 0;
	+ i__1 = *n;
	+ for (i__ = 2; i__ <= i__1; ++i__) {
	+ il[i__ - 1] = 0;
	+ ir[i__ - 1] = 0;
	+ j = 1;
	+L2:
	+ if (a[i__] > a[j]) {
	+ goto L5;
	+ }
	+/* L3: */
	+ if (il[j - 1] == 0) {
	+ goto L4;
	+ }
	+ j = il[j - 1];
	+ goto L2;
	+L4:
	+ ir[i__ - 1] = -j;
	+ il[j - 1] = i__;
	+ goto L10;
	+L5:
	+ if (ir[j - 1] <= 0) {
	+ goto L6;
	+ }
	+ j = ir[j - 1];
	+ goto L2;
	+L6:
	+ ir[i__ - 1] = ir[j - 1];
	+ ir[j - 1] = i__;
	+L10:
	+ ;
	+ }
	+ i__ = 1;
	+ j = 1;
	+ goto L8;
	+L20:
	+ j = il[j - 1];
	+L8:
	+ if (il[j - 1] > 0) {
	+ goto L20;
	+ }
	+L9:
	+ k[i__] = j;
	+ b[i__ - 1] = a[j];
	+ ++i__;
	+ if ((i__1 = ir[j - 1]) < 0) {
	+ goto L12;
	+ } else if (i__1 == 0) {
	+ goto L30;
	+ } else {
	+ goto L13;
	+ }
	+L13:
	+ j = ir[j - 1];
	+ goto L8;
	+L12:
	+ j = -ir[j - 1];
	+ goto L9;
	+L30:
	+ if ( opt == 'I') {
	+ return 0;
	+ }
	+ i__1 = *n;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+/* L31: */
	+ a[i__] = b[i__ - 1];
	+ }
	+/* L999: */
	+ return 0;
	+} /* pxsorv_ */
	+
	+/* ******************************************************************** */
	+/* CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper */
	+/* CMZ : 1.06/00 14/03/94 15.37.44 by P. Schleper */
	+/* -- Author : */
	+
	+/* +DECK,PXTRY. */
	+/* Subroutine / int pxtry_(int mode, double cosr, int ntrak,
	+ double pu, double pp, double oaxis, double naxis,
	+ double pnew, int newlis, int *ok)
	+{
	+ /* System generated locals */
	+ int i__1;
	+ double d__1, d__2, d__3;
	+
	+ /* Builtin functions */
	+ double sqrt(double);
	+
	+ /* Local variables */
	+ static double norm;
	+ static int n, mu;
	+ static double cosval;
	+ extern double pxmdpi_(double *);
	+ static double normsq;
	+ static int npp, npu;
	+
	+
	+/* +SEQ,DECLARE. */
	+/* ** Note that although PU and PP are assumed to be 2d arrays, they */
	+/* ** are used as 1d in this routine for efficiency */
	+/* ** Finds all particles in cone of size COSR about OAXIS direction. */
	+/* ** Calculates 4-momentum sum of all particles in cone (PNEW) , and */
	+/* ** returns this as new jet axis NAXIS (Both unit Vectors) */
	+
	+ /* Parameter adjustments */
	+ --newlis;
	+ --pnew;
	+ --naxis;
	+ --oaxis;
	+ --pp;
	+ --pu;
	+
	+ /* Function Body */
	+ *ok = false;
	+ for (mu = 1; mu <= 4; ++mu) {
	+ pnew[mu] = (float)0.;
	+/* L100: */
	+ }
	+ npu = -3;
	+ npp = -4;
	+ i__1 = *ntrak;
	+ for (n = 1; n <= i__1; ++n) {
	+ npu += 3;
	+ npp += 4;
	+ if (*mode != 2) {
	+ cosval = (float)0.;
	+ for (mu = 1; mu <= 3; ++mu) {
	+ cosval += oaxis[mu] * pu[mu + npu];
	+/* L120: */
	+ }
	+ } else {
	+ if ((d__1 = pu[npu + 1], abs(d__1)) >= 20. \|\| abs(oaxis[1]) >=
	+ 20.) {
	+ cosval = -1e3;
	+ } else {
	+/* Computing 2nd power */
	+ d__1 = oaxis[1] - pu[npu + 1];
	+ d__3 = oaxis[2] - pu[npu + 2];
	+/* Computing 2nd power */
	+ d__2 = pxmdpi_(&d__3);
	+ cosval = 1 - (d__1 * d__1 + d__2 * d__2);
	+ }
	+ }
	+ if (cosval >= *cosr) {
	+ newlis[n] = true;
	+ *ok = true;
	+ if (*mode != 2) {
	+ for (mu = 1; mu <= 4; ++mu) {
	+ pnew[mu] += pp[mu + npp];
	+/* L130: */
	+ }
	+ } else {
	+ pnew[1] += pp[npp + 4] / (pp[npp + 4] + pnew[4]) * (pp[npp +
	+ 1] - pnew[1]);
	+/* PNEW(2)=PNEW(2) */
	+/* + + PP(4+NPP)/(PP(4+NPP)+PNEW(4)) */
	+/* + PXMDPI(PP(2+NPP)-PNEW(2)) /
	+/* GPS 25/02/07 */
	+ d__2 = pp[npp + 2] - pnew[2];
	+ d__1 = pnew[2] + pp[npp + 4] / (pp[npp + 4] + pnew[4]) *
	+ pxmdpi_(&d__2);
	+ pnew[2] = pxmdpi_(&d__1);
	+ pnew[4] += pp[npp + 4];
	+ }
	+ } else {
	+ newlis[n] = false;
	+ }
	+/* L110: */
	+ }
	+/* ** If there are particles in the cone, calc new jet axis */
	+ if (*ok) {
	+ if (*mode != 2) {
	+ normsq = (float)0.;
	+ for (mu = 1; mu <= 3; ++mu) {
	+/* Computing 2nd power */
	+ d__1 = pnew[mu];
	+ normsq += d__1 * d__1;
	+/* L140: */
	+ }
	+ norm = sqrt(normsq);
	+ } else {
	+ norm = 1.;
	+ }
	+ for (mu = 1; mu <= 3; ++mu) {
	+ naxis[mu] = pnew[mu] / norm;
	+/* L150: */
	+ }
	+ }
	+ return 0;
	+} /* pxtry_ */
	+
	+/* ******************************************************************** */
	+/* CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper */
	+/* CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper */
	+/* -- Author : */
	+/* +DECK,PXUVEC. */
	+
	+/* Subroutine / int pxuvec_(int ntrak, double pp, double pu,
	+ int *ierr)
	+{
	+ /* System generated locals */
	+ int i__1;
	+ double d__1;
	+
	+ /* Builtin functions */
	+ double sqrt(double);
	+
	+ /* Local variables */
	+ static int n, mu;
	+ static double mag;
	+
	+
	+/* ** Routine to calculate unit vectors PU of all particles PP */
	+/* +SEQ,DECLARE. */
	+ /* Parameter adjustments */
	+ pu -= 4;
	+ pp -= 5;
	+
	+ /* Function Body */
	+ i__1 = *ntrak;
	+ for (n = 1; n <= i__1; ++n) {
	+ mag = (float)0.;
	+ for (mu = 1; mu <= 3; ++mu) {
	+/* Computing 2nd power */
	+ d__1 = pp[mu + (n << 2)];
	+ mag += d__1 * d__1;
	+/* L110: */
	+ }
	+ mag = sqrt(mag);
	+ if (mag == (float)0.) {
	+/* WRITE(6,)' PXCONE: An input particle has zero mod(p)' /
	+ printf(" PXCONE: An input particle has zero mod(p)\n");
	+ *ierr = -1;
	+ return 0;
	+ }
	+ for (mu = 1; mu <= 3; ++mu) {
	+ pu[mu + n * 3] = pp[mu + (n << 2)] / mag;
	+/* L120: */
	+ }
	+/* L100: */
	+ }
	+ return 0;
	+} /* pxuvec_ */
	+
	+/* CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper */
	+/* -- Author : */
	+/* ----------------------------------------------------------------------- */
	+/* Subroutine / int pxzeri_(int n, int *a)
	+{
	+ /* System generated locals */
	+ int i__1;
	+
	+ /* Local variables */
	+ static int i__;
	+
	+ /* Parameter adjustments */
	+ --a;
	+
	+ /* Function Body */
	+ i__1 = *n;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ a[i__] = 0;
	+/* L10: */
	+ }
	+ return 0;
	+} /* pxzeri_ */
	+
	+/* CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper */
	+/* -- Author : */
	+/* ----------------------------------------------------------------------- */
	+/* This is a set of routines written by Mike Seymour to provide the */
	+/* services presumably normally provided by standard OPAL routines */
	+/* PXZERV zeroes a vector */
	+/* PXZERI zeroes a vector of integers */
	+/* PXNORV normalizes a vector */
	+/* PXADDV adds two vectors */
	+/* PXSORV sorts a vector (copied from HERWIG) */
	+/* PXANG3 finds the angle (and its cosine) between two vectors */
	+/* PXMDPI moves its argument onto the range [-pi,pi) */
	+/* ----------------------------------------------------------------------- */
	+/* Subroutine / int pxzerv_(int n, double *a)
	+{
	+ /* System generated locals */
	+ int i__1;
	+
	+ /* Local variables */
	+ static int i__;
	+
	+ /* Parameter adjustments */
	+ --a;
	+
	+ /* Function Body */
	+ i__1 = *n;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ a[i__] = 0.;
	+/* L10: */
	+ }
	+ return 0;
	+} /* pxzerv_ */
	+
	+/* -- Author : */
	+/* ----------------------------------------------------------------------- */
	+/* Subroutine / int pxaddv_(int n, double a, double b,
	+ double c__, int iterr)
	+{
	+ /* System generated locals */
	+ int i__1;
	+
	+ /* Local variables */
	+ static int i__;
	+
	+ /* Parameter adjustments */
	+ --c__;
	+ --b;
	+ --a;
	+
	+ /* Function Body */
	+ i__1 = *n;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ c__[i__] = a[i__] + b[i__];
	+/* L10: */
	+ }
	+ return 0;
	+} /* pxaddv_ */
	+
	+/* CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper */
	+/* -- Author : */
	+/* ----------------------------------------------------------------------- */
	+/* Subroutine / int pxang3_(double a, double b, double cost,
	+ double thet, int iterr)
	+{
	+ /* System generated locals */
	+ double d__1, d__2, d__3, d__4, d__5, d__6;
	+
	+ /* Builtin functions */
	+ double sqrt(double), acos(double);
	+
	+ /* Local variables */
	+ static double c__;
	+
	+ /* Parameter adjustments */
	+ --b;
	+ --a;
	+
	+ /* Function Body */
	+/* Computing 2nd power */
	+ d__1 = a[1];
	+/* Computing 2nd power */
	+ d__2 = a[2];
	+/* Computing 2nd power */
	+ d__3 = a[3];
	+/* Computing 2nd power */
	+ d__4 = b[1];
	+/* Computing 2nd power */
	+ d__5 = b[2];
	+/* Computing 2nd power */
	+ d__6 = b[3];
	+ c__ = (d__1 * d__1 + d__2 * d__2 + d__3 * d__3) * (d__4 * d__4 + d__5 *
	+ d__5 + d__6 * d__6);
	+ if (c__ <= 0.) {
	+ return 0;
	+ }
	+ c__ = 1 / sqrt(c__);
	+ cost = (a[1] b[1] + a[2] * b[2] + a[3] * b[3]) * c__;
	+ thet = acos(cost);
	+ return 0;
	+} /* pxang3_ */
	+
	+/* CMZ : 1.06/00 14/03/94 15.41.57 by P. Schleper */
	+/* -- Author : P. Schleper 28/02/94 */
	+int pxnew_(int tstlis, int jetlis, int ntrak, int
	+ njet)
	+{
	+ /* System generated locals */
	+ int i__1, i__2;
	+ int ret_val;
	+
	+ /* Local variables */
	+ static int i__, n;
	+ static int match;
	+ static int in;
	+
	+
	+/* ** Note that although JETLIS is assumed to be a 2d array, it */
	+/* ** it is used as 1d in this routine for efficiency */
	+/* ** Checks to see if TSTLIS entries correspond to a jet already found */
	+/* ** and entered in JETLIS */
	+
	+ /* Parameter adjustments */
	+ --jetlis;
	+ --tstlis;
	+
	+ /* Function Body */
	+ ret_val = true;
	+ i__1 = *njet;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ match = true;
	+ in = i__ - 5000;
	+ i__2 = *ntrak;
	+ for (n = 1; n <= i__2; ++n) {
	+ in += 5000;
	+ if (tstlis[n] != jetlis[in]) {
	+ match = false;
	+ goto L100;
	+ }
	+/* L110: */
	+ }
	+ if (match) {
	+ ret_val = false;
	+ return ret_val;
	+ }
	+L100:
	+ ;
	+ }
	+ return ret_val;
	+} /* pxnew_ */
	+
	+/* CMZ : 1.06/00 14/03/94 15.41.57 by P. Schleper */
	+/* -- Author : P. Schleper 28/02/94 */
	+int pxsame_(int list1, int list2, int *n)
	+{
	+ /* System generated locals */
	+ int i__1;
	+ int ret_val;
	+
	+ /* Local variables */
	+ static int i__;
	+
	+
	+/* ** Returns T if the first N elements of LIST1 are the same as the */
	+/* ** first N elements of LIST2. */
	+
	+ /* Parameter adjustments */
	+ --list2;
	+ --list1;
	+
	+ /* Function Body */
	+ ret_val = true;
	+ i__1 = *n;
	+ for (i__ = 1; i__ <= i__1; ++i__) {
	+ if (list1[i__] != list2[i__]) {
	+ ret_val = false;
	+ return ret_val;
	+ }
	+/* L100: */
	+ }
	+ return ret_val;
	+} /* pxsame_ */
	+
	+/* CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper */
	+/* -- Author : */
	+/* ----------------------------------------------------------------------- */
	+double pxmdpi_(double *phi)
	+{
	+ /* System generated locals */
	+ double ret_val, d__1;
	+
	+
	+/* ---RETURNS PHI, MOVED ONTO THE RANGE [-PI,PI) */
	+ ret_val = *phi;
	+ if (ret_val <= 3.141592654) {
	+ if (ret_val > -3.141592654) {
	+ goto L100;
	+ } else if (ret_val > -9.424777961) {
	+ ret_val += 6.283185307;
	+ } else {
	+ d__1 = 3.141592654 - ret_val;
	+ ret_val = -d_mod(&d__1, &twopi) + 3.141592654;
	+ }
	+ } else if (ret_val <= 9.424777961) {
	+ ret_val += -6.283185307;
	+ } else {
	+ d__1 = ret_val + 3.141592654;
	+ ret_val = d_mod(&d__1, &twopi) - 3.141592654;
	+ }
	+L100:
	+ if (abs(ret_val) < 1e-15) {
	+ ret_val = 0.;
	+ }
	+ return ret_val;
	+} /* pxmdpi_ */
	+
	+#ifdef __cplusplus
	+// }
	+#endif
	+
	+int main() {
	+ int mode = 0;
	+ int ntrak = 0;
	+ int itkdm = 0;
	+ int mxjet = 0;
	+ int njet = 0;
	+ int ipass = 0;
	+ int ijmul = 0;
	+ int ierr = 0;
	+ double ptrak = 0.0;
	+ double coner = 0.0;
	+ double epslon = 0.0;
	+ double ovlim = 0.0;
	+ double pjet = 0.0;
	+
	+ pxcone_(&mode, &ntrak, &itkdm, &ptrak, &coner, &epslon, &ovlim,
	+ &mxjet, &njet, &pjet, &ipass, &ijmul, &ierr);
	+}
	+#ifdef __cplusplus
	+ }
	+#endif

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 5:47 PM (1 d, 17 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805449
Default Alt Text: (97 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions