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.

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) {

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).       
 
 

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

docs/part5/longexerciseanswers.md

@@ -31,7 +31,7 @@
 <details>
 <summary><b>Show answer</b></summary>
 
-<b> 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 <code> --toysH </code> </b>
+<b> 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 <code> --toysH </code> </b>
 
 </details>
   - How good is the agreement between the asymptotic and toy-based methods?
@@ -41,11 +41,11 @@
 <b> The agreement should be pretty good in this example, but will generally break down once we get to the level of 0-5 events. </b>
 
 </details>
-  - 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.
 <details>
 <summary><b>Show answer</b></summary>
 
-<b> 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.</b>
+<b> 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.</b>
 
 </details>
 

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

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**
 

interface/HybridNew.h

@@ -111,7 +111,7 @@ class HybridNew : public LimitAlgo {
 
   std::pair<double,double> 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<double,double> eval(const RooStats::HypoTestResult &hcres, const RooAbsCollection & rVals) ;
   std::pair<double,double> eval(const RooStats::HypoTestResult &hcres, double rVal) ;
 

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 &)
 

src/AsymptoticLimits.cc

@@ -54,7 +54,7 @@ LimitAlgo("AsymptoticLimits specific options") {
         //("minimizerTolerance", boost::program_options::value<float>(&minimizerTolerance_)->default_value(minimizerTolerance_),  "Tolerance for minimizer used for profiling")
         //("minimizerStrategy",  boost::program_options::value<int>(&minimizerStrategy_)->default_value(minimizerStrategy_),      "Stragegy for minimizer")
         ("qtilde", boost::program_options::value<bool>(&qtilde_)->default_value(qtilde_),  "Allow only non-negative signal strengths (default is true).")
-        ("rule",    boost::program_options::value<std::string>(&rule_)->default_value(rule_),            "Rule to use: CLs, CLsplusb")
+        ("rule",    boost::program_options::value<std::string>(&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<std::string>(&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<std::pair<float,float> > 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<std::pair<float,float> > 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();