JETAN/fastjet/fastjet/JetDefinition.hh

   1 //STARTHEADER
   2 // $Id: JetDefinition.hh 1402 2009-01-21 18:03:54Z salam $
   3 //
   4 // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
   5 //
   6 //----------------------------------------------------------------------
   7 // This file is part of FastJet.
   8 //
   9 //  FastJet is free software; you can redistribute it and/or modify
  10 //  it under the terms of the GNU General Public License as published by
  11 //  the Free Software Foundation; either version 2 of the License, or
  12 //  (at your option) any later version.
  13 //
  14 //  The algorithms that underlie FastJet have required considerable
  15 //  development and are described in hep-ph/0512210. If you use
  16 //  FastJet as part of work towards a scientific publication, please
  17 //  include a citation to the FastJet paper.
  18 //
  19 //  FastJet is distributed in the hope that it will be useful,
  20 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  21 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22 //  GNU General Public License for more details.
  23 //
  24 //  You should have received a copy of the GNU General Public License
  25 //  along with FastJet; if not, write to the Free Software
  26 //  Foundation, Inc.:
  27 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  28 //----------------------------------------------------------------------
  29 //ENDHEADER
  30
  31 #ifndef __FASTJET_JETDEFINITION_HH__
  32 #define __FASTJET_JETDEFINITION_HH__
  33
  34 #include<cassert>
  35 #include "fastjet/internal/numconsts.hh"
  36 #include "fastjet/PseudoJet.hh"
  37 #include<string>
  38 #include<memory>
  39
  40 namespace fastjet {      // defined in fastjet/internal/base.hh
  41
  42 /// return a string containing information about the release
  43 //  NB: (implemented in ClusterSequence.cc but defined here because
  44 //  this is a visible location)
  45 std::string fastjet_version_string();
  46
  47 //======================================================================
  48 /// the various options for the algorithmic strategy to adopt in
  49 /// clustering events with kt and cambridge style algorithms.
  50 enum Strategy {
  51   /// fastest form about 500..10^4
  52   N2MinHeapTiled   = -4,
  53   /// fastest from about 50..500
  54   N2Tiled     = -3,
  55   /// legacy
  56   N2PoorTiled = -2,
  57   /// fastest below 50
  58   N2Plain     = -1,
  59   /// worse even than the usual N^3 algorithms
  60   N3Dumb      =  0,
  61   /// automatic selection of the best (based on N)
  62   Best        =  1,
  63   /// best of the NlnN variants -- best overall for N>10^4
  64   NlnN        =  2,
  65   /// legacy N ln N using 3pi coverage of cylinder
  66   NlnN3pi     =  3,
  67   /// legacy N ln N using 4pi coverage of cylinder
  68   NlnN4pi     =  4,
  69   /// Chan's closest pair method (in a variant with 4pi coverage),
  70   /// for use exclusively with the Cambridge algorithm
  71   NlnNCam4pi   = 14,
  72   NlnNCam2pi2R = 13,
  73   NlnNCam      = 12, // 2piMultD
  74   /// the plugin has been used...
  75   plugin_strategy = 999
  76 };
  77
  78
  79 //======================================================================
  80 /// the various families of jet-clustering algorithm
  81 enum JetAlgorithm {
  82   /// the longitudinally invariant kt algorithm
  83   kt_algorithm=0,
  84   /// the longitudinally invariant variant of the cambridge algorithm
  85   /// (aka Aachen algoithm).
  86   cambridge_algorithm=1,
  87   /// like the k_t but with distance measures
  88   ///       dij = min(1/kti^2,1/ktj^2) Delta R_{ij}^2 / R^2
  89   ///       diB = 1/kti^2
  90   antikt_algorithm=2,
  91   /// like the k_t but with distance measures
  92   ///       dij = min(kti^{2p},ktj^{2p}) Delta R_{ij}^2 / R^2
  93   ///       diB = 1/kti^{2p}
  94   /// where p = extra_param()
  95   genkt_algorithm=3,
  96   /// a version of cambridge with a special distance measure for particles
  97   /// whose pt is < extra_param()
  98   cambridge_for_passive_algorithm=11,
  99   /// a version of genkt with a special distance measure for particles
 100   /// whose pt is < extra_param() [relevant for passive areas when p<=0]
 101   genkt_for_passive_algorithm=13,
 102   //.................................................................
 103   /// the e+e- kt algorithm
 104   ee_kt_algorithm=50,
 105   /// the e+e- genkt algorithm  (R > 2 and p=1 gives ee_kt)
 106   ee_genkt_algorithm=53,
 107   //.................................................................
 108   /// any plugin algorithm supplied by the user
 109   plugin_algorithm = 99
 110 };
 111
 112 /// make standard Les Houches nomenclature JetAlgorithm (algorithm is general
 113 /// recipe without the parameters) backward-compatible with old JetFinder
 114 typedef JetAlgorithm JetFinder;
 115
 116 /// provide other possible names for the Cambridge/Aachen algorithm?
 117 const JetAlgorithm aachen_algorithm = cambridge_algorithm;
 118 const JetAlgorithm cambridge_aachen_algorithm = cambridge_algorithm;
 119
 120 //======================================================================
 121 /// the various recombination schemes
 122 enum RecombinationScheme {
 123   /// summing the 4-momenta
 124   E_scheme=0,
 125   /// pt weighted recombination of y,phi (and summing of pt's)
 126   /// with preprocessing to make things massless by rescaling E=|\vec p|
 127   pt_scheme=1,
 128   /// pt^2 weighted recombination of y,phi (and summing of pt's)
 129   /// with preprocessing to make things massless by rescaling E=|\vec p|
 130   pt2_scheme=2,
 131   /// pt weighted recombination of y,phi (and summing of pt's)
 132   /// with preprocessing to make things massless by rescaling |\vec p|->=E
 133   Et_scheme=3,
 134   /// pt^2 weighted recombination of y,phi (and summing of pt's)
 135   /// with preprocessing to make things massless by rescaling |\vec p|->=E
 136   Et2_scheme=4,
 137   /// pt weighted recombination of y,phi (and summing of pt's), with
 138   /// no preprocessing
 139   BIpt_scheme=5,
 140   /// pt^2 weighted recombination of y,phi (and summing of pt's)
 141   /// no preprocessing
 142   BIpt2_scheme=6,
 143   /// for the user's external scheme
 144   external_scheme = 99
 145 };
 146
 147
 148
 149
 150 // forward declaration, needed in order to specify interface for the
 151 // plugin.
 152 class ClusterSequence;
 153
 154
 155
 156
 157 //======================================================================
 158 /// class that is intended to hold a full definition of the jet
 159 /// clusterer
 160 class JetDefinition {
 161
 162 public:
 163
 164   /// forward declaration of a class that allows the user to introduce
 165   /// their own plugin
 166   class Plugin;
 167
 168   // forward declaration of a class that will provide the
 169   // recombination scheme facilities and/or allow a user to
 170   // extend these facilities
 171   class Recombiner;
 172
 173
 174
 175   /// constructor with alternative ordering or arguments -- note that
 176   /// we have not provided a default jet finder, to avoid ambiguous
 177   /// JetDefinition() constructor.
 178   JetDefinition(JetAlgorithm jet_algorithm,
 179                 double R,
 180                 RecombinationScheme recomb_scheme = E_scheme,
 181                 Strategy strategy = Best) {
 182     *this = JetDefinition(jet_algorithm, R, strategy, recomb_scheme, 1);
 183   }
 184
 185   /// constructor for algorithms that have no free parameters
 186   /// (e.g. ee_kt_algorithm)
 187   JetDefinition(JetAlgorithm jet_algorithm,
 188                 RecombinationScheme recomb_scheme = E_scheme,
 189                 Strategy strategy = Best) {
 190     double dummyR = 0.0;
 191     *this = JetDefinition(jet_algorithm, dummyR, strategy, recomb_scheme, 0);
 192   }
 193
 194   /// constructor for algorithms that require R + one extra parameter to be set
 195   /// (the gen-kt series for example)
 196   JetDefinition(JetAlgorithm jet_algorithm,
 197                 double R,
 198                 double xtra_param,
 199                 RecombinationScheme recomb_scheme = E_scheme,
 200                 Strategy strategy = Best) {
 201     *this = JetDefinition(jet_algorithm, R, strategy, recomb_scheme, 2);
 202     set_extra_param(xtra_param);
 203   }
 204
 205
 206   /// constructor in a form that allows the user to provide a pointer
 207   /// to an external recombiner class (which must remain valid for the
 208   /// life of the JetDefinition object).
 209   JetDefinition(JetAlgorithm jet_algorithm,
 210                 double R,
 211                 const Recombiner * recombiner,
 212                 Strategy strategy = Best) {
 213     *this = JetDefinition(jet_algorithm, R, external_scheme, strategy);
 214     _recombiner = recombiner;
 215   }
 216
 217
 218   /// constructor for case with 0 parameters (ee_kt_algorithm) and
 219   /// and external recombiner
 220   JetDefinition(JetAlgorithm jet_algorithm,
 221                 const Recombiner * recombiner,
 222                 Strategy strategy = Best) {
 223     *this = JetDefinition(jet_algorithm, external_scheme, strategy);
 224     _recombiner = recombiner;
 225   }
 226
 227   /// constructor allowing the extra parameter to be set and a pointer to
 228   /// a recombiner
 229   JetDefinition(JetAlgorithm jet_algorithm,
 230                 double R,
 231                 double xtra_param,
 232                 const Recombiner * recombiner,
 233                 Strategy strategy = Best) {
 234     *this = JetDefinition(jet_algorithm, R, external_scheme, strategy);
 235     _recombiner = recombiner;
 236     set_extra_param(xtra_param);
 237   }
 238
 239   /// a default constructor
 240   JetDefinition() {
 241     *this = JetDefinition(kt_algorithm, 1.0);
 242   }
 243
 244   /// constructor based on a pointer to a user's plugin; the object
 245   /// pointed to must remain valid for the whole duration of existence
 246   /// of the JetDefinition and any related ClusterSequences
 247   JetDefinition(const Plugin * plugin) {
 248     _plugin = plugin;
 249     _strategy = plugin_strategy;
 250     _Rparam = _plugin->R();
 251     _jet_algorithm = plugin_algorithm;
 252     set_recombination_scheme(E_scheme);
 253   }
 254
 255
 256   /// constructor to fully specify a jet-definition (together with
 257   /// information about how algorithically to run it).
 258   ///
 259   /// the ordering of arguments here is old and deprecated (except
 260   /// as the common constructor for internal use)
 261   JetDefinition(JetAlgorithm jet_algorithm,
 262                 double R,
 263                 Strategy strategy,
 264                 RecombinationScheme recomb_scheme = E_scheme,
 265                 int nparameters = 1);
 266 // :
 267 //    _jet_algorithm(jet_algorithm), _Rparam(R), _strategy(strategy) {
 268 //    // the largest sensible value for R
 269 //    if (jet_algorithm != ee_kt_algorithm &&
 270 //      jet_algorithm != ee_genkt_algorithm) assert(_Rparam <= 0.5*pi);
 271 //    assert(_jet_algorithm != plugin_algorithm &&
 272 //           _strategy   != plugin_strategy);
 273 //    _plugin = NULL;
 274 //    set_recombination_scheme(recomb_scheme);
 275 //    set_extra_param(0.0); // make sure it's defined
 276 //  }
 277
 278
 279   /// set the recombination scheme to the one provided
 280   void set_recombination_scheme(RecombinationScheme);
 281
 282   /// set the recombiner class to the one provided
 283   void set_recombiner(const Recombiner * recomb) {
 284     _recombiner = recomb;
 285     _default_recombiner = DefaultRecombiner(external_scheme);
 286   }
 287
 288   /// return a pointer to the plugin
 289   const Plugin * plugin() const {return _plugin;};
 290
 291   /// return information about the definition...
 292   JetAlgorithm jet_algorithm  () const {return _jet_algorithm  ;}
 293   /// same as above for backward compatibility
 294   JetAlgorithm jet_finder     () const {return _jet_algorithm  ;}
 295   double    R           () const {return _Rparam      ;}
 296   // a general purpose extra parameter, whose meaning depends on
 297   // the algorithm, and may often be unused.
 298   double    extra_param () const {return _extra_param ;}
 299   Strategy  strategy    () const {return _strategy    ;}
 300   RecombinationScheme recombination_scheme() const {
 301     return _default_recombiner.scheme();}
 302
 303   /// (re)set the jet finder
 304   void set_jet_algorithm(JetAlgorithm njf) {_jet_algorithm = njf;}
 305   /// same as above for backward compatibility
 306   void set_jet_finder(JetAlgorithm njf)    {_jet_algorithm = njf;}
 307   /// (re)set the general purpose extra parameter
 308   void set_extra_param(double xtra_param) {_extra_param = xtra_param;}
 309
 310   /// return a pointer to the currently defined recombiner (it may
 311   /// be the internal one)
 312   const Recombiner * recombiner() const {
 313     return _recombiner == 0 ? & _default_recombiner : _recombiner;}
 314
 315   /// return a textual description of the current jet definition
 316   std::string description() const;
 317
 318
 319 public:
 320   //======================================================================
 321   /// An abstract base class that will provide the recombination scheme
 322   /// facilities and/or allow a user to extend these facilities
 323   class Recombiner {
 324   public:
 325     /// return a textual description of the recombination scheme
 326     /// implemented here
 327     virtual std::string description() const = 0;
 328
 329     /// recombine pa and pb and put result into pab
 330     virtual void recombine(const PseudoJet & pa, const PseudoJet & pb,
 331                            PseudoJet & pab) const = 0;
 332
 333     /// routine called to preprocess each input jet (to make all input
 334     /// jets compatible with the scheme requirements (e.g. massless).
 335     virtual void preprocess(PseudoJet & p) const {};
 336
 337     /// a destructor to be replaced if necessary in derived classes...
 338     virtual ~Recombiner() {};
 339
 340     /// pa += pb in the given recombination scheme. Not virtual -- the
 341     /// user should have no reason to want to redefine this!
 342     inline void plus_equal(PseudoJet & pa, const PseudoJet & pb) const {
 343       // put result in a temporary location in case the recombiner
 344       // does something funny (ours doesn't, but who knows about the
 345       // user's)
 346       PseudoJet pres;
 347       recombine(pa,pb,pres);
 348       pa = pres;
 349     }
 350
 351   };
 352
 353
 354   //======================================================================
 355   /// A class that will provide the recombination scheme facilities and/or
 356   /// allow a user to extend these facilities
 357   class DefaultRecombiner : public Recombiner {
 358   public:
 359     DefaultRecombiner(RecombinationScheme recomb_scheme = E_scheme) :
 360       _recomb_scheme(recomb_scheme) {}
 361
 362     virtual std::string description() const;
 363
 364     /// recombine pa and pb and put result into pab
 365     virtual void recombine(const PseudoJet & pa, const PseudoJet & pb,
 366                            PseudoJet & pab) const;
 367
 368     virtual void preprocess(PseudoJet & p) const;
 369
 370     /// return the index of the recombination scheme
 371     RecombinationScheme scheme() const {return _recomb_scheme;}
 372
 373   private:
 374     RecombinationScheme _recomb_scheme;
 375   };
 376
 377
 378   //======================================================================
 379   /// a class that allows a user to introduce their own "plugin" jet
 380   /// finder
 381   class Plugin{
 382   public:
 383     /// return a textual description of the jet-definition implemented
 384     /// in this plugin
 385     virtual std::string description() const = 0;
 386
 387     /// given a ClusterSequence that has been filled up with initial
 388     /// particles, the following function should fill up the rest of the
 389     /// ClusterSequence, using the following member functions of
 390     /// ClusterSequence:
 391     ///   - plugin_do_ij_recombination(...)
 392     ///   - plugin_do_iB_recombination(...)
 393     virtual void run_clustering(ClusterSequence &) const = 0;
 394
 395     virtual double R() const = 0;
 396
 397     /// return true if there is specific support for the measurement
 398     /// of passive areas, in the sense that areas determined from all
 399     /// particles below the ghost separation scale will be a passive
 400     /// area. [If you don't understand this, ignore it!]
 401     virtual bool supports_ghosted_passive_areas() const {return false;}
 402
 403     /// set the ghost separation scale for passive area determinations
 404     /// in future runs (strictly speaking that makes the routine
 405     /// a non const, so related internal info must be stored as a mutable)
 406     virtual void set_ghost_separation_scale(double scale) const;
 407     virtual double ghost_separation_scale() const {return 0.0;}
 408
 409     /// a destructor to be replaced if necessary in derived classes...
 410     virtual ~Plugin() {};
 411   };
 412
 413 private:
 414
 415
 416   JetAlgorithm _jet_algorithm;
 417   double    _Rparam;
 418   double    _extra_param ; ///< parameter whose meaning varies according to context
 419   Strategy  _strategy  ;
 420
 421   const Plugin * _plugin;
 422
 423   // when we use our own recombiner it's useful to point to it here
 424   // so that we don't have to worry about deleting it etc...
 425   DefaultRecombiner _default_recombiner;
 426   const Recombiner * _recombiner;
 427
 428 };
 429
 430
 431
 432
 433
 434
 435 } // fastjet namespace
 436
 437 #endif // __FASTJET_JETDEFINITION_HH__