diff --git a/QFramework/QFramework/TQHistogramUtils.h b/QFramework/QFramework/TQHistogramUtils.h index d7626009f62f5071ef9c20b1c5d31ecf186b2e53..06f223214d56457ee20419b921703b5cddaeed96 100644 --- a/QFramework/QFramework/TQHistogramUtils.h +++ b/QFramework/QFramework/TQHistogramUtils.h @@ -35,6 +35,7 @@ namespace TQHistogramUtils { kSoSqSpB, kSoB, kPoisson, + kPoissonWErr, kUndefined }; @@ -104,9 +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 getPoissonPoissonWithError(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 465c6d7cb5b8e9dc80f11b852037e04bd9e3100c..955d59a9929d95592c03b69694ce62400b61124b 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/TQDefaultPlotter.cxx b/QFramework/Root/TQDefaultPlotter.cxx index 222d9291ee56e8825eecb151cfba2d388c054ded..109d0063144110d0daf441477d4772470f9d8c94 100644 --- a/QFramework/Root/TQDefaultPlotter.cxx +++ b/QFramework/Root/TQDefaultPlotter.cxx @@ -631,7 +631,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 34a558d7b304cc8331f6db750c4aa800a85c8a85..c77d4444197a8a9975ed405a5fb6aa8c0766a47b 100644 --- a/QFramework/Root/TQHWWPlotter.cxx +++ b/QFramework/Root/TQHWWPlotter.cxx @@ -2414,7 +2414,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 cab2fa6a812164fe4ee20e45600f44aabfea8520..c26da8613c1cbf0485f64d8f4cce208ca16190be 100644 --- a/QFramework/Root/TQHistogramUtils.cxx +++ b/QFramework/Root/TQHistogramUtils.cxx @@ -3780,53 +3780,73 @@ 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) { - // Returns the Poisson significance according to <b> expected background events - // and <s> expected signal events. Returns 0 if any term becomes invalid +double TQHistogramUtils::getPoissonWithError(double b, double db, double s) { + // Returns the significance corresponding to a Poisson-Poisson model with asymptotic formula. + // Returns 0 if any term becomes invalid // - // Poisson significance = TMath::Sqrt(2. * ((s + b) * TMath::Log(1. + s / b) - s)) + // Poisson-Poisson asymptotic 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.) { + if (b <= 0. || db < 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 + 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(); } - - return TMath::Sqrt(sqrtInput); + + return TMath::Sqrt(sqrtInput) * (s > 0. ? 1. : -1.); } //__________________________________________________________________________________|___________ -double TQHistogramUtils::getPoissonPoissonWithError(double b, double berr, double s) { - // Returns the significance corresponding to a Poisson-Poisson model with asymptotic formula. - // See https://cds.cern.ch/record/2643488/files/ATL-COM-GEN-2018-026.pdf - // Returns 0 if any term becomes invalid +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. || berr < 0.) { - WARN("Invalid input to Poisson-Poisson significance calculation, returning NaN!"); + if (b <= 0. || db < 0.) { return std::numeric_limits<double>::quiet_NaN(); } + + if (std::abs(s) < 1e-5) return 0.; - double n = s+b; - double arg1 = n*TMath::Log(n*(b+TMath::Power(berr, 2))/(TMath::Power(b, 2)+n*TMath::Power(berr, 2))); - double arg2 = TMath::Power(b, 2)/TMath::Power(berr, 2)*TMath::Log(1+TMath::Power(berr, 2)*(n-b)/(b*(b+TMath::Power(berr, 2)))); - double combo = 2*(arg1-arg2); - if (combo < 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)) ); + // negative radicand check + double sqrtInput = 2. * ( ((s+b)*db*db)*LogA - b*b*LogB ); + if (sqrtInput < 0.) { return std::numeric_limits<double>::quiet_NaN(); } - return TMath::Sqrt(combo); + // 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) ); } //__________________________________________________________________________________|___________ @@ -3897,7 +3917,7 @@ double TQHistogramUtils::getSignificance(double b, double s, } else if (sgnfName.CompareTo("asymptoticPoissonPoisson", TString::kIgnoreCase) == 0) { if (sgnfTitle) *sgnfTitle = "Poisson-Poisson model Significance"; - sgnf = getPoissonPoissonWithError(b, berr, s); + sgnf = getPoissonWithError(b, berr, s); } return sgnf; @@ -6754,7 +6774,9 @@ double TQHistogramUtils::getFOM(TQHistogramUtils::FOM fom, double b, double berr return s/sqrt(b + berr*berr); case kPoisson: return TQHistogramUtils::getPoisson(b,s); - case kSoB: + case kPoissonWErr: + return TQHistogramUtils::getPoissonWithError(b,berr,s); + case kSoB: return s/b; case kSoSqSpB: return s/sqrt(s+b); @@ -6779,6 +6801,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; } @@ -6795,7 +6820,9 @@ 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::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 "Significance"; case TQHistogramUtils::kSoB: return "s/b"; case TQHistogramUtils::kSoSqSpB: return "s/#sqrt{s+b}"; default: return "<unknown>"; @@ -6810,7 +6837,8 @@ 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::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>";