From 18a74ce728066ffc6fabd275b607abc3e9237790 Mon Sep 17 00:00:00 2001 From: nucleosynthesis Date: Mon, 6 Feb 2023 17:03:03 +0100 Subject: [PATCH] Attempt to depricate CLsb and CLb where possible, code & docs now *both* uses "Pmu" and "Pb" or variants thereof instead of CLs+b and CLb. This is not possible for some method names which are used in RooStats objects so there it remains but with some code comments to explain. --- .../simple-counting/lhcClsByHand.cxx | 14 ++--- docs/part4/usefullinks.md | 2 +- docs/part5/longexercise.md | 2 +- docs/part5/longexerciseanswers.md | 6 +- docs/tutorial2020/exercise.md | 2 +- docs/tutorials-part-2.md | 4 +- interface/HybridNew.h | 2 +- macros/benchmark-simple.sh | 2 +- src/AsymptoticLimits.cc | 42 ++++++------- src/HybridNew.cc | 59 ++++++++++--------- test/makeGrid2DUsingCrab.py | 2 +- test/plotTestStatCLs.py | 4 +- test/plotting/lepLikePlot.cxx | 23 ++++---- test/plotting/qmuPlot.cxx | 48 +++++++-------- test/validation/reference.json | 12 ++-- test/validation/test_HN.py | 2 +- 16 files changed, 113 insertions(+), 113 deletions(-) diff --git a/data/benchmarks/simple-counting/lhcClsByHand.cxx b/data/benchmarks/simple-counting/lhcClsByHand.cxx index 4a5515681a1..d210e351aa0 100644 --- a/data/benchmarks/simple-counting/lhcClsByHand.cxx +++ b/data/benchmarks/simple-counting/lhcClsByHand.cxx @@ -121,7 +121,7 @@ void toyByHand(double r0 = 1.0, int toys=500, int seed=1) { double qData = (nll_S - nll_F); std::cout << "Test statistics on data: " << qData << std::endl; - // Compute CLsb + // Compute Pmu int Nabove = 0, Nbelow = 0; for (int i = 0; i < toys; ++i) { // Reset nuisances @@ -135,10 +135,10 @@ void toyByHand(double r0 = 1.0, int toys=500, int seed=1) { //std::cout << "q toy (" << i << ") = " << qToy << std::endl; if (qToy <= qData) Nbelow++; else Nabove++; } - double CLsb = Nabove/double(Nabove+Nbelow), CLsbError = sqrt(CLsb*(1.-CLsb)/double(Nabove+Nbelow)); - std::cout << "CLsb = " << Nabove << "/" << Nabove+Nbelow << " = " << CLsb << " +/- " << CLsbError << std::endl; + double Pmu = Nabove/double(Nabove+Nbelow), PmuError = sqrt(Pmu*(1.-Pmu)/double(Nabove+Nbelow)); + std::cout << "Pmu = " << Nabove << "/" << Nabove+Nbelow << " = " << Pmu << " +/- " << PmuError << std::endl; - // Compute CLb + // Compute 1-Pb int Nabove2 = 0, Nbelow2 = 0; for (int i = 0; i < toys/4; ++i) { // Reset nuisances @@ -153,10 +153,10 @@ void toyByHand(double r0 = 1.0, int toys=500, int seed=1) { //std::cout << "q toy (" << i << ") = " << qToy << std::endl; if (qToy <= qData) Nbelow2++; else Nabove2++; } - double CLb = Nabove2/double(Nabove2+Nbelow2), CLbError = sqrt(CLb*(1.-CLb)/double(Nabove2+Nbelow2)); - std::cout << "CLb = " << Nabove2 << "/" << Nabove2+Nbelow2 << " = " << CLb << " +/- " << CLbError << std::endl; + double OnemPb = Nabove2/double(Nabove2+Nbelow2), OnemPbError = sqrt(*(1.-OnemPb)/double(Nabove2+Nbelow2)); + std::cout << "1-Pb = " << Nabove2 << "/" << Nabove2+Nbelow2 << " = " << OnemPb << " +/- " << OnemPbError << std::endl; - std::cout << "CLs = "<< CLsb/CLb << " +/- " << hypot(CLsbError/CLb, CLsb*CLbError/CLb/CLb) << std::endl; + std::cout << "CLs = "<< Pmu/OnemPg << " +/- " << hypot(PmuError/OnemPg, Pmu*OnemPbError/OnemPb/OnemPb) << std::endl; } int main(int argc, char **argv) { diff --git a/docs/part4/usefullinks.md b/docs/part4/usefullinks.md index 24f7fcea639..cc7ef28d437 100644 --- a/docs/part4/usefullinks.md +++ b/docs/part4/usefullinks.md @@ -86,6 +86,6 @@ There is no document currently which can be cited for using the combine tool, ho * _My nuisances are (artificially) constrained and/or the impact plot show some strange behaviour, especially after including MC statistical uncertainties. What can I do?_ * Depending on the details of the analysis, several solutions can be adopted to mitigate these effects. We advise to run the validation tools at first, to identify possible redundant shape uncertainties that can be safely eliminated or replaced with lnN ones. Any remaining artificial constrain should be studies. Possible mitigating strategies can be to (a) smooth the templates or (b) adopt some rebinning in order to reduce statistical fluctuations in the templates. A description of possible strategies and effects can be found in [this talk by Margaret Eminizer](https://indico.cern.ch/event/788727/contributions/3401374/attachments/1831680/2999825/higgs_combine_4_17_2019_fitting_details.pdf) * _What do CLs, CLs+b and CLb in the code mean?_ - * The names CLs+b and CLb are rather outdated and should instead be referred to as p-values - $p_{\mu}$ and $1-p_{b}$, respectively. We use the CLs (which itself is not a p-value) criterion often in High energy physics as it is designed to avoid excluding a signal model when the sensitivity is low (and protects against excluding due to underfluctuations in the data). Typically, when excluding a signal model the p-value $p_{\mu}$ often refers to the p-value under the signal+background hypothesis, assuming a particular value of the signal stregth ($\mu$) while $p_{b}$ is the p-value under the background only hypothesis. You can find more details and definitions of the CLs criterion and $p_{\mu}$ and $p_{b}$ in section 39.4.2.4 the [2016 PDG review](http://pdg.lbl.gov/2016/reviews/rpp2016-rev-statistics.pdf). + * The names CLs+b and CLb what are found within some of the RooStats tools are rather outdated and should instead be referred to as p-values - $p_{\mu}$ and $1-p_{b}$, respectively. We use the CLs (which itself is not a p-value) criterion often in High energy physics as it is designed to avoid excluding a signal model when the sensitivity is low (and protects against excluding due to underfluctuations in the data). Typically, when excluding a signal model the p-value $p_{\mu}$ often refers to the p-value under the signal+background hypothesis, assuming a particular value of the signal stregth ($\mu$) while $p_{b}$ is the p-value under the background only hypothesis. You can find more details and definitions of the CLs criterion and $p_{\mu}$ and $p_{b}$ in section 39.4.2.4 of the [2016 PDG review](http://pdg.lbl.gov/2016/reviews/rpp2016-rev-statistics.pdf). diff --git a/docs/part5/longexercise.md b/docs/part5/longexercise.md index 3549809b404..30f42c5015c 100644 --- a/docs/part5/longexercise.md +++ b/docs/part5/longexercise.md @@ -141,7 +141,7 @@ Calculate the median expected limit and the 68% range. The 95% range could also **Tasks and questions:** - In contrast to `AsymptoticLimits`, with `HybridNew` each limit comes with an uncertainty. What is the origin of this uncertainty? - How good is the agreement between the asymptotic and toy-based methods? -- Why does it take longer to calculate the lower expected quantiles (e.g. 0.025, 0.16)? Think about how the statistical uncertainty on the CLs value depends on CLs+b and CLb. +- Why does it take longer to calculate the lower expected quantiles (e.g. 0.025, 0.16)? Think about how the statistical uncertainty on the CLs value depends on Pmu and Pb. Next plot the test statistic distributions stored in the output file: ```shell diff --git a/docs/part5/longexerciseanswers.md b/docs/part5/longexerciseanswers.md index d1140349439..859c8e1dd90 100644 --- a/docs/part5/longexerciseanswers.md +++ b/docs/part5/longexerciseanswers.md @@ -31,7 +31,7 @@
Show answer - The uncertainty is caused by the limited number of toys: the values of CLs+b and CLb come from counting the number of toys in the tails of the test statistic distributions. The number of toys used can be adjusted with the option --toysH + The uncertainty is caused by the limited number of toys: the values of Pmu and Pb come from counting the number of toys in the tails of the test statistic distributions. The number of toys used can be adjusted with the option --toysH
- How good is the agreement between the asymptotic and toy-based methods? @@ -41,11 +41,11 @@ The agreement should be pretty good in this example, but will generally break down once we get to the level of 0-5 events. - - Why does it take longer to calculate the lower expected quantiles (e.g. 0.025, 0.16)? Think about how the statistical uncertainty on the CLs value depends on CLs+b and CLb. + - Why does it take longer to calculate the lower expected quantiles (e.g. 0.025, 0.16)? Think about how the statistical uncertainty on the CLs value depends on Pmu and Pb.
Show answer - For this we need the definition of CLs = CLs+b / CLb. The 0.025 expected quantile is by definition where CLb = 0.025, so for a 95% CL limit we have CLs = 0.05, implying we are looking for the value of r where CLs+b = 0.00125. With 1000 s+b toys we would then only expect `1000 * 0.00125 = 1.25` toys in the tail region we have to integrate over. Contrast this to the median limit where 25 toys would be in this region. This means we have to generate a much larger numbers of toys to get the same statistical power. + For this we need the definition of CLs = Pmu / (1-Pb). The 0.025 expected quantile is by definition where Pb = 0.025, so for a 95% CL limit we have CLs = 0.05, implying we are looking for the value of r where Pmu = 0.00125. With 1000 s+b toys we would then only expect `1000 * 0.00125 = 1.25` toys in the tail region we have to integrate over. Contrast this to the median limit where 25 toys would be in this region. This means we have to generate a much larger numbers of toys to get the same statistical power.
diff --git a/docs/tutorial2020/exercise.md b/docs/tutorial2020/exercise.md index 299e9d857b8..3f7bde1fb44 100644 --- a/docs/tutorial2020/exercise.md +++ b/docs/tutorial2020/exercise.md @@ -701,7 +701,7 @@ Calculate the median expected limit and the 68% range. The 95% range could also **Tasks and questions:** - In contrast to `AsymptoticLimits`, with `HybridNew` each limit comes with an uncertainty. What is the origin of this uncertainty? - How good is the agreement between the asymptotic and toy-based methods? -- Why does it take longer to calculate the lower expected quantiles (e.g. 0.025, 0.16)? Think about how the statistical uncertainty on the CLs value depends on CLs+b and CLb. +- Why does it take longer to calculate the lower expected quantiles (e.g. 0.025, 0.16)? Think about how the statistical uncertainty on the CLs value depends on Pmu and Pb. Next plot the test statistic distributions stored in the output file: ```shell diff --git a/docs/tutorials-part-2.md b/docs/tutorials-part-2.md index 9f9924fe69e..7e69f8f7201 100644 --- a/docs/tutorials-part-2.md +++ b/docs/tutorials-part-2.md @@ -370,7 +370,7 @@ The F-C procedure for a generic model is: In combine, you can perform this test on each individual point (param1,param2,...) = (value1,value2,...) by doing - combine workspace.root -M HybridNew --freq --testStat=PL --rule=CLsplusb --singlePoint param1=value1,param2=value2,param3=value3,... [other options of HybridNew] + combine workspace.root -M HybridNew --freq --testStat=PL --rule=Pmu --singlePoint param1=value1,param2=value2,param3=value3,... [other options of HybridNew] The point belongs to your confidence region if CL<sub>s+b</sub> is larger than 1-CL (e.g. 0.3173 for a 1-sigma region, CL=0.6827). @@ -382,7 +382,7 @@ Imposing physical boundaries (such as requiring mu>0) can be achieved by sett As in general for HybridNew, you can split the task into multiple tasks and then merge the outputs, as described in the [HybridNew chapter](SWGuideHiggsAnalysisCombinedLimit#HybridNew_algorithm). -For uni-dimensional models only, and if the parameter behaves like a cross-section, the code is somewhat able to do interpolation and determine the values of your parameter on the contour (just like it does for the limits). In that case, the syntax is the same as per the CLs limits with [HybridNew chapter](HybridNew chapter) except that you want **`--testStat=PL --rule=CLsplusb`** . +For uni-dimensional models only, and if the parameter behaves like a cross-section, the code is somewhat able to do interpolation and determine the values of your parameter on the contour (just like it does for the limits). In that case, the syntax is the same as per the CLs limits with [HybridNew chapter](HybridNew chapter) except that you want **`--testStat=PL --rule=Pmu`** . **Extracting Contours** diff --git a/interface/HybridNew.h b/interface/HybridNew.h index d066ddd1c8c..650a8c46a5a 100644 --- a/interface/HybridNew.h +++ b/interface/HybridNew.h @@ -111,7 +111,7 @@ class HybridNew : public LimitAlgo { std::pair eval(RooStats::HybridCalculator &hc, const RooAbsCollection & rVals, bool adaptive=false, double clsTarget=-1) ; - // Compute CLs or CLsb from a HypoTestResult. In the case of the LHCFC test statistics, do the proper transformation to LHC or Profile (which needs to know the POI) + // Compute CLs or Pmu from a HypoTestResult. In the case of the LHCFC test statistics, do the proper transformation to LHC or Profile (which needs to know the POI) std::pair eval(const RooStats::HypoTestResult &hcres, const RooAbsCollection & rVals) ; std::pair eval(const RooStats::HypoTestResult &hcres, double rVal) ; diff --git a/macros/benchmark-simple.sh b/macros/benchmark-simple.sh index 031f5c0fab4..65da5d6ccac 100644 --- a/macros/benchmark-simple.sh +++ b/macros/benchmark-simple.sh @@ -5,5 +5,5 @@ INPUT=benchmark-simple.txt (combine $INPUT -M MarkovChainMC --proposal uniform -i 200000 -H ProfileLikelihood > $INPUT.log.MCMC.uniform 2>&1 &) (combine $INPUT -M MarkovChainMC --proposal gaus -i 200000 -H ProfileLikelihood > $INPUT.log.MCMC.gaus 2>&1 &) #(combine $INPUT -M Hybrid --rule CLs -H ProfileLikelihood > $INPUT.log.CLs.LEP 2>&1 &) -#(combine $INPUT -M Hybrid --rule CLsplusb -H ProfileLikelihood > $INPUT.log.CLsb.LEP 2>&1 &) +#(combine $INPUT -M Hybrid --rule Pmu -H ProfileLikelihood > $INPUT.log.Pmu.LEP 2>&1 &) diff --git a/src/AsymptoticLimits.cc b/src/AsymptoticLimits.cc index ba3ff391403..dfa7e8b7e62 100644 --- a/src/AsymptoticLimits.cc +++ b/src/AsymptoticLimits.cc @@ -54,7 +54,7 @@ LimitAlgo("AsymptoticLimits specific options") { //("minimizerTolerance", boost::program_options::value(&minimizerTolerance_)->default_value(minimizerTolerance_), "Tolerance for minimizer used for profiling") //("minimizerStrategy", boost::program_options::value(&minimizerStrategy_)->default_value(minimizerStrategy_), "Stragegy for minimizer") ("qtilde", boost::program_options::value(&qtilde_)->default_value(qtilde_), "Allow only non-negative signal strengths (default is true).") - ("rule", boost::program_options::value(&rule_)->default_value(rule_), "Rule to use: CLs, CLsplusb") + ("rule", boost::program_options::value(&rule_)->default_value(rule_), "Rule to use: CLs, Pmu") ("picky", "Abort on fit failures") ("noFitAsimov", "Use the pre-fit asimov dataset") ("getLimitFromGrid", boost::program_options::value(&gridFileName_), "calculates the limit from a grid of r,cls values") @@ -76,8 +76,8 @@ void AsymptoticLimits::applyOptions(const boost::program_options::variables_map noFitAsimov_ = vm.count("noFitAsimov") || vm.count("bypassFrequentistFit"); // aslo pick up base option from combine if (rule_=="CLs") doCLs_ = true; - else if (rule_=="CLsplusb") doCLs_ = false; - else throw std::invalid_argument("AsymptoticLimits: Rule must be either 'CLs' or 'CLsplusb'"); + else if (rule_=="Pmu") doCLs_ = false; + else throw std::invalid_argument("AsymptoticLimits: Rule must be either 'CLs' or 'Pmu'"); if (what_ == "blind") { what_ = "expected"; noFitAsimov_ = true; } if (noFitAsimov_) std::cout << "Will use a-priori expected background instead of a-posteriori one." << std::endl; @@ -328,34 +328,34 @@ double AsymptoticLimits::getCLs(RooRealVar &r, double rVal, bool getAlsoExpected } double qA = 2*(nllA_->getVal() - minNllA_); if (qA < 0) qA = 0; - double CLsb = ROOT::Math::normal_cdf_c(sqrt(qmu)); - double CLb = ROOT::Math::normal_cdf(sqrt(qA)-sqrt(qmu)); + double Pmu = ROOT::Math::normal_cdf_c(sqrt(qmu)); + double OnemPb = ROOT::Math::normal_cdf(sqrt(qA)-sqrt(qmu)); if (qtilde_ && qmu > qA) { // In this region, things are tricky double mos = sqrt(qA); // mu/sigma - CLsb = ROOT::Math::normal_cdf_c( (qmu + qA)/(2*mos) ); - CLb = ROOT::Math::normal_cdf_c( (qmu - qA)/(2*mos) ); + Pmu = ROOT::Math::normal_cdf_c( (qmu + qA)/(2*mos) ); + OnemPb = ROOT::Math::normal_cdf_c( (qmu - qA)/(2*mos) ); } - double CLs = (CLb == 0 ? 0 : CLsb/CLb); + double CLs = (OnemPb == 0 ? 0 : Pmu/OnemPb); sentry.clear(); if (verbose > 0) { - printf("At %s = %f:\tq_mu = %.5f\tq_A = %.5f\tCLsb = %7.5f\tCLb = %7.5f\tCLs = %7.5f\n", r.GetName(), rVal, qmu, qA, CLsb, CLb, CLs); - Logger::instance().log(std::string(Form("AsymptoticLimits.cc: %d -- At %s = %f:\tq_mu = %.5f\tq_A = %.5f\tCLsb = %7.5f\tCLb = %7.5f\tCLs = %7.5f",__LINE__,r.GetName(), rVal, qmu, qA, CLsb, CLb, CLs)),Logger::kLogLevelInfo,__func__); + printf("At %s = %f:\tq_mu = %.5f\tq_A = %.5f\tPmu = %7.5f\t1-Pb = %7.5f\tCLs = %7.5f\n", r.GetName(), rVal, qmu, qA, Pmu, OnemPb, CLs); + Logger::instance().log(std::string(Form("AsymptoticLimits.cc: %d -- At %s = %f:\tq_mu = %.5f\tq_A = %.5f\tPmu = %7.5f\t1-Pb = %7.5f\tCLs = %7.5f",__LINE__,r.GetName(), rVal, qmu, qA, Pmu, OnemPb, CLs)),Logger::kLogLevelInfo,__func__); } if (getAlsoExpected) { const double quantiles[5] = { 0.025, 0.16, 0.50, 0.84, 0.975 }; for (int iq = 0; iq < 5; ++iq) { double N = ROOT::Math::normal_quantile(quantiles[iq], 1.0); - double clb = quantiles[iq]; - double clsplusb = ROOT::Math::normal_cdf_c( sqrt(qA) - N, 1.); - if (doCLs_) { *limit = (clb != 0 ? clsplusb/clb : 0); *limitErr = 0 ; } - else { *limit = (clsplusb); *limitErr = 0; } + double pb = quantiles[iq]; // note that this is really 1-pb ! + double pmu = ROOT::Math::normal_cdf_c( sqrt(qA) - N, 1.); + if (doCLs_) { *limit = (pb != 0 ? pmu/pb : 0); *limitErr = 0 ; } + else { *limit = (pmu); *limitErr = 0; } Combine::commitPoint(true, quantiles[iq]); - if (verbose > 0) printf("Expected %4.1f%%: CLsb = %.5f CLb = %.5f CLs = %.5f\n", quantiles[iq]*100, clsplusb, clb, clsplusb/clb); + if (verbose > 0) printf("Expected %4.1f%%: Pmu = %.5f Pb = %.5f CLs = %.5f\n", quantiles[iq]*100, pmu, pb, pmu/pb); } } - return doCLs_ ? CLs : CLsb ; + return doCLs_ ? CLs : Pmu ; } std::vector > AsymptoticLimits::runLimitExpected(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint) { @@ -462,17 +462,17 @@ std::vector > AsymptoticLimits::runLimitExpected(RooWorks } -float AsymptoticLimits::findExpectedLimitFromCrossing(RooAbsReal &nll, RooRealVar *r, double rMin, double rMax, double nll0, double clb) { +float AsymptoticLimits::findExpectedLimitFromCrossing(RooAbsReal &nll, RooRealVar *r, double rMin, double rMax, double nll0, double pb) { // EQ 37 of CMS NOTE 2011-005: // mu_N = sigma * ( normal_quantile_c( (1-cl) * normal_cdf(N) ) + N ) - // --> (mu_N/sigma) = N + normal_quantile_c( (1-cl) * clb ) + // --> (mu_N/sigma) = N + normal_quantile_c( (1-cl) * (1-Pb) ) but in our code here we refer to pb=1-Pb // but qmu = (mu_N/sigma)^2 - // --> qmu = [ N + normal_quantile_c( (1-cl)*CLb ) ]^2 + // --> qmu = [ N + normal_quantile_c( (1-cl)*(1-Pb) ) ]^2 // remember that qmu = 2*nll - double N = ROOT::Math::normal_quantile(clb, 1.0); - double errorlevel = 0.5 * pow(N+ROOT::Math::normal_quantile_c((doCLs_ ? clb:1.)*(1-cl),1.0), 2); + double N = ROOT::Math::normal_quantile(pb, 1.0); + double errorlevel = 0.5 * pow(N+ROOT::Math::normal_quantile_c((doCLs_ ? pb:1.)*(1-cl),1.0), 2); int minosStat = -1; if (minosAlgo_ == "minos") { double rMax0 = r->getMax(); diff --git a/src/HybridNew.cc b/src/HybridNew.cc index 468793ff3b4..780b02784df 100644 --- a/src/HybridNew.cc +++ b/src/HybridNew.cc @@ -101,7 +101,7 @@ std::string HybridNew::mode_ = ""; HybridNew::HybridNew() : LimitAlgo("HybridNew specific options") { options_.add_options() - ("rule", boost::program_options::value(&rule_)->default_value(rule_), "Rule to use: CLs, CLsplusb") + ("rule", boost::program_options::value(&rule_)->default_value(rule_), "Rule to use: CLs, Pmu") ("testStat",boost::program_options::value(&testStat_)->default_value(testStat_), "Test statistics: LEP, TEV, LHC (previously known as Atlas), Profile.") ("singlePoint", boost::program_options::value(&rValue_)->default_value(rValue_), "Just compute CLs for the given value of the parameter of interest. In case of multiple parameters, use a syntax 'name=value,name2=value2,...'") ("onlyTestStat", "Just compute test statistic for the data (or toy if using -t N), i.e don't throw toys to calculate actual p-values (works only with --singlePoint)") @@ -109,7 +109,7 @@ LimitAlgo("HybridNew specific options") { ("generateExternalMeasurements", boost::program_options::value(&genGlobalObs_)->default_value(genGlobalObs_), "Generate external measurements for each toy, taken from the GlobalObservables of the ModelConfig") ("fitNuisances", boost::program_options::value(&fitNuisances_)->default_value(fitNuisances_), "Fit the nuisances first, before generating the toy data. Set this option to false to acheive the same results as with --bypassFrequentistFit. When not generating toys, eg as in when running with --readHybridresult, this has no effect") ("searchAlgo", boost::program_options::value(&algo_)->default_value(algo_), "Algorithm to use to search for the limit (bisection, logSecant)") - ("toysH,T", boost::program_options::value(&nToys_)->default_value(nToys_), "Number of Toy MC extractions to compute CLs+b, CLb and CLs") + ("toysH,T", boost::program_options::value(&nToys_)->default_value(nToys_), "Number of Toy MC extractions to compute Pmu, Pb and CLs") ("clsAcc", boost::program_options::value(&clsAccuracy_ )->default_value(clsAccuracy_), "Absolute accuracy on CLs to reach to terminate the scan") ("rAbsAcc", boost::program_options::value(&rAbsAccuracy_)->default_value(rAbsAccuracy_), "Absolute accuracy on r to reach to terminate the scan") ("rRelAcc", boost::program_options::value(&rRelAccuracy_)->default_value(rRelAccuracy_), "Relative accuracy on r to reach to terminate the scan") @@ -122,7 +122,7 @@ LimitAlgo("HybridNew specific options") { ("grid", boost::program_options::value(&gridFile_), "Use the specified file containing a grid of SamplingDistributions for the limit (implies readHybridResults).\n For --singlePoint or --signif use --toysFile=x.root --readHybridResult instead of this.") ("expectedFromGrid", boost::program_options::value(&quantileForExpectedFromGrid_)->default_value(0.5), "Use the grid to compute the expected limit for this quantile") ("signalForSignificance", boost::program_options::value()->default_value("1"), "Use this value of the parameter of interest when generating signal toys for expected significance (same syntax as --singlePoint)") - ("clsQuantiles", boost::program_options::value(&clsQuantiles_)->default_value(clsQuantiles_), "Compute correct quantiles of CLs or CLsplusb instead of assuming they're the same as CLb ones") + ("clsQuantiles", boost::program_options::value(&clsQuantiles_)->default_value(clsQuantiles_), "Compute correct quantiles of CLs or Pmu instead of assuming they're the same as those for Pb") //("importanceSamplingNull", boost::program_options::value(&importanceSamplingNull_)->default_value(importanceSamplingNull_), // "Enable importance sampling for null hypothesis (background only)") //("importanceSamplingAlt", boost::program_options::value(&importanceSamplingAlt_)->default_value(importanceSamplingAlt_), @@ -145,7 +145,7 @@ LimitAlgo("HybridNew specific options") { ("importantContours",boost::program_options::value(&scaleAndConfidenceSelection_)->default_value(scaleAndConfidenceSelection_), "Throw less toys far from interesting contours , format : CL_1,CL_2,..CL_N (--toysH scaled down when prob is far from any of CL_i) ") ("maxProbability", boost::program_options::value(&maxProbability_)->default_value(maxProbability_), "when point is > maxProbability countour, don't bother throwing toys") ("confidenceTolerance", boost::program_options::value(&confidenceToleranceForToyScaling_)->default_value(confidenceToleranceForToyScaling_), "Determine what 'far' means for adatptiveToys. (relative in terms of (1-cl))") - ("LHCmode", boost::program_options::value(&mode_)->default_value(mode_), "Shortcuts for LHC style running modes. --LHCmode LHC-significance: --generateNuisances=0 --generateExternalMeasurements=1 --fitNuisances=1 --testStat=LHC (Q_LHC, modified for discovery) --significance, --LHCmode LHC-limits: --generateNuisances=0 --generateExternalMeasurements=1 --fitNuisances=1 --testStat=LHC (Q_LHC, modified for upper limits) --rule CLs, --LHCmode LHC-feldman-cousins: --generateNuisances=0 --generateExternalMeasurements=1 --fitNuisances=1 --testStat=PL (Q_Profile, includes boundaries) --rule CLsplusb") + ("LHCmode", boost::program_options::value(&mode_)->default_value(mode_), "Shortcuts for LHC style running modes. --LHCmode LHC-significance: --generateNuisances=0 --generateExternalMeasurements=1 --fitNuisances=1 --testStat=LHC (Q_LHC, modified for discovery) --significance, --LHCmode LHC-limits: --generateNuisances=0 --generateExternalMeasurements=1 --fitNuisances=1 --testStat=LHC (Q_LHC, modified for upper limits) --rule CLs, --LHCmode LHC-feldman-cousins: --generateNuisances=0 --generateExternalMeasurements=1 --fitNuisances=1 --testStat=PL (Q_Profile, includes boundaries) --rule Pmu") ; } @@ -186,7 +186,7 @@ void HybridNew::applyOptions(const boost::program_options::variables_map &vm) { } else if (mode_=="LHC-feldman-cousins"){ genNuisances_ = 0; genGlobalObs_ = withSystematics; fitNuisances_ = withSystematics; testStat_ = "PL"; - rule_ = "CLsplusb"; + rule_ = "Pmu"; doFC_ = true; } else { throw std::invalid_argument(Form("HybridNew: invalid LHCmode %s, only --LHCmode 'LHC-significance', 'LHC-limits', and 'LHC-feldman-cousins' supported\n",mode_.c_str())); @@ -232,13 +232,13 @@ void HybridNew::validateOptions() { if (fork_ > 1) nToys_ /= fork_; // makes more sense if (rule_ == "CLs") { CLs_ = true; - } else if (rule_ == "CLsplusb") { + } else if (rule_ == "Pmu") { CLs_ = false; } else if (rule_ == "FC") { CLs_ = false; doFC_ = true; } else { - throw std::invalid_argument("HybridNew: Rule must be either 'CLs' or 'CLsplusb'"); + throw std::invalid_argument("HybridNew: Rule must be either 'CLs' or 'Pmu'"); } if (doFC_) noUpdateGrid_ = true; // Needed since addition of points can skew interval if (testStat_ == "SimpleLikelihoodRatio" || testStat_ == "SLR" ) { testStat_ = "LEP"; } @@ -304,7 +304,7 @@ bool HybridNew::runSignificance(RooWorkspace *w, RooStats::ModelConfig *mc_s, Ro for (unsigned int i = 1; i < iterations_; ++i) { std::unique_ptr more(evalGeneric(*hc)); hcResult->Append(more.get()); - if (verbose) std::cout << "\t1 - CLb = " << hcResult->CLb() << " +/- " << hcResult->CLbError() << std::endl; + if (verbose) std::cout << "\t1 - Pb = " << hcResult->CLb() << " +/- " << hcResult->CLbError() << std::endl; } } if (hcResult.get() == 0) { @@ -398,7 +398,7 @@ bool HybridNew::runLimit(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats: rMax = limitPlot_->GetX()[limitPlot_->GetN()-1]; for (int i = 0, n = limitPlot_->GetN(); i < n; ++i) { double x = limitPlot_->GetX()[i], y = limitPlot_->GetY()[i], ey = limitPlot_->GetErrorY(i); - if (verbose > 0) printf(" r %.2f: %s = %6.4f +/- %6.4f\n", x, CLs_ ? "CLs" : "CLsplusb", y, ey); + if (verbose > 0) printf(" r %.2f: %s = %6.4f +/- %6.4f\n", x, CLs_ ? "CLs" : "Pmu", y, ey); if (saveGrid_) { limit = x; limitErr = ey; Combine::commitPoint(false, y); } if (y-3*max(ey,0.01) >= clsTarget) { rMin = x; clsMin = CLs_t(y,ey); } if (fabs(y-clsTarget) < minDist) { nearest = x; minDist = fabs(y-clsTarget); } @@ -424,7 +424,7 @@ bool HybridNew::runLimit(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats: if (clsMax.first == 0 || clsMax.first + 3 * fabs(clsMax.second) < clsTarget ) break; rMax += rMax; if (tries == 5) { - std::cerr << "Cannot set higher limit: at " << r->GetName() << " = " << rMax << " still get " << (CLs_ ? "CLs" : "CLsplusb") << " = " << clsMax.first << std::endl; + std::cerr << "Cannot set higher limit: at " << r->GetName() << " = " << rMax << " still get " << (CLs_ ? "CLs" : "Pmu") << " = " << clsMax.first << std::endl; return false; } } @@ -441,7 +441,7 @@ bool HybridNew::runLimit(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats: if (clsMin.first == 1 || clsMin.first - 3 * fabs(clsMin.second) > clsTarget) break; rMin += rMin; if (tries == 5) { - std::cerr << "Cannot set lower limit: at " << r->GetName() << " = " << rMin << " still get " << (CLs_ ? "CLs" : "CLsplusb") << " = " << clsMin.first << std::endl; + std::cerr << "Cannot set lower limit: at " << r->GetName() << " = " << rMin << " still get " << (CLs_ ? "CLs" : "Pmu") << " = " << clsMin.first << std::endl; return false; } } @@ -632,7 +632,7 @@ bool HybridNew::runLimit(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats: bool HybridNew::runSinglePoint(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint) { std::pair result = eval(w, mc_s, mc_b, data, rValues_, clsAccuracy_ != 0); std::cout << "\n -- Hybrid New -- \n"; - std::cout << (CLs_ ? "CLs = " : "CLsplusb = ") << result.first << " +/- " << result.second << std::endl; + std::cout << (CLs_ ? "CLs = " : "Pmu = ") << result.first << " +/- " << result.second << std::endl; if (verbose > 1) std::cout << "Total toys: " << perf_totalToysRun_ << std::endl; limit = result.first; limitErr = result.second; @@ -1092,7 +1092,7 @@ HybridNew::eval(RooStats::HybridCalculator &hc, const RooAbsCollection & rVals, hcResult->SetPValueIsRightTail(!hcResult->GetPValueIsRightTail()); } std::pair cls = eval(*hcResult, rVals); - if (verbose) std::cout << (CLs_ ? "\tCLs = " : "\tCLsplusb = ") << cls.first << " +/- " << cls.second << std::endl; + if (verbose) std::cout << (CLs_ ? "\tCLs = " : "\tPmu = ") << cls.first << " +/- " << cls.second << std::endl; if (adaptive) { if (CLs_) { hc.SetToys(int(0.25*nToys_ + 1), nToys_); @@ -1111,7 +1111,7 @@ HybridNew::eval(RooStats::HybridCalculator &hc, const RooAbsCollection & rVals, hcResult->Append(more.get()); if (expectedFromGrid_) applyExpectedQuantile(*hcResult); cls = eval(*hcResult, rVals); - if (verbose) std::cout << (CLs_ ? "\tCLs = " : "\tCLsplusb = ") << cls.first << " +/- " << cls.second << std::endl; + if (verbose) std::cout << (CLs_ ? "\tCLs = " : "\tPmu = ") << cls.first << " +/- " << cls.second << std::endl; } } else if (iterations_ > 1) { for (unsigned int i = 1; i < iterations_; ++i) { @@ -1121,17 +1121,18 @@ HybridNew::eval(RooStats::HybridCalculator &hc, const RooAbsCollection & rVals, hcResult->Append(more.get()); if (expectedFromGrid_) applyExpectedQuantile(*hcResult); cls = eval(*hcResult, rVals); - if (verbose) std::cout << (CLs_ ? "\tCLs = " : "\tCLsplusb = ") << cls.first << " +/- " << cls.second << std::endl; + if (verbose) std::cout << (CLs_ ? "\tCLs = " : "\tPmu = ") << cls.first << " +/- " << cls.second << std::endl; } } if (verbose > 0) { + // Note that here, RooFit uses "CLsplusb" and "CLb", whereas we use Pmu and Pb std::cout << "\tCLs = " << hcResult->CLs() << " +/- " << hcResult->CLsError() << "\n" << - "\tCLb = " << hcResult->CLb() << " +/- " << hcResult->CLbError() << "\n" << - "\tCLsplusb = " << hcResult->CLsplusb() << " +/- " << hcResult->CLsplusbError() << "\n" << + "\t1-Pb = " << hcResult->CLb() << " +/- " << hcResult->CLbError() << "\n" << + "\tPmu = " << hcResult->CLsplusb() << " +/- " << hcResult->CLsplusbError() << "\n" << std::endl; - Logger::instance().log(std::string(Form("HybridNew.cc: %d -- CLs = %g +/- %g\n\tCLb = %g +/- %g\n\tCLsplusb = %g +/- %g",__LINE__ + Logger::instance().log(std::string(Form("HybridNew.cc: %d -- CLs = %g +/- %g\n\t1-Pb = %g +/- %g\n\tPmu = %g +/- %g",__LINE__ ,hcResult->CLs(), hcResult->CLsError() ,hcResult->CLb(), hcResult->CLbError() ,hcResult->CLsplusb(), hcResult->CLsplusbError())) @@ -1273,12 +1274,12 @@ void HybridNew::applyClsQuantile(RooStats::HypoTestResult &hcres) { std::vector >::const_iterator match = std::upper_bound(sbegin, send, std::pair(it->first, 0)); if (match == send) { //std::cout << "Did not find match, for it->first == " << it->first << ", as back = ( " << scumul.back().first << " , " << scumul.back().second << " ) " << std::endl; - double clsb = (scumul.back().second > 0.5 ? 1.0 : 0.0), clb = runningSum*binv, cls = clsb / clb; - xcumul.push_back(std::make_pair(CLs_ ? cls : clsb, *it)); + double pmu = (scumul.back().second > 0.5 ? 1.0 : 0.0), pb = runningSum*binv, cls = pmu / pb; // pb := 1-Pb + xcumul.push_back(std::make_pair(CLs_ ? cls : pmu, *it)); } else { - double clsb = match->second, clb = runningSum*binv, cls = clsb / clb; - //if ((++k) % 100 == 0) printf("At %+8.5f CLb = %6.4f, CLsplusb = %6.4f, CLs =%7.4f\n", it->first, clb, clsb, cls); - xcumul.push_back(std::make_pair(CLs_ ? cls : clsb, *it)); + double pmu = match->second, pb = runningSum*binv, cls = pmu / pb; + //if ((++k) % 100 == 0) printf("At %+8.5f 1-Pb = %6.4f, Pmu = %6.4f, CLs =%7.4f\n", it->first, pb, pmu, cls); + xcumul.push_back(std::make_pair(CLs_ ? cls : pmu, *it)); } } // sort @@ -1487,10 +1488,10 @@ RooStats::HypoTestResult * HybridNew::readToysFromFile(const RooAbsCollection & if (verbose > 0) { std::cout << "\tCLs = " << ret->CLs() << " +/- " << ret->CLsError() << "\n" << - "\tCLb = " << ret->CLb() << " +/- " << ret->CLbError() << "\n" << - "\tCLsplusb = " << ret->CLsplusb() << " +/- " << ret->CLsplusbError() << "\n" << + "\t1-Pb = " << ret->CLb() << " +/- " << ret->CLbError() << "\n" << + "\tPmu = " << ret->CLsplusb() << " +/- " << ret->CLsplusbError() << "\n" << std::endl; - Logger::instance().log(std::string(Form("HybridNew.cc: %d -- CLs = %g +/- %g\n\tCLb = %g +/- %g\n\tCLsplusb = %g +/- %g",__LINE__ + Logger::instance().log(std::string(Form("HybridNew.cc: %d -- CLs = %g +/- %g\n\t1-Pb = %g +/- %g\n\tPmu = %g +/- %g",__LINE__ ,ret->CLs(), ret->CLsError() ,ret->CLb(), ret->CLbError() ,ret->CLsplusb(), ret->CLsplusbError())) @@ -1649,12 +1650,12 @@ std::pair HybridNew::updateGridPoint(RooWorkspace *w, RooStats::M if (verbose > 1) { std::cout << "At " << r->GetName() << " = " << point->first << ":\n" << "\tCLs = " << point->second->CLs() << " +/- " << point->second->CLsError() << "\n" << - "\tCLb = " << point->second->CLb() << " +/- " << point->second->CLbError() << "\n" << - "\tCLsplusb = " << point->second->CLsplusb() << " +/- " << point->second->CLsplusbError() << "\n" << + "\t1-Pb = " << point->second->CLb() << " +/- " << point->second->CLbError() << "\n" << + "\tCPmu = " << point->second->CLsplusb() << " +/- " << point->second->CLsplusbError() << "\n" << std::endl; } if (verbose){ - Logger::instance().log(std::string(Form("HybridNew.cc: %d -- At %s = %g:\n, \tCLs = %g +/- %g\n\tCLb = %g +/- %g\n\tCLsplusb = %g +/- %g",__LINE__ + Logger::instance().log(std::string(Form("HybridNew.cc: %d -- At %s = %g:\n, \tCLs = %g +/- %g\n\t1-Pb = %g +/- %g\n\tPmu = %g +/- %g",__LINE__ ,r->GetName(), point->first ,point->second->CLs(), point->second->CLsError() ,point->second->CLb(), point->second->CLbError() diff --git a/test/makeGrid2DUsingCrab.py b/test/makeGrid2DUsingCrab.py index 433a03385d3..b9fbf38a5d7 100755 --- a/test/makeGrid2DUsingCrab.py +++ b/test/makeGrid2DUsingCrab.py @@ -233,7 +233,7 @@ toys = "$n" what = "--singlePoint %s=%g,%s=%g " % (options.POI[0], x[0], options.POI[1], x[1]) script.write( - "{cond} ./combine {wsp} -M HybridNew {opts} -m {mass} --testStat=PL --rule=CLsplusb --fork $nchild -T {T} --clsAcc 0 -v {v} -n {out} --saveHybridResult --saveToys -s {seed} -i {toys} {what} {ranges} \n".format( + "{cond} ./combine {wsp} -M HybridNew {opts} -m {mass} --testStat=PL --rule=Pmu --fork $nchild -T {T} --clsAcc 0 -v {v} -n {out} --saveHybridResult --saveToys -s {seed} -i {toys} {what} {ranges} \n".format( wsp=workspace, opts=options.options, T=options.T, diff --git a/test/plotTestStatCLs.py b/test/plotTestStatCLs.py index d5f11eaabba..797f9f41d0a 100755 --- a/test/plotTestStatCLs.py +++ b/test/plotTestStatCLs.py @@ -75,7 +75,7 @@ "--quantileExpected", default=-1, type="float", - help="Replace observation with expected quantile (under alt hyp, i.e CLb=quantileExpected)", + help="Replace observation with expected quantile (under alt hyp, i.e 1-Pb=quantileExpected)", ) parser.add_option( "", @@ -89,7 +89,7 @@ "--signif", action="store_true", default=False, - help="If significance, don't plot CLs+b, i.e make the q0 plot", + help="If significance, don't plot Pmu, i.e make the q0 plot", ) (options, args) = parser.parse_args() diff --git a/test/plotting/lepLikePlot.cxx b/test/plotting/lepLikePlot.cxx index b98251ec4da..7ec652b27d8 100644 --- a/test/plotting/lepLikePlot.cxx +++ b/test/plotting/lepLikePlot.cxx @@ -34,29 +34,28 @@ void addLepPoint(LepBand &b, double x, TString file, double r=1.0) { if (!toyDir) throw std::logic_error("Cannot use readHypoTestResult: option toysFile not specified, or input file empty"); RooStats::HypoTestResult *res = readLepFile(toyDir, r); double clsObs = res->CLs(), clsObsErr = res->CLsError(); - double clsbObs = res->CLsplusb(), clsbObsErr = res->CLsplusbError(); - //std::cout << "Observed CLs = " << res->CLs() << " +/- " << res->CLsError() << std::endl; - //std::cout << "Observed CLsb = " << res->CLsplusb() << " +/- " << res->CLsplusbError() << std::endl; + double pmuObs = res->CLsplusb(), pmuObsErr = res->CLsplusbError(); + std::vector samples = res->GetNullDistribution()->GetSamplingDistribution(); int nd = samples.size(); std::sort(samples.begin(), samples.end()); double median = (samples.size() % 2 == 0 ? 0.5*(samples[nd/2]+samples[nd/2+1]) : samples[nd/2]); double summer68 = samples[floor(nd * 0.5*(1-0.68)+0.5)], winter68 = samples[TMath::Min(int(floor(nd * 0.5*(1+0.68)+0.5)), nd-1)]; double summer95 = samples[floor(nd * 0.5*(1-0.95)+0.5)], winter95 = samples[TMath::Min(int(floor(nd * 0.5*(1+0.95)+0.5)), nd-1)]; - res->SetTestStatisticData(median+1e-6); double clsMid = res->CLs(), clsbMid = res->CLsplusb(); - res->SetTestStatisticData(summer68+1e-6); double cls68L = res->CLs(), clsb68L = res->CLsplusb(); - res->SetTestStatisticData(winter68+1e-6); double cls68H = res->CLs(), clsb68H = res->CLsplusb(); - res->SetTestStatisticData(summer95+1e-6); double cls95L = res->CLs(), clsb95L = res->CLsplusb(); - res->SetTestStatisticData(winter95+1e-6); double cls95H = res->CLs(), clsb95H = res->CLsplusb(); + res->SetTestStatisticData(median+1e-6); double clsMid = res->CLs(), pmuMid = res->CLsplusb(); + res->SetTestStatisticData(summer68+1e-6); double cls68L = res->CLs(), pmu68L = res->CLsplusb(); + res->SetTestStatisticData(winter68+1e-6); double cls68H = res->CLs(), pmu68H = res->CLsplusb(); + res->SetTestStatisticData(summer95+1e-6); double cls95L = res->CLs(), pmu95L = res->CLsplusb(); + res->SetTestStatisticData(winter95+1e-6); double cls95H = res->CLs(), pmu95H = res->CLsplusb(); /* std::cout << "Expected CLs: median " << clsMid << ", 68% band [" << cls68L << ", " << cls68H << "]" << ", 95% band [" << cls95L << ", " << cls95H << "]" << std::endl; - std::cout << "Expected CLsb: median " << clsbMid - << ", 68% band [" << clsb68L << ", " << clsb68H << "]" - << ", 95% band [" << clsb95L << ", " << clsb95H << "]" << std::endl; + std::cout << "Expected Pmu: median " << pmuMid + << ", 68% band [" << pmu68L << ", " << pmu68H << "]" + << ", 95% band [" << pmu95L << ", " << pmu95H << "]" << std::endl; */ - b.obs->Set(b.n+1); b.obs->SetPoint(b.n, x, clsObs); b.obs->SetPointError(b.n, 0, 0, clsObsErr, clsObsErr); + b.obs->Set(b.n+1); b.obs->SetPoint(b.n, x, clsObs); b.obs->SetPointError(b.n, 0, 0, clsObsErr, clsObsErr); b.exp68->Set(b.n+1); b.exp68->SetPoint(b.n, x, clsMid); b.exp68->SetPointError(b.n, 0, 0, clsMid-cls68L, cls68H-clsMid); b.exp95->Set(b.n+1); b.exp95->SetPoint(b.n, x, clsMid); b.exp95->SetPointError(b.n, 0, 0, clsMid-cls95L, cls95H-clsMid); b.n++; diff --git a/test/plotting/qmuPlot.cxx b/test/plotting/qmuPlot.cxx index c9eb1497f1c..0347684c060 100644 --- a/test/plotting/qmuPlot.cxx +++ b/test/plotting/qmuPlot.cxx @@ -137,14 +137,14 @@ TCanvas *q0Plot(float mass, std::string poinam , int rebin=0) { qO->SetLineColor(kBlack); qO->SetLineWidth(3); - double clB; - clB = tailReal(t,"qB",qObs,0); + double pB; + pB = tailReal(t,"qB",qObs,0); - double clBerr = sqrt(clB*(1-clB)/nB); - double sig = RooStats::PValueToSignificance(clB); - double sigerr = 0.5*( TMath::Abs(RooStats::PValueToSignificance(clB+clBerr)-sig) + TMath::Abs(RooStats::PValueToSignificance(clB-clBerr)-sig)); + double pBerr = sqrt(pB*(1-pB)/nB); + double sig = RooStats::PValueToSignificance(pB); + double sigerr = 0.5*( TMath::Abs(RooStats::PValueToSignificance(pB+pBerr)-sig) + TMath::Abs(RooStats::PValueToSignificance(pB-pBerr)-sig)); - printf("P-val (CLb) = %.4f +/- %.4f\n", clB , clBerr); + printf("P-val (1-Pb) = %.4f +/- %.4f\n", pB , pBerr); printf("Signif = %.1f +/- %.2f sigma\n", sig, sigerr); // Worst way to calculate ! @@ -167,7 +167,7 @@ TCanvas *q0Plot(float mass, std::string poinam , int rebin=0) { leg2->SetTextFont(42); leg2->SetTextSize(0.04); leg2->SetLineColor(1); - leg2->AddEntry(qB1, Form("CL_{b} = %.4f (%.1f#sigma)", clB,sig), "F"); + leg2->AddEntry(qB1, Form("CL_{b} = %.4f (%.1f#sigma)", pB,sig), "F"); qB->Draw(); qB1->Draw("HIST SAME"); qO->Draw(); @@ -261,24 +261,24 @@ TCanvas *qmuPlot(float mass, std::string poinam, double poival, int mode=0, int qO->SetLineColor(kBlack); qO->SetLineWidth(3); - double clSB; - double clB; + double pMU; + double pB; if (invert){ - clSB = tailReal(t,"qB",qObs,mode); - clB = tailReal(t,"qS",qObs,mode); + pMU = tailReal(t,"qB",qObs,mode); + pB = tailReal(t,"qS",qObs,mode); } else { - clSB = tailReal(t,"qS",qObs,mode); - clB = tailReal(t,"qB",qObs,mode); + pMU = tailReal(t,"qS",qObs,mode); + pB = tailReal(t,"qB",qObs,mode); } - //double clSB = qS1->Integral(), clB = qB1->Integral(), - double clS = clSB/clB; - double clSBerr = sqrt(clSB*(1-clSB)/nS); - double clBerr = sqrt(clB*(1-clB)/nB); - double clSerr = clS * TMath::Hypot(clBerr/clB, clSBerr/clSB); - printf("CLs+b = %.4f +/- %.4f\n", clSB, clSBerr); - printf("CLb = %.4f +/- %.4f\n", clB , clBerr); - printf("CLs = %.4f +/- %.4f\n", clS , clSerr); + //double pMU = qS1->Integral(), pB := qB1->Integral() (# note this is really 1-p_{b}) + double clS = pMU/pB; + double pMUerr = sqrt(pMU*(1-pMU)/nS); + double pBerr = sqrt(pB*(1-pB)/nB); + double pSerr = clS * TMath::Hypot(pBerr/pB, pMUerr/pMU); + printf("Pmu = %.4f +/- %.4f\n", pMU , pMUerr); + printf("1-Pb = %.4f +/- %.4f\n", pB , pBerr); + printf("CLs = %.4f +/- %.4f\n", pMu , pMuerr); // Worst way to calculate ! TH1F *qS1 = tail(qS, qObs,mode); @@ -313,10 +313,10 @@ TCanvas *qmuPlot(float mass, std::string poinam, double poival, int mode=0, int leg2->SetTextSize(0.04); leg2->SetLineColor(1); //if (mode==0) - leg2->AddEntry(qS1, Form("CL_{s+b} = %.4f", clSB), "F"); + leg2->AddEntry(qS1, Form("p_{#mu} = %.4f", pMU), "F"); //if (mode==0) - leg2->AddEntry(qB1, Form("CL_{b} = %.4f", clB), "F"); - leg2->AddEntry("", Form("CL_{s} = %.4f", clS), ""); + leg2->AddEntry(qB1, Form("1-p_{b} = %.4f", pB), "F"); + leg2->AddEntry("", Form("CL_{s} = %.4f", pS), ""); qB->Draw(); diff --git a/test/validation/reference.json b/test/validation/reference.json index 00fec4ec10d..56914857a27 100644 --- a/test/validation/reference.json +++ b/test/validation/reference.json @@ -654,7 +654,7 @@ }, "status": "done" }, - "Counting_pValues_CLsplusb_Hybrid": { + "Counting_pValues_Pmu_Hybrid": { "comment": "NOTE: for S = 5", "results": { "counting-B5p5-Obs6-Syst30B.txt LEP": { @@ -674,7 +674,7 @@ }, "status": "done" }, - "Counting_pValues_CLsplusb_HybridNew": { + "Counting_pValues_Pmu_HybridNew": { "comment": "NOTE: for S = 5", "results": { "counting-B5p5-Obs6-Syst30B.txt Atlas": { @@ -1055,28 +1055,28 @@ }, "status": "done" }, - "Counting_pValues_Freq_0_CLsplusb_HybridNew": { + "Counting_pValues_Freq_0_Pmu_HybridNew": { "comment": "", "results": { "counting-B5p5-Obs6-StatOnly.txt": { "comment": "note: more toys used here", "limit": 0.098489999999999994, "limitErr": 0.00094228297182958796, "status": "done", "t_real": 0.24138765037059784 } }, "status": "done" }, - "Counting_pValues_Freq_B_CLsplusb_HybridNew": { + "Counting_pValues_Freq_B_Pmu_HybridNew": { "comment": "", "results": { "counting-B5p5-Obs6-Syst30B.txt": { "comment": "note: more toys used here", "limit": 0.093240000000000003, "limitErr": 0.00091949063290497964, "status": "done", "t_real": 1.5756527185440063 } }, "status": "done" }, - "Counting_pValues_Freq_S_CLsplusb_HybridNew": { + "Counting_pValues_Freq_S_Pmu_HybridNew": { "comment": "", "results": { "counting-B5p5-Obs6-Syst30S.txt": { "comment": "note: more toys used here", "limit": 0.089819999999999997, "limitErr": 0.00090417015876437777, "status": "done", "t_real": 1.5730928182601929 } }, "status": "done" }, - "Counting_pValues_Freq_U_CLsplusb_HybridNew": { + "Counting_pValues_Freq_U_Pmu_HybridNew": { "comment": "", "results": { "counting-B5p5-Obs6-Syst30U.txt": { "comment": "note: more toys used here", "limit": 0.11123, "limitErr": 0.00099427303644421521, "status": "done", "t_real": 1.7407045364379883 } diff --git a/test/validation/test_HN.py b/test/validation/test_HN.py index 5c7c401c1a7..cd5b3129a80 100644 --- a/test/validation/test_HN.py +++ b/test/validation/test_HN.py @@ -96,7 +96,7 @@ ] ### Test the p-values -for R in ["CLs", "CLsplusb"]: +for R in ["CLs", "Pmu"]: suite += [ ( M,