[Rivet-svn] r2522 - in trunk: . include/Rivet/Math src/Analyses

[blackhole at projects.hepforge.org]

Author: buckley
Date: Wed Jun 23 16:21:14 2010
New Revision: 2522

Log:
Numerical precision improvement in calculating 4-vector masses, and a related improvement in handling jet masses in the base jet analysis

Modified:
   trunk/ChangeLog
   trunk/include/Rivet/Math/Vector4.hh
   trunk/src/Analyses/MC_JetAnalysis.cc

Modified: trunk/ChangeLog
==============================================================================
--- trunk/ChangeLog	Wed Jun 23 15:04:01 2010	(r2521)
+++ trunk/ChangeLog	Wed Jun 23 16:21:14 2010	(r2522)
@@ -1,8 +1,16 @@
 2010-06-23  Andy Buckley  <andy at insectnation.org>
 
+	* Adding protection and warning about numerical precision issues
+	in jet mass calculation/histogramming to the MC_JetAnalysis
+	analysis.
+
+	* Numerical precision improvement in calculation of
+	Vector4::mass2.
+
 	* Adding relative scale ratio plot flag to compare-histos
 
-	* Extended command completion to rivet-config, compare-histos, and make-plots.
+	* Extended command completion to rivet-config, compare-histos, and
+	make-plots.
 
 	* Providing protected log messaging macros,
 	MSG_{TRACE,DEBUG,INFO,WARNING,ERROR} cf. Athena.

Modified: trunk/include/Rivet/Math/Vector4.hh
==============================================================================
--- trunk/include/Rivet/Math/Vector4.hh	Wed Jun 23 15:04:01 2010	(r2521)
+++ trunk/include/Rivet/Math/Vector4.hh	Wed Jun 23 16:21:14 2010	(r2522)
@@ -57,7 +57,8 @@
     FourVector& setZ(const double z) { set(3, z); return *this; }
 
     double invariant() const {
-      return t()*t() - x()*x() - y()*y() - z()*z();
+      // Done this way for numerical precision
+      return (t() + z())*(t() - z()) - x()*x() - y()*y();
     }
 
     double angle(const FourVector& v) const {

Modified: trunk/src/Analyses/MC_JetAnalysis.cc
==============================================================================
--- trunk/src/Analyses/MC_JetAnalysis.cc	Wed Jun 23 15:04:01 2010	(r2521)
+++ trunk/src/Analyses/MC_JetAnalysis.cc	Wed Jun 23 16:21:14 2010	(r2522)
@@ -119,7 +119,14 @@
     for (size_t i=0; i<m_njet; ++i) {
       if (jets.size()<i+1) continue;
       _h_pT_jet[i]->fill(jets[i].momentum().pT()/GeV, weight);
-      _h_mass_jet[i]->fill(jets[i].momentum().mass()/GeV, weight);
+      // Check for numerical precision issues with jet masses
+      double m2_i = jets[i].momentum().mass2();
+      if (m2_i < 0) {
+        getLog() << Log::WARNING << "Jet mass2 is negative: " << m2_i << " GeV^2. "
+                 << "Truncating to 0.0, assuming numerical precision is to blame." << endl;
+        m2_i = 0.0;
+      }
+      _h_mass_jet[i]->fill(sqrt(m2_i)/GeV, weight);
       _h_eta_jet[i]->fill(jets[i].momentum().eta(), weight);
       _h_rap_jet[i]->fill(jets[i].momentum().rapidity(), weight);
       // cout << "Jet mass [" << i+1 << "] = " << jets[i].momentum().mass()/GeV << " GeV" << endl;