Changes in external/fastjet/ClusterSequence.cc [273e668:d69dfe4] in git

File:

• external/fastjet/ClusterSequence.cc (modified) (26 diffs)

external/fastjet/ClusterSequence.cc

@@ -1 @@
-//FJSTARTHEADER
-// $Id: ClusterSequence.cc 3685 2014-09-11 20:15:00Z salam $
+//STARTHEADER
+// $Id$
 //
-// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2013, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13 +13 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development. They are described in the original FastJet paper,
-//  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
+//  development and are described in hep-ph/0512210. If you use
 //  FastJet as part of work towards a scientific publication, please
-//  quote the version you use and include a citation to the manual and
-//  optionally also to hep-ph/0512210.
+//  include a citation to the FastJet paper.
 //
 //  FastJet is distributed in the hope that it will be useful,
@@ -27 +25 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//FJENDHEADER
+//ENDHEADER
 
 #include "fastjet/Error.hh"
@@ -34 +32 @@
 #include "fastjet/ClusterSequenceStructure.hh"
 #include "fastjet/version.hh" // stores the current version number
-#include "fastjet/internal/LazyTiling9Alt.hh"
-#include "fastjet/internal/LazyTiling9.hh"
-#include "fastjet/internal/LazyTiling25.hh"
-#ifndef __FJCORE__
-#include "fastjet/internal/LazyTiling9SeparateGhosts.hh"
-#endif  // __FJCORE__
 #include<iostream>
 #include<sstream>
@@ -195 +187 @@
 //DEP //----------------------------------------------------------------------
 //DEP void ClusterSequence::_initialise_and_run (
-//DEP                             const double R,
+//DEP                             const double & R,
 //DEP                             const Strategy & strategy,
 //DEP                             const bool & writeout_combinations) {
@@ -281 +273 @@
   //             ./fastjet_timing_plugins -kt -nhardest 30 -repeat 50000 -strategy -3 -R 0.5 -nev 1  <  ../../data/Pythia-PtMin1000-LHC-1000ev.dat
   if (_strategy == Best) {
-    _strategy = _best_strategy();
-#ifdef DROP_CGAL
-    // fall back strategy for large N when CGAL is missing
-    if (_strategy == NlnN) _strategy = N2MHTLazy25;
-#endif  // DROP_CGAL
-  } else if (_strategy == BestFJ30) {
     int N = _jets.size();
     //if (N <= 55*max(0.5,min(1.0,_Rparam))) {// old empirical scaling with R
@@ -338 +324 @@
   // run the code containing the selected strategy
   // 
-  // We order the strategies starting from the ones used by the Best
+  // We order the strategies stqrting from the ones used by the Best
   // strategy in the order of increasing N, then the remaining ones
   // again in the order of increasing N.
@@ -348 +334 @@
   } else if (_strategy == N2MinHeapTiled) {
     this->_minheap_faster_tiled_N2_cluster();
-  } else if (_strategy == N2MHTLazy9Alt) {
-    // attempt to use an external tiling routine -- it manipulates
-    // the CS history via the plugin mechanism
-    _plugin_activated = true;
-    LazyTiling9Alt tiling(*this);
-    tiling.run();
-    _plugin_activated = false;
-
-  } else if (_strategy == N2MHTLazy25) {
-    // attempt to use an external tiling routine -- it manipulates
-    // the CS history via the plugin mechanism
-    _plugin_activated = true;
-    LazyTiling25 tiling(*this);
-    tiling.run();
-    _plugin_activated = false;
-
-  } else if (_strategy == N2MHTLazy9) {
-    // attempt to use an external tiling routine -- it manipulates
-    // the CS history via the plugin mechanism
-    _plugin_activated = true;
-    LazyTiling9 tiling(*this);
-    tiling.run();
-    _plugin_activated = false;
-
-#ifndef __FJCORE__
-  } else if (_strategy == N2MHTLazy9AntiKtSeparateGhosts) {
-    // attempt to use an external tiling routine -- it manipulates
-    // the CS history via the plugin mechanism
-    _plugin_activated = true;
-    LazyTiling9SeparateGhosts tiling(*this);
-    tiling.run();
-    _plugin_activated = false;
-#else 
-    throw Error("N2MHTLazy9AntiKtSeparateGhosts strategy not supported with FJCORE");
-#endif  // __FJCORE__
-
   } else if (_strategy == NlnN) {
     this->_delaunay_cluster();
@@ -409 +359 @@
 // these needs to be defined outside the class definition.
 bool ClusterSequence::_first_time = true;
-LimitedWarning ClusterSequence::_exclusive_warnings;
+int ClusterSequence::_n_exclusive_warnings = 0;
 
 
@@ -523 +473 @@
 
 //----------------------------------------------------------------------
+// Return the component corresponding to the specified index.
+// taken from CLHEP
 string ClusterSequence::strategy_string (Strategy strategy_in)  const {
   string strategy;
@@ -540 +492 @@
   case N2PoorTiled:
     strategy = "N2PoorTiled"; break;
-  case N2MHTLazy9:
-    strategy = "N2MHTLazy9"; break;
-  case N2MHTLazy9Alt:
-    strategy = "N2MHTLazy9Alt"; break;
-  case N2MHTLazy25:
-    strategy = "N2MHTLazy25"; break;
-  case N2MHTLazy9AntiKtSeparateGhosts:
-    strategy = "N2MHTLazy9AntiKtSeparateGhosts"; break;
   case N3Dumb:
     strategy = "N3Dumb"; break;
@@ -584 +528 @@
     } else {return 1.0;}
   } else {throw Error("Unrecognised jet algorithm");}
