Changeset 35cdc46 in git for external/fastjet/tools/Subtractor.cc

Timestamp:

Sep 3, 2014, 3:18:54 PM (10 years ago)

Author:

Pavel Demin <demin@…>

Branches:

ImprovedOutputFile, Timing, dual_readout, llp, master

Children:

be2222c

Parents:

5b5a56b

Message:

upgrade FastJet to version 3.1.0-beta.1, upgrade Nsubjettiness to version 2.1.0, add SoftKiller version 1.0.0

File:

• external/fastjet/tools/Subtractor.cc (modified) (5 diffs)

external/fastjet/tools/Subtractor.cc

@@ -1 @@
-//STARTHEADER
-// $Id$
+//FJSTARTHEADER
+// $Id: Subtractor.cc 3594 2014-08-12 15:25:11Z soyez $
 //
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory 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/tools/Subtractor.hh"
@@ -38 +40 @@
 
 
+//----------------------------------------------------------------------
+// ctor
 Subtractor::Subtractor(double rho) : _bge(0), _rho(rho) {
   assert(_rho>0.0);
+  set_defaults();
 }
 
+//----------------------------------------------------------------------
+void Subtractor::set_defaults(){
+  _use_rho_m = false; // likely to change in future releases!!
+  _safe_mass = false; // likely to change in future releases!!
+
+  _sel_known_vertex = Selector();
+  _sel_leading_vertex = Selector();
+}
+
+//----------------------------------------------------------------------
+// perform the subtraction of a given jet
 PseudoJet Subtractor::result(const PseudoJet & jet) const {
   if (!jet.has_area()){
     throw Error("Trying to subtract a jet without area support");
   }
-  
+
+  PseudoJet known_lv, known_pu;
+  PseudoJet unknown = jet;
+  if (_sel_known_vertex.worker()){
+    // separate the jet constituents in 3 groups:
+    //   unknown vertex
+    //   known vertex, leading vertex
+    //   known vertex, non-leading vertex (PU)
+    vector<PseudoJet> constits_unknown, constits_known;
+    _sel_known_vertex.sift(jet.constituents(), 
+                           constits_known,
+                           constits_unknown);
+    vector<PseudoJet> constits_known_lv, constits_known_pu;
+    _sel_leading_vertex.sift(constits_known,
+                             constits_known_lv,
+                             constits_known_pu);
+
+    // For the parts related to the known vertices (LV or PU), we just
+    // sum the 4-momenta. For the unknown part, we assign it the full
+    // jet area.
+    known_lv = (constits_known_lv.size()!=0) 
+      ? SelectorIdentity().sum(constits_known_lv) : 0.0*jet;
+    known_pu = (constits_known_pu.size()!=0) 
+      ? SelectorIdentity().sum(constits_known_pu) : 0.0*jet;
+    if (constits_unknown.size()==0){
+      // no need for any form of subtraction!
+      PseudoJet subtracted_jet = jet;
+      subtracted_jet.reset_momentum(known_lv);
+      return subtracted_jet;
+    }
+    unknown = jet; // that keeps all info including area
+    unknown.reset_momentum(SelectorIdentity().sum(constits_unknown));
+  } else {
+    known_lv = jet; // ensures correct rap-phi!
+    known_lv *= 0.0;
+    known_pu = known_lv;
+  }
+
+  // prepare for the subtraction and compute the 4-vector to be
+  // subtracted
+  PseudoJet subtracted_jet = jet;
+  PseudoJet to_subtract = known_pu + _amount_to_subtract(unknown);
+
+  // sanity check for the transverse momentum
+  if (to_subtract.pt2() < jet.pt2() ) { 
+    // this subtraction should retain the jet's structural
+    // information
+    subtracted_jet -= to_subtract;
+  } else { 
+    // this sets the jet's momentum while maintaining all of the jet's
+    // structural information
+    subtracted_jet.reset_momentum(known_lv);
+    return subtracted_jet;
+  }
+
+  // make sure that in the end the pt is at least the one known to
+  // come from the leading vertex
+  if (subtracted_jet.pt2() < known_lv.pt2()){
+    subtracted_jet.reset_momentum(known_lv);
+    return subtracted_jet;
+  }
+
+  // sanity check for the mass (if needed)
+  if ((_safe_mass) && (subtracted_jet.m2() < known_lv.m2())){
+    // in this case, we keep pt and phi as obtained from the
+    // subtraction above and take rap and m from the part that comes
+    // from the leading vertex (or the original jet if nothing comes
+    // from the leading vertex)
+    subtracted_jet.reset_momentum(PtYPhiM(subtracted_jet.pt(),
+                                          known_lv.rap(),
+                                          subtracted_jet.phi(),
+                                          known_lv.m()));
+  }
+
+  return subtracted_jet;
+}
+
+//----------------------------------------------------------------------
+std::string Subtractor::description() const{
+  if (_bge != 0) {
+    string desc = "Subtractor that uses the following background estimator to determine rho: "+_bge->description();
+    if (use_rho_m()) desc += "; including the rho_m correction";
+    if (safe_mass()) desc += "; including mass safety tests";
+    if (_sel_known_vertex.worker()){
+      desc += "; using known vertex selection: "+_sel_known_vertex.description()+" and leading vertex selection: "+_sel_leading_vertex.description();
+    }
+    return desc;
+  } else if (_rho != _invalid_rho) {
+    ostringstream ostr;
+    ostr << "Subtractor that uses a fixed value of rho = " << _rho;
+    return ostr.str();
+  } else {
+    return "Uninitialised subtractor";
+  }
+}
+
+//----------------------------------------------------------------------
+// compute the 4-vector that should be subtracted from the given
+// jet
+PseudoJet Subtractor::_amount_to_subtract(const PseudoJet &jet) const{
+  // the "transverse momentum" part
   double rho;
   if (_bge != 0) {
@@ -56 +172 @@
   }
 
-  PseudoJet subtracted_jet = jet;
-  PseudoJet area4vect = jet.area_4vector();
-  // sanity check
-  if (rho*area4vect.perp() < jet.perp() ) { 
-    // this subtraction should retain the jet's structural
-    // information
-    subtracted_jet -= rho*area4vect;
-  } else { 
-    // this sets the jet's momentum to zero while
-    // maintaining all of the jet's structural information
-    subtracted_jet *= 0;
+  PseudoJet area = jet.area_4vector();
+  PseudoJet to_subtract = rho*area;
+
+  double const rho_m_warning_threshold = 1e-5;
+
+  // add an optional contribution from the unknown particles masses
+  if (_use_rho_m){
+    assert(_bge != 0); // test done in "set_use_rho_m()"
+    to_subtract += _bge->rho_m(jet) * PseudoJet(0.0, 0.0, area.pz(), area.E());
+  } else if (_bge && 
+             _bge->has_rho_m() && 
+             _bge->rho_m(jet) > rho_m_warning_threshold * rho) {
+    _unused_rho_m_warning.warn("Background estimator indicates non-zero rho_m, but use_rho_m()==false in subtractor; consider calling set_use_rho_m(true) to include the rho_m information");
   }
-  return subtracted_jet;
+
+  return to_subtract;
 }
 
-//----------------------------------------------------------------------
-std::string Subtractor::description() const{
-  if (_bge != 0) {
-    return "Subtractor that uses the following background estimator to determine rho: "+_bge->description();
-  } else if (_rho != _invalid_rho) {
-    ostringstream ostr;
-    ostr << "Subtractor that uses a fixed value of rho = " << _rho;
-    return ostr.str();
-  } else {
-    return "Uninitialised subtractor";
-  }
-}
 
 FASTJET_END_NAMESPACE