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