-}
-
-//----------------------------------------------------------------------
-// returns a suggestion for the best strategy to use on event
-// multiplicity, algorithm, R, etc.
-//
-// Some of the work to establish the best strategy is collected in
-// issue-tracker/2014-07-auto-strategy-selection;
-// transition_fit_v2.fit indicates the results of the fits that we're
-// using here. (Automatically generated by transition_fit_v2.gp).
-//
-// The transition to NlnN is always present, and it is the the
-// caller's responsibility to drop back down to N2MHTLazy25 if NlnN
-// isn't available.
-//
-// This routine should be called only if the jet alg is one of kt,
-// antikt, cam or genkt.
-Strategy ClusterSequence::_best_strategy() const {
-  int N = _jets.size();
-  // define bounded R, always above 0.1, because we don't trust any
-  // of our parametrizations below R = 0.1
-  double bounded_R = max(_Rparam, 0.1);
-
-  // the very first test thing is a quick hard-coded test to decide
-  // if we immediately opt for N2Plain
-  if (N <= 30 || N <= 39.0/(bounded_R + 0.6)) {
-    return N2Plain;
-  } 
-  
-  // Define objects that describe our various boundaries. A prefix N_
-  // indicates that boundary is for N, while L_ means it's for log(N).
-  //
-  // Hopefully having them static will ensure minimal overhead
-  // in creating them; collecting them in one place should
-  // help with updates?
-  //
-  const static _Parabola N_Tiled_to_MHT_lowR             (-45.4947,54.3528,44.6283);
-  const static _Parabola L_MHT_to_MHTLazy9_lowR          (0.677807,-1.05006,10.6994);
-  const static _Parabola L_MHTLazy9_to_MHTLazy25_akt_lowR(0.169967,-0.512589,12.1572);
-  const static _Parabola L_MHTLazy9_to_MHTLazy25_kt_lowR (0.16237,-0.484612,12.3373);
-  const static _Parabola L_MHTLazy9_to_MHTLazy25_cam_lowR = L_MHTLazy9_to_MHTLazy25_kt_lowR;
-  const static _Parabola L_MHTLazy25_to_NlnN_akt_lowR    (0.0472051,-0.22043,15.9196);
-  const static _Parabola L_MHTLazy25_to_NlnN_kt_lowR     (0.118609,-0.326811,14.8287);
-  const static _Parabola L_MHTLazy25_to_NlnN_cam_lowR    (0.10119,-0.295748,14.3924);
-
-  const static _Line     L_Tiled_to_MHTLazy9_medR         (-1.31304,7.29621);
-  const static _Parabola L_MHTLazy9_to_MHTLazy25_akt_medR = L_MHTLazy9_to_MHTLazy25_akt_lowR;
-  const static _Parabola L_MHTLazy9_to_MHTLazy25_kt_medR  = L_MHTLazy9_to_MHTLazy25_kt_lowR;
-  const static _Parabola L_MHTLazy9_to_MHTLazy25_cam_medR = L_MHTLazy9_to_MHTLazy25_cam_lowR;
-  const static _Parabola L_MHTLazy25_to_NlnN_akt_medR     = L_MHTLazy25_to_NlnN_akt_lowR;
-  const static _Parabola L_MHTLazy25_to_NlnN_kt_medR      = L_MHTLazy25_to_NlnN_kt_lowR;
-  const static _Parabola L_MHTLazy25_to_NlnN_cam_medR     = L_MHTLazy25_to_NlnN_cam_lowR;
-
-  const static double    N_Plain_to_MHTLazy9_largeR         = 75;
-  const static double    N_MHTLazy9_to_MHTLazy25_akt_largeR = 700;
-  const static double    N_MHTLazy9_to_MHTLazy25_kt_largeR  = 1000;
-  const static double    N_MHTLazy9_to_MHTLazy25_cam_largeR = 1000;
-  const static double    N_MHTLazy25_to_NlnN_akt_largeR     = 100000;
-  const static double    N_MHTLazy25_to_NlnN_kt_largeR      = 40000;
-  const static double    N_MHTLazy25_to_NlnN_cam_largeR     = 15000;
-
-  // We have timing studies only for kt, cam and antikt; for other
-  // algorithms we set the local jet_algorithm variable to the one of
-  // kt,cam,antikt that we think will be closest in behaviour to the
-  // other alg.
-  JetAlgorithm jet_algorithm;
-  if (_jet_algorithm == genkt_algorithm) {
-    // for genkt, then we set the local jet_algorithm variable (used
-    // only for strategy choice) to be either kt or antikt, depending on
-    // the p value.
-    double p   = jet_def().extra_param();
-    if (p < 0.0) jet_algorithm = antikt_algorithm;
-    else         jet_algorithm =     kt_algorithm;
-  } else if (_jet_algorithm == cambridge_for_passive_algorithm) {
-    // we assume (but haven't tested) that using the kt-alg timing
-    // transitions should be adequate for cambridge_for_passive_algorithm
-    jet_algorithm = kt_algorithm;
-  } else {
-    jet_algorithm = _jet_algorithm;
-  }
-
-  if (bounded_R < 0.65) {
-    // low R case
-    if          (N    < N_Tiled_to_MHT_lowR(bounded_R))              return N2Tiled;
-    double logN = log(double(N));
-    if          (logN < L_MHT_to_MHTLazy9_lowR(bounded_R))           return N2MinHeapTiled;
-    else {
-      if (jet_algorithm == antikt_algorithm){
-        if      (logN < L_MHTLazy9_to_MHTLazy25_akt_lowR(bounded_R)) return N2MHTLazy9;
-        else if (logN < L_MHTLazy25_to_NlnN_akt_lowR(bounded_R))     return N2MHTLazy25;
-        else                                                         return NlnN;
-      } else if (jet_algorithm == kt_algorithm){
-        if      (logN < L_MHTLazy9_to_MHTLazy25_kt_lowR(bounded_R))  return N2MHTLazy9;
-        else if (logN < L_MHTLazy25_to_NlnN_kt_lowR(bounded_R))      return N2MHTLazy25;
-        else                                                         return NlnN;
-      } else if (jet_algorithm == cambridge_algorithm)  {
-        if      (logN < L_MHTLazy9_to_MHTLazy25_cam_lowR(bounded_R)) return N2MHTLazy9;
-        else if (logN < L_MHTLazy25_to_NlnN_cam_lowR(bounded_R))     return N2MHTLazy25;
-        else                                                         return NlnNCam;
-      }
-    }
-  } else if (bounded_R < 0.5*pi) {
-    // medium R case
-    double logN = log(double(N));
-    if      (logN < L_Tiled_to_MHTLazy9_medR(bounded_R))             return N2Tiled;
-    else {
-      if (jet_algorithm == antikt_algorithm){
-        if      (logN < L_MHTLazy9_to_MHTLazy25_akt_medR(bounded_R)) return N2MHTLazy9;
-        else if (logN < L_MHTLazy25_to_NlnN_akt_medR(bounded_R))     return N2MHTLazy25;
-        else                                                         return NlnN;
-      } else if (jet_algorithm == kt_algorithm){
-        if      (logN < L_MHTLazy9_to_MHTLazy25_kt_medR(bounded_R))  return N2MHTLazy9;
-        else if (logN < L_MHTLazy25_to_NlnN_kt_medR(bounded_R))      return N2MHTLazy25;
-        else                                                         return NlnN;
-      } else if (jet_algorithm == cambridge_algorithm)  {
-        if      (logN < L_MHTLazy9_to_MHTLazy25_cam_medR(bounded_R)) return N2MHTLazy9;
-        else if (logN < L_MHTLazy25_to_NlnN_cam_medR(bounded_R))     return N2MHTLazy25;
-        else                                                         return NlnNCam;
-      }
-    }
-  } else {
-    // large R case (R > pi/2)
-    if      (N    < N_Plain_to_MHTLazy9_largeR)                      return N2Plain;
-    else {
-      if (jet_algorithm == antikt_algorithm){
-        if      (N < N_MHTLazy9_to_MHTLazy25_akt_largeR)             return N2MHTLazy9;
-        else if (N < N_MHTLazy25_to_NlnN_akt_largeR)                 return N2MHTLazy25;
-        else                                                         return NlnN;
-      } else if (jet_algorithm == kt_algorithm){
-        if      (N < N_MHTLazy9_to_MHTLazy25_kt_largeR)              return N2MHTLazy9;
-        else if (N < N_MHTLazy25_to_NlnN_kt_largeR)                  return N2MHTLazy25;
-        else                                                         return NlnN;
-      } else if (jet_algorithm == cambridge_algorithm)  {
-        if      (N < N_MHTLazy9_to_MHTLazy25_cam_largeR)             return N2MHTLazy9;
-        else if (N < N_MHTLazy25_to_NlnN_cam_largeR)                 return N2MHTLazy25;
-        else                                                         return NlnNCam;
-      }
-    }
-  }
-  
-  bool code_should_never_reach_here = false;
-  assert(code_should_never_reach_here); 
-  return N2MHTLazy9;
-
 }
 
