Skip to content

Adding option to run expected limits for mu!=0 #1128

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
nucleosynthesis wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Adding option to run expected limits for mu!=0 #1128

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0765fc4
Adding option to run expected limits for mu!=0
Oct 13, 2025
e7d161b
Add some documentation
Oct 15, 2025
a8e70c8
clang format
Oct 15, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions docs/part3/commonstatsmethods.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,11 @@ hadd limits.root higgsCombine*.AsymptoticLimits.*

combine -M AsymptoticLimits realistic-counting-experiment.txt --getLimitFromGrid limits.root
```
## Expected limits assuming non-zero signal strengths

The median expected limit and quantiles can be calculated under Hypotheses other than at $\mu=0$. The Asymptotic expected limits can be calculated by setting the value `--signalStrengthForExpected mu0`. The Asimov dataset produced for the calculation will be generated for $\mu=$`mu0` instead of the usual $\mu=0$. Note that setting ` --signalStrengthForExpected 0` will give the same results as not setting the option, but will be slower as this forces a recreation of the Asimov dataset.
Copy link
Contributor

@coderabbitai coderabbitai bot Oct 15, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟡 Minor

Fix markdown formatting: remove extra space in code span.

There's an extra space before the -- in the code span. This triggers a markdown linter warning.

Apply this diff:

-The median expected limit and quantiles can be calculated under Hypotheses other than at $\mu=0$. The Asymptotic expected limits can be calculated by setting the value `--signalStrengthForExpected mu0`. The Asimov dataset produced for the calculation will be generated for $\mu=$`mu0` instead of the usual $\mu=0$. Note that setting ` --signalStrengthForExpected 0` will give the same results as not setting the option, but will be slower as this forces a recreation of the Asimov dataset.
+The median expected limit and quantiles can be calculated under Hypotheses other than at $\mu=0$. The Asymptotic expected limits can be calculated by setting the value `--signalStrengthForExpected mu0`. The Asimov dataset produced for the calculation will be generated for $\mu=$`mu0` instead of the usual $\mu=0$. Note that setting `--signalStrengthForExpected 0` will give the same results as not setting the option, but will be slower as this forces a recreation of the Asimov dataset.

As per static analysis hints.

📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
The median expected limit and quantiles can be calculated under Hypotheses other than at $\mu=0$. The Asymptotic expected limits can be calculated by setting the value `--signalStrengthForExpected mu0`. The Asimov dataset produced for the calculation will be generated for $\mu=$`mu0` instead of the usual $\mu=0$. Note that setting ` --signalStrengthForExpected 0` will give the same results as not setting the option, but will be slower as this forces a recreation of the Asimov dataset.
The median expected limit and quantiles can be calculated under Hypotheses other than at $\mu=0$. The Asymptotic expected limits can be calculated by setting the value `--signalStrengthForExpected mu0`. The Asimov dataset produced for the calculation will be generated for $\mu=$`mu0` instead of the usual $\mu=0`. Note that setting `--signalStrengthForExpected 0` will give the same results as not setting the option, but will be slower as this forces a recreation of the Asimov dataset.
🧰 Tools
🪛 markdownlint-cli2 (0.18.1)

95-95: Spaces inside code span elements

(MD038, no-space-in-code)

