Changeset 49234af in git for external/fastjet/JetDefinition.cc

Timestamp:

Dec 9, 2014, 1:27:13 PM (10 years ago)

Author:

Michele <michele.selvaggi@…>

Branches:

ImprovedOutputFile, Timing, dual_readout, llp, master

Children:

37deb3b, 9e991f8

Parents:

f6b6ee7 (diff), e7e90df (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'TestFastJet310b1'

File:

• external/fastjet/JetDefinition.cc (modified) (11 diffs)

external/fastjet/JetDefinition.cc

@@ -1 @@
-//STARTHEADER
-// $Id$
-//
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//FJSTARTHEADER
+// $Id: JetDefinition.cc 3677 2014-09-09 22:45:25Z soyez $
+//
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13 +13 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  development. They are described in the original FastJet paper,
+//  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
 //  FastJet as part of work towards a scientific publication, please
-//  include a citation to the FastJet paper.
+//  quote the version you use and include a citation to the manual and
+//  optionally also to hep-ph/0512210.
 //
 //  FastJet is distributed in the hope that it will be useful,
@@ -25 +27 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #include "fastjet/JetDefinition.hh"
@@ -68 +70 @@
   // cross-check the number of parameters that were declared in setting up the
   // algorithm (passed internally from the public constructors)
-  switch (jet_algorithm_in) {
-  case ee_kt_algorithm:
-    if (nparameters != 0) {
-      ostringstream oss;
-      oss << "ee_kt_algorithm should be constructed with 0 parameters but was called with " 
-          << nparameters << " parameter(s)\n";
-      throw Error(oss.str()); 
-    }
-    break;
-  case genkt_algorithm: 
-  case ee_genkt_algorithm: 
-    if (nparameters != 2) {
-      ostringstream oss;
-      oss << "(ee_)genkt_algorithm should be constructed with 2 parameters but was called with " 
-          << nparameters << " parameter(s)\n";
-      throw Error(oss.str()); 
-    }
-    break;
-  default:
-    if (nparameters != 1) {
-      ostringstream oss;
-      oss << "The jet algorithm you requested ("
-          << jet_algorithm_in << ") should be constructed with 1 parameter but was called with " 
-          << nparameters << " parameter(s)\n";
-      throw Error(oss.str()); 
-    }
+  unsigned int nparameters_expected = n_parameters_for_algorithm(jet_algorithm_in);
+  if (nparameters != (int) nparameters_expected){
+    ostringstream oss;
+    oss << "The jet algorithm you requested ("
+        << jet_algorithm_in << ") should be constructed with " << nparameters_expected 
+        << " parameter(s) but was called with " << nparameters << " parameter(s)\n";
+    throw Error(oss.str()); 
   }
 
@@ -106 +89 @@
 
 //----------------------------------------------------------------------
+// returns true if the jet definition involves an algorithm
+// intended for use on a spherical geometry (e.g. e+e- algorithms,
+// as opposed to most pp algorithms, which use a cylindrical,
+// rapidity-phi geometry).
+bool JetDefinition::is_spherical() const {
+  if (jet_algorithm() == plugin_algorithm) {
+    return plugin()->is_spherical();
+  } else {
+    return (jet_algorithm() == ee_kt_algorithm ||  // as of 2013-02-14, the two
+            jet_algorithm() == ee_genkt_algorithm  // native spherical algorithms
+            );
+  }
+}
+
+//----------------------------------------------------------------------
 string JetDefinition::description() const {
+  ostringstream name;
+  
+  name << description_no_recombiner();
+
+  if ((jet_algorithm() == plugin_algorithm) || (jet_algorithm() == undefined_jet_algorithm)){
+    return name.str();
+  }
+
+  if (n_parameters_for_algorithm(jet_algorithm()) == 0)
+    name << " with ";
+  else 
+    name << " and ";
+  name << recombiner()->description();
+
+  return name.str();
+}
+
+//----------------------------------------------------------------------
+string JetDefinition::description_no_recombiner() const {
+  
   ostringstream name;
   if (jet_algorithm() == plugin_algorithm) {
     return plugin()->description();
-  } else if (jet_algorithm() == kt_algorithm) {
-    name << "Longitudinally invariant kt algorithm with R = " << R();
-    name << " and " << recombiner()->description();
-  } else if (jet_algorithm() == cambridge_algorithm) {
-    name << "Longitudinally invariant Cambridge/Aachen algorithm with R = " 
-         << R() ;
-    name << " and " << recombiner()->description();
-  } else if (jet_algorithm() == antikt_algorithm) {
-    name << "Longitudinally invariant anti-kt algorithm with R = " 
-         << R() ;
-    name << " and " << recombiner()->description();
-  } else if (jet_algorithm() == genkt_algorithm) {
-    name << "Longitudinally invariant generalised kt algorithm with R = " 
-         << R() << ", p = " << extra_param();
-    name << " and " << recombiner()->description();
-  } else if (jet_algorithm() == cambridge_for_passive_algorithm) {
-    name << "Longitudinally invariant Cambridge/Aachen algorithm with R = " 
-         << R() << "and a special hack whereby particles with kt < " 
-         << extra_param() << "are treated as passive ghosts";
-  } else if (jet_algorithm() == ee_kt_algorithm) {
-    name << "e+e- kt (Durham) algorithm (NB: no R)";
-    name << " with " << recombiner()->description();
-  } else if (jet_algorithm() == ee_genkt_algorithm) {
-    name << "e+e- generalised kt algorithm with R = " 
-         << R() << ", p = " << extra_param();
-    name << " and " << recombiner()->description();
   } else if (jet_algorithm() == undefined_jet_algorithm) {
-    name << "uninitialised JetDefinition (jet_algorithm=undefined_jet_algorithm)" ;
-  } else {
-    throw Error("JetDefinition::description(): unrecognized jet_algorithm");
-  }
+    return "uninitialised JetDefinition (jet_algorithm=undefined_jet_algorithm)" ;
+  }
+
+  name << algorithm_description(jet_algorithm());
+  switch (n_parameters_for_algorithm(jet_algorithm())){
+  case 0: name << " (NB: no R)"; break;
+  case 1: name << " with R = " << R(); break; // the parameter is always R
+  case 2: 
+    // the 1st parameter is always R
+    name << " with R = " << R();
+    // the 2nd depends on the algorithm
+    if (jet_algorithm() == cambridge_for_passive_algorithm){
+      name << "and a special hack whereby particles with kt < " 
+           << extra_param() << "are treated as passive ghosts";
+    } else {
+      name << ", p = " << extra_param();
+    }
+  };
+
   return name.str();
 }
 
-
+//----------------------------------------------------------------------
+string JetDefinition::algorithm_description(const JetAlgorithm jet_alg){
+  ostringstream name;
+  switch (jet_alg){
+  case plugin_algorithm:                return "plugin algorithm";
+  case kt_algorithm:                    return "Longitudinally invariant kt algorithm";
+  case cambridge_algorithm:             return "Longitudinally invariant Cambridge/Aachen algorithm";
+  case antikt_algorithm:                return "Longitudinally invariant anti-kt algorithm";
+  case genkt_algorithm:                 return "Longitudinally invariant generalised kt algorithm";
+  case cambridge_for_passive_algorithm: return "Longitudinally invariant Cambridge/Aachen algorithm";
+  case ee_kt_algorithm:                 return "e+e- kt (Durham) algorithm (NB: no R)";
+  case ee_genkt_algorithm:              return "e+e- generalised kt algorithm";
+  case undefined_jet_algorithm:         return "undefined jet algorithm";
+  default:
+    throw Error("JetDefinition::algorithm_description(): unrecognized jet_algorithm");
+  };
+}
+
+//----------------------------------------------------------------------
+unsigned int JetDefinition::n_parameters_for_algorithm(const JetAlgorithm jet_alg){
+  switch (jet_alg) {
+  case ee_kt_algorithm:    return 0;
+  case genkt_algorithm:
+  case ee_genkt_algorithm: return 2;
+  default:                 return 1;
+  };
+}
+
+//----------------------------------------------------------------------
 void JetDefinition::set_recombination_scheme(
                                RecombinationScheme recomb_scheme) {
@@ -150 +185 @@
 
   // do not forget to delete the existing recombiner if needed
-  if (_recombiner_shared()) _recombiner_shared.reset();
+  if (_shared_recombiner()) _shared_recombiner.reset();
 
   _recombiner = 0;
+}
+
+void JetDefinition::set_recombiner(const JetDefinition &other_jet_def){
+  // make sure the "invariants" of the other jet def are sensible
+  assert(other_jet_def._recombiner || 
+         other_jet_def.recombination_scheme() != external_scheme);
+
+  // first treat the situation where we're using the default recombiner
+  if (other_jet_def._recombiner == 0){
+    set_recombination_scheme(other_jet_def.recombination_scheme());
+    return;
+  }
+
+  // in other cases, copy the pointer to the recombiner
+  _recombiner = other_jet_def._recombiner;
+  // set the default recombiner appropriately
+  _default_recombiner = DefaultRecombiner(external_scheme);
+  // and set the _shared_recombiner to the same state
+  // as in the other_jet_def, whatever that was
+  _shared_recombiner.reset(other_jet_def._shared_recombiner);
+
+  // NB: it is tempting to go via set_recombiner and then to sort
+  // out the shared part, but this would be dangerous in the
+  // specific (rare?) case where other_jet_def is the same as this
+  // it deletes_recombiner_when_unused. In that case the shared
+  // pointer reset would delete the recombiner.
 }
 
@@ -163 +224 @@
   if (other_jd.recombination_scheme() != scheme) return false;
 
-  // if the scheme is "external", also check that they ahve the same
+  // if the scheme is "external", also check that they have the same
   // recombiner
   return (scheme != external_scheme) 
@@ -169 +230 @@
 }
 
