Merge pull request cms-sw#47 from anehrkor/CMSSW_9_1_X_tau-pog_tauID-…

…GJangle Protection against arguments to asin and acos outside of allowed range
anehrkor · Apr 25, 2017 · 1a3f320 · 1a3f320
2 parents d14f66a + f947030
commit 1a3f320
Show file tree

Hide file tree

Showing 2 changed files with 14 additions and 6 deletions.
diff --git a/RecoTauTag/RecoTau/plugins/PATTauDiscriminationByMVAIsolationRun2.cc b/RecoTauTag/RecoTau/plugins/PATTauDiscriminationByMVAIsolationRun2.cc
@@ -218,9 +218,13 @@ double PATTauDiscriminationByMVAIsolationRun2::discriminate(const TauRef& tau) c
         double mTau = 1.77682;
         double mAOne = tau->p4().M();
         double pAOneMag = tau->p();
-        double thetaGJmax = std::asin( (std::pow(mTau,2) - std::pow(mAOne,2) ) / ( 2 * mTau * pAOneMag ) );
-        double thetaGJmeasured = std::acos( ( tau->p4().px() * decayDistX + tau->p4().py() * decayDistY + tau->p4().pz() * decayDistZ ) / ( pAOneMag * decayDistMag ) );
-        gjAngleDiff = thetaGJmeasured - thetaGJmax;
+        double argumentThetaGJmax = (std::pow(mTau,2) - std::pow(mAOne,2) ) / ( 2 * mTau * pAOneMag );
+        double argumentThetaGJmeasured = ( tau->p4().px() * decayDistX + tau->p4().py() * decayDistY + tau->p4().pz() * decayDistZ ) / ( pAOneMag * decayDistMag );
+        if ( std::fabs(argumentThetaGJmax) <= 1. && std::fabs(argumentThetaGJmeasured) <= 1. ) {
+            double thetaGJmax = std::asin( argumentThetaGJmax );
+            double thetaGJmeasured = std::acos( argumentThetaGJmeasured );
+            gjAngleDiff = thetaGJmeasured - thetaGJmax;
+        }
     }
 
     if ( mvaOpt_ == kOldDMwoLT || mvaOpt_ == kNewDMwoLT ) {

diff --git a/RecoTauTag/RecoTau/plugins/PFRecoTauDiscriminationByMVAIsolationRun2.cc b/RecoTauTag/RecoTau/plugins/PFRecoTauDiscriminationByMVAIsolationRun2.cc
@@ -229,9 +229,13 @@ double PFRecoTauDiscriminationByMVAIsolationRun2::discriminate(const PFTauRef& t
         double mTau = 1.77682;
         double mAOne = tau->p4().M();
         double pAOneMag = tau->p();
-        double thetaGJmax = std::asin( (std::pow(mTau,2) - std::pow(mAOne,2) ) / ( 2 * mTau * pAOneMag ) );
-        double thetaGJmeasured = std::acos( ( tau->p4().px() * decayDistX + tau->p4().py() * decayDistY + tau->p4().pz() * decayDistZ ) / ( pAOneMag * decayDistMag ) );
-        gjAngleDiff = thetaGJmeasured - thetaGJmax;
+        double argumentThetaGJmax = (std::pow(mTau,2) - std::pow(mAOne,2) ) / ( 2 * mTau * pAOneMag );
+        double argumentThetaGJmeasured = ( tau->p4().px() * decayDistX + tau->p4().py() * decayDistY + tau->p4().pz() * decayDistZ ) / ( pAOneMag * decayDistMag );
+        if ( std::fabs(argumentThetaGJmax) <= 1. && std::fabs(argumentThetaGJmeasured) <= 1. ) {
+            double thetaGJmax = std::asin( argumentThetaGJmax );
+            double thetaGJmeasured = std::acos( argumentThetaGJmeasured );
+            gjAngleDiff = thetaGJmeasured - thetaGJmax;
+        }
     }
 
     if ( mvaOpt_ == kOldDMwoLT || mvaOpt_ == kNewDMwoLT ) {