Skip to content
52 changes: 27 additions & 25 deletions bin/rabbit_fit.py
Original file line number Diff line number Diff line change
Expand Up @@ -276,41 +276,43 @@ def fit(args, fitter, ws, dofit=True):

if args.externalPostfit is not None:
fitter.load_fitresult(args.externalPostfit, args.externalPostfitResult)
else:
if dofit:
cb = fitter.minimize()

if cb is not None:
ws.add_loss_time_hist(cb.loss_history, cb.time_history)
if fitter.binByBinStat:
fitter._profile_beta()

if not args.noHessian:
# compute the covariance matrix and estimated distance to minimum
if dofit:
cb = fitter.minimize()

val, grad, hess = fitter.loss_val_grad_hess()
edmval, cov = edmval_cov(grad, hess)
if cb is not None:
ws.add_loss_time_hist(cb.loss_history, cb.time_history)

logger.info(f"edmval: {edmval}")
if not args.noHessian:
# compute the covariance matrix and estimated distance to minimum

val, grad, hess = fitter.loss_val_grad_hess()
edmval, cov = edmval_cov(grad, hess)
logger.info(f"edmval: {edmval}")

fitter.cov.assign(cov)
fitter.cov.assign(cov)

del cov
del cov

if fitter.binByBinStat and fitter.diagnostics:
# This is the estimated distance to minimum with respect to variations of
# the implicit binByBinStat nuisances beta at fixed parameter values.
# It should be near-zero by construction as long as the analytic profiling is
# correct
_, gradbeta, hessbeta = fitter.loss_val_grad_hess_beta()
edmvalbeta, covbeta = edmval_cov(gradbeta, hessbeta)
logger.info(f"edmvalbeta: {edmvalbeta}")
if fitter.binByBinStat and fitter.diagnostics:
# This is the estimated distance to minimum with respect to variations of
# the implicit binByBinStat nuisances beta at fixed parameter values.
# It should be near-zero by construction as long as the analytic profiling is
# correct
_, gradbeta, hessbeta = fitter.loss_val_grad_hess_beta()
edmvalbeta, covbeta = edmval_cov(gradbeta, hessbeta)
logger.info(f"edmvalbeta: {edmvalbeta}")

if args.doImpacts:
ws.add_impacts_hists(*fitter.impacts_parms(hess))
if args.doImpacts:
ws.add_impacts_hists(*fitter.impacts_parms(hess))

del hess
del hess

if args.globalImpacts:
ws.add_global_impacts_hists(*fitter.global_impacts_parms())
if args.globalImpacts:
ws.add_global_impacts_hists(*fitter.global_impacts_parms())

nllvalreduced = fitter.reduced_nll().numpy()

Expand Down
63 changes: 26 additions & 37 deletions bin/rabbit_plot_likelihood_scan.py
Original file line number Diff line number Diff line change
@@ -1,9 +1,7 @@
#!/usr/bin/env python3

import argparse
import os

import matplotlib.pyplot as plt
import mplhep as hep
import numpy as np

Expand All @@ -16,32 +14,6 @@
logger = None


def writeOutput(fig, outfile, extensions=[], postfix=None, args=None, meta_info=None):
name, _ = os.path.splitext(outfile)

if postfix:
name += f"_{postfix}"

for ext in extensions:
if ext[0] != ".":
ext = "." + ext
output = name + ext
print(f"Write output file {output}")
plt.savefig(output)

output = name.rsplit("/", 1)
output[1] = os.path.splitext(output[1])[0]
if len(output) == 1:
output = (None, *output)
if args is None and meta_info is None:
return
output_tools.write_logfile(
*output,
args=args,
meta_info=meta_info,
)


def parseArgs():
parser = argparse.ArgumentParser()
parser.add_argument(
Expand Down Expand Up @@ -73,6 +45,11 @@ def parseArgs():
default="./test",
help="Folder path for output",
)
parser.add_argument(
"--eoscp",
action="store_true",
help="Override use of xrdcp and use the mount instead",
)
parser.add_argument(
"-p", "--postfix", type=str, help="Postfix for output file name"
)
Expand Down Expand Up @@ -113,6 +90,13 @@ def parseArgs():
default=None,
help="x axis limits",
)
parser.add_argument(
"--ylim",
type=float,
nargs=2,
default=None,
help="y axis limits",
)
parser.add_argument("--titlePos", type=int, default=2, help="title position")
parser.add_argument(
"--config",
Expand All @@ -133,6 +117,7 @@ def plot_scan(
subtitle=None,
titlePos=0,
xlim=None,
ylim=None,
combine=None,
ylabel=r"$-2\,\Delta \log L$",
config={},
Expand All @@ -148,13 +133,15 @@ def plot_scan(

if xlim is None:
xlim = (min(x), max(x))
if ylim is None:
ylim = (min(y), max(y))

fig, ax = plot_tools.figure(
x,
xlabel,
ylabel,
xlim=xlim,
ylim=(min(y), max(y)), # logy=args.logy
ylim=ylim, # logy=args.logy
)

ax.axhline(y=1, color="gray", linestyle="--", alpha=0.5)
Expand Down Expand Up @@ -227,6 +214,7 @@ def plot_scan(

def main():
args = parseArgs()
outdir = output_tools.make_plot_dir(args.outpath, eoscp=args.eoscp)

global logger
logger = logging.setup_logger(__file__, args.verbose, args.noColorLogger)
Expand Down Expand Up @@ -280,19 +268,20 @@ def main():
subtitle=args.subtitle,
titlePos=args.titlePos,
xlim=args.xlim,
ylim=args.ylim,
config=config,
combine=(vals, nlls) if args.combine is not None else None,
no_hessian=args.noHessian,
)
os.makedirs(args.outpath, exist_ok=True)
outfile = os.path.join(args.outpath, f"nll_scan_{param}")
writeOutput(
fig,
outfile=outfile,
extensions=["png", "pdf"],
meta_info=meta,

to_join = [f"nll_scan_{param}", args.postfix]
outfile = "_".join(filter(lambda x: x, to_join))
plot_tools.save_pdf_and_png(outdir, outfile)
output_tools.write_index_and_log(
outdir,
outfile,
analysis_meta_info=meta,
args=args,
postfix=args.postfix,
)


Expand Down
Loading