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,173 +1,172 @@
//FJSTARTHEADER
// $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 .
//----------------------------------------------------------------------
//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
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
// actual physical parameters:
//
// coner
// epsilon
// ovlim
-// extern "C" {
void pxcone_
(
int mode , // 1=>e+e-, 2=>hadron-hadron
int ntrak , // Number of particles
int itkdm , // First dimension of PTRAK array:
const double * ptrak , // Array of particle 4-momenta (Px,Py,Pz,E)
double coner , // Cone size (half angle) in radians
double epslon , // Minimum Jet energy (GeV)
double ovlim , // Maximum fraction of overlap energy in a jet
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
);
+//FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
}
#endif
diff --git a/src/Projections/PxConePlugin.cc b/src/Projections/PxConePlugin.cc
--- a/src/Projections/PxConePlugin.cc
+++ b/src/Projections/PxConePlugin.cc
@@ -1,1555 +1,1221 @@
//FJSTARTHEADER
// $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 .
//----------------------------------------------------------------------
//FJENDHEADER
#include
#include
#include
#include "Rivet/Projections/PxConePlugin.hh"
#include "fastjet/ClusterSequence.hh"
#include
// pxcone stuff
// #include "Rivet/Projections/pxcone.h"
// FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
namespace Rivet {
using namespace std;
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(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[ntrak*4+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[mxjet*5+1];
int *ipass = new int[ntrak+1];
int *ijmul = new int[mxjet+1];
int ierr;
// 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)
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
);
if (ierr != 0) throw fastjet::Error("An error occurred while running PXCONE");
// now transfer information back
valarray last_index_created(njet);
vector > 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 & 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,
fastjet::PseudoJet(pjet[5*ipxjet+0],pjet[5*ipxjet+1],
pjet[5*ipxjet+2],pjet[5*ipxjet+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 ourjets(clust_seq.inclusive_jets());
//for (vector::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 "Rivet/Projections/pxcone.h"
using namespace std;
-#ifdef __cplusplus
-// extern "C" {
-#endif
-
/* 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;
+void pxtry_(int, double *, int, double *, double *, double *, double *,
+ double *, int *, int *);
+
+void pxsorv_(int, double *, int *, char);
+
+void pxsear_(int, double *, int, double *, double *, double *, int &, int *,
+ double *, int *, int *);
+
+void pxolap_(int, int, int, int *, double *, double *, double);
+
+void pxnorv_(int *, double *, double *, int *);
+
// The standard fortran SIGN function for doubles.
-double d_sign(double a, double b) {
+inline 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) {
+inline double d_mod(double a, double p) {
return a - int(a/p)*p;
}
+/* ---RETURNS PHI, MOVED ONTO THE RANGE [-PI,PI) */
+inline double pxmdpi(double phi) {
+ while ( phi <= -M_PI ) phi += 2*M_PI;
+ while ( phi > M_PI ) phi -= 2*M_PI;
+ return abs(phi) < 1e-15? 0.0: 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;
+
+// }
+
+/* Set integer vector a to zero */
+inline void pxzeri(int n, int *a){
+ for (int i = 0; i < n; ++i) a[i] = 0;
+}
+
+/* Set vector a to zero */
+inline void pxzerv(int n, double *a) {
+ for (int i = 0; i < n; ++i) a[i] = 0.;
+}
+
+/* add vectors c = a + b */
+inline void pxaddv(int n, double *a, double *b, double *c) {
+ for (int i = 0; i < n; ++i) c[i] = a[i] + b[i];
+}
+
+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;
+}
+
+/* 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. */
+bool pxsame(int *list1, int *list2, int n) {
+ for (int i = 0; i < n; ++i)
+ if (list1[i] != list2[i])
+ return false;
+ return true;
+}
+
+/* ** Routine to put jets into order and eliminate tose less than EPSLON */
+/* ** Puts jets in order of energy: 1 = highest energy etc. */
+/* ** Then Eliminate jets with energy below EPSLON */
+void pxord(double epslon, int & njet, int ntrak,
+ int *jetlis, double *pj)
+{
+ /* Local variables */
+ static int index[5000];
+ static double elist[5000], ptemp[20000] /* was [4][5000] */;
+ static int logtmp[25000000] /* was [5000][5000] */;
+
+
+
+/* ** Copy input arrays. */
+ /* Parameter adjustments */
+ pj -= 5;
+ jetlis -= 5001;
+
+ /* Function Body */
+ for (int i = 1; i <= njet; ++i) {
+ for (int j = 1; j <= 4; ++j) {
+ ptemp[j + (i << 2) - 5] = pj[j + (i << 2)];
+ }
+ for (int j = 1; j <= ntrak; ++j) {
+ logtmp[i + j * 5000 - 5001] = jetlis[i + j * 5000];
+ }
+ }
+ for (int i = 1; i <= njet; ++i) {
+ elist[i - 1] = pj[(i << 2) + 4];
+ }
+
+/* ** Sort the energies... */
+ pxsorv_(njet, elist, index, 'I');
+/* ** Fill PJ and JETLIS according to sort ( sort is in ascending order!!) */
+ for (int i = 1; i <= njet; ++i) {
+ for (int j = 1; j <= 4; ++j) {
+ pj[j + (i << 2)] = ptemp[j + (index[njet + 1 - i - 1] << 2) - 5];
+ }
+ for (int j = 1; j <= ntrak; ++j) {
+ jetlis[i + j * 5000] =
+ logtmp[index[njet + 1 - i - 1] + j * 5000 - 5001];
+ }
+ }
+/* * Jets are now in order */
+/* ** Now eliminate jets with less than Epsilon energy */
+ int nold = njet;
+ for (int i = 1; i <= nold; ++i) {
+ if (pj[(i << 2) + 4] < epslon) {
+ --njet;
+ pj[(i << 2) + 4] = 0.0;
+ }
+ }
+}
+
// The main PXCONE function.
void pxcone_(int mode, int ntrak, int itkdm,
const 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.;
+ static set rold;
+ static set epsold;
+ static set ovold;
/* 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,,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) {
+ if ((rold.find(coner) == rold.end() ||
+ epsold.find(epslon) == epsold.end() ||
+ ovold.find(ovlim) == ovold.end()) ) {
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%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;
+
+ rold.insert(coner);
+ epsold.insert(epslon);
+ ovold.insert(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]);
+ 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-R**2 */
/* Computing 2nd power */
d__1 = coner;
cosr = 1 - d__1 * d__1;
/* MW -- select Rsep: 1-(Rsep*CONER)**2 */
/* Computing 2nd power */
d__1 = rsep * coner;
cos2r = 1 - d__1 * d__1;
/* ORIGINAL COS2R = 1-(2*CONER)**2 */
}
unstbl = false;
if (mode != 2) {
- pxuvec_(ntrak, pp, pu, ierr);
- if (*ierr != 0) {
- return;
- }
+ 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, &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);
+ 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);
+ 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);
+ 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);
+ 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)
+
+void pxnorv_(int *n, double *a, double *b, int *iterr)
{
/* System generated locals */
int i__1;
double d__1;
/* Builtin functions */
/* 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;
+ return;
}
c__ = 1 / sqrt(c__);
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
b[i__] = a[i__] * c__;
-/* L20: */
}
- return 0;
+ return;
} /* 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,
+void 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;
+ return;
}
/* ** 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);
+ 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);
+ 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);
+ 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;
+ return;
} /* 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)
+void 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;
+ return;
}
- if (pxsame_(newlis, oldlis, ntrak)) {
+ if (pxsame(newlis, oldlis, ntrak)) {
/* ** We have a stable jet. */
- if (pxnew_(newlis, &jetlis[5001], ntrak, njet)) {
+ 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;
+ return;
}
++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;
+ return;
}
/* ** 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;
+ return;
} /* 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)
+
+void 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;
+ return;
}
i__1 = n;
for (i__ = 1; i__ <= i__1; ++i__) {
-/* L31: */
a[i__] = b[i__ - 1];
}
-/* L999: */
- return 0;
+ return;
} /* 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,
+void 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 */
/* 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);
+ 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);
+ 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 */
-
- /* 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 */
-
- /* 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_ */
-
-/* ---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;
}
-
-
-#ifdef __cplusplus
-// }
-#endif
-}