From b596f5c6089aec508614e6a77a7a9797549995e6 Mon Sep 17 00:00:00 2001
From: Alexander Gavrilyuk <alexander.gavrilyuk@cern.ch>
Date: Wed, 20 Nov 2019 12:16:30 +0100
Subject: [PATCH 1/4] Correct plotted formula of poisson significance

Compare https://gitlab.cern.ch/atlas-caf/CAFCore/blob/master/QFramework/Root/TQHistogramUtils.cxx#L3532
---
 QFramework/Root/TQHistogramUtils.cxx | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/QFramework/Root/TQHistogramUtils.cxx b/QFramework/Root/TQHistogramUtils.cxx
index 7ae456032..2dd6dfffd 100644
--- a/QFramework/Root/TQHistogramUtils.cxx
+++ b/QFramework/Root/TQHistogramUtils.cxx
@@ -6500,7 +6500,7 @@ TString TQHistogramUtils::getFOMTitleROOT(TQHistogramUtils::FOM fom){
   switch(fom){
   case TQHistogramUtils::kSoSqB: return "s/#sqrt{b}";
   case TQHistogramUtils::kSoSqBpdB: return "s/#sqrt{b+db^{2}}";
-  case TQHistogramUtils::kPoisson: return "#sqrt{2(s+b)ln(1+s/b)-s}";
+  case TQHistogramUtils::kPoisson: return "#sqrt{2((s+b)ln(1+s/b)-s)}";
   case TQHistogramUtils::kSoB: return "s/b";
   case TQHistogramUtils::kSoSqSpB: return "s/#sqrt{s+b}";
   default: return "<unknown>";
@@ -6515,7 +6515,7 @@ TString TQHistogramUtils::getFOMTitleLaTeX(TQHistogramUtils::FOM fom){
   switch(fom){
   case TQHistogramUtils::kSoSqB: return "$s/\\sqrt{b}$";
   case TQHistogramUtils::kSoSqBpdB: return "$s/\\sqrt{b+db^{2}}$";
-  case TQHistogramUtils::kPoisson: return "$\\sqrt{2(s+b)\\ln(1+s/b)-s}$";
+  case TQHistogramUtils::kPoisson: return "$\\sqrt{2((s+b)\\ln(1+s/b)-s)}$";
   case TQHistogramUtils::kSoB: return "$s/b$";
   case TQHistogramUtils::kSoSqSpB: return "$s/\\sqrt{s+b}$";
   default: return "<unknown>";
-- 
GitLab


From 8b9d7425368cc4c7b9dd84de56a0786a419f476f Mon Sep 17 00:00:00 2001
From: Alexander Gavrilyuk <alexander.gavrilyuk@cern.ch>
Date: Wed, 20 Nov 2019 16:56:39 +0100
Subject: [PATCH 2/4] Fix error in PoissonWithError formula and make it
 available as a figure-of-merit for a subplot

---
 QFramework/QFramework/TQHistogramUtils.h |  1 +
 QFramework/Root/TQDefaultPlotter.cxx     |  2 +-
 QFramework/Root/TQHWWPlotter.cxx         |  2 +-
 QFramework/Root/TQHistogramUtils.cxx     | 26 ++++++++++++++++++------
 4 files changed, 23 insertions(+), 8 deletions(-)

diff --git a/QFramework/QFramework/TQHistogramUtils.h b/QFramework/QFramework/TQHistogramUtils.h
index 57922f307..64c629aed 100644
--- a/QFramework/QFramework/TQHistogramUtils.h
+++ b/QFramework/QFramework/TQHistogramUtils.h
@@ -35,6 +35,7 @@ namespace TQHistogramUtils {
     kSoSqSpB,
     kSoB,
     kPoisson,
+    kPoissonWErr,
     kUndefined
   };
 
diff --git a/QFramework/Root/TQDefaultPlotter.cxx b/QFramework/Root/TQDefaultPlotter.cxx
index 98931d484..41823fbc6 100644
--- a/QFramework/Root/TQDefaultPlotter.cxx
+++ b/QFramework/Root/TQDefaultPlotter.cxx
@@ -564,7 +564,7 @@ void TQDefaultPlotter::drawSub_CutOptimization(TQTaggable& tags){
   }
   if(!hSig) return;
 
-  //@tag:optScan.FOMmode: figure of merit to be used. currently available are: s/sqrt(s+b),s/b,poisson,s/sqrt(b),s/sqrt(b+db2)
+  //@tag:optScan.FOMmode: figure of merit to be used. currently available are: s/sqrt(s+b),s/b,poisson,poissonwerr,s/sqrt(b),s/sqrt(b+db2)
   //@tag:optScan.FOMbbb: evaluate the figure-of-merit bin-by-bin instead of integrated left and right (default:false)
   //@tag:style.FOMmode: deprecated, use optScan.FOMmode
   //@tag:style.FOMbbb: deprecated, use optScan.FOMbbb
diff --git a/QFramework/Root/TQHWWPlotter.cxx b/QFramework/Root/TQHWWPlotter.cxx
index 9d1360d01..d1dd39c71 100644
--- a/QFramework/Root/TQHWWPlotter.cxx
+++ b/QFramework/Root/TQHWWPlotter.cxx
@@ -2411,7 +2411,7 @@ void TQHWWPlotter::drawOptScan(TQTaggable& tags){
   TH1* hFOMr = 0;
   TH1* hFOM = 0;
 
-  //@tag:optScan.FOMmode: figure of merit to be used. currently available are: s/sqrt(s+b),s/b,poisson,s/sqrt(b),s/sqrt(b+db2)
+  //@tag:optScan.FOMmode: figure of merit to be used. currently available are: s/sqrt(s+b),s/b,poisson,poissonwerr,s/sqrt(b),s/sqrt(b+db2)
   //@tag:optScan.FOMbbb: evaluate the figure-of-merit bin-by-bin instead of integrated left and right (default:false)
   //@tag:style.FOMmode: deprecated, use optScan.FOMmode
   //@tag:style.FOMbbb: deprecated, use optScan.FOMbbb
diff --git a/QFramework/Root/TQHistogramUtils.cxx b/QFramework/Root/TQHistogramUtils.cxx
index 2dd6dfffd..74a7956cf 100644
--- a/QFramework/Root/TQHistogramUtils.cxx
+++ b/QFramework/Root/TQHistogramUtils.cxx
@@ -3561,26 +3561,32 @@ double TQHistogramUtils::getPoissonError(double b, double s, double db, double d
 
   double logTerm = TMath::Log(1 + s / b);
 
-  // dP/ds, dP/db
-  double dPds = ( logTerm - 1 +     (b + s) / (b *     (1 + s / b)) ) / P;
-  double dPdb = ( logTerm     - s * (b + s) / (b * b * (1 + s / b)) ) / P;
+  // dP/ds = ( logTerm - 1 +     (b + s) / (b *     (1 + s / b)) ) / P
+  // dP/db = ( logTerm     - s * (b + s) / (b * b * (1 + s / b)) ) / P
+  double dPds = ( logTerm       ) / P;
+  double dPdb = ( logTerm - s/b ) / P;
 
   return TMath::Sqrt( TMath::Power((dPds * ds), 2) + TMath::Power((dPdb * db), 2) );
 }
 
 //__________________________________________________________________________________|___________
 
-double TQHistogramUtils::getPoissonWithError(double b, double berr, double s) {
+double TQHistogramUtils::getPoissonWithError(double b, double db, double s) {
   // Returns the Poisson significance according to <b> expected background events
   // and <s> expected signal events. Returns 0 if any term becomes invalid
   //
-  // Poisson significance = TMath::Sqrt(2. * ((s + b) * TMath::Log(1. + s / b) - s))
+  // Poisson significance = TMath::Sqrt(2. *(
+  //         (s+b) * TMath::Log((s+b)*(b+db*db)/(b*b+(s+b)*db*db) ) - 
+  //   b*b/(db*db) * TMath::Log(1. +  (db*db*s)/(     b*(b+db*db)))    ) )
+  // See slide 9 https://indico.cern.ch/event/768968/contributions/3195413/attachments/1743878/2822499/statsCommitteeAtlasWeeklyOct2018.pdf
 
   if (b <= 0. || s < 0.) {
     return std::numeric_limits<double>::quiet_NaN();
   }
 
-  double sqrtInput = 2. * ((s + b) * TMath::Log( (s + b) * (b + berr*berr) / ( b*b + (s + b) * berr * berr) ) - b*b / berr*berr * TMath::Log(1. + berr*berr*s / (b * (b + berr*berr) ) ) );
+  double sqrtInput = 2. * (
+          (s + b) * TMath::Log( (s + b) * (b + db*db) / ( b*b + (s + b) * db * db) ) -
+    b*b / (db*db) * TMath::Log( 1. +          db*db*s / (b * (b + db*db) )    )        );
   if (sqrtInput < 0.) {
     return std::numeric_limits<double>::quiet_NaN();
   }
@@ -6459,6 +6465,8 @@ double TQHistogramUtils::getFOM(TQHistogramUtils::FOM fom, double b, double berr
     return s/sqrt(b + berr*berr);
   case kPoisson: 
     return TQHistogramUtils::getPoisson(b,s);
+  case kPoissonWErr: 
+    return TQHistogramUtils::getPoissonWithError(b,berr,s);
   case kSoB: 
     return s/b;
   case kSoSqSpB: 
@@ -6484,6 +6492,9 @@ TQHistogramUtils::FOM TQHistogramUtils::readFOM(TString fom){
   if(TQStringUtils::matches(fom,"poisson")){
     return TQHistogramUtils::kPoisson;
   }
+  if(TQStringUtils::matches(fom,"poissonwerr")){
+    return TQHistogramUtils::kPoissonWErr;
+  }
   if(TQStringUtils::matches(fom,"s/sqrt(b)")){
     return TQHistogramUtils::kSoSqB;
   }
@@ -6501,6 +6512,8 @@ TString TQHistogramUtils::getFOMTitleROOT(TQHistogramUtils::FOM fom){
   case TQHistogramUtils::kSoSqB: return "s/#sqrt{b}";
   case TQHistogramUtils::kSoSqBpdB: return "s/#sqrt{b+db^{2}}";
   case TQHistogramUtils::kPoisson: return "#sqrt{2((s+b)ln(1+s/b)-s)}";
+  //case TQHistogramUtils::kPoissonWErr: return "#sqrt{2{(s+b)ln[#frac{(s+b)(b+#sigma^{2})}{b^{2}+(s+b)#sigma^{2}}] - #frac{b^{2}}{#sigma^{2}}ln[1+#frac{#sigma^{2}s}{b(b+#sigma^{2})}]}}";
+  case TQHistogramUtils::kPoissonWErr: return "#sqrt{2{(s+b)ln[A]-#frac{b^{2}}{#sigma^{2}}ln[B]}}";
   case TQHistogramUtils::kSoB: return "s/b";
   case TQHistogramUtils::kSoSqSpB: return "s/#sqrt{s+b}";
   default: return "<unknown>";
@@ -6516,6 +6529,7 @@ TString TQHistogramUtils::getFOMTitleLaTeX(TQHistogramUtils::FOM fom){
   case TQHistogramUtils::kSoSqB: return "$s/\\sqrt{b}$";
   case TQHistogramUtils::kSoSqBpdB: return "$s/\\sqrt{b+db^{2}}$";
   case TQHistogramUtils::kPoisson: return "$\\sqrt{2((s+b)\\ln(1+s/b)-s)}$";
+  case TQHistogramUtils::kPoissonWErr: return "$\\sqrt{2{(s+b)\\ln[\\frac{(s+b)(b+\\sigma^{2})}{b^{2}+(s+b)\\sigma^{2}}] - \\frac{b^{2}}{\\sigma^{2}}\\ln[1+\\frac{\\sigma^{2}s}{b(b+\\sigma^{2})}]}}$";
   case TQHistogramUtils::kSoB: return "$s/b$";
   case TQHistogramUtils::kSoSqSpB: return "$s/\\sqrt{s+b}$";
   default: return "<unknown>";
-- 
GitLab


From d2037940cb6895bf12dda0878b42dc465a4e58b3 Mon Sep 17 00:00:00 2001
From: Alexander Gavrilyuk <alexander.gavrilyuk@cern.ch>
Date: Wed, 20 Nov 2019 19:57:08 +0100
Subject: [PATCH 3/4] Add error calculating method for PoissonWithErrors and
 switch cutflow poisson significances to include uncertainties

---
 QFramework/QFramework/TQHistogramUtils.h |  3 ++-
 QFramework/Root/TQCutflowPrinter.cxx     |  6 ++---
 QFramework/Root/TQHistogramUtils.cxx     | 31 ++++++++++++++++++++++++
 3 files changed, 36 insertions(+), 4 deletions(-)

diff --git a/QFramework/QFramework/TQHistogramUtils.h b/QFramework/QFramework/TQHistogramUtils.h
index 64c629aed..5cad61032 100644
--- a/QFramework/QFramework/TQHistogramUtils.h
+++ b/QFramework/QFramework/TQHistogramUtils.h
@@ -105,7 +105,8 @@ namespace TQHistogramUtils {
   double getPoisson(double b, double s);
   double getPoissonError(double b, double s, double db, double ds);
 
-  double getPoissonWithError(double b, double berr, double s);
+  double getPoissonWithError(double b, double db, double s);
+  double getPoissonWithErrorError(double b, double s, double db, double ds);
 
   double getSoverSqrtB(double b, double s);
   double getSoverSqrtSplusB(double b, double s);
diff --git a/QFramework/Root/TQCutflowPrinter.cxx b/QFramework/Root/TQCutflowPrinter.cxx
index 465c6d7cb..955d59a99 100644
--- a/QFramework/Root/TQCutflowPrinter.cxx
+++ b/QFramework/Root/TQCutflowPrinter.cxx
@@ -365,10 +365,10 @@ bool TQCutflowPrinter::getValues(TQTaggable& tags, TQNamedTaggable* cut, TQNamed
       }
       DEBUG("making significance");
       showNumbers = true;
-      number = TQHistogramUtils::getPoisson(number2, number1);
-      statErr = TQHistogramUtils::getPoissonError(number2, number1, statErr2, statErr1);
+      number = TQHistogramUtils::getPoissonWithError(number2, statErr2, number1);
+      statErr = TQHistogramUtils::getPoissonWithErrorError(number2, number1, statErr2, statErr1);
       if (includeScaleUncertainty) {
-        theoSysErr = TQHistogramUtils::getPoissonError(number2, number1, theoSysErr2, theoSysErr1);
+        theoSysErr = TQHistogramUtils::getPoissonWithErrorError(number2, number1, theoSysErr2, theoSysErr1);
       }
       raw = -1;
       if (includeScaleUncertainty) {
diff --git a/QFramework/Root/TQHistogramUtils.cxx b/QFramework/Root/TQHistogramUtils.cxx
index 74a7956cf..4348f965c 100644
--- a/QFramework/Root/TQHistogramUtils.cxx
+++ b/QFramework/Root/TQHistogramUtils.cxx
@@ -3594,6 +3594,37 @@ double TQHistogramUtils::getPoissonWithError(double b, double db, double s) {
   return TMath::Sqrt(sqrtInput);
 }
 
+//__________________________________________________________________________________|___________
+
+double TQHistogramUtils::getPoissonWithErrorError(double b, double s, double db, double ds) {
+  // Returns the error on the Poisson significance according to <b> expected background events
+  // and <s> expected signal events. Returns 0 if any term becomes invalid
+  //
+  // Poisson significance = TMath::Sqrt(2. *(
+  //         (s+b) * TMath::Log((s+b)*(b+db*db)/(b*b+(s+b)*db*db) ) - 
+  //   b*b/(db*db) * TMath::Log(1. +  (db*db*s)/(     b*(b+db*db)))    ) )
+  // The error is propagated from errors in b and s according to the most general formula
+  // for independent error propagation.
+
+  if (b <= 0. || s < 0.) {
+    return std::numeric_limits<double>::quiet_NaN();
+  }
+
+  double LogA=TMath::Log(  (s+b)*(b + db*db) /((s+b)*db*db + b*b) );
+  double LogB=TMath::Log( ((s+b)*db*db + b*b)/(b*(b+db*db)) );
+  // negative radicand check
+  double sqrtInput = 2. * ( ((s+b)*db*db)*LogA - b*b*LogB );
+  if (sqrtInput < 0.) {
+    return std::numeric_limits<double>::quiet_NaN();
+  }
+  // Poisson significance
+  //double P = TMath::Sqrt(sqrtInput);
+
+  double dPds = db*LogA / TMath::Sqrt(sqrtInput);
+  double dPdb = ( (b+db*db)*(db*db)*LogA - 2*b*(b+db*db)*LogB + s*db*db )  / (db*(b+db*db)*TMath::Sqrt(sqrtInput));
+
+  return TMath::Sqrt( TMath::Power((dPds * ds), 2) + TMath::Power((dPdb * db), 2) );
+}
 
 //__________________________________________________________________________________|___________
 
-- 
GitLab


From 00bb2731f5407b72227183cda5d52d6c10b83fb7 Mon Sep 17 00:00:00 2001
From: Alexander Gavrilyuk <alexander.gavrilyuk@cern.ch>
Date: Sun, 17 May 2020 15:21:00 +0200
Subject: [PATCH 4/4] Apply asymptotic workaround against division by zero

---
 QFramework/Root/TQHistogramUtils.cxx | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/QFramework/Root/TQHistogramUtils.cxx b/QFramework/Root/TQHistogramUtils.cxx
index 4348f965c..daa9ceab3 100644
--- a/QFramework/Root/TQHistogramUtils.cxx
+++ b/QFramework/Root/TQHistogramUtils.cxx
@@ -3580,9 +3580,11 @@ double TQHistogramUtils::getPoissonWithError(double b, double db, double s) {
   //   b*b/(db*db) * TMath::Log(1. +  (db*db*s)/(     b*(b+db*db)))    ) )
   // See slide 9 https://indico.cern.ch/event/768968/contributions/3195413/attachments/1743878/2822499/statsCommitteeAtlasWeeklyOct2018.pdf
 
-  if (b <= 0. || s < 0.) {
+  if (b <= 0.) {
     return std::numeric_limits<double>::quiet_NaN();
   }
+  
+  if (std::abs(s)<1e-5) return 0;
 
   double sqrtInput = 2. * (
           (s + b) * TMath::Log( (s + b) * (b + db*db) / ( b*b + (s + b) * db * db) ) -
@@ -3591,7 +3593,7 @@ double TQHistogramUtils::getPoissonWithError(double b, double db, double s) {
     return std::numeric_limits<double>::quiet_NaN();
   }
  
-  return TMath::Sqrt(sqrtInput);
+  return TMath::Sqrt(sqrtInput) * (s>0 ? 1: -1);
 }
 
 //__________________________________________________________________________________|___________
@@ -3606,9 +3608,11 @@ double TQHistogramUtils::getPoissonWithErrorError(double b, double s, double db,
   // The error is propagated from errors in b and s according to the most general formula
   // for independent error propagation.
 
-  if (b <= 0. || s < 0.) {
+  if (b <= 0.) {
     return std::numeric_limits<double>::quiet_NaN();
   }
+  
+  if (std::abs(s)<1e-5) return 0;
 
   double LogA=TMath::Log(  (s+b)*(b + db*db) /((s+b)*db*db + b*b) );
   double LogB=TMath::Log( ((s+b)*db*db + b*b)/(b*(b+db*db)) );
@@ -6544,7 +6548,7 @@ TString TQHistogramUtils::getFOMTitleROOT(TQHistogramUtils::FOM fom){
   case TQHistogramUtils::kSoSqBpdB: return "s/#sqrt{b+db^{2}}";
   case TQHistogramUtils::kPoisson: return "#sqrt{2((s+b)ln(1+s/b)-s)}";
   //case TQHistogramUtils::kPoissonWErr: return "#sqrt{2{(s+b)ln[#frac{(s+b)(b+#sigma^{2})}{b^{2}+(s+b)#sigma^{2}}] - #frac{b^{2}}{#sigma^{2}}ln[1+#frac{#sigma^{2}s}{b(b+#sigma^{2})}]}}";
-  case TQHistogramUtils::kPoissonWErr: return "#sqrt{2{(s+b)ln[A]-#frac{b^{2}}{#sigma^{2}}ln[B]}}";
+  case TQHistogramUtils::kPoissonWErr: return "Poisson Significance";
   case TQHistogramUtils::kSoB: return "s/b";
   case TQHistogramUtils::kSoSqSpB: return "s/#sqrt{s+b}";
   default: return "<unknown>";
-- 
GitLab