🤖 Prompt for AI Agents 
In docs/part3/commonstatsmethods.md around line 95, there's an extra space
inside the inline code span before the option (` ` --signalStrengthForExpected
0` ), which triggers a markdown linter warning; remove the leading space so the
span reads `--signalStrengthForExpected 0` (update the single inline code
occurrence at that line to eliminate the extra space).


Be aware that while expected values are derived from the Asimov created at the specified signal strength, setting this option *does not* alter the definition of $CL_s$ in that the denominator is still defined with respect to the no signal ($\mu=0$) hypothesis.

## Asymptotic Significances

Expand Down
5 changes: 4 additions & 1 deletion interface/AsymptoticLimits.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@ class AsymptoticLimits : public LimitAlgo {
//static int minimizerStrategy_;

static double rValue_;
static double signalStrengthForExpected_;
Copy link
Contributor

@github-actions github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: invalid case style for variable 'signalStrengthForExpected_' [readability-identifier-naming]

Suggested change
static double signalStrengthForExpected_;
static double signalStrengthForExpected;


static bool strictBounds_;

Expand All @@ -66,7 +67,7 @@ class AsymptoticLimits : public LimitAlgo {

float calculateLimitFromGrid(RooRealVar *, double, double);

RooAbsData *asimovDataset(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data);
RooAbsData *asimovDataset(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, bool overwrite=false);
double getCLs(RooRealVar &r, double rVal, bool getAlsoExpected=false, double *limit=0, double *limitErr=0);

TFile *gridFile_;
Expand All @@ -75,6 +76,8 @@ class AsymptoticLimits : public LimitAlgo {
double readMU_;
bool doCLs_;

bool doNonStandardAsimov_;
Copy link
Contributor

@github-actions github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: invalid case style for private member 'doNonStandardAsimov_' [readability-identifier-naming]

Suggested change
bool doNonStandardAsimov_;
bool m_doNonStandardAsimov;


};

#endif
55 changes: 45 additions & 10 deletions src/AsymptoticLimits.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@ std::string AsymptoticLimits::minosAlgo_ = "stepping";
//float AsymptoticLimits::minimizerTolerance_ = 0.01;
//int AsymptoticLimits::minimizerStrategy_ = 0;
double AsymptoticLimits::rValue_ = 1.0;
double AsymptoticLimits::signalStrengthForExpected_ = 0.0;
bool AsymptoticLimits::strictBounds_ = false;

RooAbsData * AsymptoticLimits::asimovDataset_ = nullptr;
Expand All @@ -46,6 +47,7 @@ LimitAlgo("AsymptoticLimits specific options") {
options_.add_options()
("rAbsAcc", boost::program_options::value<double>(&rAbsAccuracy_)->default_value(rAbsAccuracy_), "Absolute accuracy on r to reach to terminate the scan")
("rRelAcc", boost::program_options::value<double>(&rRelAccuracy_)->default_value(rRelAccuracy_), "Relative accuracy on r to reach to terminate the scan")
("signalStrengthForExpected", boost::program_options::value<double>(&signalStrengthForExpected_)->default_value(signalStrengthForExpected_), "Signal strength for expected limits (0=background only, default is 0)")
("run", boost::program_options::value<std::string>(&what_)->default_value(what_), "What to run: both (default), observed, expected, blind.")
("singlePoint", boost::program_options::value<double>(&rValue_), "Just compute CLs for the given value of r")
//("minimizerAlgo", boost::program_options::value<std::string>(&minimizerAlgo_)->default_value(minimizerAlgo_), "Choice of minimizer used for profiling (Minuit vs Minuit2)")
Expand Down Expand Up @@ -89,7 +91,8 @@ void AsymptoticLimits::applyOptions(const boost::program_options::variables_map
limitsTree_->SetBranchAddress("limit",&readCL_);
limitsTree_->SetBranchAddress("r",&readMU_);
}


doNonStandardAsimov_ = vm.count("signalStrengthForExpected") && !vm["signalStrengthForExpected"].defaulted();
}

void AsymptoticLimits::applyDefaultOptions() {
Expand Down Expand Up @@ -151,7 +154,12 @@ bool AsymptoticLimits::runLimit(RooWorkspace *w, RooStats::ModelConfig *mc_s, Ro
}

w->loadSnapshot("clean");
RooAbsData &asimov = *asimovDataset(w, mc_s, mc_b, data);

// Bit of a waste of time if we are not using a non-standard value
double tmpsexp = signalStrengthForExpected_;
signalStrengthForExpected_ = 0.0;
RooAbsData &asimov = *asimovDataset(w, mc_s, mc_b, data, /*overwrite=*/doNonStandardAsimov_);
signalStrengthForExpected_ = tmpsexp;
w->loadSnapshot("clean");

r->setConstant(false);
Expand Down Expand Up @@ -381,7 +389,7 @@ std::vector<std::pair<float,float> > AsymptoticLimits::runLimitExpected(RooWorks
}

// 2) get asimov dataset
RooAbsData *asimov = asimovDataset(w, mc_s, mc_b, data);
RooAbsData *asimov = asimovDataset(w, mc_s, mc_b, data, /*overwrite=*/ doNonStandardAsimov_);

// 2b) load asimov global observables
if (params_.get() == 0) params_.reset(mc_s->GetPdf()->getParameters(data));
Expand All @@ -408,7 +416,7 @@ std::vector<std::pair<float,float> > AsymptoticLimits::runLimitExpected(RooWorks
std::unique_ptr<RooFitResult> res(minim.save());
res->Print("V");
}
if (r->getVal()/r->getMax() > 1e-3) {
if (r->getVal()/r->getMax() > 1e-3 && !doNonStandardAsimov_) {
if (verbose) {
CombineLogger::instance().log("AsymptoticLimits.cc",__LINE__,std::string(Form("[WARNING] Best fit of asimov dataset is at %s = %f (%f times %sMax), while it should be at zero",r->GetName(), r->getVal(), r->getVal()/r->getMax(), r->GetName())),__func__);
}
Expand Down Expand Up @@ -463,8 +471,33 @@ float AsymptoticLimits::findExpectedLimitFromCrossing(RooAbsReal &nll, RooRealVa
// crossing value of mu that gives the specified qmu in the above.
// note that as in CCGV the asymptotic formula for upper limits in qmu and qmutilde are identical so can use qmu here.

// only need to modify the value of pb compared to the typical case where mu'=0 for the asimov dataset
double pb_expected = pb;
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);

// Things get tricker here so first, we use the Asimov value of the test stat and plug it into
if (doNonStandardAsimov_) {

// std::cout << "Using non-standard asimov dataset with signal strength " << signalStrengthForExpected_ << std::endl;
// Need to find q(0)
r->setVal(0); r->setConstant(true);
CascadeMinimizer minim2(nll, CascadeMinimizer::Constrained);
if (hasDiscreteParams_) minim2.minimize(verbose-2);
else minim2.improve(verbose-2);
double q_At_0 = 2*nll.getVal();
r->setVal(signalStrengthForExpected_);
if (hasDiscreteParams_) minim2.minimize(verbose-2);
else minim2.improve(verbose-2);
double q_At_muA = 2*nll.getVal();
r->setConstant(false);

double qMuAsimov = q_At_0-q_At_muA;
double N_for_expected = N+ROOT::Math::sqrt(qMuAsimov);
Copy link
Contributor

@github-actions github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "ROOT::Math::sqrt" is directly included [misc-include-cleaner]

src/AsymptoticLimits.cc:3:

- #include <RooRealVar.h>
+ #include <Math/Math.h>
+ #include <RooRealVar.h>

pb_expected = ROOT::Math::normal_cdf(N_for_expected, 1.0);
// std::cout << " --> this gives pb = " << pb_expected << " (N=" << N_for_expected << ")" << std::endl;
}

double errorlevel = 0.5 * pow(N+ROOT::Math::normal_quantile_c((doCLs_ ? pb_expected:1.)*(1-cl),1.0), 2);
Copy link
Contributor

@github-actions github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no header providing "ROOT::Math::normal_quantile_c" is directly included [misc-include-cleaner]

    double errorlevel = 0.5 * pow(N+ROOT::Math::normal_quantile_c((doCLs_ ? pb_expected:1.)*(1-cl),1.0), 2);
                                                ^

int minosStat = -1;
if (minosAlgo_ == "minos") {
double rMax0 = r->getMax();
Expand Down Expand Up @@ -730,13 +763,15 @@ float AsymptoticLimits::calculateLimitFromGrid(RooRealVar *r , double quantile,
return rlim;
}

RooAbsData * AsymptoticLimits::asimovDataset(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data) {
RooAbsData * AsymptoticLimits::asimovDataset(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, bool overwrite) {
// Do this only once
// if (w->data("_Asymptotic_asimovDataset_") != 0) {
// return w->data("_Asymptotic_asimovDataset_");
// }

if (asimovDataset_) {

if (asimovDataset_ && !overwrite) {
//std::cerr << "Reusing asimov dataset" << std::endl;
return asimovDataset_;
}
// snapshot data global observables
Expand All @@ -748,9 +783,9 @@ RooAbsData * AsymptoticLimits::asimovDataset(RooWorkspace *w, RooStats::ModelCon
gobs.snapshot(snapGlobalObsData);
}
// get asimov dataset and global observables
asimovDataset_ = (noFitAsimov_ ? asimovutils::asimovDatasetNominal(mc_s, 0.0, verbose) :
asimovutils::asimovDatasetWithFit(mc_s, data, snapGlobalObsAsimov,!bypassFrequentistFit_, 0.0, verbose));
asimovDataset_->SetName("_Asymptotic_asimovDataset_");
asimovDataset_ = (noFitAsimov_ ? asimovutils::asimovDatasetNominal(mc_s, signalStrengthForExpected_, verbose) :
asimovutils::asimovDatasetWithFit(mc_s, data, snapGlobalObsAsimov,!bypassFrequentistFit_, signalStrengthForExpected_, verbose));
asimovDataset_->SetName(Form("_Asymptotic_asimovDataset_%d_%g",doNonStandardAsimov_,signalStrengthForExpected_)); // in case we want to keep multiple asimov datasets in the same workspace
Copy link
Contributor

@github-actions github-actions bot Oct 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: do not call c-style vararg functions [cppcoreguidelines-pro-type-vararg]

    asimovDataset_->SetName(Form("_Asymptotic_asimovDataset_%d_%g",doNonStandardAsimov_,signalStrengthForExpected_)); // in case we want to keep multiple asimov datasets in the same workspace
                            ^

// w->import(*asimovData); // I'm assuming the Workspace takes ownership. Might be false.
// delete asimovData; // ^^^^^^^^----- now assuming that the workspace clones.
return asimovDataset_;
Expand Down
Loading