-/// allows to let the JetDefinition handle the deletion of the
+/// causes the JetDefinition to handle the deletion of the
 /// recombiner when it is no longer used
 void JetDefinition::delete_recombiner_when_unused(){
   if (_recombiner == 0){
     throw Error("tried to call JetDefinition::delete_recombiner_when_unused() for a JetDefinition without a user-defined recombination scheme");
-  }
-
-  _recombiner_shared.reset(_recombiner);
+  } else if (_shared_recombiner.get()) {
+    throw Error("Error in JetDefinition::delete_recombiner_when_unused: the recombiner is already scheduled for deletion when unused (or was already set as shared)");
+  }
+
+  _shared_recombiner.reset(_recombiner);
 }
 
@@ -207 +270 @@
   case BIpt2_scheme:
     return "boost-invariant pt2 scheme recombination";
+  case WTA_pt_scheme:
+    return "pt-ordered Winner-Takes-All recombination";
+  // Energy-ordering can lead to dangerous situations with particles at
+  // rest. We instead implement the WTA_modp_scheme
+  //
+  //   case WTA_E_scheme:
+  //     return "energy-ordered Winner-Takes-All recombination";
+  case WTA_modp_scheme:
+    return "|3-momentum|-ordered Winner-Takes-All recombination";
   default:
     ostringstream err;
