Fork me on GitHub

Changeset 273e668 in git for external/fastjet/plugins


Ignore:
Timestamp:
Oct 15, 2014, 10:55:55 AM (10 years ago)
Author:
Pavel Demin <pavel.demin@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
35b9204, b25d4cf
Parents:
f14bd6a
git-author:
Pavel Demin <pavel.demin@…> (10/10/14 08:56:40)
git-committer:
Pavel Demin <pavel.demin@…> (10/15/14 10:55:55)
Message:

upgrade FastJet to version 3.1.0

Location:
external/fastjet/plugins
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • external/fastjet/plugins/D0RunICone/fastjet/D0RunIConePlugin.hh

    rf14bd6a r273e668  
    3232//FJENDHEADER
    3333
     34#include "fastjet/internal/base.hh" // namespace macros (include explicitly to help Doxygen)
    3435#include "D0RunIBaseConePlugin.hh"
    3536
  • external/fastjet/plugins/D0RunICone/fastjet/D0RunIpre96ConePlugin.hh

    rf14bd6a r273e668  
    3232//FJENDHEADER
    3333
     34#include "fastjet/internal/base.hh" // namespace macros (include explicitly to help Doxygen)
    3435#include "D0RunIBaseConePlugin.hh"
    3536
  • external/fastjet/plugins/SISCone/SISConePlugin.cc

    rf14bd6a r273e668  
    2121                      four_vector.E);
    2222}
     23
     24//======================================================================
     25// wrap-up around siscone's user-defined scales
     26namespace siscone_plugin_internal{
     27  /// @ingroup internal
     28  /// \class SISConeUserScale
     29  /// class that makes the transition between the internal SISCone
     30  /// user-defined scale choice (using SISCone's Cjet) and
     31  /// user-defined scale choices in the plugn above (using FastJet's
     32  /// PseudoJets)
     33  class SISConeUserScale  : public siscone::Csplit_merge::Cuser_scale_base{
     34  public:
     35    /// ctor takes the "fastjet-style" user-defined scale as well as a
     36    /// reference to the current cluster sequence (to access the
     37    /// particles if needed)
     38    SISConeUserScale(const SISConePlugin::UserScaleBase *user_scale,
     39                     const ClusterSequence &cs)
     40      : _user_scale(user_scale), _cs(cs){}
     41
     42    /// returns the scale associated to a given jet
     43    virtual double operator()(const siscone::Cjet &jet) const{
     44      return _user_scale->result(_build_PJ_from_Cjet(jet));
     45    }
     46
     47    /// returns true id the scasle associated to jet a is larger than
     48    /// the scale associated to jet b
     49    virtual bool is_larger(const siscone::Cjet &a, const siscone::Cjet &b) const{
     50      return _user_scale->is_larger(_build_PJ_from_Cjet(a), _build_PJ_from_Cjet(b));
     51    }
     52
     53  private:
     54    /// constructs a PseudoJet from a siscone::Cjet
     55    ///
     56    /// Note that it is tempting to overload the PseudoJet ctor. This
     57    /// would not work because down the line we need to access the
     58    /// original PseudoJet through the ClusterSequence and therefore
     59    /// the PseudoJet structure needs to be aware of the
     60    /// ClusterSequence.
     61    PseudoJet _build_PJ_from_Cjet(const siscone::Cjet &jet) const{
     62      PseudoJet j(jet.v.px, jet.v.py, jet.v.pz, jet.v.E);
     63      j.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(
     64                                   new SISConePlugin::UserScaleBaseStructureType<siscone::Cjet>(jet,_cs)));     
     65      return j;
     66    }
     67
     68    const SISConePlugin::UserScaleBase *_user_scale;
     69    const ClusterSequence & _cs;
     70  };
     71}
     72
     73// end of the internal material
     74//======================================================================
    2375
    2476
     
    4597  desc << "SISCone jet algorithm with " ;
    4698  desc << "cone_radius = "       << cone_radius        () << ", ";
    47   desc << "overlap_threshold = " << overlap_threshold  () << ", ";
     99  if (_progressive_removal)
     100    desc << "progressive-removal mode, ";
     101  else 
     102    desc << "overlap_threshold = " << overlap_threshold  () << ", ";
    48103  desc << "n_pass_max = "        << n_pass_max         () << ", ";
    49104  desc << "protojet_ptmin = "    << protojet_ptmin()      << ", ";
    50   desc <<  sm_scale_string                                << ", ";
    51   desc << "caching turned "      << (caching() ? on : off);
    52   desc << ", SM stop scale = "     << _split_merge_stopping_scale;
     105  if (_progressive_removal && _user_scale) {
     106    desc << "using a user-defined scale for ordering of stable cones";
     107    string user_scale_desc = _user_scale->description();
     108    if (user_scale_desc != "") desc << " (" << user_scale_desc << ")";
     109  } else {
     110    desc <<  sm_scale_string;
     111  }
     112  if (!_progressive_removal){
     113    desc << ", caching turned "      << (caching() ? on : off);
     114    desc << ", SM stop scale = "     << _split_merge_stopping_scale;
     115  }
    53116
    54117  // add a note to the description if we use the pt-weighted splitting
     
    85148  bool new_siscone = true; // by default we'll be running it
    86149
    87   if (caching()) {
     150  if (caching() && !_progressive_removal) {
    88151
    89152    // Establish if we have a cached run with the same R, npass and
     
    138201    // run the jet finding
    139202    //cout << "plg sms: " << split_merge_scale() << endl;
    140     siscone->compute_jets(siscone_momenta, cone_radius(), overlap_threshold(),
    141                           n_pass_max(), protojet_or_ghost_ptmin(),
    142                           Esplit_merge_scale(split_merge_scale()));
     203    if (_progressive_removal){
     204      // handle the optional user-defined scale choice
     205      SharedPtr<siscone_plugin_internal::SISConeUserScale> internal_scale;
     206      if (_user_scale){
     207        internal_scale.reset(new siscone_plugin_internal::SISConeUserScale(_user_scale, clust_seq));
     208        siscone->set_user_scale(internal_scale.get());
     209      }
     210      siscone->compute_jets_progressive_removal(siscone_momenta, cone_radius(),
     211                                                n_pass_max(), protojet_or_ghost_ptmin(),
     212                                                Esplit_merge_scale(split_merge_scale()));
     213    } else {
     214      siscone->compute_jets(siscone_momenta, cone_radius(), overlap_threshold(),
     215                            n_pass_max(), protojet_or_ghost_ptmin(),
     216                            Esplit_merge_scale(split_merge_scale()));
     217    }
    143218  } else {
    144219    // rerun the jet finding
     
    155230  SISConeExtras * extras = new SISConeExtras(n);
    156231
     232  // the ordering in which the inclusive jets are transfered here is
     233  // deliberate and ensures that when a user asks for
     234  // inclusive_jets(), they are provided in the order in which SISCone
     235  // created them.
    157236  for (int ijet = njet-1; ijet >= 0; ijet--) {
    158237    const Cjet & jet = siscone->jets[ijet]; // shorthand
     
    220299}
    221300
     301// //======================================================================
     302// // material to handle user-defined scales
     303//
     304// //--------------------------------------------------
     305// // SISCone structure type
     306//
     307// // the textual descripotion
     308// std::string SISConePlugin::UserScaleBase::StructureType::description() const{
     309//   return "PseudoJet wrapping a siscone::Cjet stable cone";
     310// }
     311//
     312// // retrieve the constituents
     313// //
     314// // if you simply need to iterate over the constituents, it will be
     315// // faster to access them via constituent(i)
     316// vector<PseudoJet> SISConePlugin::UserScaleBase::StructureType::constituents(const PseudoJet &) const{
     317//   vector<PseudoJet> constits;
     318//   constits.reserve(_jet.n);
     319//   for (unsigned int i=0; i<(unsigned int)_jet.n;i++)
     320//     constits.push_back(constituent(i));
     321//   return constits;
     322// }
     323//
     324// // returns the number of constituents
     325// unsigned int SISConePlugin::UserScaleBase::StructureType::size() const{
     326//   return _jet.n;
     327// }
     328//
     329// // returns the index (in the original particle list) of the ith
     330// // constituent
     331// int SISConePlugin::UserScaleBase::StructureType::constituent_index(unsigned int i) const{
     332//   return _jet.contents[i];
     333// }
     334//
     335// // returns the ith constituent (as a PseusoJet)
     336// const PseudoJet & SISConePlugin::UserScaleBase::StructureType::constituent(unsigned int i) const{
     337//   return _cs.jets()[_jet.contents[i]];
     338// }
     339//
     340// // returns the scalar pt of this stable cone
     341// double SISConePlugin::UserScaleBase::StructureType::pt_tilde() const{
     342//   return _jet.pt_tilde;
     343// }
     344//
     345// // returns the sm_var2 (signed ordering variable squared) for this stable cone
     346// double SISConePlugin::UserScaleBase::StructureType::ordering_var2() const{
     347//   return _jet.sm_var2;
     348// }
     349
     350
    222351FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
  • external/fastjet/plugins/SISCone/config.h

    rf14bd6a r273e668  
    4949
    5050/* Define to the full name and version of this package. */
    51 #define PACKAGE_STRING "SISCone 2.0.6"
     51#define PACKAGE_STRING "SISCone 3.0.0"
    5252
    5353/* Define to the one symbol short name of this package. */
    5454#define PACKAGE_TARNAME "siscone"
    5555
     56/* Define to the home page for this package. */
     57/* #undef PACKAGE_URL */
     58
    5659/* Define to the version of this package. */
    57 #define PACKAGE_VERSION "2.0.6"
     60#define PACKAGE_VERSION "3.0.0"
    5861
    5962/* Define to 1 if you have the ANSI C header files. */
     
    6164
    6265/* Version number of package */
    63 #define VERSION "2.0.6"
     66#define VERSION "3.0.0"
  • external/fastjet/plugins/SISCone/fastjet/SISConeBasePlugin.hh

    rf14bd6a r273e668  
    4343  SISConeBasePlugin (){
    4444    _use_jet_def_recombiner = false;
     45    set_progressive_removal(false);
    4546  }
    4647
     
    4950    *this = plugin;
    5051  }
     52
     53  /// set whether to use SISCone with progressive removal instead of
     54  /// the default split_merge step.
     55  ///
     56  /// If progressive removal is enabled, the following SISCone
     57  /// variables are not used:
     58  ///
     59  /// - overlap_threshold
     60  /// - caching
     61  /// - split_merge_stopping_scale
     62  ///
     63  /// The split_merge_scale choice is reinterpreted as the ordering
     64  /// variable for progressive removal. It is also possible for the
     65  /// user to supply his/her own function for the scale that orders
     66  /// progressive removal, with set_user_scale(...)
     67  void set_progressive_removal(bool progressive_removal_in=true){
     68    _progressive_removal = progressive_removal_in;
     69  }
     70
     71  /// returns true if progressive_removal is enabled
     72  bool progressive_removal() const{ return _progressive_removal;}
    5173
    5274  /// the cone radius
     
    105127  }
    106128
     129  // user-defined scale for progressive removal
     130  //------------------------------------------------------------
     131
     132  /// \class UserScaleBase
     133  /// base class for user-defined ordering of stable cones (used for
     134  /// prorgessive removal)
     135  ///
     136  /// derived classes have to implement the () operator that returns
     137  /// the scale associated with a given jet.
     138  ///
     139  /// It is also highly recommended to implement the is_larger()
     140  /// method whenever possible, in order to avoid rounding issues
     141  /// known to lead to possible infrared unsafeties.
     142  ///
     143  /// The jets that are passed to this class will carry the structure
     144  /// of type SISConePlugin::StructureType which allows to retreive
     145  /// easily the following information:
     146  ///
     147  ///   vector<PseudoJet> constituents = jet.constituents();
     148  ///   unsigned int n_constituents = jet.structure_of<SISConePlugin::UserScaleBase>().size();
     149  ///   int index = jet.structure_of<SISConePlugin::UserScaleBase>().constituent_index(index i);
     150  ///   const PseudoJet & p = jet.structure_of<SISConePlugin::UserScaleBase>().constituent(index i);
     151  ///   double scalar_pt = jet.structure_of<SISConePlugin::UserScaleBase>().pt_tilde();
     152  ///
     153  /// see SISConePlugin::StructureType below for further details
     154  class UserScaleBase : public FunctionOfPseudoJet<double>{
     155  public:
     156    /// empty virtual dtor
     157    virtual ~UserScaleBase(){}
     158
     159    /// returns the scale associated with a given jet
     160    ///
     161    /// "progressive removal" iteratively removes the stable cone with
     162    /// the largest scale
     163    virtual double result(const PseudoJet & jet) const = 0;
     164
     165    /// returns true when the scale associated with jet a is larger than
     166    /// the scale associated with jet b
     167    ///
     168    /// By default this does a simple direct comparison but it can be
     169    /// overloaded for higher precision [recommended if possible]
     170    virtual bool is_larger(const PseudoJet & a, const PseudoJet & b) const;
     171
     172    class StructureType; // defined below
     173  };
     174
     175  // template class derived from UserScaleBase::StryctureType that
     176  // works for both SISCone jet classes
     177  // implemented below
     178  template<class Tjet>
     179  class UserScaleBaseStructureType;
     180
     181  /// set a user-defined scale for stable-cone ordering in
     182  /// progressive removal
     183  void set_user_scale(const UserScaleBase *user_scale_in){ _user_scale = user_scale_in;}
     184
     185  /// returns the user-defined scale in use (0 if none)
     186  const UserScaleBase * user_scale() const{ return _user_scale;}
     187
     188
    107189  // the things that one MUST overload required by base class
    108190  //---------------------------------------------------------
     
    120202  double _split_merge_stopping_scale;
    121203  bool   _use_jet_def_recombiner;
     204  bool   _progressive_removal;
    122205
    123206  mutable double _ghost_sep_scale;
     
    126209  /// call the re-clustering itself
    127210  virtual void reset_stored_plugin() const =0;
     211
     212  const UserScaleBase * _user_scale;
    128213
    129214};
     
    185270inline SISConeBaseExtras::~SISConeBaseExtras(){}
    186271
     272//----------------------------------------------------------------------
     273// implementation of the structure type associated with the UserScaleBase class
     274
     275/// \class SISConeBasePlugin::UserScaleBase::StructureType
     276/// the structure that allows to store the information contained
     277/// into a siscone::Cjet (built internally in SISCone from a stable
     278/// cone) into a PseudoJet
     279class SISConeBasePlugin::UserScaleBase::StructureType : public PseudoJetStructureBase {
     280public:
     281  /// base ctor (constructed from a ClusterSequence tin order to have
     282  /// access to the initial particles
     283  StructureType(const ClusterSequence &cs)
     284    : _cs(cs){}
     285
     286  /// empty virtual dtor
     287  virtual ~StructureType(){}
     288 
     289  //--------------------------------------------------
     290  // members inherited from the base class
     291  /// the textual descripotion
     292  virtual std::string description() const{
     293    return "PseudoJet wrapping a siscone jet from a stable cone";
     294  }
     295
     296  /// this structure has constituents
     297  virtual bool has_constituents() const {return true;}
     298
     299  /// retrieve the constituents
     300  ///
     301  /// if you simply need to iterate over the constituents, it will be
     302  /// faster to access them via constituent(i)
     303  virtual std::vector<PseudoJet> constituents(const PseudoJet & /*reference*/) const{
     304    std::vector<PseudoJet> constits;
     305    constits.reserve(size());
     306    for (unsigned int i=0; i<size();i++)
     307      constits.push_back(constituent(i));
     308    return constits;
     309  }
     310 
     311  //--------------------------------------------------
     312  // additional information relevant for this structure
     313
     314  /// returns the number of constituents
     315  virtual unsigned int size() const = 0;
     316
     317  /// returns the index (in the original particle list) of the ith
     318  /// constituent
     319  virtual int constituent_index(unsigned int i) const = 0;
     320
     321  /// returns the ith constituent (as a PseusoJet)
     322  const PseudoJet & constituent(unsigned int i) const{
     323    return _cs.jets()[constituent_index(i)];
     324  }
     325
     326  // /// returns the scalar pt of this stable cone
     327  // virtual double pt_tilde() const = 0;
     328
     329  /// returns the sm_var2 (signed ordering variable squared) for this stable cone
     330  virtual double ordering_var2() const = 0;
     331
     332protected:
     333  const ClusterSequence &_cs; ///< a reference to the CS (for access to the particles)
     334};
     335
     336
     337///@ingroup internal
     338/// template class derived from UserScaleBase::StryctureType that
     339/// works for both SISCone jet classes
     340/// implemented below
     341template<class Tjet>
     342class SISConeBasePlugin::UserScaleBaseStructureType : public UserScaleBase::StructureType{
     343public:
     344  UserScaleBaseStructureType(const Tjet &jet, const ClusterSequence &cs)
     345    : UserScaleBase::StructureType(cs), _jet(jet){}
     346
     347  /// empty virtual dtor
     348  virtual ~UserScaleBaseStructureType(){}
     349
     350  //--------------------------------------------------
     351  // additional information relevant for this structure
     352
     353  /// returns the number of constituents
     354  virtual unsigned int size() const{
     355    return _jet.n;
     356  }
     357
     358  /// returns the index (in the original particle list) of the ith
     359  /// constituent
     360  virtual int constituent_index(unsigned int i) const{
     361    return _jet.contents[i];
     362  }
     363
     364  // /// returns the scalar pt of this stable cone
     365  // virtual double pt_tilde() const{
     366  //   return _jet.pt_tilde;
     367  // }
     368
     369  /// returns the sm_var2 (signed ordering variable squared) for this stable cone
     370  virtual double ordering_var2() const{
     371    return _jet.sm_var2;
     372  }
     373
     374protected:
     375  const Tjet &_jet; ///< a reference to the internal jet in SISCone
     376};
     377
    187378
    188379FASTJET_END_NAMESPACE        // defined in fastjet/internal/base.hh
  • external/fastjet/plugins/SISCone/fastjet/SISConePlugin.hh

    rf14bd6a r273e668  
    77namespace siscone{
    88  class Csiscone;
     9  class Cjet;
    910}
    1011
     
    7475  /// enum for the different split-merge scale choices;
    7576  /// Note that order _must_ be the same as in siscone
    76   enum SplitMergeScale {SM_pt,     ///< transverse momentum (E-scheme), IR unsafe
    77                         SM_Et,     ///< transverse energy (E-scheme), not long. boost invariant
    78                                    ///< original run-II choice [may not be implemented]
    79                         SM_mt,     ///< transverse mass (E-scheme), IR safe except
    80                                    ///< in decays of two identical narrow heavy particles
    81                         SM_pttilde ///< pt-scheme pt = \sum_{i in jet} |p_{ti}|, should
    82                                    ///< be IR safe in all cases
     77  enum SplitMergeScale {
     78    SM_pt,     ///< transverse momentum (E-scheme), IR unsafe
     79    SM_Et,     ///< transverse energy (E-scheme), not long. boost invariant
     80               ///< original run-II choice [may not be implemented]
     81    SM_mt,     ///< transverse mass (E-scheme), IR safe except
     82               ///< in decays of two identical narrow heavy particles
     83    SM_pttilde ///< pt-scheme pt = \sum_{i in jet} |p_{ti}|, should
     84               ///< be IR safe in all cases
    8385  };
    8486
     
    108110    _split_merge_stopping_scale = split_merge_stopping_scale_in;
    109111    _ghost_sep_scale       = 0.0;
    110     _use_pt_weighted_splitting = false;}
     112    _use_pt_weighted_splitting = false;
     113    _user_scale = 0;}
    111114
    112115
     
    125128    _split_merge_stopping_scale = 0.0;
    126129    _split_merge_scale     = split_merge_on_transverse_mass_in ? SM_mt : SM_pttilde;
    127     _ghost_sep_scale       = 0.0;}
     130    _ghost_sep_scale       = 0.0;
     131    _user_scale = 0;}
    128132 
    129133  /// backwards compatible constructor for the SISCone Plugin class
     
    141145    _split_merge_stopping_scale = 0.0;
    142146    _ghost_sep_scale       = 0.0;
    143     _use_pt_weighted_splitting = false;}
     147    _use_pt_weighted_splitting = false;
     148    _user_scale = 0;}
    144149
    145150  /// minimum pt for a protojet to be considered in the split-merge step
     
    191196};
    192197
     198
     199/////\class SISConePlugin::UserScaleBase::StructureType
     200///// the structure that allows to store the information contained
     201///// into a siscone::Cjet (built internally in SISCone from a stable
     202///// cone) into a PseudoJet
     203//class SISConePlugin::UserScaleBase::StructureType : public PseudoJetStructureBase {
     204//public:
     205//  StructureType(const siscone::Cjet & jet, const ClusterSequence &cs)
     206//    : _jet(jet), _cs(cs){}
     207//
     208//  //--------------------------------------------------
     209//  // members inherited from the base class
     210//  /// the textual descripotion
     211//  virtual std::string description() const;
     212//
     213//  /// this structure has constituents
     214//  virtual bool has_constituents() const {return true;}
     215//
     216//  /// retrieve the constituents
     217//  ///
     218//  /// if you simply need to iterate over the constituents, it will be
     219//  /// faster to access them via constituent(i)
     220//  virtual std::vector<PseudoJet> constituents(const PseudoJet & /*reference*/) const;
     221//
     222//  //--------------------------------------------------
     223//  // additional information relevant for this structure
     224//
     225//  /// returns the number of constituents
     226//  unsigned int size() const;
     227//
     228//  /// returns the index (in the original particle list) of the ith
     229//  /// constituent
     230//  int constituent_index(unsigned int i) const;
     231//
     232//  /// returns the ith constituent (as a PseusoJet)
     233//  const PseudoJet & constituent(unsigned int i) const;
     234//
     235//  /// returns the scalar pt of this stable cone
     236//  double pt_tilde() const;
     237//
     238//  /// returns the sm_var2 (signed ordering variable squared) for this stable cone
     239//  double ordering_var2() const;
     240//
     241//protected:
     242//  const siscone::Cjet &_jet;  ///< a dreference to the internal jet in SISCone
     243//  const ClusterSequence &_cs; ///< a reference to the CS (for access to the particles)
     244//};
    193245
    194246//======================================================================
  • external/fastjet/plugins/SISCone/siscone.cc

    rf14bd6a r273e668  
    2121// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
    2222//                                                                           //
    23 // $Revision:: 341                                                          $//
    24 // $Date:: 2013-04-08 12:21:06 +0200 (Mon, 08 Apr 2013)                     $//
     23// $Revision:: 371                                                          $//
     24// $Date:: 2014-09-09 10:05:32 +0200 (Tue, 09 Sep 2014)                     $//
    2525///////////////////////////////////////////////////////////////////////////////
    2626
     
    2929//#else
    3030//#define PACKAGE_NAME "SISCone"
    31 //#define VERSION "2.0.6"
     31//#define VERSION "3.0.0"
    3232//#warning "No config.h file available, using preset values"
    3333//#endif
     
    8787                           int _n_pass_max, double _ptmin,
    8888                           Esplit_merge_scale _split_merge_scale){
    89   // initialise random number generator
    90   if (!init_done){
    91     // initialise random number generator
    92     ranlux_init();
    93 
    94     // do not do this again
    95     init_done=true;
    96 
    97     // print the banner
    98     if (_banner_ostr != 0){
    99       (*_banner_ostr) << "#ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << endl;
    100       (*_banner_ostr) << "#                    SISCone   version " << setw(28) << left << siscone_version() << "o" << endl;
    101       (*_banner_ostr) << "#              http://projects.hepforge.org/siscone                o" << endl;
    102       (*_banner_ostr) << "#                                                                  o" << endl;
    103       (*_banner_ostr) << "# This is SISCone: the Seedless Infrared Safe Cone Jet Algorithm   o" << endl;
    104       (*_banner_ostr) << "# SISCone was written by Gavin Salam and Gregory Soyez             o" << endl;
    105       (*_banner_ostr) << "# It is released under the terms of the GNU General Public License o" << endl;
    106       (*_banner_ostr) << "#                                                                  o" << endl;
    107       (*_banner_ostr) << "# A description of the algorithm is available in the publication   o" << endl;
    108       (*_banner_ostr) << "# JHEP 05 (2007) 086 [arXiv:0704.0292 (hep-ph)].                   o" << endl;
    109       (*_banner_ostr) << "# Please cite it if you use SISCone.                               o" << endl;
    110       (*_banner_ostr) << "#ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << endl;
    111       (*_banner_ostr) << endl;
    112 
    113       _banner_ostr->flush();
    114     }
    115   }
     89  _initialise_if_needed();
    11690
    11791  // run some general safety tests (NB: f will be checked in split-merge)
     
    174148}
    175149
     150
     151/*
     152 * compute the jets from a given particle set doing multiple passes
     153 * such pass N looks for jets among all particles not put into jets
     154 * during previous passes.
     155 *  - _particles   list of particles
     156 *  - _radius      cone radius
     157 *  - _n_pass_max  maximum number of runs
     158 *  - _ptmin       minimum pT of the protojets
     159 *  - _ordering_scale    the ordering scale to decide which stable
     160 *                       cone is removed
     161 * return the number of jets found.
     162 **********************************************************************/
     163int Csiscone::compute_jets_progressive_removal(vector<Cmomentum> &_particles, double _radius,
     164                                               int _n_pass_max, double _ptmin,
     165                                               Esplit_merge_scale _ordering_scale){
     166  _initialise_if_needed();
     167
     168  // run some general safety tests (NB: f will be checked in split-merge)
     169  if (_radius <= 0.0 || _radius >= 0.5*M_PI) {
     170    ostringstream message;
     171    message << "Illegal value for cone radius, R = " << _radius
     172            << " (legal values are 0<R<pi/2)";
     173    throw Csiscone_error(message.str());
     174  }
     175
     176  ptcomparison.split_merge_scale = _ordering_scale;
     177  partial_clear(); // make sure some things are initialised properly
     178
     179  // init the split_merge algorithm with the initial list of particles
     180  // this initialises particle list p_left of remaining particles to deal with
     181  //
     182  // this stores the "processed" particles in p_uncol_hard
     183  init_particles(_particles);
     184  jets.clear();
     185
     186  bool unclustered_left;
     187  rerun_allowed = false;
     188  protocones_list.clear();
     189
     190  do{
     191    //cout << n_left << " particle left" << endl;
     192
     193    // initialise stable_cone finder
     194    // here we use the list of remaining particles
     195    // AFTER COLLINEAR CLUSTERING !!!!!!
     196    Cstable_cones::init(p_uncol_hard);
     197
     198    // get stable cones (stored in 'protocones')
     199    unclustered_left = get_stable_cones(_radius);
     200
     201    // add the hardest stable cone to the list of jets
     202    if (add_hardest_protocone_to_jets(&protocones, R2, _ptmin)) break;
     203 
     204    _n_pass_max--;
     205  } while ((unclustered_left) && (n_left>0) && (_n_pass_max!=0));
     206
     207  // split & merge
     208  return jets.size();
     209}
     210
     211
    176212/*
    177213 * recompute the jets with a different overlap parameter.
     
    205241
    206242
     243// ensure things are initialised
     244void Csiscone::_initialise_if_needed(){
     245  // initialise random number generator
     246  if (init_done) return;
     247
     248  // initialise random number generator
     249  ranlux_init();
     250
     251  // do not do this again
     252  init_done=true;
     253
     254  // print the banner
     255  if (_banner_ostr != 0){
     256    (*_banner_ostr) << "#ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << endl;
     257    (*_banner_ostr) << "#                    SISCone   version " << setw(28) << left << siscone_version() << "o" << endl;
     258    (*_banner_ostr) << "#              http://projects.hepforge.org/siscone                o" << endl;
     259    (*_banner_ostr) << "#                                                                  o" << endl;
     260    (*_banner_ostr) << "# This is SISCone: the Seedless Infrared Safe Cone Jet Algorithm   o" << endl;
     261    (*_banner_ostr) << "# SISCone was written by Gavin Salam and Gregory Soyez             o" << endl;
     262    (*_banner_ostr) << "# It is released under the terms of the GNU General Public License o" << endl;
     263    (*_banner_ostr) << "#                                                                  o" << endl;
     264    (*_banner_ostr) << "# A description of the algorithm is available in the publication   o" << endl;
     265    (*_banner_ostr) << "# JHEP 05 (2007) 086 [arXiv:0704.0292 (hep-ph)].                   o" << endl;
     266    (*_banner_ostr) << "# Please cite it if you use SISCone.                               o" << endl;
     267    (*_banner_ostr) << "#ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << endl;
     268    (*_banner_ostr) << endl;
     269
     270    _banner_ostr->flush();
     271  }
     272}
    207273
    208274// finally, a bunch of functions to access to
  • external/fastjet/plugins/SISCone/siscone.h

    rf14bd6a r273e668  
    2222// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
    2323//                                                                           //
    24 // $Revision:: 325                                                          $//
    25 // $Date:: 2011-11-25 12:41:17 +0100 (Fri, 25 Nov 2011)                     $//
     24// $Revision:: 369                                                          $//
     25// $Date:: 2014-09-04 16:57:55 +0200 (Thu, 04 Sep 2014)                     $//
    2626///////////////////////////////////////////////////////////////////////////////
    2727
     
    7979
    8080  /**
     81   * compute the jets from a given particle set.
     82   * We are doing multiple passes such pass n_pass looks for jets among
     83   * all particles not put into jets during previous passes.
     84   * By default the number of passes is infinite (0).
     85   * \param _particles   list of particles
     86   * \param _radius      cone radius
     87   * \param _n_pass_max  maximum number of passes (0=full search)
     88   * \param _ptmin       minimum pT of the protojets
     89   * \param _ordering_scale    the ordering scale to decide which stable
     90   *                           cone is removed
     91   *
     92   * Note that the Csplit_merge::SM_var2_hardest_cut_off cut is not
     93   * used in the progressive removal variant.
     94   *
     95   * \return the number of jets found.
     96   */
     97  int compute_jets_progressive_removal(std::vector<Cmomentum> &_particles, double _radius,
     98                                       int _n_pass_max=0, double _ptmin=0.0,
     99                                       Esplit_merge_scale _ordering_scale=SM_pttilde);
     100
     101  /**
    81102   * recompute the jets with a different overlap parameter.
    82103   * we use the same particles and R as in the preceeding call.
     
    93114                     Esplit_merge_scale _split_merge_scale=SM_pttilde);
    94115
    95   /// list of protocones found pass-by-pass
     116  /// list of protocones found pass-by-pass (not filled by compute_jets_progressive_removal)
    96117  std::vector<std::vector<Cmomentum> > protocones_list;
    97118
     
    125146  bool rerun_allowed;         ///< is recompute_jets allowed ?
    126147  static std::ostream * _banner_ostr; ///< stream to use for banners
     148
     149  /// ensure things are initialised
     150  void _initialise_if_needed();
     151
    127152};
    128153
  • external/fastjet/plugins/SISCone/split_merge.cc

    rf14bd6a r273e668  
    2121// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
    2222//                                                                           //
    23 // $Revision:: 322                                                          $//
    24 // $Date:: 2011-11-15 10:12:36 +0100 (Tue, 15 Nov 2011)                     $//
     23// $Revision:: 370                                                          $//
     24// $Date:: 2014-09-04 17:03:15 +0200 (Thu, 04 Sep 2014)                     $//
    2525///////////////////////////////////////////////////////////////////////////////
    2626
     
    2828#include "siscone_error.h"
    2929#include "momentum.h"
    30 #include <math.h>
    3130#include <limits>   // for max
    3231#include <iostream>
     
    3433#include <sstream>
    3534#include <cassert>
     35#include <cmath>
    3636
    3737namespace siscone{
     
    229229#endif
    230230#endif
     231  _user_scale = NULL;
    231232  indices = NULL;
    232233
     
    237238
    238239  // no hardest cut (col-unsafe)
    239   SM_var2_hardest_cut_off = -1.0;
     240  SM_var2_hardest_cut_off = -numeric_limits<double>::max();
    240241
    241242  // no pt cutoff for the particles to put in p_uncol_hard
     
    555556}
    556557
     558
     559/*
     560 * remove the hardest protocone and declare it a a jet
     561 *  - protocones  list of protocones (initial jet candidates)
     562 *  - R2          cone radius (squared)
     563 *  - ptmin       minimal pT allowed for jets
     564 * return 0 on success, 1 on error
     565 *
     566 * The list of remaining particles (and the uncollinear-hard ones)
     567 * is updated.
     568 */
     569int Csplit_merge::add_hardest_protocone_to_jets(std::vector<Cmomentum> *protocones, double R2, double ptmin){
     570
     571  int i;
     572  Cmomentum *c;
     573  Cmomentum *v;
     574  double eta, phi;
     575  double dx, dy;
     576  double R;
     577  Cjet jet, jet_candidate;
     578  bool found_jet = false;
     579
     580  if (protocones->size()==0)
     581    return 1;
     582
     583  pt_min2 = ptmin*ptmin;
     584  R = sqrt(R2);
     585
     586  // browse protocones
     587  // for each of them, build the list of particles in them
     588  for (vector<Cmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){
     589    // initialise variables
     590    c = &(*p_it);
     591
     592    // note: cones have been tested => their (eta,phi) coordinates are computed
     593    eta = c->eta;
     594    phi = c->phi;
     595
     596    // NOTE: this is common to this method and add_protocones, so it
     597    // could be moved into a 'build_jet_from_protocone' method
     598    //
     599    // browse particles to create cone contents
     600    jet_candidate.v = Cmomentum();
     601    jet_candidate.pt_tilde=0;
     602    jet_candidate.contents.clear();
     603    for (i=0;i<n_left;i++){
     604      v = &(p_remain[i]);
     605      // for security, avoid including particles with infinite rapidity)
     606      // NO NEEDED ANYMORE SINCE REMOVED FROM p_remain FROM THE BEGINNING
     607      //if (fabs(v->pz)!=v->E){
     608      dx = eta - v->eta;
     609      dy = fabs(phi - v->phi);
     610      if (dy>M_PI)
     611        dy -= twopi;
     612      if (dx*dx+dy*dy<R2){
     613        jet_candidate.contents.push_back(v->parent_index);
     614        jet_candidate.v+= *v;
     615        jet_candidate.pt_tilde+= pt[v->parent_index];
     616        v->index=0;
     617      }
     618    }
     619    jet_candidate.n=jet_candidate.contents.size();
     620
     621    // set the momentum in protocones
     622    // (it was only known through eta and phi up to now)
     623    *c = jet_candidate.v;
     624    c->eta = eta; // restore exact original coords
     625    c->phi = phi; // to avoid rounding error inconsistencies
     626
     627    // set the jet range
     628    jet_candidate.range=Ceta_phi_range(eta,phi,R);
     629
     630    // check that the protojet has large enough pt
     631    if (jet_candidate.v.perp2()<pt_min2)
     632      continue;
     633
     634    // assign the split-merge (or progressive-removal) squared scale variable
     635    if (_user_scale) {
     636      // sm_var2 is the signed square of the user scale returned
     637      // for the jet candidate
     638      jet_candidate.sm_var2 = (*_user_scale)(jet_candidate);
     639      jet_candidate.sm_var2 *= abs(jet_candidate.sm_var2);
     640    } else {
     641      jet_candidate.sm_var2 = get_sm_var2(jet_candidate.v, jet_candidate.pt_tilde);
     642    }
     643
     644    // now check if it is possibly the hardest
     645    if ((! found_jet) ||
     646        (_user_scale ? _user_scale->is_larger(jet_candidate, jet)
     647                     : ptcomparison(jet_candidate, jet))){
     648      jet = jet_candidate;
     649      found_jet = true;
     650    }
     651  }
     652
     653  // make sure at least one of the jets has passed the selection
     654  if (!found_jet) return 1; 
     655 
     656  // add the jet to the list of jets
     657  jets.push_back(jet);
     658  jets[jets.size()-1].v.build_etaphi();
     659
     660#ifdef DEBUG_SPLIT_MERGE
     661  cout << "PR-Jet " << jets.size() << " [size " << next_jet.contents.size() << "]:";
     662#endif
     663   
     664  // update the list of what particles are left
     665  int p_remain_index = 0;
     666  int contents_index = 0;
     667  //sort(next_jet.contents.begin(),next_jet.contents.end());
     668  for (int index=0;index<n_left;index++){
     669    if ((contents_index<(int) jet.contents.size()) &&
     670        (p_remain[index].parent_index == jet.contents[contents_index])){
     671      // this particle belongs to the newly found jet
     672      // set pass in initial list
     673      particles[p_remain[index].parent_index].index = n_pass;
     674#ifdef DEBUG_SPLIT_MERGE
     675      cout << " " << jet.contents[contents_index];
     676#endif
     677      contents_index++;
     678    } else {
     679      // this particle still has to be clustered
     680      p_remain[p_remain_index] = p_remain[index];
     681      p_remain[p_remain_index].parent_index = p_remain[index].parent_index;
     682      p_remain[p_remain_index].index=1;
     683      p_remain_index++;
     684    }
     685  }
     686  p_remain.resize(n_left-jet.contents.size());
     687  n_left = p_remain.size();
     688  jets[jets.size()-1].pass = particles[jet.contents[0]].index;
     689
     690  // increase the pass index
     691  n_pass++;
     692
     693#ifdef DEBUG_SPLIT_MERGE
     694  cout << endl;
     695#endif
     696
     697  // male sure the list of uncol_hard particles (used for the next
     698  // stable cone finding) is updated [could probably be more
     699  // efficient]
     700  merge_collinear_and_remove_soft();
     701 
     702  return 0;
     703}
    557704
    558705/*
  • external/fastjet/plugins/SISCone/split_merge.h

    rf14bd6a r273e668  
    2222// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
    2323//                                                                           //
    24 // $Revision:: 268                                                          $//
    25 // $Date:: 2009-03-12 21:24:16 +0100 (Thu, 12 Mar 2009)                     $//
     24// $Revision:: 367                                                          $//
     25// $Date:: 2014-09-04 15:57:37 +0200 (Thu, 04 Sep 2014)                     $//
    2626///////////////////////////////////////////////////////////////////////////////
    2727
     
    141141  /// the split-merge process i.e. the variable we use for
    142142  ///  1. ordering jet candidates;
    143   ///  2. computing te overlap fraction of two candidates.
     143  ///  2. computing the overlap fraction of two candidates.
    144144  /// The default value uses pttile (p-scheme pt). Other alternatives are
    145145  /// pt, mt=sqrt(pt^2+m^2)=sqrt(E^2-pz^2) or Et.
     
    151151  ///   the default value i.e.  to use pt only for the purpose of
    152152  ///   investigating the IR issue
    153   /// - using Et is safe but do not respect boost invariance
     153  /// - using Et is safe but does not respect boost invariance
    154154  /// - using mt solves the IR unsafety issues with the pt variable
    155155  ///   for QCD jets but the IR unsafety remains for nack-to-back
     
    234234  int full_clear();
    235235
     236  ///////////////////////////////////////
     237  // user-defined stable-cone ordering //
     238  ///////////////////////////////////////
     239
     240  /// \class Cuser_scale_base
     241  /// base class for user-defined ordering of stable cones
     242  ///
     243  /// derived classes have to implement the () operator that returns
     244  /// the scale associated with a given jet.
     245  class Cuser_scale_base{
     246  public:
     247    /// empty virtual dtor
     248    virtual ~Cuser_scale_base(){}
     249
     250    /// the scale associated with a given jet
     251    ///
     252    /// "progressive removal" iteratively removes the stable cone with
     253    /// the largest scale
     254    virtual double operator()(const Cjet & jet) const = 0;
     255
     256    /// returns true when the scale associated with jet a is larger than
     257    /// the scale associated with jet b
     258    ///
     259    /// By default this does a simple direct comparison but it can be
     260    /// overloaded for higher precision [recommended if possible]
     261    ///
     262    /// This function assumes that a.sm_var2 and b.sm_var2 have been
     263    /// correctly initialised with the signed squared output of
     264    /// operator(), as is by default the case when is_larger is called
     265    /// from within siscone.
     266    virtual bool is_larger(const Cjet & a, const Cjet & b) const{
     267      return (a.sm_var2 > b.sm_var2);
     268    }
     269  };
     270
     271  /// associate a user-defined scale to order the stable cones
     272  ///
     273  /// Note that this is only used in "progressive-removal mode",
     274  /// e.g. in add_hardest_protocone_to_jets().
     275  void set_user_scale(const Cuser_scale_base * user_scale_in){
     276    _user_scale = user_scale_in;
     277  }
     278
     279  /// return the user-defined scale (NULL if none)
     280  const Cuser_scale_base * user_scale() const { return _user_scale; }
     281
    236282
    237283  /////////////////////////////////
     
    256302   */
    257303  int add_protocones(std::vector<Cmomentum> *protocones, double R2, double ptmin=0.0);
     304
     305  /**
     306   * remove the hardest protocone and declare it a jet
     307   * \param protocones  list of protocones (initial jet candidates)
     308   * \param R2          cone radius (squared)
     309   * \param ptmin       minimal pT allowed for jets
     310   * \return 0 on success, 1 on error
     311   *
     312   * The list of remaining particles (and the uncollinear-hard ones)
     313   * is updated.
     314   */
     315  int add_hardest_protocone_to_jets(std::vector<Cmomentum> *protocones, double R2, double ptmin=0.0);
    258316
    259317  /**
     
    314372  Csplit_merge_ptcomparison ptcomparison;
    315373
    316   /// stop split--merge when the SM_var of the hardest protojet
    317   /// is below this cut-off.
     374  /// stop split--merge or progressive-removal when the squared SM_var
     375  /// of the hardest protojet is below this cut-off. Note that this is
     376  /// a signed square (ie SM_var*|SM_var|) to be able to handle
     377  /// negative values.
     378  ///
     379  /// Note that the cut-off is set on the variable squared.
     380  double SM_var2_hardest_cut_off;
     381
     382  /// pt cutoff for the particles to put in p_uncol_hard
     383  /// this is meant to allow removing soft particles in the
     384  /// stable-cone search.
     385  ///
    318386  /// This is not collinear-safe so you should not use this
    319387  /// variable unless you really know what you are doing
    320388  /// Note that the cut-off is set on the variable squared.
    321   double SM_var2_hardest_cut_off;
    322 
    323   /// pt cutoff for the particles to put in p_uncol_hard
    324   /// this is meant to allow removing soft particles in the
    325   /// stable-cone search.
    326389  double stable_cone_soft_pt2_cutoff;
    327390
     
    390453  bool use_pt_weighted_splitting;
    391454
     455  /// use a user-defined scale to order the stable cones and jet
     456  /// candidates
     457  const Cuser_scale_base *_user_scale;
     458
    392459#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
    393460  /// checkxor for the candidates (to avoid having twice the same contents)
Note: See TracChangeset for help on using the changeset viewer.