diff --git a/include/Rivet/Makefile.am b/include/Rivet/Makefile.am --- a/include/Rivet/Makefile.am +++ b/include/Rivet/Makefile.am @@ -1,186 +1,187 @@ ## 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/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/FastJets.hh b/include/Rivet/Projections/FastJets.hh --- a/include/Rivet/Projections/FastJets.hh +++ b/include/Rivet/Projections/FastJets.hh @@ -1,292 +1,292 @@ // -*- C++ -*- #ifndef RIVET_FastJets_HH #define RIVET_FastJets_HH #include "Rivet/Jet.hh" #include "Rivet/Particle.hh" #include "Rivet/Projection.hh" #include "Rivet/Projections/JetAlg.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Tools/RivetFastJet.hh" #include "fastjet/SISConePlugin.hh" #include "fastjet/ATLASConePlugin.hh" #include "fastjet/CMSIterativeConePlugin.hh" #include "fastjet/CDFJetCluPlugin.hh" #include "fastjet/CDFMidPointPlugin.hh" #include "fastjet/D0RunIIConePlugin.hh" #include "fastjet/TrackJetPlugin.hh" #include "fastjet/JadePlugin.hh" +#include "Rivet/Projections/PxConePlugin.hh" //#include "fastjet/PxConePlugin.hh" namespace Rivet { /// Project out jets found using the FastJet package jet algorithms. class FastJets : public JetAlg { public: /// Wrapper enum for selected FastJet jet algorithms. /// @todo Move to JetAlg and alias here? - enum JetAlgName { KT, CAM, SISCONE, ANTIKT, - // PXCONE, + enum JetAlgName { KT, CAM, SISCONE, ANTIKT, PXCONE, ATLASCONE, CMSCONE, CDFJETCLU, CDFMIDPOINT, D0ILCONE, JADE, DURHAM, TRACKJET, GENKTEE }; /// @name Constructors etc. //@{ /// Constructor from a FastJet JetDefinition /// /// @warning The AreaDefinition pointer must be heap-allocated: it will be stored/deleted via a shared_ptr. FastJets(const FinalState& fsp, const fastjet::JetDefinition& jdef, JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS, JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES, fastjet::AreaDefinition* adef=nullptr) : JetAlg(fsp, usemuons, useinvis), _jdef(jdef), _adef(adef) { _initBase(); } /// JetDefinition-based constructor with reordered args for easier specification of jet area definition /// /// @warning The AreaDefinition pointer must be heap-allocated: it will be stored/deleted via a shared_ptr. FastJets(const FinalState& fsp, const fastjet::JetDefinition& jdef, fastjet::AreaDefinition* adef, JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS, JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES) : FastJets(fsp, jdef, usemuons, useinvis, adef) { } /// Native argument constructor, using FastJet alg/scheme enums. /// /// @warning The AreaDefinition pointer must be heap-allocated: it will be stored/deleted via a shared_ptr. FastJets(const FinalState& fsp, fastjet::JetAlgorithm type, fastjet::RecombinationScheme recom, double rparameter, JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS, JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES, fastjet::AreaDefinition* adef=nullptr) : FastJets(fsp, fastjet::JetDefinition(type, rparameter, recom), usemuons, useinvis, adef) { } /// Native argument constructor with reordered args for easier specification of jet area definition /// /// @warning The AreaDefinition pointer must be heap-allocated: it will be stored/deleted via a shared_ptr. FastJets(const FinalState& fsp, fastjet::JetAlgorithm type, fastjet::RecombinationScheme recom, double rparameter, fastjet::AreaDefinition* adef, JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS, JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES) : FastJets(fsp, type, recom, rparameter, usemuons, useinvis, adef) { } /// @brief Explicitly pass in an externally-constructed plugin /// /// @warning Provided plugin and area definition pointers must be heap-allocated; Rivet will store/delete via a shared_ptr FastJets(const FinalState& fsp, fastjet::JetDefinition::Plugin* plugin, JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS, JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES, fastjet::AreaDefinition* adef=nullptr) : FastJets(fsp, fastjet::JetDefinition(plugin), usemuons, useinvis, adef) { _plugin.reset(plugin); } /// @brief Explicitly pass in an externally-constructed plugin, with reordered args for easier specification of jet area definition /// /// @warning Provided plugin and area definition pointers must be heap-allocated; Rivet will store/delete via a shared_ptr FastJets(const FinalState& fsp, fastjet::JetDefinition::Plugin* plugin, fastjet::AreaDefinition* adef, JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS, JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES) : FastJets(fsp, plugin, usemuons, useinvis, adef) { } /// @brief Convenience constructor using Rivet enums for most common jet algs (including some plugins). /// /// For the built-in algs, E-scheme recombination is used. For full control /// of FastJet built-in jet algs, use the constructors from native-args or a /// plugin pointer. /// /// @warning Provided area definition pointer must be heap-allocated; Rivet will store/delete via a shared_ptr FastJets(const FinalState& fsp, JetAlgName alg, double rparameter, JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS, JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES, fastjet::AreaDefinition* adef=nullptr, double seed_threshold=1.0) : JetAlg(fsp, usemuons, useinvis) { _initBase(); _initJdef(alg, rparameter, seed_threshold); } // /// Same thing as above, but without an FS (for when we want to pass the particles directly to the calc method) // /// @todo Does this work properly, without internal HeavyQuarks etc.? // FastJets(JetAlgName alg, double rparameter, double seed_threshold=1.0) { _initJdef(alg, rparameter, seed_threshold); } // /// Same thing as above, but without an FS (for when we want to pass the particles directly to the calc method) // /// @todo Does this work properly, without internal HeavyQuarks etc.? // FastJets(fastjet::JetAlgorithm type, fastjet::RecombinationScheme recom, double rparameter) { _initJdef(type, recom, rparameter); } // /// Same thing as above, but without an FS (for when we want to pass the particles directly to the calc method) // /// @todo Does this work properly, without internal HeavyQuarks etc.? // FastJets(fastjet::JetDefinition::Plugin* plugin) : _jdef(plugin), _plugin(plugin) { // // _plugin.reset(plugin); // // _jdef = fastjet::JetDefinition(plugin); // } /// Clone on the heap. DEFAULT_RIVET_PROJ_CLONE(FastJets); //@} /// @name Static helper functions for FastJet interaction, with tagging //@{ /// Make PseudoJets for input to a ClusterSequence, with user_index codes for constituent- and tag-particle linking static PseudoJets mkClusterInputs(const Particles& fsparticles, const Particles& tagparticles=Particles()); /// Make a Rivet Jet from a PseudoJet holding a user_index code for lookup of Rivet fsparticle or tagparticle links static Jet mkJet(const PseudoJet& pj, const Particles& fsparticles, const Particles& tagparticles=Particles()); /// Convert a whole list of PseudoJets to a list of Jets, with mkJet-style unpacking static Jets mkJets(const PseudoJets& pjs, const Particles& fsparticles, const Particles& tagparticles=Particles()); //@} /// Reset the projection. Jet def, etc. are unchanged. void reset(); /// @brief Use provided jet area definition /// /// @warning The provided pointer must be heap-allocated: it will be stored/deleted via a shared_ptr. /// @note Provide an adef null pointer to re-disable jet area calculation void useJetArea(fastjet::AreaDefinition* adef) { _adef.reset(adef); } /// @name Access to the jets //@{ /// Get the jets (unordered) with pT > ptmin. Jets _jets() const; /// Get the pseudo jets (unordered). PseudoJets pseudoJets(double ptmin=0.0) const; /// Alias PseudoJets pseudojets(double ptmin=0.0) const { return pseudoJets(ptmin); } /// Get the pseudo jets, ordered by \f$ p_T \f$. PseudoJets pseudoJetsByPt(double ptmin=0.0) const { return sorted_by_pt(pseudoJets(ptmin)); } /// Alias PseudoJets pseudojetsByPt(double ptmin=0.0) const { return pseudoJetsByPt(ptmin); } /// Get the pseudo jets, ordered by \f$ E \f$. PseudoJets pseudoJetsByE(double ptmin=0.0) const { return sorted_by_E(pseudoJets(ptmin)); } /// Alias PseudoJets pseudojetsByE(double ptmin=0.0) const { return pseudoJetsByE(ptmin); } /// Get the pseudo jets, ordered by rapidity. PseudoJets pseudoJetsByRapidity(double ptmin=0.0) const { return sorted_by_rapidity(pseudoJets(ptmin)); } /// Alias PseudoJets pseudojetsByRapidity(double ptmin=0.0) const { return pseudoJetsByRapidity(ptmin); } /// Trim (filter) a jet, keeping tag and constituent info in the resulting jet Jet trimJet(const Jet& input, const fastjet::Filter& trimmer) const; //@} /// @name Access to the FastJet clustering objects such as jet def, area def, and cluster //@{ /// Return the cluster sequence. /// @todo Care needed re. const shared_ptr vs. shared_ptr const shared_ptr clusterSeq() const { return _cseq; } /// Return the area-enabled cluster sequence (if an area defn exists, otherwise returns a null ptr). /// @todo Care needed re. const shared_ptr vs. shared_ptr const shared_ptr clusterSeqArea() const { return areaDef() ? dynamic_pointer_cast(_cseq) : nullptr; } /// Return the jet definition. const fastjet::JetDefinition& jetDef() const { return _jdef; } /// @brief Return the area definition. /// /// @warning May be null! /// @todo Care needed re. const shared_ptr vs. shared_ptr const shared_ptr areaDef() const { return _adef; } //@} private: /// Shared utility functions to implement constructor behaviour void _initBase(); void _initJdef(JetAlgName alg, double rparameter, double seed_threshold); protected: /// Perform the projection on the Event. void project(const Event& e); /// Compare projections. int compare(const Projection& p) const; public: /// Do the calculation locally (no caching). void calc(const Particles& fsparticles, const Particles& tagparticles=Particles()); private: /// Jet definition fastjet::JetDefinition _jdef; /// Pointer to user-handled area definition std::shared_ptr _adef; /// Cluster sequence std::shared_ptr _cseq; /// FastJet external plugin std::shared_ptr _plugin; /// Map of vectors of y scales. This is mutable so we can use caching/lazy evaluation. mutable std::map > _yscales; /// Particles used for constituent and tag lookup Particles _fsparticles, _tagparticles; }; } #endif diff --git a/include/Rivet/Projections/PxConePlugin.hh b/include/Rivet/Projections/PxConePlugin.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Projections/PxConePlugin.hh @@ -0,0 +1,165 @@ +//FJSTARTHEADER +// $Id: PxConePlugin.hh,v 1.1 2019/01/04 09:47:31 leif Exp $ +// +// 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 +} + +//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/src/Projections/FastJets.cc b/src/Projections/FastJets.cc --- a/src/Projections/FastJets.cc +++ b/src/Projections/FastJets.cc @@ -1,215 +1,216 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/HeavyHadrons.hh" #include "Rivet/Projections/TauFinder.hh" namespace Rivet { void FastJets::_initBase() { setName("FastJets"); addProjection(HeavyHadrons(), "HFHadrons"); addProjection(TauFinder(TauFinder::HADRONIC), "Taus"); } void FastJets::_initJdef(JetAlgName alg, double rparameter, double seed_threshold) { MSG_DEBUG("JetAlg = " << alg); MSG_DEBUG("R parameter = " << rparameter); MSG_DEBUG("Seed threshold = " << seed_threshold); if (alg == KT) { _jdef = fastjet::JetDefinition(fastjet::kt_algorithm, rparameter, fastjet::E_scheme); } else if (alg == CAM) { _jdef = fastjet::JetDefinition(fastjet::cambridge_algorithm, rparameter, fastjet::E_scheme); } else if (alg == ANTIKT) { _jdef = fastjet::JetDefinition(fastjet::antikt_algorithm, rparameter, fastjet::E_scheme); } else if (alg == DURHAM) { _jdef = fastjet::JetDefinition(fastjet::ee_kt_algorithm, fastjet::E_scheme); } else if (alg == GENKTEE) { _jdef = fastjet::JetDefinition(fastjet::ee_genkt_algorithm, rparameter, -1); } else { // Plugins: if (alg == SISCONE) { const double OVERLAP_THRESHOLD = 0.75; _plugin.reset(new fastjet::SISConePlugin(rparameter, OVERLAP_THRESHOLD)); - // } else if (alg == PXCONE) { - // string msg = "PxCone currently not supported, since FastJet doesn't install it by default. "; - // msg += "Please notify the Rivet authors if this behaviour should be changed."; - // throw Error(msg); - // _plugin.reset(new fastjet::PxConePlugin(rparameter)); + } else if (alg == PXCONE) { + string msg = "Using own c++ version of PxCone, since FastJet doesn't install it by default. "; + msg += "Please notify the Rivet authors if this behaviour should be changed."; + MSG_WARNING(msg); + // _plugin.reset(new fastjet::PxConePlugin(rparameter)); + _plugin.reset(new Rivet::PxConePlugin(rparameter)); } else if (alg == ATLASCONE) { const double OVERLAP_THRESHOLD = 0.5; _plugin.reset(new fastjet::ATLASConePlugin(rparameter, seed_threshold, OVERLAP_THRESHOLD)); } else if (alg == CMSCONE) { _plugin.reset(new fastjet::CMSIterativeConePlugin(rparameter, seed_threshold)); } else if (alg == CDFJETCLU) { const double OVERLAP_THRESHOLD = 0.75; _plugin.reset(new fastjet::CDFJetCluPlugin(rparameter, OVERLAP_THRESHOLD, seed_threshold)); } else if (alg == CDFMIDPOINT) { const double OVERLAP_THRESHOLD = 0.5; _plugin.reset(new fastjet::CDFMidPointPlugin(rparameter, OVERLAP_THRESHOLD, seed_threshold)); } else if (alg == D0ILCONE) { const double min_jet_Et = 6.0; _plugin.reset(new fastjet::D0RunIIConePlugin(rparameter, min_jet_Et)); } else if (alg == JADE) { _plugin.reset(new fastjet::JadePlugin()); } else if (alg == TRACKJET) { _plugin.reset(new fastjet::TrackJetPlugin(rparameter)); } _jdef = fastjet::JetDefinition(_plugin.get()); } } int FastJets::compare(const Projection& p) const { const FastJets& other = dynamic_cast(p); return \ cmp(_useMuons, other._useMuons) || cmp(_useInvisibles, other._useInvisibles) || mkNamedPCmp(other, "FS") || cmp(_jdef.jet_algorithm(), other._jdef.jet_algorithm()) || cmp(_jdef.recombination_scheme(), other._jdef.recombination_scheme()) || cmp(_jdef.plugin(), other._jdef.plugin()) || cmp(_jdef.R(), other._jdef.R()) || cmp(_adef, other._adef); } // STATIC PseudoJets FastJets::mkClusterInputs(const Particles& fsparticles, const Particles& tagparticles) { PseudoJets pjs; /// @todo Use FastJet3's UserInfo system to store Particle pointers directly? // Store 4 vector data about each particle into FastJet's PseudoJets for (size_t i = 0; i < fsparticles.size(); ++i) { fastjet::PseudoJet pj = fsparticles[i]; pj.set_user_index(i+1); pjs.push_back(pj); } // And the same for ghost tagging particles (with negative user indices) for (size_t i = 0; i < tagparticles.size(); ++i) { fastjet::PseudoJet pj = tagparticles[i]; pj *= 1e-20; ///< Ghostify the momentum pj.set_user_index(-i-1); pjs.push_back(pj); } return pjs; } // STATIC Jet FastJets::mkJet(const PseudoJet& pj, const Particles& fsparticles, const Particles& tagparticles) { const PseudoJets pjconstituents = pj.constituents(); Particles constituents, tags; constituents.reserve(pjconstituents.size()); for (const fastjet::PseudoJet& pjc : pjconstituents) { // Pure ghosts don't have corresponding particles if (pjc.has_area() && pjc.is_pure_ghost()) continue; // Default user index = 0 doesn't give valid particle lookup if (pjc.user_index() == 0) continue; // Split by index sign into constituent & tag lookup if (pjc.user_index() > 0) { // Find constituents if index > 0 const size_t i = pjc.user_index() - 1; if (i >= fsparticles.size()) throw RangeError("FS particle lookup failed in jet construction"); constituents.push_back(fsparticles.at(i)); } else if (!tagparticles.empty()) { // Find tags if index < 0 const size_t i = abs(pjc.user_index()) - 1; if (i >= tagparticles.size()) throw RangeError("Tag particle lookup failed in jet construction"); tags.push_back(tagparticles.at(i)); } } return Jet(pj, constituents, tags); } // STATIC Jets FastJets::mkJets(const PseudoJets& pjs, const Particles& fsparticles, const Particles& tagparticles) { Jets rtn; rtn.reserve(pjs.size()); for (const PseudoJet pj : pjs) { rtn.push_back(FastJets::mkJet(pj, fsparticles, tagparticles)); } return rtn; } void FastJets::project(const Event& e) { // Assemble final state particles const string fskey = (_useInvisibles == JetAlg::NO_INVISIBLES) ? "VFS" : "FS"; Particles fsparticles = applyProjection(e, fskey).particles(); // Remove prompt invisibles if needed (already done by VFS if using NO_INVISIBLES) if (_useInvisibles == JetAlg::DECAY_INVISIBLES) { ifilter_discard(fsparticles, [](const Particle& p) { return !(p.isVisible() || p.fromDecay()); }); } // Remove prompt/all muons if needed if (_useMuons == JetAlg::DECAY_MUONS) { ifilter_discard(fsparticles, [](const Particle& p) { return isMuon(p) && !p.fromDecay(); }); } else if (_useMuons == JetAlg::NO_MUONS) { ifilter_discard(fsparticles, isMuon); } // Tagging particles const Particles chadrons = applyProjection(e, "HFHadrons").cHadrons(); const Particles bhadrons = applyProjection(e, "HFHadrons").bHadrons(); const Particles taus = applyProjection(e, "Taus").particles(); calc(fsparticles, chadrons+bhadrons+taus); } void FastJets::calc(const Particles& fsparticles, const Particles& tagparticles) { MSG_DEBUG("Finding jets from " << fsparticles.size() << " input particles + " << tagparticles.size() << " tagging particles"); _fsparticles = fsparticles; _tagparticles = tagparticles; // Make pseudojets, with mapping info to Rivet FS and tag particles PseudoJets pjs = mkClusterInputs(_fsparticles, _tagparticles); // Run either basic or area-calculating cluster sequence as reqd. if (_adef) { _cseq.reset(new fastjet::ClusterSequenceArea(pjs, _jdef, *_adef)); } else { _cseq.reset(new fastjet::ClusterSequence(pjs, _jdef)); } MSG_DEBUG("ClusterSequence constructed; Njets_tot = " << _cseq->inclusive_jets().size() << ", Njets(pT > 10 GeV) = " << _cseq->inclusive_jets(10*GeV).size()); } void FastJets::reset() { _yscales.clear(); _fsparticles.clear(); _tagparticles.clear(); /// @todo _cseq = fastjet::ClusterSequence(); } Jets FastJets::_jets() const { /// @todo Cache? return mkJets(pseudojets(), _fsparticles, _tagparticles); } /// @todo "Automate" trimming as part of project() with pre-registered Filters Jet FastJets::trimJet(const Jet& input, const fastjet::Filter& trimmer) const { if (input.pseudojet().associated_cluster_sequence() != clusterSeq().get()) throw Error("To trim a Rivet::Jet, its associated PseudoJet must have come from this FastJets' ClusterSequence"); PseudoJet pj = trimmer(input); return mkJet(pj, _fsparticles, _tagparticles); } PseudoJets FastJets::pseudoJets(double ptmin) const { return clusterSeq() ? clusterSeq()->inclusive_jets(ptmin) : PseudoJets(); } } diff --git a/src/Projections/Makefile.am b/src/Projections/Makefile.am --- a/src/Projections/Makefile.am +++ b/src/Projections/Makefile.am @@ -1,46 +1,47 @@ 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 \ 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 new file mode 100644 --- /dev/null +++ b/src/Projections/PxConePlugin.cc @@ -0,0 +1,1387 @@ +//FJSTARTHEADER +// $Id: PxConePlugin.cc,v 1.1 2019/01/04 09:47:31 leif Exp $ +// +// 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 "Rivet/Projections/PxConePlugin.hh" + +#include "fastjet/ClusterSequence.hh" +#include + +// 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 { + // 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; + 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 + 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 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, + 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 +#include +#include + +/* 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,,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-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) { + 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; + +} +