@@ -246 +318 @@
     weightb = pb.perp2();
     break;
+  case WTA_pt_scheme:{
+    const PseudoJet & phard = (pa.pt2() >= pb.pt2()) ? pa : pb;
+    /// keep y,phi and m from the hardest, sum pt
+    pab.reset_PtYPhiM(pa.pt()+pb.pt(), 
+                      phard.rap(), phard.phi(), phard.m());
+    return;}
+  // Energy-ordering can lead to dangerous situations with particles at
+  // rest. We instead implement the WTA_modp_scheme
+  //
+  //   case WTA_E_scheme:{
+  //     const PseudoJet & phard = (pa.E() >= pb.E()) ? pa : pb;
+  //     /// keep 3-momentum direction and mass from the hardest, sum energies
+  //     ///
+  //     /// If the particle with the largest energy is at rest, the sum
+  //     /// remains at rest, implying that the mass of the sum is larger
+  //     /// than the mass of pa.
+  //     double Eab = pa.E() + pb.E();
+  //     double scale = (phard.modp2()==0.0)
+  //       ? 0.0
+  //       : sqrt((Eab*Eab - phard.m2())/phard.modp2());
+  //     pab.reset(phard.px()*scale, phard.py()*scale, phard.pz()*scale, Eab);
+  //     return;}
+  case WTA_modp_scheme:{
+    // Note: we need to compute both a and b modp. And we need pthard
+    // and its modp. If we want to avoid repeating the test and do
+    // only 2 modp calculations, we'd have to duplicate the code (or
+    // use a pair<const PJ&>). An alternative is to write modp_soft as
+    // modp_ab-modp_hard but this could suffer from larger rounding
+    // errors
+    bool a_hardest = (pa.modp2() >= pb.modp2());
+    const PseudoJet & phard = a_hardest ? pa : pb;
+    const PseudoJet & psoft = a_hardest ? pb : pa;
+    /// keep 3-momentum direction and mass from the hardest, sum modp
+    ///
+    /// If the hardest particle is at rest, the sum remains at rest
+    /// (the energy of the sum is therefore the mass of pa)
+    double modp_hard = phard.modp();
+    double modp_ab = modp_hard + psoft.modp();
+    if (phard.modp2()==0.0){
+      pab.reset(0.0, 0.0, 0.0, phard.m());
+    } else {
+      double scale = modp_ab/modp_hard;
+      pab.reset(phard.px()*scale, phard.py()*scale, phard.pz()*scale,
+                sqrt(modp_ab*modp_ab + phard.m2()));
+    }
+    return;}
   default:
     ostringstream err;
@@ -281 +399 @@
   case BIpt_scheme:
   case BIpt2_scheme:
+  case WTA_pt_scheme:
+  //case WTA_E_scheme:
+  case WTA_modp_scheme:
     break;
   case pt_scheme: