Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add protection against sqrt(negative value) due to numerical rounding ef... #1752

Merged
merged 1 commit into from Dec 10, 2013
Merged

Add protection against sqrt(negative value) due to numerical rounding ef... #1752

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
8 changes: 4 additions & 4 deletions RecoJets/JetProducers/src/VirtualJetProducerHelper.cc
View file
Open in desktop
@@ -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)
Expand All @@ -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;
}
Expand All @@ -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;
}
Expand Down