Merge pull request #1752 from Martin-Grunewald/FixJetAreaCalculation

Reco fixes -- Add protection against sqrt(negative value) due to numerical rounding effects.
cms-sw · Dec 10, 2013 · a4b6ffa · a4b6ffa
2 parents 099fd5f + 203d07e
commit a4b6ffa
Showing 1 changed file with 4 additions and 4 deletions.
diff --git a/RecoJets/JetProducers/src/VirtualJetProducerHelper.cc b/RecoJets/JetProducers/src/VirtualJetProducerHelper.cc
@@ -1,13 +1,13 @@
 #include "RecoJets/JetProducers/interface/VirtualJetProducerHelper.h"
 
 #include <cmath>
-
+#include <algorithm>
 
 double reco::helper::VirtualJetProducerHelper::intersection(double r12) 
 {
   if (r12 == 0)         return M_PI;
   if (r12 >= 2)         return 0;
-  return 2 * acos(r12/2) - 0.5*r12*sqrt(4 - r12*r12);
+  return 2 * acos(r12/2) - 0.5*r12*sqrt(std::max(0.0 , 4 - r12*r12));
 }
 
 double reco::helper::VirtualJetProducerHelper::intersection(double r12, double r23, double r13) 
@@ -16,7 +16,7 @@ double reco::helper::VirtualJetProducerHelper::intersection(double r12, double r
   const double          r12_2   = r12*r12;
   const double          r13_2   = r13*r13;
   const double          temp    = (r12_2 + r13_2 - r23*r23);
-  const double          T2      = 4*r12_2*r13_2 - temp*temp;
+  const double          T2      = std::max(0.0 , 4*r12_2*r13_2 - temp*temp);
   const double          common  = 0.5*( intersection(r12) + intersection(r13) + intersection(r23) - M_PI + sqrt(T2)/2 );
   return common;
 }
@@ -27,7 +27,7 @@ double reco::helper::VirtualJetProducerHelper::intersection(double r12, double r
   const double          r12_2   = r12*r12;
   const double          r13_2   = r13*r13;
   const double          temp    = (r12_2 + r13_2 - r23*r23);
-  const double          T2	= 4*r12_2*r13_2 - temp*temp;
+  const double          T2	= std::max(0.0 , 4*r12_2*r13_2 - temp*temp);
   const double          common  = 0.5*( a12 + a13 + a23 - M_PI + sqrt(T2)/2 );
   return common;
 }