@@ -862 +662 @@
 //----------------------------------------------------------------------
 // return all inclusive jets with pt > ptmin
-vector<PseudoJet> ClusterSequence::inclusive_jets (const double ptmin) const{
+vector<PseudoJet> ClusterSequence::inclusive_jets (const double & ptmin) const{
   double dcut = ptmin*ptmin;
   int i = _history.size() - 1; // last jet
@@ -914 +714 @@
 // return the number of exclusive jets that would have been obtained
 // running the algorithm in exclusive mode with the given dcut
-int ClusterSequence::n_exclusive_jets (const double dcut) const {
+int ClusterSequence::n_exclusive_jets (const double & dcut) const {
 
   // first locate the point where clustering would have stopped (i.e. the
@@ -933 +733 @@
 // return all exclusive jets that would have been obtained running
 // the algorithm in exclusive mode with the given dcut
-vector<PseudoJet> ClusterSequence::exclusive_jets (const double dcut) const {
+vector<PseudoJet> ClusterSequence::exclusive_jets (const double & dcut) const {
   int njets = n_exclusive_jets(dcut);
   return exclusive_jets(njets);
@@ -942 +742 @@
 // return the jets obtained by clustering the event to n jets.
 // Throw an error if there are fewer than n particles.
-vector<PseudoJet> ClusterSequence::exclusive_jets (const int njets) const {
+vector<PseudoJet> ClusterSequence::exclusive_jets (const int & njets) const {
 
   // make sure the user does not ask for more than jets than there
@@ -959 +759 @@
 // return the jets obtained by clustering the event to n jets.
 // If there are fewer than n particles, simply return all particles
-vector<PseudoJet> ClusterSequence::exclusive_jets_up_to (const int njets) const {
+vector<PseudoJet> ClusterSequence::exclusive_jets_up_to (const int & njets) const {
 
   // provide a warning when extracting exclusive jets for algorithms 
   // that does not support it explicitly.
-  // Native algorithm that support it are: kt, ee_kt, Cambridge/Aachen, 
+  // Native algorithm that support it are: kt, ee_kt, cambridge, 
   //   genkt and ee_genkt (both with p>=0)
   // For plugins, we check Plugin::exclusive_sequence_meaningful()
@@ -973 +773 @@
        (_jet_def.extra_param() <0)) &&
       ((_jet_def.jet_algorithm() != plugin_algorithm) ||
-       (!_jet_def.plugin()->exclusive_sequence_meaningful()))) {
-    _exclusive_warnings.warn("dcut and exclusive jets for jet-finders other than kt, C/A or genkt with p>=0 should be interpreted with care.");
+       (!_jet_def.plugin()->exclusive_sequence_meaningful())) &&
+      (_n_exclusive_warnings < 5)) {
+    _n_exclusive_warnings++;
+    cerr << "FastJet WARNING: dcut and exclusive jets for jet-finders other than kt should be interpreted with care." << endl;
   }
 
@@ -1027 +829 @@
 /// return the dmin corresponding to the recombination that went from
 /// n+1 to n jets
-double ClusterSequence::exclusive_dmerge (const int njets) const {
+double ClusterSequence::exclusive_dmerge (const int & njets) const {
   assert(njets >= 0);
   if (njets >= _initial_n) {return 0.0;}
@@ -1039 +841 @@
 /// exclusive_dmerge, except in cases where the dmin do not increase
 /// monotonically.
-double ClusterSequence::exclusive_dmerge_max (const int njets) const {
+double ClusterSequence::exclusive_dmerge_max (const int & njets) const {
   assert(njets >= 0);
   if (njets >= _initial_n) {return 0.0;}
@@ -1051 +853 @@
 /// the algorithm with the given dcut.
 std::vector<PseudoJet> ClusterSequence::exclusive_subjets 
-   (const PseudoJet & jet, const double dcut) const {
+   (const PseudoJet & jet, const double & dcut) const {
 
   set<const history_element*> subhist;
@@ -1074 +876 @@
 /// exclusive_subjets.size()
 int ClusterSequence::n_exclusive_subjets(const PseudoJet & jet, 
-                        const double dcut) const {
+                        const double & dcut) const {
   set<const history_element*> subhist;
   // get the set of history elements that correspond to subjets at
@@ -1379 +1181 @@
 // //----------------------------------------------------------------------
 // // print out all inclusive jets with pt > ptmin
-// void ClusterSequence::print_jets (const double ptmin) const{
+// void ClusterSequence::print_jets (const double & ptmin) const{
 //     vector<PseudoJet> jets = sorted_by_pt(inclusive_jets(ptmin));
 // 
@@ -1455 +1257 @@
 // initialise the history in a standard way
 void ClusterSequence::_add_step_to_history (
-               const int step_number, const int parent1, 
-               const int parent2, const int jetp_index,
-               const double dij) {
+               const int & step_number, const int & parent1, 
+               const int & parent2, const int & jetp_index,
+               const double & dij) {
 
   history_element element;
@@ -1627 +1429 @@
 /// of the recombined jet, newjet_k.
 void ClusterSequence::_do_ij_recombination_step(
-                               const int jet_i, const int jet_j, 
-                               const double dij, 
+                               const int & jet_i, const int & jet_j, 
+                               const double & dij, 
                                int & newjet_k) {
 
@@ -1664 +1466 @@
 /// jet_i with the beam
 void ClusterSequence::_do_iB_recombination_step(
-                                  const int jet_i, const double diB) {
+                                  const int & jet_i, const double & diB) {
   // get history index
   int newstep_k = _history.size();