include/fastjet/GhostedAreaSpec.hh

0001 //FJSTARTHEADER
0002 // $Id$
0003 //
0004 // Copyright (c) 2005-2021, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
0005 //
0006 //----------------------------------------------------------------------
0007 // This file is part of FastJet.
0008 //
0009 //  FastJet is free software; you can redistribute it and/or modify
0010 //  it under the terms of the GNU General Public License as published by
0011 //  the Free Software Foundation; either version 2 of the License, or
0012 //  (at your option) any later version.
0013 //
0014 //  The algorithms that underlie FastJet have required considerable
0015 //  development. They are described in the original FastJet paper,
0016 //  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
0017 //  FastJet as part of work towards a scientific publication, please
0018 //  quote the version you use and include a citation to the manual and
0019 //  optionally also to hep-ph/0512210.
0020 //
0021 //  FastJet is distributed in the hope that it will be useful,
0022 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
0023 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0024 //  GNU General Public License for more details.
0025 //
0026 //  You should have received a copy of the GNU General Public License
0027 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
0028 //----------------------------------------------------------------------
0029 //FJENDHEADER
0030
0031
0032 #ifndef __FASTJET_GHOSTEDAREASPEC_HH__
0033 #define __FASTJET_GHOSTEDAREASPEC_HH__
0034
0035 #include<vector>
0036 #include<string>
0037 #include "fastjet/PseudoJet.hh"
0038 #include "fastjet/internal/BasicRandom.hh"
0039 #include "fastjet/Selector.hh"
0040 #include "fastjet/LimitedWarning.hh"
0041 #include "fastjet/internal/deprecated.hh"
0042 #include "fastjet/SharedPtr.hh"
0043
0044 //
0045 #define STATIC_GENERATOR 1
0046
0047 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0048
0049 /// namespace to hold default parameters for the active area spec
0050 namespace gas {
0051   const double def_ghost_maxrap  = 6.0;
0052   const int    def_repeat        = 1;
0053   const double def_ghost_area    = 0.01;
0054   const double def_grid_scatter  = 1.0;
0055   const double def_pt_scatter    = 0.1;
0056   const double def_mean_ghost_pt = 1e-100;
0057 }
0058
0059 //----------------------------------------------------------------------
0060 /// @ingroup area_classes
0061 /// \class GhostedAreaSpec
0062 /// Parameters to configure the computation of jet areas using ghosts
0063 ///
0064 /// Class that defines the parameters that go into the measurement
0065 /// of active jet areas.
0066 ///
0067 /// Notes about thread-safety.
0068 /// --------------------------
0069 ///
0070 /// Ghosts are generated randomly, using by default a static random
0071 /// number generator.
0072 ///
0073 /// By default, we will lock the number generator during the period
0074 /// over which we generate the required random numbers.  The procedure
0075 /// will keep track of the seeds that have been used to generate a
0076 /// particular set of ghosts and, ultimately, these seeds will be
0077 /// made available from ClusterSequenceArea via
0078 ///
0079 ///   ClusterSequenceArea::area_def().ghost_spec().get_last_seed(vector<int>);
0080 ///
0081 /// To use user-specified seeds in a thread-safe way, the end-user
0082 /// should use
0083 ///
0084 ///   ClusterSequenceArea csa(particles, jet_def,
0085 ///                           area_def.with_fixed_seed(user_defined_seed));
0086 ///
0087 /// or explicitly make a copy of the AreaDefinition before doing
0088 /// the clustering:
0089 ///
0090 ///   AreaDefinition local_area_def
0091 ///     = area_def.with_fixed_seed(user_defined_seed);
0092 ///   ClusterSequenceArea csa(particles, jet_def, area_def,local_area_def);
0093 ///
0094 /// This will use a local random generator to compute the ghosts (in
0095 /// particular, it will not affect the static global generator)
0096 ///
0097 /// Note that each clustering done with the GhostedAreaSpec obtained
0098 /// through area_def.with_seed(user_defined_seed) will use exactly the
0099 /// same set of ghosts.  Using
0100 /// area_def.with_fixed_seed(user_defined_seed), with an empty vector
0101 /// passed as argument, will return to using the common static random
0102 /// generator.
0103 class GhostedAreaSpec {
0104 public:
0105   /// default constructor
0106   GhostedAreaSpec():
0107     _ghost_maxrap (gas::def_ghost_maxrap),
0108     _ghost_rap_offset(0.0),
0109     _repeat       (gas::def_repeat),
0110     _ghost_area   (gas::def_ghost_area),
0111     _grid_scatter (gas::def_grid_scatter),
0112     _pt_scatter   (gas::def_pt_scatter),
0113     _mean_ghost_pt(gas::def_mean_ghost_pt),
0114     _fj2_placement(false),
0115     _user_random_generator(){_initialize();}
0116
0117   /// explicit constructor
0118   ///
0119   /// It takes as parameters the maximal (abs) rapidity for the ghosts
0120   /// and an optional user-specified random number generator.
0121   ///
0122   /// For the latter, ownership is transferred to the GhostedAreaSpec
0123   /// class (i.e. it is stored internally as a shared pointer)
0124   explicit GhostedAreaSpec(double ghost_maxrap_in,
0125                            BasicRandom<double> *user_random_generator):
0126     _ghost_maxrap(ghost_maxrap_in),
0127     _ghost_rap_offset(0.0),
0128     _repeat       (gas::def_repeat),
0129     _ghost_area   (gas::def_ghost_area),
0130     _grid_scatter (gas::def_grid_scatter),
0131     _pt_scatter   (gas::def_pt_scatter),
0132     _mean_ghost_pt(gas::def_mean_ghost_pt),
0133     _fj2_placement(false),
0134     _user_random_generator(user_random_generator) {_initialize();}
0135
0136   /// explicit constructor
0137   explicit GhostedAreaSpec(double ghost_maxrap_in,
0138                            int    repeat_in        = gas::def_repeat,
0139                            double ghost_area_in    = gas::def_ghost_area,
0140                            double grid_scatter_in  = gas::def_grid_scatter,
0141                            double pt_scatter_in    = gas::def_pt_scatter,
0142                            double mean_ghost_pt_in = gas::def_mean_ghost_pt,
0143                            BasicRandom<double> *user_random_generator=NULL):
0144     _ghost_maxrap(ghost_maxrap_in),
0145     _ghost_rap_offset(0.0),
0146     _repeat(repeat_in),
0147     _ghost_area(ghost_area_in),
0148     _grid_scatter(grid_scatter_in),
0149     _pt_scatter(pt_scatter_in),
0150     _mean_ghost_pt(mean_ghost_pt_in),
0151     _fj2_placement(false),
0152     _user_random_generator(user_random_generator) {_initialize();}
0153
0154   /// explicit constructor
0155   explicit GhostedAreaSpec(double ghost_minrap_in,
0156                double ghost_maxrap_in,
0157                            int    repeat_in        = gas::def_repeat,
0158                            double ghost_area_in    = gas::def_ghost_area,
0159                            double grid_scatter_in  = gas::def_grid_scatter,
0160                            double pt_scatter_in    = gas::def_pt_scatter,
0161                            double mean_ghost_pt_in = gas::def_mean_ghost_pt,
0162                            BasicRandom<double> *user_random_generator=NULL):
0163     _ghost_maxrap    (0.5*(ghost_maxrap_in - ghost_minrap_in)),
0164     _ghost_rap_offset(0.5*(ghost_maxrap_in + ghost_minrap_in)),
0165     _repeat(repeat_in),
0166     _ghost_area(ghost_area_in),
0167     _grid_scatter(grid_scatter_in),
0168     _pt_scatter(pt_scatter_in),
0169     _mean_ghost_pt(mean_ghost_pt_in),
0170     _fj2_placement(false),
0171     _user_random_generator(user_random_generator) {_initialize();}
0172
0173
0174   /// constructor based on a Selector
0175   explicit GhostedAreaSpec(const Selector & selector,
0176                            int    repeat_in        = gas::def_repeat,
0177                            double ghost_area_in    = gas::def_ghost_area,
0178                            double grid_scatter_in  = gas::def_grid_scatter,
0179                            double pt_scatter_in    = gas::def_pt_scatter,
0180                            double mean_ghost_pt_in = gas::def_mean_ghost_pt,
0181                            BasicRandom<double> *user_random_generator=NULL);
0182
0183
0184   /// does the initialization of actual ghost parameters
0185   void _initialize();
0186
0187   // for accessing values set by the user
0188   inline double ghost_rapmax () const {return _ghost_maxrap;}
0189   inline double ghost_maxrap () const {return _ghost_maxrap;}
0190   inline double ghost_etamax () const {return _ghost_maxrap;}
0191   inline double ghost_maxeta () const {return _ghost_maxrap;}
0192   inline double ghost_area   () const {return _ghost_area   ;}
0193   inline double grid_scatter () const {return _grid_scatter;}
0194   inline double pt_scatter   () const {return _pt_scatter  ;}
0195   inline double mean_ghost_pt() const {return _mean_ghost_pt  ;}
0196   inline int    repeat       () const {return _repeat      ;}
0197   inline bool   fj2_placement() const {return _fj2_placement;}
0198
0199   inline double kt_scatter   () const {return _pt_scatter  ;}
0200   inline double mean_ghost_kt() const {return _mean_ghost_pt  ;}
0201
0202   // for accessing values
0203   inline double actual_ghost_area() const {return _actual_ghost_area;}
0204   inline int    n_ghosts()          const {return _n_ghosts;}
0205
0206   // when explicitly modifying values, sometimes call the initializer
0207   inline void set_ghost_area   (double val) {_ghost_area    = val; _initialize();}
0208   inline void set_ghost_rapmax (double val) {_ghost_maxrap = val; _initialize();}
0209   inline void set_ghost_maxrap (double val) {_ghost_maxrap = val; _initialize();}
0210   inline void set_ghost_etamax (double val) {_ghost_maxrap = val; _initialize();}
0211   inline void set_ghost_maxeta (double val) {_ghost_maxrap = val; _initialize();}
0212   inline void set_grid_scatter (double val) {_grid_scatter   = val; }
0213   inline void set_pt_scatter   (double val) {_pt_scatter     = val; }
0214   inline void set_mean_ghost_pt(double val) {_mean_ghost_pt  = val; }
0215   inline void set_repeat       (int    val) {_repeat         = val; }
0216
0217   inline void set_kt_scatter   (double val) {_pt_scatter     = val; }
0218   inline void set_mean_ghost_kt(double val) {_mean_ghost_pt  = val; }
0219
0220   /// if val is true, set ghost placement as it was in FastJet 2.X. The
0221   /// main differences between FJ2 and FJ3 ghost placement are
0222   ///
0223   ///  - in FJ2 the rapidity spacing was
0224   ///    ceil((maxrap-minrap)/sqrt(area)), while in FJ3 it is
0225   ///    int((maxrap-minrap)/sqrt(area) + 0.5) [similarly for phi].
0226   ///    The FJ3 option offers more stability when trying to specify a
0227   ///    spacing that exactly fits the extent.
0228   ///
0229   /// - in FJ2, the ghosts are placed at the corners of grid cells
0230   ///   (i.e. extending up to maxrap), while in FJ3 they are placed at
0231   ///   the centres of grid cells (i.e. extending roughly up to
0232   ///   maxrap-sqrt(area)). The FJ2 behaviour effectively skews the
0233   ///   total area coverage when maxrap is small, by an amount
0234   ///   sqrt(area)/(2*maxrap).
0235   ///
0236   /// FJ2 placement is now deprecated.
0237   FASTJET_DEPRECATED_MSG("This is deprecated since we strongly recommend use of the new ghost placement instead",
0238   void set_fj2_placement(bool  val));
0239
0240   /// return nphi (ghosts layed out (-nrap, 0..nphi-1), (-nrap+1,0..nphi-1),
0241   /// ... (nrap,0..nphi-1)
0242   inline int nphi() const {return _nphi;}
0243   inline int nrap() const {return _nrap;}
0244
0245   /// get all relevant information about the status of the
0246   /// random number generator, so that it can be reset subsequently
0247   /// with set_random_status.
0248   ///
0249   ///
0250   inline void get_random_status(std::vector<int> & __iseed) const {
0251     if (_user_random_generator){
0252       _user_random_generator->get_status(__iseed);
0253     } else {
0254       _random_generator.get_status(__iseed);
0255     }
0256   }
0257
0258   /// set the status of the random number generator, as obtained
0259   /// previously with get_random_status. Note that the random
0260   /// generator is a static member of the class, i.e. common to all
0261   /// instances of the class --- so if you modify the random for this
0262   /// instance, you modify it for all instances.
0263   inline void set_random_status(const std::vector<int> & __iseed) {
0264     if (_user_random_generator){
0265       _user_random_generator->set_status(__iseed);
0266     } else {
0267       _random_generator.set_status(__iseed);
0268     }
0269   }
0270
0271   /// allows to return a copy of this GhostedAreaSpec with a local set
0272   /// of seeds
0273   GhostedAreaSpec with_fixed_seed(const std::vector<int> & __iseed) const {
0274     GhostedAreaSpec new_spec = (*this);
0275     new_spec._fixed_seed = __iseed;
0276     return new_spec;
0277   }
0278
0279   /// returns the current fixed seed
0280   void get_fixed_seed(std::vector<int> & __iseed) const {
0281     __iseed = _fixed_seed;
0282   }
0283
0284   /// allows the user to get the seed that was used at the start of the
0285   /// last generation of ghosts.
0286   ///
0287   /// This should typically be access through the area definition held
0288   /// by the ClusterSequenceArea, because the CSA class takes a copy of the
0289   /// AreaDefinition and it is that copy that stored the
0290   void get_last_seed(std::vector<int> & __iseed) const {
0291     if (_repeat > 1) _warn_fixed_last_seeds_nrepeat_gt_1
0292                       .warn("Using fixed seeds (or accessing last used seeds) not sensible with repeat>1");
0293     __iseed = _last_used_seed;
0294   }
0295
0296   inline void checkpoint_random() {get_random_status(_random_checkpoint);}
0297   inline void restore_checkpoint_random() {set_random_status(_random_checkpoint);}
0298
0299   /// for a summary
0300   std::string description() const;
0301
0302   /// push a set of ghost 4-momenta onto the back of the vector of
0303   /// PseudoJets
0304   void add_ghosts(std::vector<PseudoJet> & ) const;
0305
0306   /// very deprecated public access to a random number
0307   /// from the internal generator
0308   inline double random_at_own_risk() const {return _our_rand();}
0309   /// very deprecated public access to the generator itself
0310   inline BasicRandom<double> & generator_at_own_risk() const {
0311     return _user_random_generator ? *_user_random_generator : _random_generator;}
0312   /// access to the user-defined random-number generator. Will be empty if not set.
0313   inline SharedPtr<BasicRandom<double> > & user_random_generator_at_own_risk(){
0314     return _user_random_generator;}
0315
0316 private:
0317
0318   // quantities that determine nature and distribution of ghosts
0319   double _ghost_maxrap;
0320   double _ghost_rap_offset;
0321   int    _repeat      ;
0322   double _ghost_area   ;
0323   double _grid_scatter;
0324   double _pt_scatter  ;
0325   double _mean_ghost_pt;
0326   bool   _fj2_placement;
0327
0328   Selector _selector;
0329
0330   // derived quantities
0331   double _actual_ghost_area, _dphi, _drap;
0332   int    _n_ghosts, _nphi, _nrap;
0333
0334
0335   std::vector<int> _random_checkpoint;
0336
0337   // optional fixed seeds
0338   std::vector<int> _fixed_seed;
0339
0340   // access to the seeds used the very last time
0341   mutable std::vector<int> _last_used_seed;
0342
0343   // in order to keep thread-safety, have an independent random
0344   // generator for each thread
0345   //
0346   // 2015-09-24: thread_local is not supported by Apple's version of
0347   // clang! So we'll rely on something different here (see comment at
0348   // the top of the class)
0349   //#ifdef FASTJET_HAVE_LIMITED_THREAD_SAFETY
0350   //  static thread_local BasicRandom<double> _random_generator;
0351   //#else
0352   static BasicRandom<double> _random_generator;
0353   //#endif
0354   //mutable BasicRandom<double> _random_generator;
0355
0356   /// a set of seeds as defined by the end-user
0357   std::vector<int> _user_defined_seeds;
0358
0359   // allow for a user-defined random generator
0360   SharedPtr<BasicRandom<double> > _user_random_generator;
0361
0362   static LimitedWarning _warn_fj2_placement_deprecated;
0363   static LimitedWarning _warn_fixed_last_seeds_nrepeat_gt_1;
0364
0365   inline double _our_rand() const {
0366     return _user_random_generator ? (*_user_random_generator)() : _random_generator();}
0367
0368   inline void _our_rand(unsigned int npoints, double *pointer,
0369                         std::vector<int> & used_init_seed) const {
0370     return _user_random_generator
0371       ? (*_user_random_generator)(npoints, pointer, used_init_seed)
0372       : _random_generator(npoints, pointer, used_init_seed);
0373   }
0374
0375 };
0376
0377 /// just provide a typedef for backwards compatibility with programs
0378 /// based on versions 2.0 and 2.1 of fastjet. Since there is no
0379 /// easy way of telling people this is deprecated at compile or run
0380 /// time, we should be careful before removing this in the future.
0381 typedef GhostedAreaSpec ActiveAreaSpec;
0382
0383
0384 FASTJET_END_NAMESPACE
0385
0386 #endif // __FASTJET_GHOSTEDAREASPEC_HH__