Index: /trunk/SISCone/area.cc =================================================================== --- /trunk/SISCone/area.cc (revision 20) +++ /trunk/SISCone/area.cc (revision 20) @@ -0,0 +1,339 @@ +// -*- C++ -*- +/////////////////////////////////////////////////////////////////////////////// +// File: area.h // +// Description: header file for the computation of jet area // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:24 $// +/////////////////////////////////////////////////////////////////////////////// + +#include "area.h" +#include "momentum.h" +#include +#include + +namespace siscone{ +using namespace std; + +/******************************************************* + * Cjet_area implementation * + * real Jet information, including its area(s) * + * * + * This class contains information for one single jet. * + * That is, first, its momentum carrying information * + * about its centre and pT, and second, its particle * + * contents (from CJeT). * + * Compared to the Cjet class, it also includes the * + * passive and active areas of the jet computed using * + * the Carea class. * + *******************************************************/ + +// default ctor +//-------------- +Cjet_area::Cjet_area(){ + active_area = passive_area = 0.0; +} + +// jet-initiated ctor +//------------------- +Cjet_area::Cjet_area(Cjet &j){ + v = j.v; + n = j.n; + contents = j.contents; + + pass = j.pass; + + pt_tilde = j.pt_tilde; + sm_var2 = j.sm_var2; + + active_area = passive_area = 0.0; +} + +// default dtor +//-------------- +Cjet_area::~Cjet_area(){ + +} + + +/****************************************************************** + * Csiscone_area implementation * + * class for the computation of jet areas. * + * * + * This is the class user should use whenever you want to compute * + * the jet area (passive and active). * + * It uses the SISCone algorithm to perform the jet analysis. * + ******************************************************************/ + +// default ctor +//------------- +Carea::Carea(){ + grid_size = 60; // 3600 particles added + grid_eta_max = 6.0; // maybe having twice more points in eta than in phi should be nice + grid_shift = 0.5; // 50% "shacking" + + pt_soft = 1e-100; + pt_shift = 0.05; + pt_soft_min = 1e-90; +} + +// default dtor +//------------- +Carea::~Carea(){ + +} + +/* + * compute the jet areas from a given particle set. + * The parameters of this method are the ones which control the jet clustering alghorithm. + * Note that the pt_min is not allowed here soince the jet-area determination involves soft + * particles/jets and thus is used internally. + * - _particles list of particles + * - _radius cone radius + * - _f shared energy threshold for splitting&merging + * - _n_pass_max maximum number of passes (0=full search, the default) + * - _split_merge_scale the scale choice for the split-merge procedure + * NOTE: SM_pt leads to IR unsafety for some events with momentum conservation. + * SM_Et is IR safe but not boost invariant and not implemented(!) + * SM_mt is IR safe for hadronic events, but not for decays of two + * back-to-back particles of identical mass + * SM_pttilde + * is always IR safe, and also boost invariant (default) + * - _hard_only when this is set on, only hard jets are computed + * and not the purely ghosted jets (default: false) + * return the jets together with their areas + ********************************************************************************************/ +int Carea::compute_areas(std::vector &_particles, double _radius, double _f, + int _n_pass_max, Esplit_merge_scale _split_merge_scale, + bool _hard_only){ + + vector all_particles; + + // put "hardest cut-off" if needed + // this avoids computation of ghosted jets when not required and + // significantly shortens the SM. + if (_hard_only){ + SM_var2_hardest_cut_off = pt_soft_min*pt_soft_min; + } + + // clear potential previous runs + jet_areas.clear(); + + // put initial set of particles in the list + int n_hard = _particles.size(); + all_particles = _particles; + + // build the set of ghost particles and add them to the particle list + int i,j; + double eta,phi,pt; + + for (i=0;i= n_hard + // and deduce the number of ghosts in the jet, hence the area. + double area_factor = (2.0*grid_eta_max/grid_size)*(twopi/grid_size); + + for (i=0;i<(int) jets.size();i++){ + jet_areas.push_back(jets[i]); + j=0; + while ((j &_particles, double _radius, double _f, + int _n_pass_max, Esplit_merge_scale _split_merge_scale){ + + vector all_particles; + + // clear potential previous runs + jet_areas.clear(); + + // put initial set of particles in the list + int n_hard = _particles.size(); + all_particles = _particles; + + // build the set of ghost particles and add them to the particle list + int i,j; + double eta,phi,pt; + + for (i=0;i= n_hard + // and deduce the number of ghosts in the jet, hence the area. + double area_factor = (2.0*grid_eta_max/grid_size)*(twopi/grid_size); + for (i=0;i<(int) jets.size();i++){ + j=0; + while ((j &_particles, double _radius, double _f, + int _n_pass_max, Esplit_merge_scale _split_merge_scale, + bool _hard_only){ + + vector all_particles; + + // put "hardest cut-off" if needed + // this avoids computation of ghosted jets when not required and + // significantly shortens the SM. + if (_hard_only){ + SM_var2_hardest_cut_off = pt_soft_min*pt_soft_min; + } + + // clear potential previous runs + jet_areas.clear(); + + // put initial set of particles in the list + int n_hard = _particles.size(); + all_particles = _particles; + + // build the set of ghost particles and add them to the particle list + int i,j; + double eta,phi,pt; + + for (i=0;i= n_hard + // and deduce the number of ghosts in the jet, hence the area. + double area_factor = (2.0*grid_eta_max/grid_size)*(twopi/grid_size); + + for (i=0;i<(int) jets.size();i++){ + jet_areas.push_back(jets[i]); + j=0; + while ((j &_particles, double _radius, double _f, + int _n_pass_max=0, Esplit_merge_scale _split_merge_scale=SM_pttilde, + bool _hard_only=false); + + /** + * compute the jet active areas from a given particle set. + * The parameters of this method are the ones which control the jet clustering alghorithn. + * Note that the pt_min is not allowed here soince the jet-area determination involves soft + * particles/jets and thus is used internally. + * See compute_areas for paramters definition. + * \return the jets together with their active areas + */ + int compute_active_areas(std::vector &_particles, double _radius, double _f, + int _n_pass_max=0, Esplit_merge_scale _split_merge_scale=SM_pttilde, + bool _hard_only=false); + + /** + * compute the jet passive areas from a given particle set. + * The parameters of this method are the ones which control the jet clustering alghorithn. + * Note that the pt_min is not allowed here soince the jet-area determination involves soft + * particles/jets and thus is used internally. + * See compute_areas for paramters definition. + * \return the jets together with their passive areas + */ + int compute_passive_areas(std::vector &_particles, double _radius, double _f, + int _n_pass_max=0, Esplit_merge_scale _split_merge_scale=SM_pttilde); + + int grid_size; ///< size of the grid we add soft particles on (N_soft=(grid_size^2)) + double grid_eta_max; ///< maximal value of eta we add soft particles on + double grid_shift; ///< fractional (random) displacement of the points om the grid + + double pt_soft; ///< pt of the soft particles added + double pt_shift; ///< amplitude of the pt random shift + double pt_soft_min; ///< pt_min used in SM to compute passive areas + + /// jets with their areas + std::vector jet_areas; +}; + +} +#endif Index: /trunk/SISCone/circulator.h =================================================================== --- /trunk/SISCone/circulator.h (revision 20) +++ /trunk/SISCone/circulator.h (revision 20) @@ -0,0 +1,71 @@ +// -*- C++ -*- +/////////////////////////////////////////////////////////////////////////////// +// File: circulator.h // +// Description: header file for circulator (circulator class) // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:24 $// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef __CIRCULATOR_H__ +#define __CIRCULATOR_H__ + +namespace siscone{ + +/// a circulator that is foreseen to take as template member either a +/// pointer or an iterator; +template class circulator { + +public: + inline circulator(T here, T begin, T end) : m_here(here), m_begin(begin), m_end(end) {} + + inline circulator(const circulator & other) : m_here(other.m_here), m_begin(other.m_begin), m_end(other.m_end) {} + + /// set just the position without resetting the begin and end elements + void set_position(const circulator & other) {m_here = other.m_here;} + void set_position(T pointer) {m_here = pointer;} + + T operator()() {return m_here;} + + inline circulator & operator++() { + ++m_here; + if (m_here == m_end) m_here = m_begin; + return *this; + } + + inline circulator & operator--() { + if (m_here == m_begin) m_here = m_end; + --m_here; + return *this; + } + + /// NB: for efficiency, this checks only the here element + bool operator==(const circulator & other) const {return m_here == other.m_here;} + /// NB: for efficiency, this checks only the here element + bool operator!=(const circulator & other) const {return m_here != other.m_here;} + +private: + T m_here, m_begin, m_end; +}; + +} + +#endif // __CIRCULATOR_H__ Index: /trunk/SISCone/config.h =================================================================== --- /trunk/SISCone/config.h (revision 20) +++ /trunk/SISCone/config.h (revision 20) @@ -0,0 +1,59 @@ +/* siscone/config.h. Generated from config.h.in by configure. */ +/* config.h.in. Generated from configure.ac by autoheader. */ + +/* Define to 1 if you have the header file. */ +#define HAVE_DLFCN_H 1 + +/* Define to 1 if you have the header file. */ +#define HAVE_INTTYPES_H 1 + +/* Define to 1 if you have the `m' library (-lm). */ +#define HAVE_LIBM 1 + +/* Define to 1 if you have the header file. */ +#define HAVE_MEMORY_H 1 + +/* Define to 1 if you have the header file. */ +#define HAVE_STDINT_H 1 + +/* Define to 1 if you have the header file. */ +#define HAVE_STDLIB_H 1 + +/* Define to 1 if you have the header file. */ +#define HAVE_STRINGS_H 1 + +/* Define to 1 if you have the header file. */ +#define HAVE_STRING_H 1 + +/* Define to 1 if you have the header file. */ +#define HAVE_SYS_STAT_H 1 + +/* Define to 1 if you have the header file. */ +#define HAVE_SYS_TYPES_H 1 + +/* Define to 1 if you have the header file. */ +#define HAVE_UNISTD_H 1 + +/* Name of package */ +#define PACKAGE "siscone" + +/* Define to the address where bug reports for this package should be sent. */ +#define PACKAGE_BUGREPORT "" + +/* Define to the full name of this package. */ +#define PACKAGE_NAME "SISCone" + +/* Define to the full name and version of this package. */ +#define PACKAGE_STRING "SISCone 1.3.3" + +/* Define to the one symbol short name of this package. */ +#define PACKAGE_TARNAME "siscone" + +/* Define to the version of this package. */ +#define PACKAGE_VERSION "1.3.3" + +/* Define to 1 if you have the ANSI C header files. */ +#define STDC_HEADERS 1 + +/* Version number of package */ +#define VERSION "1.3.3" Index: /trunk/SISCone/defines.h =================================================================== --- /trunk/SISCone/defines.h (revision 20) +++ /trunk/SISCone/defines.h (revision 20) @@ -0,0 +1,128 @@ +// -*- C++ -*- +/////////////////////////////////////////////////////////////////////////////// +// File: defines.h // +// Description: header file for generic parameters definitions // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:25 $// +/////////////////////////////////////////////////////////////////////////////// + +//! \file defines.h +#ifndef __DEFINES_H__ +#define __DEFINES_H__ + +/// program name +// we get "SISCone" by calling +// siscone::siscone_package_name +// defined in siscone.h +// Otherwise, config.h +// It is also defined as "PACKAGE_NAME" in config.h but this method +// might lead to conflicts +//#define PROGRAM PACKAGE_NAME + +// program version +// we get it from +// siscone::siscone_version +// defined in siscone.h +// It is also defined as "VERSION" in config.h but this method +// might lead to conflicts + +/// perform final stability check using the quadtree +/// With the following define enabled, the final check for stability +/// is performed using the quadtree rather than the explicit list of +/// particles (see Cstable_cone::proceed_with_stability()) +//#define USE_QUADTREE_FOR_STABILITY_TEST + + +/// threshold for recomoutation of the cone (see Cstable_cones::update_cone()) +/// When traversing cone candidates along the angular ordering, +/// the momentum of the protojet candidate is computed incrementally +/// from the particles that enter and leave the cone. +/// When the cumulative change in "|px|+|py|" exceeds the cone "|px|+|py|" +/// we explicitely recompute the cone contents +#define PT_TSHOLD 1000.0 + + +/// The following parameter controls collinear safety. For the set of +/// particles used in the search of stable cones, we gather particles +/// if their distance in eta and phi is smaller than EPSILON_COLLINEAR. +/// +/// NB: for things to behave sensibly one requires +/// 1e-15 << EPSILON_COCIRCULAR << EPSILON_COLLINEAR << 1 +/// +/// among the scales that appear in practice (e.g. in deciding to use +/// special strategies), we have EPSILON_COLLINEAR, EPSILON_COCIRCULAR, +/// sqrt(EPSILON_COCIRCULAR) and EPSILON_COLLINEAR / EPSILON_COCIRCULAR. +/// +#define EPSILON_COLLINEAR 1e-8 + + +/// The following parameter controls cocircular situations. +/// When consecutive particles in the ordered vicinity list are separated +/// (in angle) by less that that limit, we consider that we face a situation +/// of cocircularity. +#define EPSILON_COCIRCULAR 1e-12 + + +/// The following define enables you to allow for identical protocones +/// to be merged automatically after each split-merge step before +/// anything else happens. Whether this happens depends on teh value +/// of the merge_identical_protocones flag in Csplit_merge. +/// +/// It we allow such a merging and define allow +/// MERGE_IDENTICAL_PROTOCONES_DEFAULT_TRUE then the +/// 'merge_identical_protocones' flag in Csplit_merge to be set to +/// 'true'. It may be manually reset to false in which case the +/// merging of identical protocones (protojets) will be turned off. +/// +/// Note that this merging identical protocones makes the algorithm +/// infrared-unsafe, so it should be disabled except for testing +/// purposes. +//#define ALLOW_MERGE_IDENTICAL_PROTOCONES +//#define MERGE_IDENTICAL_PROTOCONES_DEFAULT_TRUE + + +/// if EPSILON_SPLITMERGE is defined then, during the split-merge +/// step, when two jets are found with PTs that are identical to +/// within a relative difference of EPSILON_SPLITMERGE they are +/// compared element-by-element to see where the differences are, and +/// one then uses pt1^2-pt2^2 = (pt1-pt2).(pt1+pt2) as an estimator of +/// which is harder. NB: in unfortunate cases, this can take the split +/// merge step up to N n * ln N time, though on normal events there +/// don't seem to have been any major problems yet. +#define EPSILON_SPLITMERGE 1e-12 + +/// definition of 2*M_PI which is useful a bit everyhere! +const double twopi = 6.283185307179586476925286766559005768394; + +/// debugging information +//#define DEBUG_STABLE_CONES ///< debug messages in stable cones search (not implemented yet) +//#define DEBUG_SPLIT_MERGE ///< debug messages in split-merge +//#define DEBUG ///< all debug messages ! + +// in case all debug massages allowed, allow them in practice ! +#ifdef DEBUG +#define DEBUG_STABLE_CONES +#define DEBUG_SPLIT_MERGE +#endif + +#endif // __DEFINES_H__ + Index: /trunk/SISCone/geom_2d.cc =================================================================== --- /trunk/SISCone/geom_2d.cc (revision 20) +++ /trunk/SISCone/geom_2d.cc (revision 20) @@ -0,0 +1,150 @@ +/////////////////////////////////////////////////////////////////////////////// +// File: geom_2d.cpp // +// Description: source file for two-dimensional geometry tools // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:25 $// +/////////////////////////////////////////////////////////////////////////////// + +#include "geom_2d.h" +#include + +namespace siscone{ + +#define PHI_RANGE_MASK 0xFFFFFFFF + +/********************************************************* + * class Ceta_phi_range implementation * + * class for holding a covering range in eta-phi * + * * + * This class deals with ranges in the eta-phi plane. It * + * implements methods to test if two ranges overlap and * + * to take the union of two overlapping intervals. * + *********************************************************/ + +using namespace std; + +// static member default init +//---------------------------- +double Ceta_phi_range::eta_min = -100.0; +double Ceta_phi_range::eta_max = 100.0; + +// default ctor +//-------------- +Ceta_phi_range::Ceta_phi_range(){ + eta_range = 0; + phi_range = 0; +} + +// ctor with initialisation +// we initialise with a centre (in eta,phi) and a radius +// - c_eta eta coordinate of the centre +// - c_phi phi coordinate of the centre +// - R radius +//------------------------------------------------------- +Ceta_phi_range::Ceta_phi_range(double c_eta, double c_phi, double R){ + // determination of the eta range + //------------------------------- + double xmin = max(c_eta-R,eta_min+0.0001); + double xmax = min(c_eta+R,eta_max-0.0001); + + unsigned int cell_min = get_eta_cell(xmin); + unsigned int cell_max = get_eta_cell(xmax); + + // warning: if cell_max==2^31, 2*cell_max==0 hence, + // even if the next formula is formally (2*cell_max-cell_min), + // expressing it as (cell_max-cell_min)+cell_max is safe. + eta_range = (cell_max-cell_min)+cell_max; + + // determination of the phi range + // !! taking care of periodicity !! + //--------------------------------- + xmin = phi_in_range(c_phi-R); + xmax = phi_in_range(c_phi+R); + + cell_min = get_phi_cell(xmin); + cell_max = get_phi_cell(xmax); + + // Also, if the interval goes through pi, inversion is needed + if (xmax>xmin) + phi_range = (cell_max-cell_min)+cell_max; + else { + phi_range = (cell_min==cell_max) + ? PHI_RANGE_MASK + : ((PHI_RANGE_MASK^(cell_min-cell_max)) + cell_max); + } +} + +// assignment of range +// - r range to assign to current one +//--------------------------------------- +Ceta_phi_range& Ceta_phi_range::operator = (const Ceta_phi_range &r){ + eta_range = r.eta_range; + phi_range = r.phi_range; + + return *this; +} + +// add a particle to the range +// - eta eta coordinate of the particle +// - phi phi coordinate of the particle +// \return 0 on success, 1 on error +//---------------------------------------- +int Ceta_phi_range::add_particle(const double eta, const double phi){ + // deal with the eta coordinate + eta_range |= get_eta_cell(eta); + + // deal with the phi coordinate + phi_range |= get_phi_cell(phi); + + return 0; +} + + +// test overlap +// - r1 first range +// - r2 second range +// return true if overlap, false otherwise. +//------------------------------------------ +bool is_range_overlap(const Ceta_phi_range &r1, const Ceta_phi_range &r2){ + // check overlap in eta AND phi + return ((r1.eta_range & r2.eta_range) && (r1.phi_range & r2.phi_range)); +} + +// compute union +// Note: we assume that the two intervals overlap +// - r1 first range +// - r2 second range +// \return union of the two ranges +//------------------------------------------ +const Ceta_phi_range range_union (const Ceta_phi_range &r1, const Ceta_phi_range &r2){ + Ceta_phi_range tmp; + + // compute union in eta + tmp.eta_range = r1.eta_range | r2.eta_range; + + // compute union in phi + tmp.phi_range = r1.phi_range | r2.phi_range; + + return tmp; +} + +} Index: /trunk/SISCone/geom_2d.h =================================================================== --- /trunk/SISCone/geom_2d.h (revision 20) +++ /trunk/SISCone/geom_2d.h (revision 20) @@ -0,0 +1,179 @@ +// -*- C++ -*- +/////////////////////////////////////////////////////////////////////////////// +// File: geom_2d.h // +// Description: header file for two-dimensional geometry tools // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:25 $// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef __GEOM_2D_H__ +#define __GEOM_2D_H__ + +#include +#include +#include "defines.h" + +#ifndef M_PI +#define M_PI 3.141592653589793238462643383279502884197 +#endif + +namespace siscone{ + +/// return a result that corresponds to phi, but in the +/// range (-pi..pi]; the result is only correct if -3pi < phi <= 3pi +inline double phi_in_range(double phi) { + if (phi <= -M_PI) phi += twopi; + else if (phi > M_PI) phi -= twopi; + return phi; +} + +/// return the difference between the two phi values, +/// placed in the correct range (-pi..pi], , assuming that phi1,phi2 +/// are already in the correct range. +inline double dphi(double phi1, double phi2) { + return phi_in_range(phi1-phi2); +} + + +/// return the absolute difference between the two phi values, +/// placed in the correct range, assuming that phi1,phi2 are already +/// in the correct range. +inline double abs_dphi(double phi1, double phi2) { + double delta = fabs(phi1-phi2); + return delta > M_PI ? twopi-delta : delta; +} + +/// return the square of the argument +inline double pow2(double x) {return x*x;} + + +/** + * \class Ctwovect + * \brief class for holding a two-vector + */ +class Ctwovect { +public: + /// default ctor + Ctwovect() : x(0.0), y(0.0) {} + + /// ctor with initialisation + /// \param _x first coordinate + /// \param _y second coordinate + Ctwovect(double _x, double _y) : x(_x), y(_y) {} + + /// vector coordinates + double x, y; + + /// norm (modulud square) of the vector + inline double mod2() const {return pow2(x)+pow2(y);} + + /// modulus of the vector + inline double modulus() const {return sqrt(mod2());} +}; + + +/// dot product of two 2-vectors +/// \param a first 2-vect +/// \param b second 2-vect +/// \return a.b is returned +inline double dot_product(const Ctwovect & a, const Ctwovect & b) { + return a.x*b.x + a.y*b.y; +} + + +/// cross product of two 2-vectors +/// \param a first 2-vect +/// \param b second 2-vect +/// \return a x b is returned +inline double cross_product(const Ctwovect & a, const Ctwovect & b) { + return a.x*b.y - a.y*b.x; +} + + +/** + * \class Ceta_phi_range + * \brief class for holding a covering range in eta-phi + * + * This class deals with ranges in the eta-phi plane. It + * implements methods to test if two ranges overlap and + * to take the union of two overlapping intervals. + */ +class Ceta_phi_range{ +public: + /// default ctor + Ceta_phi_range(); + + /// ctor with initialisation + /// we initialise with a centre (in eta,phi) and a radius + /// \param c_eta eta coordinate of the centre + /// \param c_phi phi coordinate of the centre + /// \param R radius + Ceta_phi_range(double c_eta, double c_phi, double R); + + /// assignment of range + /// \param r range to assign to current one + Ceta_phi_range& operator = (const Ceta_phi_range &r); + + /// add a particle to the range + /// \param eta eta coordinate of the particle + /// \param phi phi coordinate of the particle + /// \return 0 on success, 1 on error + int add_particle(const double eta, const double phi); + + /// eta range as a binary coding of covered cells + unsigned int eta_range; + + /// phi range as a binary coding of covered cells + unsigned int phi_range; + + /// maximal value for eta + static double eta_min; + static double eta_max; + +private: + /// return the cell index corrsponding to an eta value + inline unsigned int get_eta_cell(double eta){ + return (unsigned int) (1 << ((int) (32*((eta-eta_min)/(eta_max-eta_min))))); + } + + /// return the cell index corrsponding to a phi value + inline unsigned int get_phi_cell(double phi){ + return (unsigned int) (1 << ((int) (32*phi/twopi+16)%32)); + } +}; + +/// test overlap +/// \param r1 first range +/// \param r2 second range +/// \return true if overlap, false otherwise. +bool is_range_overlap(const Ceta_phi_range &r1, const Ceta_phi_range &r2); + +/// compute union +/// Note: we assume that the two intervals overlap +/// \param r1 first range +/// \param r2 second range +/// \return union of the two ranges +const Ceta_phi_range range_union(const Ceta_phi_range &r1, const Ceta_phi_range &r2); + +} + +#endif Index: /trunk/SISCone/hash.cc =================================================================== --- /trunk/SISCone/hash.cc (revision 20) +++ /trunk/SISCone/hash.cc (revision 20) @@ -0,0 +1,217 @@ +/////////////////////////////////////////////////////////////////////////////// +// File: hash.cpp // +// Description: source file for classes hash_element and hash_cones // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:26 $// +/////////////////////////////////////////////////////////////////////////////// + +#include +#include +#include "hash.h" +#include + +namespace siscone{ + +using namespace std; + +/************************************************************** + * implementation of hash_cones * + * list of cones candidates. * + * We store in this class all the hash_elements and give * + * functions to manipulate them. * + **************************************************************/ + +// constructor with initialisation +// - _Np number of particles +// - _R2 cone radius (squared) +//----------------------------------- +hash_cones::hash_cones(int _Np, double _R2){ + int i; + + n_cones = 0; + + // determine hash size + mask = 1 << (int) (2*log(double(_Np))/log(2.0)); + if (mask<=1) mask=2; + + // create hash + hash_array = new hash_element*[mask]; + mask--; + + // set the array to 0 + //? needed ? + for (i=0;inext; + delete elm; + } + } + + delete[] hash_array; +} + + +/* + * insert a new candidate into the hash. + * - v 4-momentum of the cone to add + * - parent parent particle defining the cone + * - child child particle defining the cone + * - p_io whether the parent has to belong to the cone or not + * - c_io whether the child has to belong to the cone or not + * return 0 on success, 1 on error + ***********************************************************************/ +int hash_cones::insert(Cmomentum *v, Cmomentum *parent, Cmomentum *child, bool p_io, bool c_io){ + hash_element *elm; + int index = (v->ref.ref[0]) & mask; + + // check the array cell corresponding to our reference + elm = hash_array[index]; + do{ + // if it is not present, add it + if (elm==NULL){ + // create element + elm = new hash_element; + + // set its varibles + // Note: at this level, eta and phi have already been computed + // through Cmomentum::build_etaphi. + elm->ref = v->ref; + + //compute vectors centre + v->build_etaphi(); + elm->eta = v->eta; + elm->phi = v->phi; + // if at least one of the two is_inside tests gives a result != from the expected, + // the || will be true hence !(...) false as wanted + elm->is_stable = !((is_inside(v, parent)^p_io)||(is_inside(v, child)^c_io)); + //cout << "-- new status of " << v->ref[0] << ":" << elm->is_stable << endl; + + // update hash + elm->next = hash_array[index]; + hash_array[index] = elm; + + n_cones++; + return 0; + } + + // if the cone is already there, simply update stability status + if (v->ref == elm->ref){ + // there is only an update to perform to see if the cone is still stable + if (elm->is_stable){ + v->build_etaphi(); + elm->is_stable = !((is_inside(v, parent)^p_io)||(is_inside(v, child)^c_io)); + //cout << " parent/child: " + // << parent->ref[0] << ":" << is_inside(v, parent) << ":" << p_io << " " + // << child->ref[0] << ":" << is_inside(v, child) << ":" << c_io << endl; + //cout << "-- rep status of " << v->ref[0] << ":" << elm->is_stable << endl; + //cout << v->eta << " " << v->phi << endl; + //cout << (child->eta) << " " << child->phi << endl; + } + return 0; + } + + elm = elm->next; + } while (1); + + return 1; +} + +/* + * insert a new candidate into the hash. + * - v 4-momentum of te cone to add + * Note, in this case, we assume stability. We also assume + * that eta and phi are computed for v + * return 0 on success, 1 on error + ***********************************************************************/ +int hash_cones::insert(Cmomentum *v){ + hash_element *elm; + int index = (v->ref.ref[0]) & mask; + //cout << "-- stable candidate: " << v->ref[0] << ":" << endl; + + // check the array cell corresponding to our reference + elm = hash_array[index]; + do{ + // if it is not present, add it + if (elm==NULL){ + // create element + elm = new hash_element; + + // set its varibles + // Note: at this level, eta and phi have already been computed + // through Cmomentum::build_etaphi. + elm->ref = v->ref; + elm->eta = v->eta; + elm->phi = v->phi; + elm->is_stable = true; + + // update hash + elm->next = hash_array[index]; + hash_array[index] = elm; + + n_cones++; + return 0; + } + + // if the cone is already there, we have nothing to do + if (v->ref == elm->ref){ + return 0; + } + + elm = elm->next; + } while (1); + + return 1; +} + +/* + * test if a particle is inside a cone of given centre. + * check if the particle of coordinates 'v' is inside the circle of radius R + * centered at 'centre'. + * - centre centre of the circle + * - v particle to test + * return true if inside, false if outside + ******************************************************************************/ +inline bool hash_cones::is_inside(Cmomentum *centre, Cmomentum *v){ + double dx, dy; + + dx = centre->eta - v->eta; + dy = fabs(centre->phi - v->phi); + if (dy>M_PI) + dy -= 2.0*M_PI; + + return dx*dx+dy*dy +#include + +namespace siscone{ + +/************************************************************************* + * class Cmomentum * + * This class contains the information for particle or group of * + * particles management. * + * It includes all Lorentz properties as well as tools for summing them. * + *************************************************************************/ + +// default ctor +//-------------- +Cmomentum::Cmomentum(){ + eta = 0.0; + phi = 0.0; + px = py = pz = E = 0.0; + ref = Creference(); + index = -1; +} + +// ctor with initialisation +//-------------------------- +Cmomentum::Cmomentum(double _px, double _py, double _pz, double _E){ + px = _px; + py = _py; + pz = _pz; + E = _E; + + // compute eta and phi + build_etaphi(); + ref = Creference(); +} + +// ctor with detailed initialisation +//----------------------------------- +Cmomentum::Cmomentum(double _eta, double _phi, Creference _ref){ + eta = _eta; + phi = _phi; + + ref = _ref; +} + +// default dtor +//-------------- +Cmomentum::~Cmomentum(){ + +} + +// assignment of vectors +//----------------------- +Cmomentum& Cmomentum::operator = (const Cmomentum &v){ + px = v.px; + py = v.py; + pz = v.pz; + E = v.E; + + eta = v.eta; + phi = v.phi; + + ref = v.ref; + return *this; +} + +// addition of vectors +// !!! WARNING !!! no updating of eta and phi !!! +//------------------------------------------------ +const Cmomentum Cmomentum::operator + (const Cmomentum &v){ + Cmomentum tmp = *this; + return tmp+=v; +} + +// incrementation of vectors +// !!! WARNING !!! no updating of eta and phi !!! +//------------------------------------------------ +Cmomentum& Cmomentum::operator += (const Cmomentum &v){ + px+=v.px; + py+=v.py; + pz+=v.pz; + E +=v.E; + + ref+=v.ref; + + return *this; +} + +// incrementation of vectors +// !!! WARNING !!! no updating of eta and phi !!! +//------------------------------------------------ +Cmomentum& Cmomentum::operator -= (const Cmomentum &v){ + px-=v.px; + py-=v.py; + pz-=v.pz; + E -=v.E; + + ref-=v.ref; + return *this; +} + +// build eta-phi from 4-momentum info +// !!! WARNING !!! +// !!! computing eta and phi is time-consuming !!! +// !!! use this whenever you need eta or phi !!! +// !!! automatically called for single-particle !!! +//-------------------------------------------------- +void Cmomentum::build_etaphi(){ + // note: the factor n (ref.nb) cancels in all expressions !! + eta = 0.5*log((E+pz)/(E-pz)); + phi = atan2(py,px); +} + + +// ordering of two vectors +// the default ordering is w.r.t. their references +//------------------------------------------------- +bool operator < (const Cmomentum &v1, const Cmomentum &v2){ + return v1.ref < v2.ref; +} + +// ordering of vectors in eta (e.g. used in collinear tests) +//----------------------------------------------------------- +bool momentum_eta_less(const Cmomentum &v1, const Cmomentum &v2){ + return v1.eta < v2.eta; +} + +// ordering of vectors in pt +//--------------------------- +bool momentum_pt_less(const Cmomentum &v1, const Cmomentum &v2){ + return v1.perp2() < v2.perp2(); +} + +} + Index: /trunk/SISCone/momentum.h =================================================================== --- /trunk/SISCone/momentum.h (revision 20) +++ /trunk/SISCone/momentum.h (revision 20) @@ -0,0 +1,157 @@ +// -*- C++ -*- +/////////////////////////////////////////////////////////////////////////////// +// File: momentum.h // +// Description: header file for 4-momentum class Cmomentum // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:26 $// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef __VECTOR_H__ +#define __VECTOR_H__ + +#include +#include +#include "reference.h" +#include "geom_2d.h" +#include "defines.h" + +namespace siscone{ + +/** + * \class Cmomentum + * \brief base class for dynamic coordinates management + * + * This class contains the information for particle or group of + * particles management. + * It includes all Lorentz properties as well as tools for summing them. + * Note: 'sums' over phi angles are indeed averages. This allows to + * deal with periodicity at each step + */ +class Cmomentum{ + public: + /// default ctor + Cmomentum(); + + /// ctor with initialisation + Cmomentum(double _px, double _py, double _pz, double _E); + + /// ctor with detailed initialisation + Cmomentum(double _eta, double _phi, Creference _ref); + + /// default dtor + ~Cmomentum(); + + /// computes pT + inline double perp() const {return sqrt(perp2());} + + /// computes pT^2 + inline double perp2() const {return px*px+py*py;} + + /// computes m + inline double mass() const {return sqrt(mass2());} + + /// computes m^2 + inline double mass2() const {return perpmass2()-perp2();} + + /// transverse mass, mt = sqrt(pt^2+m^2) = sqrt(E^2 - pz^2) + inline double perpmass() const {return sqrt((E-pz)*(E+pz));} + + /// transverse mass squared, mt^2 = pt^2+m^2 = E^2 - pz^2 + inline double perpmass2() const {return (E-pz)*(E+pz);} + + /// computes transverse energy + inline double Et() const {return E/sqrt(1.0+pz*pz/perp2());} + + /// computes transverse energy (squared) + inline double Et2() const {return E*E/(1.0+pz*pz/perp2());} + + /// assignment of vectors + Cmomentum& operator = (const Cmomentum &v); + + /// addition of vectors + /// !!! WARNING !!! no updating of eta and phi !!! + const Cmomentum operator + (const Cmomentum &v); + + /// incrementation of vectors + /// !!! WARNING !!! no updating of eta and phi !!! + Cmomentum& operator += (const Cmomentum &v); + + /// decrementation of vectors + /// !!! WARNING !!! no updating of eta and phi !!! + Cmomentum& operator -= (const Cmomentum &v); + + /// build eta-phi from 4-momentum info + /// !!! WARNING !!! + /// !!! computing eta and phi is time-consuming !!! + /// !!! use this whenever you need eta or phi !!! + /// !!! automatically called for single-particle !!! + void build_etaphi(); + + double px; ///< x-momentum + double py; ///< y-momentum + double pz; ///< z-momentum + double E; ///< energy + + double eta; ///< particle pseudo-rapidity + double phi; ///< particle azimuthal angle + int parent_index; ///< particle number in the parent list + int index; ///< internal particle number + + ////////////////////////////////////////////// + // the following part is used for checksums // + ////////////////////////////////////////////// + Creference ref; ///< reference number for the vector +}; + +/// ordering of two vectors +/// this is by default done w.r.t. their references +bool operator < (const Cmomentum &v1, const Cmomentum &v2); + +/// ordering of vectors in eta (e.g. used in collinear tests) +bool momentum_eta_less(const Cmomentum &v1, const Cmomentum &v2); + +/// ordering of vectors in pt +bool momentum_pt_less(const Cmomentum &v1, const Cmomentum &v2); + + +////////////////////////// +// some handy utilities // +////////////////////////// + +/// get distance between to eta-phi points +/// \param eta eta coordinate of first point +/// \param phi phi coordinate of first point +/// \param v vector defining the second point +inline double get_distance(double eta, double phi, Cmomentum *v){ + double dx, dy; + + dx = eta - v->eta; + dy = fabs(phi - v->phi); + if (dy>M_PI) + dy -= twopi; + + return dx*dx+dy*dy; +} + +} + +#endif Index: /trunk/SISCone/protocones.cc =================================================================== --- /trunk/SISCone/protocones.cc (revision 20) +++ /trunk/SISCone/protocones.cc (revision 20) @@ -0,0 +1,838 @@ +/////////////////////////////////////////////////////////////////////////////// +// File: protocones.cpp // +// Description: source file for stable cones determination (Cstable_cones) // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:27 $// +/////////////////////////////////////////////////////////////////////////////// + +/******************************************************* + * Introductory note: * + * Since this file has many member functions, we have * + * structured them in categories: * + * INITIALISATION FUNCTIONS * + * - ctor() * + * - ctor(particle_list) * + * - dtor() * + * - init(particle_list) * + * ALGORITHM MAIN ENTRY * + * - get_stable_cone(radius) * + * ALGORITHM MAIN STEPS * + * - init_cone() * + * - test_cone() * + * - update_cone() * + * - proceed_with_stability() * + * ALGORITHM MAIN STEPS FOR COCIRCULAR SITUATIONS * + * - cocircular_pt_less(v1, v2) * + * - prepare_cocircular_list() * + * - test_cone_cocircular() * + * - test_stability(candidate, border_list) * + * - updat_cone_cocircular() * + * RECOMPUTATION OF CONE CONTENTS * + * - compute_cone_contents() * + * - recompute_cone_contents() * + * - recompute_cone_contents_if_needed() * + * VARIOUS TOOLS * + * - circle_intersect() * + * - is_inside() * + * - abs_dangle() * + *******************************************************/ + +#include "protocones.h" +#include "siscone_error.h" +#include "defines.h" +#include +#include +#include "circulator.h" +#include + +namespace siscone{ + +using namespace std; + +/********************************************************************** + * Cstable_cones implementation * + * Computes the list of stable comes from a particle list. * + * This class does the first fundamental task of te cone algorithm: * + * it is used to compute the list of stable cones given a list * + * of particles. * + **********************************************************************/ + +//////////////////////////////////////////////////////// +// INITIALISATION FUNCTIONS // +// - ctor() // +// - ctor(particle_list) // +// - dtor() // +// - init(particle_list) // +//////////////////////////////////////////////////////// + +// default ctor +//-------------- +Cstable_cones::Cstable_cones(){ + nb_tot = 0; + hc = NULL; +} + +// ctor with initialisation +//-------------------------- +Cstable_cones::Cstable_cones(vector &_particle_list) + : Cvicinity(_particle_list){ + + nb_tot = 0; + hc = NULL; +} + +// default dtor +//-------------- +Cstable_cones::~Cstable_cones(){ + if (hc!=NULL) delete hc; +} + +/* + * initialisation + * - _particle_list list of particles + * - _n number of particles + *********************************************************************/ +void Cstable_cones::init(vector &_particle_list){ + // check already allocated mem + if (hc!=NULL){ + delete hc; + } + if (protocones.size()!=0) + protocones.clear(); + + multiple_centre_done.clear(); + + // initialisation + set_particle_list(_particle_list); +} + + +//////////////////////////////////////////////////////// +// ALGORITHM MAIN ENTRY // +// - get_stable_cone(radius) // +//////////////////////////////////////////////////////// + +/* + * compute stable cones. + * This function really does the job i.e. computes + * the list of stable cones (in a seedless way) + * - _radius: radius of the cones + * The number of stable cones found is returned + *********************************************************************/ +int Cstable_cones::get_stable_cones(double _radius){ + int p_idx; + + // check if everything is correctly initialised + if (n_part==0){ + return 0; + } + + R = _radius; + R2 = R*R; + + // allow hash for cones candidates + hc = new hash_cones(n_part, R2); + + // browse all particles + for (p_idx=0;p_idxv; + centre_idx = first_cone; + + // build the initial cone (nodist: avoids calculating distances -- + // just deduces contents by circulating around all in/out operations) + // this function also sets the list of included particles + compute_cone_contents(); + + return 0; +} + + +/* + * test cones. + * We check if the cone(s) built with the present parent and child + * are stable + * return 0 on success 1 on error + *********************************************************************/ +int Cstable_cones::test_cone(){ + Creference weighted_cone_ref; + + // depending on the side we are taking the child particle, + // we test different configuration. + // Each time, two configurations are tested in such a way that + // all 4 possible cases (parent or child in or out the cone) + // are tested when taking the pair of particle parent+child + // and child+parent. + + // here are the tests entering the first series: + // 1. check if the cone is already inserted + // 2. check cone stability for the parent and child particles + + if (centre->side){ + // test when both particles are not in the cone + // or when both are in. + // Note: for the totally exclusive case, test emptyness before + cone_candidate = cone; + if (cone.ref.not_empty()){ + hc->insert(&cone_candidate, parent, child, false, false); + } + + cone_candidate = cone; + cone_candidate+= *parent + *child; + hc->insert(&cone_candidate, parent, child, true, true); + } else { + // test when 1! of the particles is in the cone + cone_candidate = cone + *parent; + hc->insert(&cone_candidate, parent, child, true, false); + + cone_candidate = cone + *child; + hc->insert(&cone_candidate, parent, child, false, true); + } + + nb_tot+=2; + + return 0; +} + + +/* + * update the cone + * go to the next child for that parent and update 'cone' appropriately + * return 0 if update candidate found, 1 otherwise + ***********************************************************************/ +int Cstable_cones::update_cone(){ + // get the next child and centre + centre_idx++; + if (centre_idx==vicinity_size) + centre_idx=0; + if (centre_idx==first_cone) + return 1; + + // update the cone w.r.t. the old child + // only required if the old child is entering inside in which + // case we need to add it. We also know that the child is + // inside iff its side is -. + if (!centre->side){ + // update cone + cone += (*child); + + // update info on particles inside + centre->is_inside->cone = true; + + // update stability check quantities + dpt += fabs(child->px)+fabs(child->py); + } + + // update centre and child to correspond to the new position + centre = vicinity[centre_idx]; + child = centre->v; + + // check cocircularity + // note that if cocirculaity is detected (i.e. if we receive 1 + // in the next test), we need to recall 'update_cone' directly + // since tests and remaining part of te update has been performed + //if (cocircular_check()) + if (cocircular_check()) + return update_cone(); + + + // update the cone w.r.t. the new child + // only required if the new child was already inside in which + // case we need to remove it. We also know that the child is + // inside iff its side is +. + if ((centre->side) && (cone.ref.not_empty())){ + // update cone + cone -= (*child); + + // update info on particles inside + centre->is_inside->cone = false; + + // update stability check quantities + dpt += fabs(child->px)+fabs(child->py); //child->perp2(); + } + + // check that the addition and subtraction of vectors does + // not lead to too much rounding error + // for that, we compute the sum of pt modifications and of |pt| + // since last recomputation and once the ratio overpasses a threshold + // we recompute vicinity. + if ((dpt>PT_TSHOLD*(fabs(cone.px)+fabs(cone.py))) && (cone.ref.not_empty())){ + recompute_cone_contents(); + } + if (cone.ref.is_empty()){ + cone = Cmomentum(); + dpt=0.0; + } + + return 0; +} + + +/* + * compute stability of all enumerated candidates. + * For all candidate cones which are stable w.r.t. their border particles, + * pass the last test: stability with quadtree intersection + ************************************************************************/ +int Cstable_cones::proceed_with_stability(){ + int i,n; + hash_element *elm; + + n=0; + for (i=0;i<=hc->mask;i++){ + // test ith cell of the hash array + elm = hc->hash_array[i]; + + // browse elements therein + while (elm!=NULL){ + // test stability + if (elm->is_stable){ + // stability is not ensured by all pairs of "edges" already browsed +#ifdef USE_QUADTREE_FOR_STABILITY_TEST + // => testing stability with quadtree intersection + if (quadtree->circle_intersect(elm->eta, elm->phi, R2)==elm->ref) +#else + // => testing stability with the particle-list intersection + if (circle_intersect(elm->eta, elm->phi)==elm->ref) +#endif + protocones.push_back(Cmomentum(elm->eta, elm->phi, elm->ref)); + } + + // jump to the next one + elm = elm->next; + } + } + + // free hash + // we do that at this level because hash eats rather a lot of memory + // we want to free it before running the split/merge algorithm + delete hc; + hc=NULL; + + return protocones.size(); +} + + +//////////////////////////////////////////////////////// +// ALGORITHM MAIN STEPS FOR COCIRCULAR SITUATIONS // +// - cocircular_pt_less(v1, v2) // +// - prepare_cocircular_list() // +// - test_cone_cocircular() // +// - test_stability(candidate, border_vect) // +// - updat_cone_cocircular() // +//////////////////////////////////////////////////////// + +/// pt-ordering of momenta used for the cocircular case +bool cocircular_pt_less(Cmomentum *v1, Cmomentum *v2){ + return v1->perp2() < v2->perp2(); +} + +/* + * run through the vicinity of the current parent and for each child + * establish which other members are cocircular... Note that the list + * associated with each child contains references to vicinity + * elements: thus two vicinity elements each associated with one given + * particle may appear in a list -- this needs to be watched out for + * later on... + **********************************************************************/ +void Cstable_cones::prepare_cocircular_lists() { + circulator::iterator > here(vicinity.begin(), + vicinity.begin(), + vicinity.end()); + + circulator::iterator > search(here); + + do { + Cvicinity_elm* here_pntr = *here(); + search.set_position(here); + + // search forwards for things that should have "here" included in + // their cocircularity list + while (true) { + ++search; + if ( abs_dphi((*search())->angle, here_pntr->angle) < + here_pntr->cocircular_range + && search() != here()) { + (*search())->cocircular.push_back(here_pntr); + } else { + break; + } + } + + // search backwards + search.set_position(here); + while (true) { + --search; + if ( abs_dphi((*search())->angle, here_pntr->angle) < + here_pntr->cocircular_range + && search() != here()) { + (*search())->cocircular.push_back(here_pntr); + } else { + break; + } + } + + ++here; + } while (here() != vicinity.begin()); + +} + +/* + * Testing cocircular configurations in p^3 time, + * rather than 2^p time; we will test all contiguous subsets of points + * on the border --- note that this is till probably overkill, since + * in principle we only have to test situations where up to a + * half-circle is filled (but going to a full circle is simpler) + ******************************************************************/ +void Cstable_cones::test_cone_cocircular(Cmomentum & borderless_cone, + list & border_list) { + vector border_vect; + + border_vect.reserve(border_list.size()); + for (list::iterator it = border_list.begin(); + it != border_list.end(); it++) { + border_vect.push_back(Cborder_store(*it, centre->eta, centre->phi)); + } + + // get them into order of angle + sort(border_vect.begin(), border_vect.end()); + + // set up some circulators, since these will help us go around the + // circle easily + circulator::iterator > + start(border_vect.begin(), border_vect.begin(),border_vect.end()); + circulator::iterator > mid(start), end(start); + + // test the borderless cone + Cmomentum candidate = borderless_cone; + candidate.build_etaphi(); + if (candidate.ref.not_empty()) + test_stability(candidate, border_vect); + + do { + // reset status wrt inclusion in the cone + mid = start; + do { + mid()->is_in = false; + } while (++mid != start); + + // now run over all inclusion possibilities with this starting point + candidate = borderless_cone; + while (++mid != start) { + // will begin with start+1 and go up to start-1 + mid()->is_in = true; + candidate += *(mid()->mom); + test_stability(candidate, border_vect); + } + + } while (++start != end); + + // mid corresponds to momentum that we need to include to get the + // full cone + mid()->is_in = true; + candidate += *(mid()->mom); + test_stability(candidate, border_vect); +} + + +/** + * carry out the computations needed for the stability check of the + * candidate, using the border_vect to indicate which particles + * should / should not be in the stable cone; if the cone is stable + * insert it into the hash. + **********************************************************************/ +void Cstable_cones::test_stability(Cmomentum & candidate, const vector & border_vect) { + + // this almost certainly has not been done... + candidate.build_etaphi(); + + bool stable = true; + for (unsigned i = 0; i < border_vect.size(); i++) { + if (is_inside(&candidate, border_vect[i].mom) ^ (border_vect[i].is_in)) { + stable = false; + break; // it's unstable so there's no point continuing + } + } + + if (stable) hc->insert(&candidate); +} + +/* + * check if we are in a situation of cocircularity. + * if it is the case, update and test in the corresponding way + * return 'false' if no cocircularity detected, 'true' otherwise + * Note that if cocircularity is detected, we need to + * recall 'update' from 'update' !!! + ***************************************************************/ +bool Cstable_cones::cocircular_check(){ + // check if many configurations have the same centre. + // if this is the case, branch on the algorithm for this + // special case. + // Note that those situation, being considered separately in + // test_cone_multiple, must only be considered here if all + // angles are on the same side (this avoid multiple counting) + + if (centre->cocircular.empty()) return false; + + // first get cone into status required at end... + if ((centre->side) && (cone.ref.not_empty())){ + // update cone + cone -= (*child); + + // update info on particles inside + centre->is_inside->cone = false; + + // update stability check quantities + dpt += fabs(child->px)+fabs(child->py); //child->perp2(); + } + + + // now establish the list of unique children in the list + // first make sure parent and child are in! + + list removed_from_cone; + list put_in_border; + list border_list; + + Cmomentum cone_removal; + Cmomentum border = *parent; + border_list.push_back(parent); + + // make sure child appears in the border region + centre->cocircular.push_back(centre); + + // now establish the full contents of the cone minus the cocircular + // region and of the cocircular region itself + for(list::iterator it = centre->cocircular.begin(); + it != centre->cocircular.end(); it++) { + + if ((*it)->is_inside->cone) { + cone_removal += *((*it)->v); + (*it)->is_inside->cone = false; + removed_from_cone.push_back((*it)->is_inside); + } + + // if a point appears twice (i.e. with + and - sign) in the list of + // points on the border, we take care not to include it twice. + // Note that this situation may appear when a point is at a distance + // close to 2R from the parent + if (!(*it)->is_inside->cocirc) { + border += *((*it)->v); + (*it)->is_inside->cocirc = true; + put_in_border.push_back((*it)->is_inside); + border_list.push_back((*it)->v); + } + } + + + // figure out whether this pairing has been observed before + Cmomentum borderless_cone = cone; + borderless_cone -= cone_removal; + bool consider = true; + for (unsigned int i=0;i(borderless_cone.ref, + border.ref)); + + // first figure out whether our cone momentum is good + double local_dpt = fabs(cone_removal.px) + fabs(cone_removal.py); + double total_dpt = dpt + local_dpt; + + recompute_cone_contents_if_needed(borderless_cone, total_dpt); + if (total_dpt == 0) { + // a recomputation has taken place -- so take advantage of this + // and update the member cone momentum + cone = borderless_cone + cone_removal; + dpt = local_dpt; + } + + test_cone_cocircular(borderless_cone, border_list); + } + + + // relabel things that were in the cone but got removed + for(list::iterator is_in = removed_from_cone.begin(); + is_in != removed_from_cone.end(); is_in++) { + (*is_in)->cone = true; + } + + // relabel things that got put into the border + for(list::iterator is_in = put_in_border.begin(); + is_in != put_in_border.end(); is_in++) { + (*is_in)->cocirc = false; + } + + // we're done with everything -- return true to signal to user that we've + // been through the co-circularity rigmarole + return true; +} + + +//////////////////////////////////////////////////////// +// RECOMPUTATION OF CONE CONTENTS // +// - compute_cone_contents() // +// - recompute_cone_contents() // +// - recompute_cone_contents_if_needed() // +//////////////////////////////////////////////////////// + +/** + * compute the cone contents by going once around the full set of + * circles and tracking the entry/exit status each time + * given parent, child and centre compute the momentum + * of the particle inside the cone + * This sets up the inclusion information, which can then be directly + * used to calculate the cone momentum. + **********************************************************************/ +void Cstable_cones::compute_cone_contents() { + circulator::iterator > + start(vicinity.begin()+first_cone, vicinity.begin(), vicinity.end()); + + circulator::iterator > here(start); + + // note that in the following algorithm, the cone contents never includes + // the child. Indeed, if it has positive sign, then it will be set as + // outside at the last step in the loop. If it has negative sign, then the + // loop will at some point go to the corresponding situation with positive + // sign and set the inclusion status to 0. + + do { + // as we leave this position a particle enters if its side is + // negative (i.e. the centre is the one at -ve angle wrt to the + // parent-child line + if (!(*here())->side) ((*here())->is_inside->cone) = 1; + + // move on to the next position + ++here; + + // as we arrive at this position a particle leaves if its side is positive + if ((*here())->side) ((*here())->is_inside->cone) = 0; + } while (here != start); + + // once we've reached the start the 'is_inside' information should be + // 100% complete, so we can use it to calculate the cone contents + // and then exit + recompute_cone_contents(); + return; + +} + + +/* + * compute the cone momentum from particle list. + * in this version, we use the 'pincluded' information + * from the Cvicinity class + */ +void Cstable_cones::recompute_cone_contents(){ + unsigned int i; + + // set momentum to 0 + cone = Cmomentum(); + + // Important note: we can browse only the particles + // in vicinity since all particles in the cone are + // withing a distance 2R w.r.t. parent hence in vicinity. + // Among those, we only add the particles for which 'is_inside' is true ! + // This methos rather than a direct comparison avoids rounding errors + for (i=0;iside) && (vicinity[i]->is_inside->cone)) + cone += *vicinity[i]->v; + } + + // set check variables back to 0 + dpt = 0.0; +} + + +/* + * if we have gone beyond the acceptable threshold of change, compute + * the cone momentum from particle list. in this version, we use the + * 'pincluded' information from the Cvicinity class, but we don't + * change the member cone, only the locally supplied one + */ +void Cstable_cones::recompute_cone_contents_if_needed(Cmomentum & this_cone, + double & this_dpt){ + + if (this_dpt > PT_TSHOLD*(fabs(this_cone.px)+fabs(this_cone.py))) { + if (cone.ref.is_empty()) { + this_cone = Cmomentum(); + } else { + // set momentum to 0 + this_cone = Cmomentum(); + + // Important note: we can browse only the particles + // in vicinity since all particles in the this_cone are + // withing a distance 2R w.r.t. parent hence in vicinity. + // Among those, we only add the particles for which 'is_inside' is true ! + // This methos rather than a direct comparison avoids rounding errors + for (unsigned int i=0;iside) && (vicinity[i]->is_inside->cone)) + this_cone += *vicinity[i]->v; + } + + } + // set check variables back to 0 + this_dpt = 0.0; + } + +} + + +//////////////////////////////////////////////////////// +// VARIOUS TOOLS // +// - circle_intersect() // +// - is_inside() // +// - abs_dangle() // +//////////////////////////////////////////////////////// + + +/* + * circle intersection. + * computes the intersection with a circle of given centre and radius. + * The output takes the form of a checkxor of the intersection's particles + * - cx circle centre x coordinate + * - cy circle centre y coordinate + * return the checkxor for the intersection + ******************************************************************/ +Creference Cstable_cones::circle_intersect(double cx, double cy){ + Creference intersection; + int i; + double dx, dy; + + for (i=0;iM_PI) + dy -= twopi; + + // really check if the distance is less than VR + if (dx*dx+dy*dyeta - v->eta; + dy = fabs(centre->phi - v->phi); + if (dy>M_PI) + dy -= twopi; + + return dx*dx+dy*dyM_PI) + dphi = dphi-twopi; + + return dphi; +} + +} Index: /trunk/SISCone/protocones.h =================================================================== --- /trunk/SISCone/protocones.h (revision 20) +++ /trunk/SISCone/protocones.h (revision 20) @@ -0,0 +1,268 @@ +// -*- C++ -*- +/////////////////////////////////////////////////////////////////////////////// +// File: protocones.h // +// Description: header file for stable cones determination (Cstable_cones) // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:27 $// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef __PROTOCONES_H__ +#define __PROTOCONES_H__ + +#include "momentum.h" +#include "vicinity.h" +#include +#include +#include +#include "hash.h" + + + +namespace siscone{ + +/** + * \class Cborder_store + * + * class for storing a border momentum (in context of co-circularity + * checks). + + * This class essentially calculates angle of border point w.r.t. + * circle center (eta & phi), and provides a store of information + * about whether we are currently including this point in the + * candidate + */ +class Cborder_store{ +public: + /// default ctor + Cborder_store(Cmomentum * momentum, double centre_eta, double centre_phi) : + mom(momentum), is_in(false) { + angle = atan2(mom->phi - centre_phi, mom->eta - centre_eta); + } + + Cmomentum * mom; ///< particle momentum + double angle; ///< angle w.r.t. circle centre + bool is_in; ///< inclusion status of the particle +}; + + +/// allows easy sorting of Cborder_store objects (which need to be +/// ordered in angle). +inline bool operator<(const Cborder_store & a, const Cborder_store & b) { + return a.angle < b.angle; +} + + +/** + * \class Cstable_cones + * \brief Computes the list of stable comes from a particle list. + * + * This class does the first fundamental task of te cone algorithm: + * it is used to compute the list of stable cones given a list + * of particles. + */ +class Cstable_cones : public Cvicinity{ + public: + /// default ctor + Cstable_cones(); + + /// ctor with initialisation (sse init for details) + Cstable_cones(std::vector &_particle_list); + + /// default dtor + ~Cstable_cones(); + + /** + * initialisation + * \param _particle_list list of particles + */ + void init(std::vector &_particle_list); + + /** + * compute stable cones. + * This function really does the job i.e. computes + * the list of stable cones (in a seedless way) + * \param _radius radius of the cones + * \return The number of stable cones found is returned + */ + int get_stable_cones(double _radius); + + /// list of stable cones + std::vector protocones; + + /// list of candidates + hash_cones *hc; + + /// total number of tested cones + int nb_tot; + + protected: + /// cone radius + double R; + + /// cone radius SQUARED + double R2; + + private: + /// cone with a given particle as parent + /// this reduction to a single vector assumes we trust the checksums + Cmomentum cone; + + /// child particle, taken in the 'vicinity' list + Cmomentum *child; + + /// centre of the tested cone + Cvicinity_elm *centre; + + /// index in the particle list; + unsigned int centre_idx; + + /// first cone used in the vicinity list + unsigned int first_cone; + + /** + * initialise the cone. + * We take the first particle in the angular ordering to compute this one + * \return 0 on success, 1 on error + */ + int init_cone(); + + /** + * test cones. + * We check if the cone(s) build with the present parent and child + * are stable + * \return 0 on success 1 on error + */ + int test_cone(); + + /** + * update the cone + * go to the next child for that parent and update 'cone' appropriately + * \return 0 if update candidate found, 1 otherwise + */ + int update_cone(); + + /* + * run through the vicinity of the current parent and for each child + * indicate which members are cocircular... + */ + void prepare_cocircular_lists(); + + /** + * check if we are in a situation of cocircularity. + * if it is the case, update and test in the corresponding way + * \return 'false' if no cocircularity detected, 'true' otherwise + * Note that if cocircularity is detected, we need to + * recall 'update' from 'update' !!! + */ + bool cocircular_check(); + + /** + * Routine for testing cocircular configurations in p^3 time, + * rather than 2^p time; + */ + void test_cone_cocircular(Cmomentum & borderless_cone, + std::list & border_list); + + /** + * carry out the computations needed for the stability check of the + * candidate, using the border_vect to indicate which particles + * should / should not be in the stable cone; if the cone is stable + * insert it into the hash. + */ + void test_stability(Cmomentum & candidate, + const std::vector & border_vect); + + /** + * compute the cone contents by going once around the full set of + * circles and tracking the entry/exit status each time -- this sets + * up the inclusion information, which can then be directly used to + * calculate the cone momentum. + */ + void compute_cone_contents(); + + /** + * compute the cone momentum from particle list. + * in this version, we use the 'pincluded' information + * from the Cviinity class + */ + void recompute_cone_contents(); + + /* + * if we have gone beyond the acceptable threshold of change, compute + * the cone momentum from particle list. in this version, we use the + * 'pincluded' information from the Cvicinity class, but we don't + * change the member cone, only the locally supplied one + */ + void recompute_cone_contents_if_needed(Cmomentum & this_cone, double & this_dpt); + + /** + * compute stability of all enumerated candidates. + * For all candidate cones which are stable w.r.t. their border particles, + * pass the last test: stability with quadtree intersection + */ + int proceed_with_stability(); + + /* + * circle intersection. + * computes the intersection with a circle of given centre and radius. + * The output takes the form of a checkxor of the intersection's particles + * - cx circle centre x coordinate + * - cy circle centre y coordinate + * return the checkxor for the intersection + ******************************************************************/ + Creference circle_intersect(double cx, double cy); + + /// present candidate cone + Cmomentum cone_candidate; + + /// in case of co-circular points, vector for them + std::vector child_list; + + /// list of cocircular enclusures already studied + /// first element if cone contents, second is cone border + std::vector< std::pair > multiple_centre_done; + + // information for updating cone contents to avoid rounding errors + double dpt; ///< sums of Delta P_t + + /** + * test if a particle is inside a cone of given centre. + * check if the particle of coordinates 'v' is inside the circle of radius R + * centered at 'centre'. + * \param centre centre of the circle + * \param v particle to test + * \return true if inside, false if outside + */ + inline bool is_inside(Cmomentum *centre, Cmomentum *v); +}; + +/* + * compute the absolute value of the difference between 2 angles. + * We take care of the 2pi periodicity + * \param angle1 first angle + * \param angle2 second angle + * \return the absolute value of the difference between the angles + *****************************************************************/ +inline double abs_dangle(double &angle1, double &angle2); + +} +#endif Index: /trunk/SISCone/quadtree.cc =================================================================== --- /trunk/SISCone/quadtree.cc (revision 20) +++ /trunk/SISCone/quadtree.cc (revision 20) @@ -0,0 +1,304 @@ +/////////////////////////////////////////////////////////////////////////////// +// File: quadtree.cpp // +// Description: source file for quadtree management (Cquadtree class) // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:27 $// +/////////////////////////////////////////////////////////////////////////////// + +#include "quadtree.h" +#include +#include +#include + +namespace siscone{ + +using namespace std; + +/******************************************************************* + * Cquadtree implementation * + * Implementation of a 2D quadtree. * + * This class implements the traditional two-dimensional quadtree. * + * The elements at each node are of 'Cmomentum' type. * + *******************************************************************/ + +// default ctor +//-------------- +Cquadtree::Cquadtree(){ + v = NULL; + + children[0][0] = children[0][1] = children[1][0] = children[1][1] = NULL; + has_child = false; +} + + +// ctor with initialisation (see init for details) +//-------------------------- +Cquadtree::Cquadtree(double _x, double _y, double _half_size_x, double _half_size_y){ + v = NULL; + + children[0][0] = children[0][1] = children[1][0] = children[1][1] = NULL; + has_child = false; + + init(_x, _y, _half_size_x, _half_size_y); +} + + +// default destructor +// at destruction, everything is destroyed except +// physical values at the leaves +//------------------------------------------------ +Cquadtree::~Cquadtree(){ + if (has_child){ + if (v!=NULL) delete v; + delete children[0][0]; + delete children[0][1]; + delete children[1][0]; + delete children[1][1]; + } +} + + +/* + * init the tree. + * By initializing the tree, we mean setting the cell parameters + * and preparing the object to act as a seed for a new tree. + * - _x x-position of the center + * - _y y-position of the center + * - half_size_x half x-size of the cell + * - half_size_y half y-size of the cell + * return 0 on success, 1 on error. Note that if the cell + * is already filled, we return an error. + ******************************************************************/ +int Cquadtree::init(double _x, double _y, double _half_size_x, double _half_size_y){ + if (v!=NULL) + return 1; + + centre_x = _x; + centre_y = _y; + half_size_x = _half_size_x; + half_size_y = _half_size_y; + + return 0; +} + + +/* + * adding a particle to the tree. + * This method adds one vector to the quadtree structure which + * is updated consequently. + * - v vector to add + * return 0 on success 1 on error + ******************************************************************/ +int Cquadtree::add(Cmomentum *v_add){ + // Description of the method: + // -------------------------- + // the addition process goes as follows: + // 1. check if the cell is empty, in which case, add the particle + // here and leave. + // 2. If there is a unique particle already inside, + // (a) create children + // (b) forward the existing particle to the appropriate child + // 3. Add current particle to this cell and forward to the + // adequate child + // NOTE: we assume in the whole procedure that the particle is + // indeed inside the cell ! + + // step 1: the case of empty cells + if (v==NULL){ + v = v_add; + return 0; + } + + // step 2: additional work if 1! particle already present + // we use the fact that only 1-particle systems have no child + if (!has_child){ + double new_half_size_x = 0.5*half_size_x; + double new_half_size_y = 0.5*half_size_y; + // create children + children[0][0] = new Cquadtree(centre_x-new_half_size_x, centre_y-new_half_size_y, + new_half_size_x, new_half_size_y); + children[0][1] = new Cquadtree(centre_x-new_half_size_x, centre_y+new_half_size_y, + new_half_size_x, new_half_size_y); + children[1][0] = new Cquadtree(centre_x+new_half_size_x, centre_y-new_half_size_y, + new_half_size_x, new_half_size_y); + children[1][1] = new Cquadtree(centre_x+new_half_size_x, centre_y+new_half_size_y, + new_half_size_x, new_half_size_y); + + has_child = true; + + // forward to child + //? The following line assumes 'true'==1 and 'false'==0 + // Note: v being a single particle, eta and phi are correct + children[v->eta>centre_x][v->phi>centre_y]->add(v); + + // copy physical params + v = new Cmomentum(*v); + } + + // step 3: add new particle + // Note: v_add being a single particle, eta and phi are correct + children[v_add->eta>centre_x][v_add->phi>centre_y]->add(v_add); + *v+=*v_add; + + return 0; +} + + +/* + * circle intersection. + * computes the intersection with a circle of given centre and radius. + * The output takes the form of a quadtree with all squares included + * in the circle. + * - cx circle centre x coordinate + * - cy circle centre y coordinate + * - cR2 circle radius SQUARED + * return the checksum for the intersection + ******************************************************************/ +Creference Cquadtree::circle_intersect(double cx, double cy, double cR2){ + // Description of the method: + // -------------------------- + // 1. check if cell is empty => no intersection + // 2. if cell has 1! particle, check if it is inside the circle. + // If yes, add it and return, if not simply return. + // 3. check if the circle intersects the square. If not, return. + // 4. check if the square is inside the circle. + // If yes, add it to qt and return. + // 5. check intersections with children. + + // step 1: if there is no particle inside te square, no reason to go further + if (v==NULL) + return Creference(); + + double dx, dy; + + // step 2: if there is only one particle inside the square, test if it is in + // the circle, in which case return associated reference + if (!has_child){ + // compute the distance + // Note: v has only one particle => eta and phi are defined + dx = cx - v->eta; + dy = fabs(cy - v->phi); + if (dy>M_PI) + dy -= 2.0*M_PI; + + // test distance + if (dx*dx+dy*dyref; + } + + return Creference(); + } + + // step 3: check if there is an intersection + //double ryp, rym; + double dx_c, dy_c; + + // store distance with the centre of the square + dx_c = fabs(cx-centre_x); + dy_c = fabs(cy-centre_y); + if (dy_c>M_PI) dy_c = 2.0*M_PI-dy_c; + + // compute (minimal) the distance (pay attention to the periodicity in phi). + dx = dx_c-half_size_x; + if (dx<0) dx=0; + dy = dy_c-half_size_y; + if (dy<0) dy=0; + + // check the distance + if (dx*dx+dy*dy>=cR2){ + // no intersection + return Creference(); + } + + // step 4: check if included + + // compute the (maximal) distance + dx = dx_c+half_size_x; + dy = dy_c+half_size_y; + if (dy>M_PI) dy = M_PI; + + // compute the distance + if (dx*dx+dy*dyref; + } + + // step 5: the square is not fully in. Recurse to children + return children[0][0]->circle_intersect(cx, cy, cR2) + + children[0][1]->circle_intersect(cx, cy, cR2) + + children[1][0]->circle_intersect(cx, cy, cR2) + + children[1][1]->circle_intersect(cx, cy, cR2); +} + + +/* + * output a data file for drawing the grid. + * This can be used to output a data file containing all the + * grid subdivisions. The file contents is as follows: + * first and second columns give center of the cell, the third + * gives the size. + * - flux opened stream to write to + * return 0 on success, 1 on error + ******************************************************************/ +int Cquadtree::save(FILE *flux){ + + if (flux==NULL) + return 1; + + if (has_child){ + fprintf(flux, "%le\t%le\t%le\t%le\n", centre_x, centre_y, half_size_x, half_size_y); + children[0][0]->save(flux); + children[0][1]->save(flux); + children[1][0]->save(flux); + children[1][1]->save(flux); + } + + return 0; +} + + +/* + * output a data file for drawing the tree leaves. + * This can be used to output a data file containing all the + * tree leaves. The file contents is as follows: + * first and second columns give center of the cell, the third + * gives the size. + * - flux opened stream to write to + * return 0 on success, 1 on error + ******************************************************************/ +int Cquadtree::save_leaves(FILE *flux){ + + if (flux==NULL) + return 1; + + if (has_child){ + if (children[0][0]!=NULL) children[0][0]->save_leaves(flux); + if (children[0][1]!=NULL) children[0][1]->save_leaves(flux); + if (children[1][0]!=NULL) children[1][0]->save_leaves(flux); + if (children[1][1]!=NULL) children[1][1]->save_leaves(flux); + } else { + fprintf(flux, "%le\t%le\t%le\t%le\n", centre_x, centre_y, half_size_x, half_size_y); + } + + return 0; +} + +} Index: /trunk/SISCone/quadtree.h =================================================================== --- /trunk/SISCone/quadtree.h (revision 20) +++ /trunk/SISCone/quadtree.h (revision 20) @@ -0,0 +1,124 @@ +// -*- C++ -*- +/////////////////////////////////////////////////////////////////////////////// +// File: quadtree.h // +// Description: header file for quadtree management (Cquadtree class) // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:27 $// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef __QUADTREE_H__ +#define __QUADTREE_H__ + +#include "momentum.h" +#include + +namespace siscone{ + +/** + * \class Cquadtree + * \brief Implementation of a 2D quadtree. + * + * This class implements the traditional two-dimensional quadtree. + * The elements at each node are of 'Cmomentum' type. + */ +class Cquadtree{ + public: + /// default ctor + Cquadtree(); + + /// ctor with initialisation (see init for details) + Cquadtree(double _x, double _y, double _half_size_x, double _half_size_y); + + /// default destructor + /// at destruction, everything is destroyed except + /// physical values at the leaves + ~Cquadtree(); + + /** + * init the tree. + * By initializing the tree, we mean setting the cell parameters + * and preparing the object to act as a seed for a new tree. + * \param _x x-position of the center + * \param _y y-position of the center + * \param _half_size_x x-size of the cell + * \param _half_size_y y-size of the cell + * \return 0 on success, 1 on error. Note that if the cell or its + * parent is already filled, we return an error. + */ + int init(double _x, double _y, double _half_size_x, double _half_size_y); + + /** + * adding a particle to the tree. + * This method adds one vector to the quadtree structure which + * is updated consequently. + * \param v_add vector to add + * \return 0 on success 1 on error + */ + int add(Cmomentum *v_add); + + /** + * circle intersection. + * computes the intersection with a circle of given centre and radius. + * The output takes the form of a quadtree with all squares included + * in the circle. + * \param cx circle centre x coordinate + * \param cy circle centre y coordinate + * \param cR2 circle radius SQUARED + * \return the checksum for that intersection + */ + Creference circle_intersect(double cx, double cy, double cR2); + + /** + * output a data file for drawing the grid. + * This can be used to output a data file containing all the + * grid subdivisions. The file contents is as follows: + * first and second columns give center of the cell, the third + * gives the size. + * \param flux opened stream to write to + * \return 0 on success, 1 on error + */ + int save(FILE *flux); + + /** + * output a data file for drawing the tree leaves. + * This can be used to output a data file containing all the + * tree leaves. The file contents is as follows: + * first and second columns give center of the cell, the third + * gives the size. + * \param flux opened stream to write to + * \return 0 on success, 1 on error + */ + int save_leaves(FILE *flux); + + double centre_x; ///< x-position of the centre of the cell + double centre_y; ///< y-position of the centre of the cell + double half_size_x; ///< HALF size of the cell + double half_size_y; ///< HALF size of the cell + + Cmomentum *v; ///< physical contents + + Cquadtree* children[2][2]; ///< sub-cells ( 0,1->left-right; 0,1->bottom,top) + bool has_child; ///< true if not a leaf +}; + +} +#endif Index: /trunk/SISCone/ranlux.cc =================================================================== --- /trunk/SISCone/ranlux.cc (revision 20) +++ /trunk/SISCone/ranlux.cc (revision 20) @@ -0,0 +1,171 @@ +// file: ranlux.xpp +#include "ranlux.h" +#include +#include + +/* This is a lagged fibonacci generator with skipping developed by Luescher. + The sequence is a series of 24-bit integers, x_n, + + x_n = d_n + b_n + + where d_n = x_{n-10} - x_{n-24} - c_{n-1}, b_n = 0 if d_n >= 0 and + b_n = 2^24 if d_n < 0, c_n = 0 if d_n >= 0 and c_n = 1 if d_n < 0, + where after 24 samples a group of p integers are "skipped", to + reduce correlations. By default p = 199, but can be increased to + 365. + + The period of the generator is around 10^171. + + From: M. Luescher, "A portable high-quality random number generator + for lattice field theory calculations", Computer Physics + Communications, 79 (1994) 100-110. + + Available on the net as hep-lat/9309020 at http://xxx.lanl.gov/ + + See also, + + F. James, "RANLUX: A Fortran implementation of the high-quality + pseudo-random number generator of Luscher", Computer Physics + Communications, 79 (1994) 111-114 + + Kenneth G. Hamilton, F. James, "Acceleration of RANLUX", Computer + Physics Communications, 101 (1997) 241-248 + + Kenneth G. Hamilton, "Assembler RANLUX for PCs", Computer Physics + Communications, 101 (1997) 249-253 */ + +namespace siscone{ + +static const unsigned long int mask_lo = 0x00ffffffUL; // 2^24 - 1 +static const unsigned long int mask_hi = ~0x00ffffffUL; +static const unsigned long int two24 = 16777216; // 2^24 + + +// internal generator structure +//------------------------------ +typedef struct { + unsigned int i; + unsigned int j; + unsigned int n; + unsigned int skip; + unsigned int carry; + unsigned long int u[24]; +} ranlux_state_t; + + +// internal generator state +//-------------------------- +ranlux_state_t local_ranlux_state; + + +// incrementation of the generator state +//--------------------------------------- +static inline unsigned long int increment_state(){ + unsigned int i = local_ranlux_state.i; + unsigned int j = local_ranlux_state.j; + long int delta = local_ranlux_state.u[j] - local_ranlux_state.u[i] + - local_ranlux_state.carry; + + if (delta & mask_hi){ + local_ranlux_state.carry = 1; + delta &= mask_lo; + } else { + local_ranlux_state.carry = 0; + } + + local_ranlux_state.u[i] = delta; + + if (i==0) + i = 23; + else + i--; + + local_ranlux_state.i = i; + + if (j == 0) + j = 23; + else + j--; + + local_ranlux_state.j = j; + + return delta; +} + + +// set generator state +//--------------------- +static void ranlux_set(unsigned long int s){ + int i; + long int seed; + + if (s==0) + s = 314159265; /* default seed is 314159265 */ + + seed = s; + + /* This is the initialization algorithm of F. James, widely in use + for RANLUX. */ + + for (i=0;i<24;i++){ + unsigned long int k = seed/53668; + seed = 40014*(seed-k*53668)-k*12211; + if (seed<0){ + seed += 2147483563; + } + local_ranlux_state.u[i] = seed%two24; + } + + local_ranlux_state.i = 23; + local_ranlux_state.j = 9; + local_ranlux_state.n = 0; + local_ranlux_state.skip = 389-24; // 389 => best decorrelation + + if (local_ranlux_state.u[23]&mask_hi){ + local_ranlux_state.carry = 1; + } else { + local_ranlux_state.carry = 0; + } +} + + +// generator initialization +//-------------------------- +void ranlux_init(){ + // seed the generator + ranlux_set(0); +} + + +// get random number +//------------------- +unsigned long int ranlux_get(){ + const unsigned int skip = local_ranlux_state.skip; + unsigned long int r = increment_state(); + + local_ranlux_state.n++; + + if (local_ranlux_state.n == 24){ + unsigned int i; + local_ranlux_state.n = 0; + for (i = 0; i < skip; i++) + increment_state(); + } + + return r; +} + +// print generator state +//----------------------- +void ranlux_print_state(){ + size_t i; + unsigned char *p = (unsigned char *) (&local_ranlux_state); + const size_t n = sizeof (ranlux_state_t); + + for (i=0;i= 0 and + b_n = 2^24 if d_n < 0, c_n = 0 if d_n >= 0 and c_n = 1 if d_n < 0, + where after 24 samples a group of p integers are "skipped", to + reduce correlations. By default p = 199, but can be increased to + 365. + + The period of the generator is around 10^171. + + From: M. Luescher, "A portable high-quality random number generator + for lattice field theory calculations", Computer Physics + Communications, 79 (1994) 100-110. + + Available on the net as hep-lat/9309020 at http://xxx.lanl.gov/ + + See also, + + F. James, "RANLUX: A Fortran implementation of the high-quality + pseudo-random number generator of Luscher", Computer Physics + Communications, 79 (1994) 111-114 + + Kenneth G. Hamilton, F. James, "Acceleration of RANLUX", Computer + Physics Communications, 101 (1997) 241-248 + + Kenneth G. Hamilton, "Assembler RANLUX for PCs", Computer Physics + Communications, 101 (1997) 249-253 */ + +namespace siscone{ + +/// initialize 'ranlux' generator +void ranlux_init(); + +/// generate random value (24 bits) +unsigned long int ranlux_get(); + +/// save state of the generator +void ranlux_print_state(); + +} +#endif Index: /trunk/SISCone/reference.cc =================================================================== --- /trunk/SISCone/reference.cc (revision 20) +++ /trunk/SISCone/reference.cc (revision 20) @@ -0,0 +1,127 @@ +/////////////////////////////////////////////////////////////////////////////// +// File: reference.cpp // +// Description: source file for checkxor management (Creference class) // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:28 $// +/////////////////////////////////////////////////////////////////////////////// + +#include "reference.h" +#include "ranlux.h" +#include + +namespace siscone{ + +/******************************************************* + * Creference implementation * + * references used for checksums. * + * * + * This class implements some reference variable * + * that can be used for checksums. Those checksums * + * are useful to disentengle between contents of two * + * cones without looking into their explicit particle * + * contents. * + *******************************************************/ + +// default constructor +////////////////////// +Creference::Creference(){ + ref[0] = ref[1] = ref[2] = 0; +} + + //static unsigned int reference_bit = 1; + +// create a random reference +//--------------------------- +void Creference::randomize(){ +// ref[0] = reference_bit; +// ref[1] = 0; +// ref[2] = 0; +// reference_bit <<= 1; + + unsigned int r1 = ranlux_get(); + unsigned int r2 = ranlux_get(); + unsigned int r3 = ranlux_get(); + unsigned int r4 = ranlux_get(); + // since ranlux only produces 24 bits, take r4 and add 8 bits + // from it to each of r1,r2, r3 to get 3*32 bits. + ref[0] = r1+((r4 & 0x00ff0000) << 8); + ref[1] = r2+((r4 & 0x0000ff00) << 16); + ref[2] = r3+((r4 & 0x000000ff) << 24); + + if (is_empty()) randomize(); +} + +// test emptyness +//---------------- +bool Creference::is_empty(){ + return (ref[0]==0) && (ref[1]==0) && (ref[2]==0); +} + +// test non-emptyness +//-------------------- +bool Creference::not_empty(){ + return (ref[0]!=0) || (ref[1]!=0) || (ref[2]!=0); +} + +// assignment of reference +//------------------------- +Creference& Creference::operator = (const Creference &r){ + ref[0] = r.ref[0]; + ref[1] = r.ref[1]; + ref[2] = r.ref[2]; + return *this; +} + +// addition of reference +//----------------------- +Creference Creference::operator + (const Creference &r){ + Creference tmp = *this; + return tmp+=r; +} + +// incrementation of reference +//----------------------------- +Creference& Creference::operator += (const Creference &r){ + ref[0] ^= r.ref[0]; + ref[1] ^= r.ref[1]; + ref[2] ^= r.ref[2]; + return *this; +} + +// decrementation of reference +//----------------------------- +Creference& Creference::operator -= (const Creference &r){ + ref[0] ^= r.ref[0]; + ref[1] ^= r.ref[1]; + ref[2] ^= r.ref[2]; + return *this; +} + +// addition with 2 references +//---------------------------- +Creference operator + (Creference &r1, Creference &r2){ + Creference tmp = r1; + return r1+=r2; +} + +} + Index: /trunk/SISCone/reference.h =================================================================== --- /trunk/SISCone/reference.h (revision 20) +++ /trunk/SISCone/reference.h (revision 20) @@ -0,0 +1,111 @@ +// -*- C++ -*- +/////////////////////////////////////////////////////////////////////////////// +// File: reference.h // +// Description: header file for checkxor management (Creference class) // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:28 $// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef __REFERENCE_H__ +#define __REFERENCE_H__ + +namespace siscone{ + +/** + * \class Creference + * \brief references used for checksums. + * + * This class implements some reference variable + * that can be used for checksums. Those checksums + * are useful to disentengle between contents of two + * cones without looking into their explicit particle + * contents. + */ +class Creference{ + public: + /// default constructor + Creference(); + + /// create a random reference + void randomize(); + + /// test emptyness + bool is_empty(); + + /// test non-emptyness + bool not_empty(); + + /// assignment of reference + Creference& operator = (const Creference &r); + + /// addition of reference + Creference operator + (const Creference &r); + + /// incrementation of reference + Creference& operator += (const Creference &r); + + /// decrementation of reference + Creference& operator -= (const Creference &r); + + /// accessing the reference + inline unsigned int operator[] (int i) {return ref[i];} + + unsigned int ref[3]; ///< actual data for the reference +}; + +/// addition of two references +Creference operator + (Creference &r1, Creference &r2); + +/// equality test of two references +bool operator == (const Creference &r1, const Creference &r2); + +/// difference test of two references +bool operator != (const Creference &r1, const Creference &r2); + +/// ordering of two references +bool operator < (const Creference &r1, const Creference &r2); + + +//=============== inline material ================ + +// equality test for two references +//---------------------------------- +inline bool operator == (const Creference &r1, const Creference &r2){ + return (r1.ref[0]==r2.ref[0]) && (r1.ref[1]==r2.ref[1]) && (r1.ref[2]==r2.ref[2]); +} + +// difference test for two references +//---------------------------------- +inline bool operator != (const Creference &r1, const Creference &r2){ + return (r1.ref[0]!=r2.ref[0]) || (r1.ref[1]!=r2.ref[1]) || (r1.ref[2]!=r2.ref[2]); +} + +// difference test for two references +//---------------------------------- +inline bool operator < (const Creference &r1, const Creference &r2){ + return (r1.ref[0] +#include +#include + +namespace siscone{ +using namespace std; + +/*************************************************************** + * Csiscone implementation * + * final class: gather everything to compute the jet contents. * + * * + * This is the class user should use. * + * It computes the jet contents of a list of particles * + * given a cone radius and a threshold for splitting/merging. * + ***************************************************************/ + +// default ctor +//-------------- +Csiscone::Csiscone(){ + rerun_allowed = false; +} + +// default dtor +//-------------- +Csiscone::~Csiscone(){ + rerun_allowed = false; +} + +bool Csiscone::init_done=false; + +/* + * compute the jets from a given particle set doing multiple passes + * such pass N looks for jets among all particles not put into jets + * during previous passes. + * - _particles list of particles + * - _radius cone radius + * - _f shared energy threshold for splitting&merging + * - _n_pass_max maximum number of runs + * - _ptmin minimum pT of the protojets + * - _split_merge_scale the scale choice for the split-merge procedure + * NOTE: using pt leads to IR unsafety for some events with momentum + * conservation. So we strongly advise not to change the default + * value. + * return the number of jets found. + **********************************************************************/ +int Csiscone::compute_jets(vector &_particles, double _radius, double _f, + int _n_pass_max, double _ptmin, + Esplit_merge_scale _split_merge_scale){ + // initialise random number generator + if (!init_done){ + // initialise random number generator + ranlux_init(); + + // print the banner + cout << "#ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << endl; + cout << "# SISCone version " << setw(28) << left << siscone_version() << "o" << endl; + cout << "# http://projects.hepforge.org/siscone o" << endl; + cout << "# o" << endl; + cout << "# This is SISCone: the Seedless Infrared Safe Cone Jet Algorithm o" << endl; + cout << "# SISCone was written by Gavin Salam and Gregory Soyez o" << endl; + cout << "# It is released under the terms of the GNU General Public License o" << endl; + cout << "# o" << endl; + cout << "# A description of the algorithm is available in the publication o" << endl; + cout << "# JHEP 05 (2007) 086 [arXiv:0704.0292 (hep-ph)]. o" << endl; + cout << "# Please cite it if you use SISCone. o" << endl; + cout << "#ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << endl; + cout << endl; + // do not do this again + init_done=true; + } + + // run some general safety tests (NB: f will be checked in split-merge) + if (_radius <= 0.0 || _radius >= 0.5*M_PI) { + ostringstream message; + message << "Illegal value for cone radius, R = " << _radius + << " (legal values are 00) && (_n_pass_max!=0)); + + rerun_allowed = true; + + // split & merge + return perform(_f, _ptmin); +} + +/* + * recompute the jets with a different overlap parameter. + * we use the same particles and R as in the preceeding call. + * - _f shared energy threshold for splitting&merging + * - _ptmin minimum pT of the protojets + * - _split_merge_scale the scale choice for the split-merge procedure + * NOTE: using pt leads to IR unsafety for some events with momentum + * conservation. So we strongly advise not to change the default + * value. + * return the number of jets found, -1 if recomputation not allowed. + ********************************************************************/ +int Csiscone::recompute_jets(double _f, double _ptmin, + Esplit_merge_scale _split_merge_scale){ + if (!rerun_allowed) + return -1; + + ptcomparison.split_merge_scale = _split_merge_scale; + + // restore particle list + partial_clear(); + init_pleft(); + + // initialise split/merge algorithm + unsigned int i; + for (i=0;i &_particles, double _radius, double _f, + int _n_pass_max=0, double _ptmin=0.0, + Esplit_merge_scale _split_merge_scale=SM_pttilde); + + /** + * recompute the jets with a different overlap parameter. + * we use the same particles and R as in the preceeding call. + * \param _f shared energy threshold for splitting&merging + * \param _ptmin minimum pT of the protojets + * \param _split_merge_scale the scale choice for the split-merge procedure + * split--merge variable + * NOTE: using pt leads to IR unsafety for some events with momentum + * conservation. So we strongly advise not to change the default + * value. + * \return the number of jets found, -1 if recomputation not allowed. + */ + int recompute_jets(double _f, double _ptmin = 0.0, + Esplit_merge_scale _split_merge_scale=SM_pttilde); + + /// list of protocones found pass-by-pass + std::vector > protocones_list; + + // random number initialisation + static bool init_done; ///< check random generator initialisation + + private: + bool rerun_allowed; ///< is recompute_jets allowed ? +}; + + +// finally, a bunch of functions to access to +// basic information (package name, version) +//--------------------------------------------- + +/** + * return SISCone package name. + * This is nothing but "SISCone", it is a replacement to the + * PACKAGE_NAME string defined in config.h and which is not + * public by default. + * \return the SISCone name as a string + */ +std::string siscone_package_name(); + +/** + * return SISCone version number. + * \return a string of the form "X.Y.Z" with possible additional tag + * (alpha, beta, devel) to mention stability status + */ +std::string siscone_version(); + +} +#endif Index: /trunk/SISCone/siscone_error.cc =================================================================== --- /trunk/SISCone/siscone_error.cc (revision 20) +++ /trunk/SISCone/siscone_error.cc (revision 20) @@ -0,0 +1,33 @@ +/////////////////////////////////////////////////////////////////////////////// +// File: siscone_error.cpp // +// Description: source file for SISCone error messages (Csiscone_error) // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:29 $// +/////////////////////////////////////////////////////////////////////////////// + +#include "siscone_error.h" + +namespace siscone{ + +bool Csiscone_error::m_print_errors = true; + +} Index: /trunk/SISCone/siscone_error.h =================================================================== --- /trunk/SISCone/siscone_error.h (revision 20) +++ /trunk/SISCone/siscone_error.h (revision 20) @@ -0,0 +1,57 @@ +// -*- C++ -*- +/////////////////////////////////////////////////////////////////////////////// +// File: siscone_error.h // +// Description: header file for SISCone error messages (Csiscone_error) // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:29 $// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef __SISCONE_ERROR_H__ +#define __SISCONE_ERROR_H__ + +#include +#include + +namespace siscone{ + +/// class corresponding to errors that will be thrown by siscone +class Csiscone_error { +public: + // constructors + Csiscone_error() {;}; + Csiscone_error(const std::string & message) { + m_message = message; + if (m_print_errors) std::cerr << "siscone::Csiscone_error: "< +#include // for max +#include +#include +#include +#include + +namespace siscone{ + +using namespace std; + +/******************************************************** + * class Cjet implementation * + * real Jet information. * + * This class contains information for one single jet. * + * That is, first, its momentum carrying information * + * about its centre and pT, and second, its particle * + * contents * + ********************************************************/ +// default ctor +//-------------- +Cjet::Cjet(){ + n = 0; + v = Cmomentum(); + pt_tilde = 0.0; + sm_var2 = 0.0; +} + +// default dtor +//-------------- +Cjet::~Cjet(){ + +} + +// ordering of jets in pt (e.g. used in final jets ordering) +//----------------------------------------------------------- +bool jets_pt_less(const Cjet &j1, const Cjet &j2){ + return j1.v.perp2() > j2.v.perp2(); +} + + +/******************************************************** + * Csplit_merge_ptcomparison implementation * + * This deals with the ordering of the jets candidates * + ********************************************************/ + +// odering of two jets +// The variable of the ordering is pt or mt +// depending on 'split_merge_scale' choice +// +// with EPSILON_SPLITMERGE defined, this attempts to identify +// delicate cases where two jets have identical momenta except for +// softish particles -- the difference of pt's may not be correctly +// identified normally and so lead to problems for the fate of the +// softish particle. +// +// NB: there is a potential issue in momentum-conserving events, +// whereby the harder of two jets becomes ill-defined when a soft +// particle is emitted --- this may have a knock-on effect on +// subsequent split-merge steps in cases with sufficiently large R +// (but we don't know what the limit is...) +//------------------------------------------------------------------ +bool Csplit_merge_ptcomparison::operator ()(const Cjet &jet1, const Cjet &jet2) const{ + double q1, q2; + + // compute the value for comparison for both jets + // This depends on the choice of variable (mt is the default) + q1 = jet1.sm_var2; + q2 = jet2.sm_var2; + + bool res = q1 > q2; + + // if we enable the refined version of the comparison (see defines.h), + // we compute the difference more precisely when the two jets are very + // close in the ordering variable. +#ifdef EPSILON_SPLITMERGE + if ( (fabs(q1-q2) < EPSILON_SPLITMERGE*max(q1,q2)) && + (jet1.v.ref != jet2.v.ref) ) { + // get the momentum of the difference + Cmomentum difference; + double pt_tilde_difference; + get_difference(jet1,jet2,&difference,&pt_tilde_difference); + + // use the following relation: pt1^2 - pt2^2 = (pt1+pt2)*(pt1-pt2) + double qdiff; + Cmomentum sum = jet1.v ; + sum += jet2.v; + double pt_tilde_sum = jet1.pt_tilde + jet2.pt_tilde; + + // depending on the choice of ordering variable, set the result + switch (split_merge_scale){ + case SM_mt: + qdiff = sum.E*difference.E - sum.pz*difference.pz; + break; + case SM_pt: + qdiff = sum.px*difference.px + sum.py*difference.py; + break; + case SM_pttilde: + qdiff = pt_tilde_sum*pt_tilde_difference; + break; + case SM_Et: + // diff = E^2 (dpt^2 pz^2- pt^2 dpz^2) + // + dE^2 (pt^2+pz^2) pt2^2 + // where, unless explicitely specified the variables + // refer to the first jet or differences jet1-jet2. + qdiff = jet1.v.E*jet1.v.E* + ((sum.px*difference.px + sum.py*difference.py)*jet1.v.pz*jet1.v.pz + -jet1.v.perp2()*sum.pz*difference.pz) + +sum.E*difference.E*(jet1.v.perp2()+jet1.v.pz*jet1.v.pz)*jet2.v.perp2(); + break; + default: + throw Csiscone_error("Unsupported split-merge scale choice: " + + SM_scale_name()); + } + res = qdiff > 0; + } +#endif // EPSILON_SPLITMERGE + + return res; +} + + +/// return a name for the sm scale choice +/// NB: used internally and also by fastjet +std::string split_merge_scale_name(Esplit_merge_scale sms) { + switch(sms) { + case SM_pt: + return "pt (IR unsafe)"; + case SM_Et: + return "Et (boost dep.)"; + case SM_mt: + return "mt (IR safe except for pairs of identical decayed heavy particles)"; + case SM_pttilde: + return "pttilde (scalar sum of pt's)"; + default: + return "[SM scale without a name]"; + } +} + + +// get the difference between 2 jets +// - j1 first jet +// - j2 second jet +// - v jet1-jet2 +// - pt_tilde jet1-jet2 pt_tilde +// return true if overlapping, false if disjoint +//----------------------------------------------- +void Csplit_merge_ptcomparison::get_difference(const Cjet &j1, const Cjet &j2, Cmomentum *v, double *pt_tilde) const { + int i1,i2; + + // initialise + i1=i2=0; + *v = Cmomentum(); + *pt_tilde = 0.0; + + // compute overlap + // at the same time, we store union in indices + do{ + if (j1.contents[i1]==j2.contents[i2]) { + i1++; + i2++; + } else if (j1.contents[i1]j2.contents[i2]){ + (*v) -= (*particles)[j2.contents[i2]]; + (*pt_tilde) -= (*pt)[j2.contents[i2]]; + i2++; + } else { + throw Csiscone_error("get_non_overlap reached part it should never have seen..."); + } + } while ((i1(ptcomparison)); + + // no hardest cut (col-unsafe) + SM_var2_hardest_cut_off = -1.0; + + // no pt cutoff for the particles to put in p_uncol_hard + stable_cone_soft_pt2_cutoff = -1.0; +} + + +// default dtor +//-------------- +Csplit_merge::~Csplit_merge(){ + full_clear(); +} + + +// initialisation function +// - _particles list of particles +// - protocones list of protocones (initial jet candidates) +// - R2 cone radius (squared) +// - ptmin minimal pT allowed for jets +//------------------------------------------------------------- +int Csplit_merge::init(vector &_particles, vector *protocones, double R2, double ptmin){ + // browse protocones + return add_protocones(protocones, R2, ptmin); +} + + +// initialisation function for particle list +// - _particles list of particles +//------------------------------------------------------------- +int Csplit_merge::init_particles(vector &_particles){ + full_clear(); + + // compute the list of particles + // here, we do not need to throw away particles + // with infinite rapidity (colinear with the beam) + particles = _particles; + n = particles.size(); + + // build the vector of particles' pt and determine min,max eta + double eta_min=0.0; /// for the Ceta_phi_range static member! + double eta_max=0.0; /// for the Ceta_phi_range static member! + pt.resize(n); + for (int i=0;i(ptcomparison)); + + // start off with huge number + most_ambiguous_split = numeric_limits::max(); + + jets.clear(); +#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES + if (merge_identical_protocones) + cand_refs.clear(); +#endif + + p_remain.clear(); + + return 0; +} + + +// full clearance +//---------------- +int Csplit_merge::full_clear(){ + partial_clear(); + + // clear previously allocated memory + if (indices != NULL){ + delete[] indices; + } + particles.clear(); + + return 0; +} + + +// build the list 'p_uncol_hard' from p_remain by clustering collinear particles +// note that thins in only used for stable-cone detection +// so the parent_index field is unnecessary +//------------------------------------------------------------------------- +int Csplit_merge::merge_collinear_and_remove_soft(){ + int i,j; + vector p_sorted; + bool collinear; + double dphi; + + p_uncol_hard.clear(); + + // we first sort the particles according to their rapidity + for (i=0;iM_PI) dphi = twopi-dphi; + if (dphi *protocones, double R2, double ptmin){ + int i; + Cmomentum *c; + Cmomentum *v; + double eta, phi; + double dx, dy; + double R; + Cjet jet; + + if (protocones->size()==0) + return 1; + + pt_min2 = ptmin*ptmin; + R = sqrt(R2); + + // browse protocones + // for each of them, build the list of particles in them + for (vector::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){ + // initialise variables + c = &(*p_it); + + // note: cones have been tested => their (eta,phi) coordinates are computed + eta = c->eta; + phi = c->phi; + + // browse particles to create cone contents + // note that jet is always initialised with default values at this level + jet.v = Cmomentum(); + jet.pt_tilde=0; + jet.contents.clear(); + for (i=0;ipz)!=v->E){ + dx = eta - v->eta; + dy = fabs(phi - v->phi); + if (dy>M_PI) + dy -= twopi; + if (dx*dx+dy*dyparent_index); + jet.v+= *v; + jet.pt_tilde+= pt[v->parent_index]; + v->index=0; + } + } + jet.n=jet.contents.size(); + + // set the momentum in protocones + // (it was only known through eta and phi up to now) + *c = jet.v; + c->eta = eta; // restore exact original coords + c->phi = phi; // to avoid rounding error inconsistencies + + // set the jet range + jet.range=Ceta_phi_range(eta,phi,R); + +#ifdef DEBUG_SPLIT_MERGE + cout << "adding jet: "; + for (int i2=0;i2size()==0) + return 0; + + if (overlap_tshold>=1.0 || overlap_tshold <= 0) { + ostringstream message; + message << "Illegal value for overlap_tshold, f = " << overlap_tshold; + message << " (legal values are 0size()>0){ + // browse for the first jet + j1 = candidates->begin(); + + // if hardest jet does not pass threshold then nothing else will + // either so one stops the split merge. + if (j1->sm_var2end()){ +#ifdef DEBUG_SPLIT_MERGE + show(); +#endif + // check overlapping + if (get_overlap(*j1, *j2, &overlap2)){ + // check if overlapping energy passes threshold + // Note that this depends on the ordering variable +#ifdef DEBUG_SPLIT_MERGE + cout << "overlap between cdt 1 and cdt " << j2_relindex+1 << " with overlap " + << sqrt(overlap2/j2->sm_var2) << endl<sm_var2){ + // split jets + split(j1, j2); + + // update iterators + j2 = j1 = candidates->begin(); + j2_relindex = 0; + } else { + // merge jets + merge(j1, j2); + + // update iterators + j2 = j1 = candidates->begin(); + j2_relindex = 0; + } + } + // watch out: split/merge might have caused new jets with pt < + // ptmin to disappear, so the total number of jets may + // have changed by more than expected and j2 might already by + // the end of the candidates list... + j2_relindex++; + if (j2 != candidates->end()) j2++; + } // end of loop on the second jet + + if (j1 != candidates->end()) { + // all "second jet" passed without overlapping + // (otherwise we won't leave the j2 loop) + // => store jet 1 as real jet + jets.push_back(*j1); + jets[jets.size()-1].v.build_etaphi(); + // a bug where the contents has non-zero size has been cropping + // up in many contexts -- so catch it! + assert(j1->contents.size() > 0); + jets[jets.size()-1].pass = particles[j1->contents[0]].index; +#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES + cand_refs.erase(j1->v.ref); +#endif + candidates->erase(j1); + + //// test that the hardest jet pass the potential cut-off + //if ((candidates->size()!=0) && + // (candidates->begin()->sm_var2clear(); + //} + } + } + } while (candidates->size()>0); + + // sort jets by pT + sort(jets.begin(), jets.end(), jets_pt_less); +#ifdef DEBUG_SPLIT_MERGE + show(); +#endif + + return jets.size(); +} + + + +// save the event on disk +// - flux stream used to save jet contents +//-------------------------------------------- +int Csplit_merge::save_contents(FILE *flux){ + jet_iterator it_j; + Cjet *j1; + int i1, i2; + + fprintf(flux, "# %d jets found\n", (int) jets.size()); + fprintf(flux, "# columns are: eta, phi, pt and number of particles for each jet\n"); + for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){ + j1 = &(*it_j); + j1->v.build_etaphi(); + fprintf(flux, "%lf\t%lf\t%le\t%d\n", + j1->v.eta, j1->v.phi, j1->v.perp(), j1->n); + } + + fprintf(flux, "# jet contents\n"); + fprintf(flux, "# columns are: eta, phi, pt, particle index and jet number\n"); + for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){ + j1 = &(*it_j); + for (i2=0;i2n;i2++) + fprintf(flux, "%lf\t%lf\t%le\t%d\t%d\n", + particles[j1->contents[i2]].eta, particles[j1->contents[i2]].phi, + particles[j1->contents[i2]].perp(), j1->contents[i2], i1); + } + + return 0; +} + + +// show current jets/candidate status +//------------------------------------ +int Csplit_merge::show(){ + jet_iterator it_j; + cjet_iterator it_c; + Cjet *j; + const Cjet *c; + int i1, i2; + + for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){ + j = &(*it_j); + fprintf(stdout, "jet %2d: %le\t%le\t%le\t%le\t", i1+1, + j->v.px, j->v.py, j->v.pz, j->v.E); + for (i2=0;i2n;i2++) + fprintf(stdout, "%d ", j->contents[i2]); + fprintf(stdout, "\n"); + } + + for (it_c = candidates->begin(), i1=0 ; it_c != candidates->end() ; it_c++, i1++){ + c = &(*it_c); + fprintf(stdout, "cdt %2d: %le\t%le\t%le\t%le\t%le\t", i1+1, + c->v.px, c->v.py, c->v.pz, c->v.E, sqrt(c->sm_var2)); + for (i2=0;i2n;i2++) + fprintf(stdout, "%d ", c->contents[i2]); + fprintf(stdout, "\n"); + } + + fprintf(stdout, "\n"); + return 0; +} + + +// get the overlap between 2 jets +// - j1 first jet +// - j2 second jet +// - overlap2 returned overlap^2 (determined by the choice of SM variable) +// return true if overlapping, false if disjoint +//--------------------------------------------------------------------- +bool Csplit_merge::get_overlap(const Cjet &j1, const Cjet &j2, double *overlap2){ + // check if ranges overlap + if (!is_range_overlap(j1.range,j2.range)) + return false; + + int i1,i2; + bool is_overlap; + + // initialise + i1=i2=idx_size=0; + is_overlap = false; + Cmomentum v; + double pt_tilde=0.0; + + // compute overlap + // at the same time, we store union in indices + do{ + if (j1.contents[i1]j2.contents[i2]){ + indices[idx_size] = j2.contents[i2]; + i2++; + } else { // (j1.contents[i1]==j2.contents[i2]) + v += particles[j1.contents[i1]]; + pt_tilde += pt[j1.contents[i1]]; + indices[idx_size] = j1.contents[i1]; + i1++; + i2++; + is_overlap = true; + } + idx_size++; + } while ((i1" all over the place + const Cjet & j1 = * it_j1; + const Cjet & j2 = * it_j2; + + i1=i2=0; + jet2.v = jet1.v = Cmomentum(); + jet2.pt_tilde = jet1.pt_tilde = 0.0; + + // compute centroids + tmp = j1.v; + tmp.build_etaphi(); + eta1 = tmp.eta; + phi1 = tmp.phi; + + tmp = j2.v; + tmp.build_etaphi(); + eta2 = tmp.eta; + phi2 = tmp.phi; + + jet1.v = jet2.v = Cmomentum(); + + // compute jet splitting + do{ + if (j1.contents[i1]eta,v->phi); + } else if (j1.contents[i1]>j2.contents[i2]){ + // particle i2 belong only to jet 2 + v = &(particles[j2.contents[i2]]); + jet2.contents.push_back(j2.contents[i2]); + jet2.v += *v; + jet2.pt_tilde += pt[j2.contents[i2]]; + i2++; + jet2.range.add_particle(v->eta,v->phi); + } else { // (j1.contents[i1]==j2.contents[i2]) + // common particle, decide which is the closest centre + v = &(particles[j1.contents[i1]]); + + // distance w.r.t. centroid 1 + dx1 = eta1 - v->eta; + dy1 = fabs(phi1 - v->phi); + if (dy1>M_PI) + dy1 -= twopi; + + // distance w.r.t. centroid 2 + dx2 = eta2 - v->eta; + dy2 = fabs(phi2 - v->phi); + if (dy2>M_PI) + dy2 -= twopi; + + //? what when == ? + double d1sq = dx1*dx1+dy1*dy1; + double d2sq = dx2*dx2+dy2*dy2; + // do bookkeeping on most ambiguous split + if (fabs(d1sq-d2sq) < most_ambiguous_split) + most_ambiguous_split = fabs(d1sq-d2sq); + + if (d1sqeta,v->phi); + } else { + // particle i2 belong only to jet 2 + jet2.contents.push_back(j2.contents[i2]); + jet2.v += *v; + jet2.pt_tilde += pt[j2.contents[i2]]; + jet2.range.add_particle(v->eta,v->phi); + } + + i1++; + i2++; + } + } while ((i1eta,v->phi); + } + while (i2eta,v->phi); + } + + // finalise jets + jet1.n = jet1.contents.size(); + jet2.n = jet2.contents.size(); + + //jet1.range = j1.range; + //jet2.range = j2.range; + + // remove previous jets +#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES + cand_refs.erase(j1.v.ref); + cand_refs.erase(j2.v.ref); +#endif + candidates->erase(it_j1); + candidates->erase(it_j2); + + // reinsert new ones + insert(jet1); + insert(jet2); + + return true; +} + +// merge the two given jet. +// during this procedure, the jets j1 & j2 are replaced +// by 1 single jets containing both of them. +// - it_j1 iterator of the first jet in 'candidates' +// - it_j2 iterator of the second jet in 'candidates' +// return true on success, false on error +//////////////////////////////////////////////////////////////// +bool Csplit_merge::merge(cjet_iterator &it_j1, cjet_iterator &it_j2){ + Cjet jet; + int i; + + // build new jet + // note: particles within j1 & j2 have already been stored in indices + for (i=0;irange, it_j2->range); + + // remove old candidates +#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES + if (merge_identical_protocones){ + cand_refs.erase(it_j1->v.ref); + cand_refs.erase(it_j2->v.ref); + } +#endif + candidates->erase(it_j1); + candidates->erase(it_j2); + + // reinsert new candidate + insert(jet); + + return true; +} + +/** + * Check whether or not a jet has to be inserted in the + * list of protojets. If it has, set its sm_variable and + * insert it to the list of protojets. + */ +bool Csplit_merge::insert(Cjet &jet){ + + // eventually check that no other candidate are present with the + // same cone contents. We recall that this automatic merging of + // identical protocones can lead to infrared-unsafe situations. +#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES + if ((merge_identical_protocones) && (!cand_refs.insert(jet.v.ref).second)) + return false; +#endif + + // check that the protojet has large enough pt + if (jet.v.perp2()insert(jet); + + return true; +} + +/** + * given a 4-momentum and its associated pT, return the + * variable that has to be used for SM + * \param v 4 momentum of the protojet + * \param pt_tilde pt_tilde of the protojet + */ +double Csplit_merge::get_sm_var2(Cmomentum &v, double &pt_tilde){ + switch(ptcomparison.split_merge_scale) { + case SM_pt: return v.perp2(); + case SM_mt: return v.perpmass2(); + case SM_pttilde: return pt_tilde*pt_tilde; + case SM_Et: return v.Et2(); + default: + throw Csiscone_error("Unsupported split-merge scale choice: " + + ptcomparison.SM_scale_name()); + } + + return 0.0; +} + +} Index: /trunk/SISCone/split_merge.h =================================================================== --- /trunk/SISCone/split_merge.h (revision 20) +++ /trunk/SISCone/split_merge.h (revision 20) @@ -0,0 +1,375 @@ +// -*- C++ -*- +/////////////////////////////////////////////////////////////////////////////// +// File: split_merge.h // +// Description: header file for splitting/merging (contains the CJet class) // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:30 $// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef __SPLIT_MERGE_H__ +#define __SPLIT_MERGE_H__ + +#include "defines.h" +#include "geom_2d.h" +#include "momentum.h" +#include +#include +#include +#include +#include + +namespace siscone{ + +/** + * \class Cjet + * real Jet information. + * + * This class contains information for one single jet. + * That is, first, its momentum carrying information + * about its centre and pT, and second, its particle + * contents + */ +class Cjet{ + public: + /// default ctor + Cjet(); + + /// default dtor + ~Cjet(); + + Cmomentum v; ///< jet momentum + double pt_tilde; ///< p-scheme pt + int n; ///< number of particles inside + std::vector contents; ///< particle contents (list of indices) + + /// ordering variable used for ordering and overlap in the + /// split--merge. This variable is automatically set either to + /// pt_tilde, or to mt or to pt, depending on the siscone + /// parameter. Note that the default behaviour is pt_tilde and that + /// other chices may lead to infrared unsafe situations. + /// Note: we use the square of the varible rather than the variable itself + double sm_var2; + + /// covered range in eta-phi + Ceta_phi_range range; + + /// pass at which the jet has been found + /// It starts at 0 (first pass), -1 means infinite rapidity + int pass; +}; + +/// ordering of jets in pt (e.g. used in final jets ordering) +bool jets_pt_less(const Cjet &j1, const Cjet &j2); + + +/// the choices of scale variable that can be used in the split-merge +/// step, both for ordering the protojets and for measuing their +/// overlap; pt, Et and mt=sqrt(pt^2+m^2) are all defined in E-scheme +/// (4-momentum) recombination; pttilde = \sum_{i\in jet} |p_{t,i}| +/// +/// NB: if one changes the order here, one _MUST_ also change the order +/// in the SISCone plugin +enum Esplit_merge_scale { + SM_pt, ///< transverse momentum (E-scheme), IR unsafe + SM_Et, ///< transverse energy (E-scheme), not long. boost inv. + ///< original run-II choice [may not be implemented] + SM_mt, ///< transverse mass (E-scheme), IR safe except + ///< in decays of two identical narrow heavy particles + SM_pttilde ///< pt-scheme pt = \sum_{i in jet} |p_{ti}|, should + ///< be IR safe in all cases +}; + +/// return the name of the split-merge scale choice +std::string split_merge_scale_name(Esplit_merge_scale sms); + +/** + * \class Csplit_merge_comparison + * a class that allows us to carry out comparisons of pt of jets, using + * information from exact particle contents where necessary. + * + */ +class Csplit_merge_ptcomparison{ +public: + Csplit_merge_ptcomparison() : + particles(0), split_merge_scale(SM_pttilde){}; + + // return the name corresponding to the SM scale variable + std::string SM_scale_name() const { + return split_merge_scale_name(split_merge_scale);} + + std::vector * particles; + std::vector * pt; + bool operator()(const Cjet &jet1, const Cjet &jet2) const; + + /** + * get the difference between 2 jets, calculated such that rounding + * errors will not affect the result even if the two jets have + * almost the same content (so that the difference is below the + * rounding errors) + * + * \param j1 first jet + * \param j2 second jet + * \param v jet1-jet2 + * \param pt_tilde jet1-jet2 pt_tilde + */ + void get_difference(const Cjet &j1, const Cjet &j2, Cmomentum *v, double *pt_tilde) const; + + /// the following parameter controls the variable we're using for + /// the split-merge process i.e. the variable we use for + /// 1. ordering jet candidates; + /// 2. computing te overlap fraction of two candidates. + /// The default value uses pttile (p-scheme pt). Other alternatives are + /// pt, mt=sqrt(pt^2+m^2)=sqrt(E^2-pz^2) or Et. + /// NOTE: Modifying the default choice can have nasty effects: + /// - using pt leads to some IR unsafety when we have two jets, + /// e.g. back-to-back, with the same pt. In that case, their ordering + /// in pt is random and can be affected by the addition of a + /// soft particle. Hence, we highly recommand to keep this to + /// the default value i.e. to use pt only for the purpose of + /// investigating the IR issue + /// - using Et is safe but do not respect boost invariance + /// - using mt solves the IR unsafety issues with the pt variable + /// for QCD jets but the IR unsafety remains for nack-to-back + /// jets of unstable narrow-width particles (e.g. Higgs). + /// Therefore, keeping the default value is strongly advised. + Esplit_merge_scale split_merge_scale; +}; + + +// iterator types +/// iterator definition for the jet candidates structure +typedef std::multiset::iterator cjet_iterator; + +/// iterator definition for the jet structure +typedef std::vector::iterator jet_iterator; + + + +/** + * \class Csplit_merge + * Class used to split and merge jets. + */ +class Csplit_merge{ + public: + /// default ctor + Csplit_merge(); + + /// default dtor + ~Csplit_merge(); + + + ////////////////////////////// + // initialisation functions // + ////////////////////////////// + + /** + * initialisation function + * \param _particles list of particles + * \param protocones list of protocones (initial jet candidates) + * \param R2 cone radius (squared) + * \param ptmin minimal pT allowed for jets + * \return 0 on success, 1 on error + */ + int init(std::vector &_particles, std::vector *protocones, double R2, double ptmin=0.0); + + /** + * initialisation function for particle list + * \param _particles list of particles + * \return 0 on success, 1 on error + */ + int init_particles(std::vector &_particles); + + /** + * build initial list of left particles + */ + int init_pleft(); + + + //////////////////////// + // cleaning functions // + //////////////////////// + + /// partial clearance + int partial_clear(); + + /// full clearance + int full_clear(); + + + ///////////////////////////////// + // main parts of the algorithm // + ///////////////////////////////// + + /** + * build the list 'p_uncol_hard' from p_remain by clustering + * collinear particles and removing particles softer than + * stable_cone_soft_pt2_cutoff + * note that thins in only used for stable-cone detection + * so the parent_index field is unnecessary + */ + int merge_collinear_and_remove_soft(); + + /** + * add a list of protocones + * \param protocones list of protocones (initial jet candidates) + * \param R2 cone radius (squared) + * \param ptmin minimal pT allowed for jets + * \return 0 on success, 1 on error + */ + int add_protocones(std::vector *protocones, double R2, double ptmin=0.0); + + /** + * really do the splitting and merging + * At the end, the vector jets is filled with the jets found. + * the 'contents' field of each jets contains the indices + * of the particles included in that jet. + * \param overlap_tshold threshold for splitting/merging transition + * \param ptmin minimal pT allowed for jets + * \return the number of jets is returned + */ + int perform(double overlap_tshold, double ptmin=0.0); + + + ////////////////////////////// + // save and debug functions // + ////////////////////////////// + + /// save final jets + /// \param flux stream to save the jet contentss + int save_contents(FILE *flux); + + /// show jets/candidates status + int show(); + + // particle information + int n; ///< number of particles + std::vector particles; ///< list of particles + std::vector pt; ///< list of particles' pt + int n_left; ///< numer of particles that does not belong to any jet + std::vector p_remain; ///< list of particles remaining to deal with + std::vector p_uncol_hard; ///< list of particles remaining with collinear clustering + int n_pass; ///< index of the run + + /// minimal difference in squared distance between a particle and + /// two overlapping protojets when doing a split (useful when + /// testing approx. collinear safety) + double most_ambiguous_split; + + // jets information + std::vector jets; ///< list of jets + + // working entries + int *indices; ///< maximal size array for indices works + int idx_size; ///< number of elements in indices1 + + /// The following flag indicates that identical protocones + /// are to be merged automatically each time around the split-merge + /// loop and before anything else happens. + /// + /// This flag is only effective if ALLOW_MERGE_IDENTICAL_PROTOCONES + /// is set in 'defines.h' + /// Note that this lead to infrared-unsafety so it is disabled + /// by default + bool merge_identical_protocones; + + /// member used for detailed comparisons of pt's + Csplit_merge_ptcomparison ptcomparison; + + /// stop split--merge when the SM_var of the hardest protojet + /// is below this cut-off. + /// This is not collinear-safe so you should not use this + /// variable unless you really know what you are doing + /// Note that the cut-off is set on the variable squared. + double SM_var2_hardest_cut_off; + + /// pt cutoff for the particles to put in p_uncol_hard + /// this is meant to allow removing soft particles in the + /// stable-cone search. + double stable_cone_soft_pt2_cutoff; + + private: + /** + * get the overlap between 2 jets + * \param j1 first jet + * \param j2 second jet + * \param v returned overlap^2 (determined by the choice of SM variable) + * \return true if overlapping, false if disjoint + */ + bool get_overlap(const Cjet &j1, const Cjet &j2, double *v); + + + /** + * split the two given jets. + * during this procedure, the jets j1 & j2 are replaced + * by 2 new jets. Common particles are associted to the + * closest initial jet. + * \param it_j1 iterator of the first jet in 'candidates' + * \param it_j2 iterator of the second jet in 'candidates' + * \param j1 first jet (Cjet instance) + * \param j2 second jet (Cjet instance) + * \return true on success, false on error + */ + bool split(cjet_iterator &it_j1, cjet_iterator &it_j2); + + /** + * merge the two given jet. + * during this procedure, the jets j1 & j2 are replaced + * by 1 single jets containing both of them. + * \param it_j1 iterator of the first jet in 'candidates' + * \param it_j2 iterator of the second jet in 'candidates' + * \return true on success, false on error + */ + bool merge(cjet_iterator &it_j1, cjet_iterator &it_j2); + + /** + * Check whether or not a jet has to be inserted in the + * list of protojets. If it has, set its sm_variable and + * insert it to the list of protojets. + * \param jet jet to insert + */ + bool insert(Cjet &jet); + + /** + * given a 4-momentum and its associated pT, return the + * variable tht has to be used for SM + * \param v 4 momentum of the protojet + * \param pt_tilde pt_tilde of the protojet + */ + double get_sm_var2(Cmomentum &v, double &pt_tilde); + + // jet information + /// list of jet candidates + std::auto_ptr > candidates; + + /// minimal pt2 + double pt_min2; + +#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES + /// checkxor for the candidates (to avoid having twice the same contents) + std::set cand_refs; +#endif +}; + +} + + +#endif Index: /trunk/SISCone/vicinity.cc =================================================================== --- /trunk/SISCone/vicinity.cc (revision 20) +++ /trunk/SISCone/vicinity.cc (revision 20) @@ -0,0 +1,295 @@ +/////////////////////////////////////////////////////////////////////////////// +// File: vicinity.cpp // +// Description: source file for particle vicinity (Cvicinity class) // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:30 $// +/////////////////////////////////////////////////////////////////////////////// + +#include "vicinity.h" +#include +#include +#include + +namespace siscone{ + +using namespace std; + +/************************************************************* + * Cvicinity_elm implementation * + * element in the vicinity of a parent. * + * class used to manage one points in the vicinity * + * of a parent point. * + *************************************************************/ + +// ordering pointers to Cvicinity_elm +//------------------------------------ +bool ve_less(Cvicinity_elm *ve1, Cvicinity_elm *ve2){ + return ve1->angle < ve2->angle; +} + + +/************************************************************* + * Cvicinity implementation * + * list of element in the vicinity of a parent. * + * class used to manage the points which are in the vicinity * + * of a parent point. The construction of the list can be * + * made from a list of points or from a quadtree. * + *************************************************************/ + +// default constructor +//--------------------- +Cvicinity::Cvicinity(){ + n_part = 0; + + ve_list = NULL; +#ifdef USE_QUADTREE_FOR_STABILITY_TEST + quadtree = NULL; +#endif + + parent = NULL; + VR2 = VR = 0.0; + +} + +// constructor with initialisation +//--------------------------------- +Cvicinity::Cvicinity(vector &_particle_list){ + parent = NULL; +#ifdef USE_QUADTREE_FOR_STABILITY_TEST + quadtree = NULL; +#endif + VR2 = VR = 0.0; + + set_particle_list(_particle_list); +} + +// default destructor +//-------------------- +Cvicinity::~Cvicinity(){ + if (ve_list!=NULL) + delete[] ve_list; + +#ifdef USE_QUADTREE_FOR_STABILITY_TEST + if (quadtree!=NULL) + delete quadtree; +#endif +} + +/* + * set the particle_list + * - particle_list list of particles (type Cmomentum) + * - n number of particles in the list + ************************************************************/ +void Cvicinity::set_particle_list(vector &_particle_list){ + int i,j; +#ifdef USE_QUADTREE_FOR_STABILITY_TEST + double eta_max=0.0; +#endif + + // if the particle list is not empty, destroy it ! + if (ve_list!=NULL){ + delete[] ve_list; + } + vicinity.clear(); +#ifdef USE_QUADTREE_FOR_STABILITY_TEST + if (quadtree!=NULL) + delete quadtree; +#endif + + // allocate memory array for particles + // Note: - we compute max for |eta| + // - we allocate indices to particles + n_part = 0; + plist.clear(); + pincluded.clear(); + for (i=0;i<(int) _particle_list.size();i++){ + // if a particle is colinear with the beam (infinite rapidity) + // we do not take it into account + if (fabs(_particle_list[i].pz)!=_particle_list[i].E){ + plist.push_back(_particle_list[i]); + pincluded.push_back(Cvicinity_inclusion()); // zero inclusion status + + // the parent_index is handled in the split_merge because + // of our multiple-pass procedure. + // Hence, it is not required here any longer. + // plist[n_part].parent_index = i; + plist[n_part].index = n_part; + + // make sure the reference is randomly created + plist[n_part].ref.randomize(); + +#ifdef USE_QUADTREE_FOR_STABILITY_TEST + if (fabs(plist[n_part].eta)>eta_max) eta_max=fabs(plist[n_part].eta); +#endif + + n_part++; + } + } + + // allocate quadtree and vicinity_elm list + // note: we set phi in [-pi:pi] as it is the natural range for atan2! + ve_list = new Cvicinity_elm[2*n_part]; +#ifdef USE_QUADTREE_FOR_STABILITY_TEST + eta_max+=0.1; + quadtree = new Cquadtree(0.0, 0.0, eta_max, M_PI); +#endif + + // append particle to the vicinity_elm list + j = 0; + for (i=0;iadd(&plist[i]); +#endif + ve_list[j].v = ve_list[j+1].v = &plist[i]; + ve_list[j].is_inside = ve_list[j+1].is_inside = &(pincluded[i]); + j+=2; + } + +} + + +/* + * build the vicinity list from a list of points. + * - _parent reference particle + * - _VR vicinity radius + ************************************************************/ +void Cvicinity::build(Cmomentum *_parent, double _VR){ + int i; + + // set parent and radius + parent = _parent; + VR = _VR; + VR2 = VR*VR; + R2 = 0.25*VR2; + R = 0.5*VR; + inv_R_EPS_COCIRC = 1.0 / R / EPSILON_COCIRCULAR; + inv_R_2EPS_COCIRC = 0.5 / R / EPSILON_COCIRCULAR; + + // clear vicinity + vicinity.clear(); + + // init parent variables + pcx = parent->eta; + pcy = parent->phi; + + // really browse the particle list + for (i=0;i0) ? 0.0 : 2.0; + double t=c/s; + return (s>0) ? 1-t/(1+fabs(t)) : 3-t/(1+fabs(t)); +} + + +/* + * append a particle to the 'vicinity' list after + * having computed the angular-ordering quantities + * - v vector to test + **********************************************************/ +void Cvicinity::append_to_vicinity(Cmomentum *v){ + double dx, dy, d2; + + // skip the particle itself) + if (v==parent) + return; + + int i=2*(v->index); + + // compute the distance of the i-th particle with the parent + dx = v->eta - pcx; + dy = v->phi - pcy; + + // pay attention to the periodicity in phi ! + if (dy>M_PI) + dy -= twopi; + else if (dy<-M_PI) + dy += twopi; + + d2 = dx*dx+dy*dy; + + // really check if the distance is less than VR + if (d2eta - ve_list[i+1].eta, + phi_in_range(v->phi-ve_list[i+1].phi)); + + // two sources of error are (GPS CCN29-19) epsilon/(R sin theta) + // and sqrt(2*epsilon/(R (1-cos theta))) and the way things work + // out, it is the _smaller_ of the two that is relevant [NB have + // changed definition of theta here relative to that used in + // CCN29] [NB2: write things so as to avoid zero denominators and + // to minimize the multiplications, divisions and above all sqrts + // -- that means that c & s are defined including a factor of VR2] + c = dot_product(OP,OC); + s = fabs(cross_product(OP,OC)); + double inv_err1 = s * inv_R_EPS_COCIRC; + double inv_err2_sq = (R2-c) * inv_R_2EPS_COCIRC; + ve_list[i].cocircular_range = pow2(inv_err1) > inv_err2_sq ? + 1.0/inv_err1 : + sqrt(1.0/inv_err2_sq); + ve_list[i+1].cocircular_range = ve_list[i].cocircular_range; + } +} + +} Index: /trunk/SISCone/vicinity.h =================================================================== --- /trunk/SISCone/vicinity.h (revision 20) +++ /trunk/SISCone/vicinity.h (revision 20) @@ -0,0 +1,156 @@ +// -*- C++ -*- +/////////////////////////////////////////////////////////////////////////////// +// File: vicinity.h // +// Description: header file for particle vicinity (Cvicinity class) // +// This file is part of the SISCone project. // +// For more details, see http://projects.hepforge.org/siscone // +// // +// Copyright (c) 2006 Gavin Salam and Gregory Soyez // +// // +// This program 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. // +// // +// This program 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 this program; if not, write to the Free Software // +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // +// // +// $Revision: 1.1 $// +// $Date: 2008-10-02 15:20:30 $// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef __VICINITY_H__ +#define __VICINITY_H__ + +#include +#include +#include "momentum.h" +#include "defines.h" +#include "quadtree.h" + +namespace siscone{ + + + +/** + * \class Cvicinity_inclusion + * \brief a class to keep track of inclusion status in cone and in cocircular region + * while using minimal resources + */ +class Cvicinity_inclusion { +public: + /// default ctor + Cvicinity_inclusion() : cone(false), cocirc(false) {} + + bool cone; ///< flag for particle inclusion in the cone + bool cocirc; ///< flag for particle inclusion in the border +}; + + +/** + * \class Cvicinity_elm + * \brief element in the vicinity of a parent. + * + * class used to manage one points in the vicinity + * of a parent point. + */ +class Cvicinity_elm{ + public: + /// pointer to the second borderline particle + Cmomentum *v; + + /// variable to tell if the particle is inside or outside the cone + Cvicinity_inclusion *is_inside; + + // centre variables + double eta; ///< eta coordinate of the center + double phi; ///< phi coordinate of the center + double angle; ///< angle with parent + bool side; ///< true if angle on the positive side, false otherwise + double cocircular_range; ///< amount by which the angle can be varied while + ///< maintaining this point within co-circularity margin + + /// list of elements co-circular with this one + /// NB: empty list uses less mem than vector + std::list cocircular; +}; + +/// ordering pointers to Cvicinity_elm +bool ve_less(Cvicinity_elm *ve1, Cvicinity_elm *ve2); + + +/** + * \class Cvicinity + * \brief list of element in the vicinity of a parent. + * + * class used to manage the points which are in the vicinity + * of a parent point. + */ +class Cvicinity{ + public: + /// default constructor + Cvicinity(); + + /// constructor with initialisation (see set_particle_list) + Cvicinity(std::vector &_particle_list); + + /// default destructor + ~Cvicinity(); + + /** + * set the particle_list + * \param _particle_list list of particles (type Cmomentum) + */ + void set_particle_list(std::vector &_particle_list); + + /** + * build the vicinity list from the list of points. + * \param _parent reference particle + * \param _VR vicinity radius + */ + void build(Cmomentum *_parent, double _VR); + + // cone kinematical information + Cmomentum *parent; ///< parent vector + double VR; ///< radius of the vicinity + double VR2; ///< squared radius of the vicinity + double R; ///< normal radius + double R2; ///< squared normal radius + double inv_R_EPS_COCIRC; ///< R / EPSILON_COCIRCULAR + double inv_R_2EPS_COCIRC; ///< R / (2*EPSILON_COCIRCULAR) + + // particle list information + int n_part; ///< number of particles + std::vector plist; ///< the list of particles + std::vector pincluded; ///< the inclusion state of particles + Cvicinity_elm *ve_list; ///< list of vicinity elements built from particle list (size=2*n) +#ifdef USE_QUADTREE_FOR_STABILITY_TEST + Cquadtree *quadtree; ///< quadtree used for final stability tests +#endif + + // vicinity information + std::vector vicinity; ///< list of points in parent's vicinity + unsigned int vicinity_size; ///< number of elements in vicinity + + protected: + /** + * append a particle to the 'vicinity' list after + * having tested it and computed the angular-ordering quantities + * \param v vector to test + */ + void append_to_vicinity(Cmomentum *v); + + // internal variables + double pcx; ///< parent centre (eta) + double pcy; ///< parent centre (phi) +}; + +} + +#endif