siscone/spherical/protocones.h

0001 // -*- C++ -*-
0002 ///////////////////////////////////////////////////////////////////////////////
0003 // File: protocones.h                                                        //
0004 // Description: header file for stable cones determination (Cstable_cones)   //
0005 // This file is part of the SISCone project.                                 //
0006 // WARNING: this is not the main SISCone trunk but                           //
0007 //          an adaptation to spherical coordinates                           //
0008 // For more details, see http://projects.hepforge.org/siscone                //
0009 //                                                                           //
0010 // Copyright (c) 2006-2008 Gavin Salam and Gregory Soyez                          //
0011 //                                                                           //
0012 // This program is free software; you can redistribute it and/or modify      //
0013 // it under the terms of the GNU General Public License as published by      //
0014 // the Free Software Foundation; either version 2 of the License, or         //
0015 // (at your option) any later version.                                       //
0016 //                                                                           //
0017 // This program is distributed in the hope that it will be useful,           //
0018 // but WITHOUT ANY WARRANTY; without even the implied warranty of            //
0019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             //
0020 // GNU General Public License for more details.                              //
0021 //                                                                           //
0022 // You should have received a copy of the GNU General Public License         //
0023 // along with this program; if not, write to the Free Software               //
0024 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
0025 //                                                                           //
0026 // $Revision::                                                              $//
0027 // $Date::                                                                  $//
0028 ///////////////////////////////////////////////////////////////////////////////
0029
0030 #ifndef __SPH_PROTOCONES_H__
0031 #define __SPH_PROTOCONES_H__
0032
0033 #include "momentum.h"
0034 #include "vicinity.h"
0035 #include <stdio.h>
0036 #include <vector>
0037 #include <list>
0038 #include "hash.h"
0039
0040 #include <siscone/defines.h>
0041
0042 namespace siscone_spherical{
0043
0044 /**
0045  * \class CSphborder_store
0046  *
0047  * class for storing a border momentum (in context of co-circularity
0048  * checks).
0049
0050  * This class essentially calculates angle of border point w.r.t.
0051  * circle center (eta & phi), and provides a store of information
0052  * about whether we are currently including this point in the
0053  * candidate
0054  */
0055 class CSphborder_store{
0056 public:
0057   /// default ctor
0058   CSphborder_store(CSphmomentum * momentum, CSph3vector &centre, CSph3vector &angl_dir1, CSph3vector &angl_dir2) :
0059     mom(momentum),  is_in(false) {
0060     CSph3vector diff = (*momentum) - centre;
0061     angle = atan2(dot_product3(diff, angl_dir2), dot_product3(diff, angl_dir1));
0062 #ifdef DEBUG_STABLE_CONES
0063     std::cout << "  adding point " << momentum->_theta << ", " << momentum->_phi
0064           << " at an angle of " << angle << std::endl;
0065 #endif
0066   }
0067
0068   CSphmomentum * mom;  ///< particle momentum
0069   double angle;     ///< angle w.r.t. circle centre
0070   bool   is_in;     ///< inclusion status of the particle
0071 };
0072
0073
0074 /// allows easy sorting of CSphborder_store objects (which need to be
0075 /// ordered in angle).
0076 inline bool operator<(const CSphborder_store & a, const CSphborder_store & b) {
0077   return a.angle < b.angle;
0078 }
0079
0080
0081 /**
0082  * \class CSphstable_cones
0083  * \brief Computes the list of stable comes from a particle list.
0084  *
0085  * This class does the first fundamental task of te cone algorithm:
0086  * it is used to compute the list of stable cones given a list
0087  * of particles.
0088  */
0089 class CSphstable_cones : public CSphvicinity{
0090  public:
0091   /// default ctor
0092   CSphstable_cones();
0093
0094   /// ctor with initialisation (sse init for details)
0095   CSphstable_cones(std::vector<CSphmomentum> &_particle_list);
0096
0097   /// default dtor
0098   ~CSphstable_cones();
0099
0100   /**
0101    * initialisation
0102    * \param _particle_list  list of particles
0103    */
0104   void init(std::vector<CSphmomentum> &_particle_list);
0105
0106   /**
0107    * compute stable cones.
0108    * This function really does the job i.e. computes
0109    * the list of stable cones (in a seedless way)
0110    * \param _radius   radius of the cones
0111    * \return The number of stable cones found is returned
0112    */
0113   int get_stable_cones(double _radius);
0114
0115   /// list of stable cones
0116   std::vector<CSphmomentum> protocones;
0117
0118   /// list of candidates
0119   sph_hash_cones *hc;
0120
0121   /// total number of tested cones
0122   int nb_tot;
0123 #ifdef DEBUG_STABLE_CONES
0124   int nb_hash_cones, nb_hash_occupied;
0125 #endif
0126
0127  protected:
0128   /// cone radius
0129   double R;
0130
0131   /// cone radius SQUARED
0132   double R2;
0133
0134   /// squared tangent of the cone radius
0135   double tan2R;
0136
0137  private:
0138   /// cone with a given particle as parent
0139   /// this reduction to a single vector assumes we trust the checksums
0140   CSphmomentum cone;
0141
0142   /// child particle, taken in the 'vicinity' list
0143   CSphmomentum *child;
0144
0145   /// centre of the tested cone
0146   CSphvicinity_elm *centre;
0147
0148   /// index in the particle list;
0149   unsigned int centre_idx;
0150
0151   /// first cone used in the vicinity list
0152   unsigned int first_cone;
0153
0154   /**
0155    * initialise the cone.
0156    * We take the first particle in the angular ordering to compute this one
0157    * \return 0 on success, 1 on error
0158    */
0159   int init_cone();
0160
0161   /**
0162    * test cones.
0163    * We check if the cone(s) build with the present parent and child
0164    * are stable
0165    * \return 0 on success 1 on error
0166    */
0167   int test_cone();
0168
0169   /**
0170    * update the cone
0171    * go to the next child for that parent and update 'cone' appropriately
0172    * \return 0 if update candidate found, 1 otherwise
0173    */
0174   int update_cone();
0175
0176   /*
0177    * run through the vicinity of the current parent and for each child
0178    * indicate which members are cocircular...
0179    */
0180   void prepare_cocircular_lists();
0181
0182   /**
0183    * check if we are in a situation of cocircularity.
0184    * if it is the case, update and test in the corresponding way
0185    * \return 'false' if no cocircularity detected, 'true' otherwise
0186    * Note that if cocircularity is detected, we need to
0187    * recall 'update' from 'update' !!!
0188    */
0189   bool cocircular_check();
0190
0191   /**
0192    * Routine for testing cocircular configurations in p^3 time,
0193    * rather than 2^p time;
0194    */
0195   void test_cone_cocircular(CSphmomentum & borderless_cone,
0196                 std::list<CSphmomentum *> & border_list);
0197
0198   /**
0199    * carry out the computations needed for the stability check of the
0200    * candidate, using the border_vect to indicate which particles
0201    * should / should not be in the stable cone; if the cone is stable
0202    * insert it into the hash.
0203    */
0204   void test_stability(CSphmomentum & candidate,
0205                       const std::vector<CSphborder_store> & border_vect);
0206
0207   /**
0208    * compute the cone contents by going once around the full set of
0209    * circles and tracking the entry/exit status each time -- this sets
0210    * up the inclusion information, which can then be directly used to
0211    * calculate the cone momentum.
0212    */
0213   void compute_cone_contents();
0214
0215   /**
0216    * compute the cone momentum from particle list.
0217    * in this version, we use the 'pincluded' information
0218    * from the CSphvicinity class
0219    */
0220   void recompute_cone_contents();
0221
0222   /*
0223    * if we have gone beyond the acceptable threshold of change, compute
0224    * the cone momentum from particle list.  in this version, we use the
0225    * 'pincluded' information from the CSphvicinity class, but we don't
0226    * change the member cone, only the locally supplied one
0227    */
0228   void recompute_cone_contents_if_needed(CSphmomentum & this_cone, double & this_dpt);
0229
0230   /**
0231    * compute stability of all enumerated candidates.
0232    * For all candidate cones which are stable w.r.t. their border particles,
0233    * pass the last test: stability with quadtree intersection
0234    */
0235   int proceed_with_stability();
0236
0237   /*
0238    * circle intersection.
0239    * computes the intersection with a circle of given centre and radius.
0240    * The output takes the form of a checkxor of the intersection's particles
0241    *  - cx    circle centre x coordinate
0242    *  - cy    circle centre y coordinate
0243    * return the checkxor for the intersection
0244    ******************************************************************/
0245   siscone::Creference circle_intersect(CSph3vector &cone_centre);
0246
0247   /// present candidate cone
0248   CSphmomentum cone_candidate;
0249
0250   /// in case of co-circular points, vector for them
0251   std::vector<CSphmomentum*> child_list;
0252
0253   /// list of cocircular enclusures already studied
0254   /// first element if cone contents, second is cone border
0255   std::vector< std::pair<siscone::Creference,siscone::Creference> > multiple_centre_done;
0256
0257   // information for updating cone contents to avoid rounding errors
0258   double dpt;          ///< sums of Delta P_t
0259 };
0260
0261 }
0262 #endif