From ad209f26a19b937fed0955001401f0e45f603c74 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Thu, 5 Dec 2024 14:24:12 -0500 Subject: [PATCH 01/11] backwards_ecal: produce file lists --- benchmarks/backwards_ecal/Snakefile | 11 +++++++---- benchmarks/backwards_ecal/backwards_ecal.org | 10 +++++++--- benchmarks/backwards_ecal/config.yml | 2 +- 3 files changed, 15 insertions(+), 8 deletions(-) diff --git a/benchmarks/backwards_ecal/Snakefile b/benchmarks/backwards_ecal/Snakefile index 9b0b1559..8b94d911 100644 --- a/benchmarks/backwards_ecal/Snakefile +++ b/benchmarks/backwards_ecal/Snakefile @@ -64,14 +64,17 @@ exec env DETECTOR_CONFIG={params.DETECTOR_CONFIG} \ """ -rule backwards_ecal_recon_many: +rule backwards_ecal_local_sim_list: input: expand( "sim_output/backwards_ecal/{{DETECTOR_CONFIG}}/{{PARTICLE}}/{{ENERGY}}/{{PHASE_SPACE}}/{{PARTICLE}}_{{ENERGY}}_{{PHASE_SPACE}}.{INDEX:04d}.eicrecon.tree.edm4eic.root", INDEX=range(20), ), output: - touch("sim_output/backwards_ecal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/flag"), + "listing/backwards_ecal/local/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}.lst", + run: + with open(output[0], "wt") as fp: + fp.write("\n".join(input)) DETECTOR_CONFIG=os.environ["DETECTOR_CONFIG"] @@ -79,7 +82,7 @@ DETECTOR_CONFIG=os.environ["DETECTOR_CONFIG"] rule backwards_ecal: input: expand( - "sim_output/backwards_ecal/" + DETECTOR_CONFIG + "/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/flag", + "listing/backwards_ecal/local/" + DETECTOR_CONFIG + "/{PARTICLE}/{ENERGY}/{PHASE_SPACE}.lst", PARTICLE=["pi-", "e-"], ENERGY=[ "100MeV", @@ -103,7 +106,7 @@ env \ MATPLOTLIBRC={input.matplotlibrc} \ DETECTOR_CONFIG=""" + DETECTOR_CONFIG + """ \ PLOT_TITLE=""" + DETECTOR_CONFIG + """ \ -INPUT_PATH_FORMAT=sim_output/backwards_ecal/""" + DETECTOR_CONFIG + """/{{particle}}/{{energy}}/130to177deg/{{particle}}_{{energy}}_130to177deg.{{ix:04d}}.eicrecon.tree.edm4eic.root \ +INPUT_PATH_FORMAT=listing/backwards_ecal/local/""" + DETECTOR_CONFIG + """/{{particle}}/{{energy}}/130to177deg.lst \ OUTPUT_DIR={output} \ python {input.script} """ diff --git a/benchmarks/backwards_ecal/backwards_ecal.org b/benchmarks/backwards_ecal/backwards_ecal.org index 27f1094c..1028e318 100644 --- a/benchmarks/backwards_ecal/backwards_ecal.org +++ b/benchmarks/backwards_ecal/backwards_ecal.org @@ -57,7 +57,6 @@ setup_presentation_style() DETECTOR_CONFIG=os.environ.get("DETECTOR_CONFIG") PLOT_TITLE=os.environ.get("PLOT_TITLE") INPUT_PATH_FORMAT=os.environ.get("INPUT_PATH_FORMAT", "EPIC/RECO/24.04.0/epic_craterlake/SINGLE/{particle}/{energy}/130to177deg/{particle}_{energy}_130to177deg.{ix:04d}.eicrecon.tree.edm4eic.root") -INPUT_PATH_INDEX_RANGE=list(range(20)) output_dir=Path(os.environ.get("OUTPUT_DIR", "./")) output_dir.mkdir(parents=True, exist_ok=True) @@ -102,13 +101,18 @@ filter_name = [ pi_eval = {} e_eval = {} +def readlist(path): + with open(path, "rt") as fp: + paths = [line.rstrip() for line in fp.readlines()] + return paths + for energy in energies: pi_eval[energy] = filter_pointing(uproot.concatenate( - {INPUT_PATH_FORMAT.format(particle="pi-", energy=energy, ix=ix): "events" for ix in INPUT_PATH_INDEX_RANGE}, + {path: "events" for path in readlist(INPUT_PATH_FORMAT.format(particle="pi-", energy=energy))}, filter_name=filter_name, )) e_eval[energy] = filter_pointing(uproot.concatenate( - {INPUT_PATH_FORMAT.format(particle="e-", energy=energy, ix=ix): "events" for ix in INPUT_PATH_INDEX_RANGE}, + {path: "events" for path in readlist(INPUT_PATH_FORMAT.format(particle="e-", energy=energy))}, filter_name=filter_name, )) #+end_src diff --git a/benchmarks/backwards_ecal/config.yml b/benchmarks/backwards_ecal/config.yml index 7048f8e7..8f2967f6 100644 --- a/benchmarks/backwards_ecal/config.yml +++ b/benchmarks/backwards_ecal/config.yml @@ -16,7 +16,7 @@ sim:backwards_ecal: ] script: - | - snakemake $SNAKEMAKE_FLAGS --cores 5 sim_output/backwards_ecal/${DETECTOR_CONFIG}/${PARTICLE}/${MOMENTUM}/130to177deg/flag + snakemake $SNAKEMAKE_FLAGS --cores 5 listing/backwards_ecal/local/${DETECTOR_CONFIG}/${PARTICLE}/${MOMENTUM}/130to177deg.lst bench:backwards_ecal: extends: .det_benchmark From d2b9416d003ef7edd8097c7c4bdcf2c46972667b Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 7 Dec 2024 00:54:07 -0500 Subject: [PATCH 02/11] backwards_ecal: add support for campaign processing --- benchmarks/backwards_ecal/Snakefile | 58 ++++++++++++++++++++++++++-- benchmarks/backwards_ecal/config.yml | 4 +- 2 files changed, 57 insertions(+), 5 deletions(-) diff --git a/benchmarks/backwards_ecal/Snakefile b/benchmarks/backwards_ecal/Snakefile index 8b94d911..302cd1b6 100644 --- a/benchmarks/backwards_ecal/Snakefile +++ b/benchmarks/backwards_ecal/Snakefile @@ -77,12 +77,64 @@ rule backwards_ecal_local_sim_list: fp.write("\n".join(input)) +if False: + rule backwards_ecal_campaign_sim_list: + output: + "listing/backwards_ecal/{CAMPAIGN}/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}.lst", + params: + search_path=lambda wildcards: f"EPIC/RECO/{wildcards.CAMPAIGN}/epic_craterlake/SINGLE/{wildcards.PARTICLE}/{wildcards.ENERGY}/{wildcards.PHASE_SPACE}/", + shell: """ + xrdfs root://dtn-eic.jlab.org/ ls /work/eic2/{params.search_path} \ + | awk '{{ print "root://dtn-eic.jlab.org/"$1; }}' \ + | sort \ + > {output} + if [ ! -s {output} ]; then + echo "Got an empty file listing for path \"\"" + exit 1 + fi + """ +else: + checkpoint backwards_ecal_campaign_sim_list_checkpoint: + output: + "listing/backwards_ecal/{CAMPAIGN}/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}.lst.orig", + params: + search_path=lambda wildcards: f"EPIC/RECO/{wildcards.CAMPAIGN}/epic_craterlake/SINGLE/{wildcards.PARTICLE}/{wildcards.ENERGY}/{wildcards.PHASE_SPACE}/", + shell: """ + xrdfs root://dtn-eic.jlab.org/ ls /work/eic2/{params.search_path} \ + | sed -e 's#^/work/eic2/##' \ + | sort \ + > {output} + if [ ! -s {output} ]; then + echo "Got an empty file listing for path \"\"" + exit 1 + fi + """ + + def get_backwards_ecal_campaign_sim_list(wildcards): + with checkpoints.backwards_ecal_campaign_sim_list_checkpoint.get(**wildcards).output[0].open() as fp: + return [line.rstrip() for line in fp.readlines()] + + rule backwards_ecal_campaign_sim_list: + input: + # depend on paths from the file list + get_backwards_ecal_campaign_sim_list, + orig_list="listing/backwards_ecal/{CAMPAIGN}/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}.lst.orig", + output: + "listing/backwards_ecal/{CAMPAIGN}/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}.lst", + shell: """ + cp {input.orig_list} {output} + """ + + +ruleorder: backwards_ecal_local_sim_list > backwards_ecal_campaign_sim_list + + DETECTOR_CONFIG=os.environ["DETECTOR_CONFIG"] rule backwards_ecal: input: expand( - "listing/backwards_ecal/local/" + DETECTOR_CONFIG + "/{PARTICLE}/{ENERGY}/{PHASE_SPACE}.lst", + "listing/backwards_ecal/{{CAMPAIGN}}/" + DETECTOR_CONFIG + "/{PARTICLE}/{ENERGY}/{PHASE_SPACE}.lst", PARTICLE=["pi-", "e-"], ENERGY=[ "100MeV", @@ -99,14 +151,14 @@ rule backwards_ecal: matplotlibrc=".matplotlibrc", script="benchmarks/backwards_ecal/backwards_ecal.py", output: - directory("results/backwards_ecal") + directory("results/backwards_ecal/{CAMPAIGN}/") shell: """ env \ MATPLOTLIBRC={input.matplotlibrc} \ DETECTOR_CONFIG=""" + DETECTOR_CONFIG + """ \ PLOT_TITLE=""" + DETECTOR_CONFIG + """ \ -INPUT_PATH_FORMAT=listing/backwards_ecal/local/""" + DETECTOR_CONFIG + """/{{particle}}/{{energy}}/130to177deg.lst \ +INPUT_PATH_FORMAT=listing/backwards_ecal/{wildcards.CAMPAIGN}/""" + DETECTOR_CONFIG + """/{{particle}}/{{energy}}/130to177deg.lst \ OUTPUT_DIR={output} \ python {input.script} """ diff --git a/benchmarks/backwards_ecal/config.yml b/benchmarks/backwards_ecal/config.yml index 8f2967f6..75d8f464 100644 --- a/benchmarks/backwards_ecal/config.yml +++ b/benchmarks/backwards_ecal/config.yml @@ -26,7 +26,7 @@ bench:backwards_ecal: script: - export PYTHONUSERBASE=$LOCAL_DATA_PATH/deps - pip install -r benchmarks/backwards_ecal/requirements.txt - - snakemake $SNAKEMAKE_FLAGS --cores 1 backwards_ecal + - snakemake $SNAKEMAKE_FLAGS --cores 1 results/backwards_ecal/local collect_results:backwards_ecal: extends: .det_benchmark @@ -36,5 +36,5 @@ collect_results:backwards_ecal: script: - ls -lrht - mv results{,_save}/ # move results directory out of the way to preserve it - - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output backwards_ecal + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/backwards_ecal/local - mv results{_save,}/ From 0bc6a46718b9fd1b197e5cb1d0453f2bfcc5da7d Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 7 Dec 2024 01:18:48 -0500 Subject: [PATCH 03/11] backwards_ecal: update YR requirement based on correction from Carlos --- benchmarks/backwards_ecal/backwards_ecal.org | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/benchmarks/backwards_ecal/backwards_ecal.org b/benchmarks/backwards_ecal/backwards_ecal.org index 1028e318..c917ce1f 100644 --- a/benchmarks/backwards_ecal/backwards_ecal.org +++ b/benchmarks/backwards_ecal/backwards_ecal.org @@ -221,10 +221,11 @@ for clf_label, sigma_rel_FWHM_cb in sigmas_rel_FWHM_cb.items(): lw=0.5, label=f"{clf_label}, ${np.ceil(stochastic * 10) / 10:.1f}\% / \sqrt{{E}} \oplus {np.ceil(constant * 10) / 10:.1f}\%$", ) -plt.plot( +plt.fill_between( energy_values, - np.sqrt((1 / energy_values) ** 2 + (1 / np.sqrt(energy_values)) ** 2 + 1 ** 2), - color="black", label=r"YR requirement $1\% / E \oplus 2.5\% / \sqrt{E} \oplus 1\%$", + np.sqrt((2 / np.sqrt(energy_values)) ** 2 + 1 ** 2), + np.sqrt((2 / np.sqrt(energy_values)) ** 2 + 3 ** 2), + alpha=0.2, color="black", label=r"YR requirement $2\% / \sqrt{E} \oplus (1-3)\%$", ) plt.title(INPUT_PATH_FORMAT) plt.legend() From d601cb05fda967f37d35476e89153233b2717537 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 7 Dec 2024 01:21:15 -0500 Subject: [PATCH 04/11] backwards_ecal: drop above plot titles --- benchmarks/backwards_ecal/Snakefile | 1 - benchmarks/backwards_ecal/backwards_ecal.org | 8 -------- 2 files changed, 9 deletions(-) diff --git a/benchmarks/backwards_ecal/Snakefile b/benchmarks/backwards_ecal/Snakefile index 302cd1b6..9bcbb9d6 100644 --- a/benchmarks/backwards_ecal/Snakefile +++ b/benchmarks/backwards_ecal/Snakefile @@ -157,7 +157,6 @@ rule backwards_ecal: env \ MATPLOTLIBRC={input.matplotlibrc} \ DETECTOR_CONFIG=""" + DETECTOR_CONFIG + """ \ -PLOT_TITLE=""" + DETECTOR_CONFIG + """ \ INPUT_PATH_FORMAT=listing/backwards_ecal/{wildcards.CAMPAIGN}/""" + DETECTOR_CONFIG + """/{{particle}}/{{energy}}/130to177deg.lst \ OUTPUT_DIR={output} \ python {input.script} diff --git a/benchmarks/backwards_ecal/backwards_ecal.org b/benchmarks/backwards_ecal/backwards_ecal.org index c917ce1f..34be9b7b 100644 --- a/benchmarks/backwards_ecal/backwards_ecal.org +++ b/benchmarks/backwards_ecal/backwards_ecal.org @@ -122,8 +122,6 @@ for energy in energies: #+begin_src jupyter-python fig, axs = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6)) -fig.suptitle(PLOT_TITLE) - axs = np.ravel(np.array(axs)) sigmas_rel_FWHM_cb = {} @@ -227,7 +225,6 @@ plt.fill_between( np.sqrt((2 / np.sqrt(energy_values)) ** 2 + 3 ** 2), alpha=0.2, color="black", label=r"YR requirement $2\% / \sqrt{E} \oplus (1-3)\%$", ) -plt.title(INPUT_PATH_FORMAT) plt.legend() plt.xlabel("Energy, GeV", loc="right") plt.ylabel(r"$\sigma_{E} / E$ derived from FWHM, %", loc="top") @@ -243,10 +240,6 @@ fig, axs = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6)) fig_log, axs_log = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6)) fig_roc, axs_roc = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6)) -fig.suptitle(PLOT_TITLE) -fig_log.suptitle(PLOT_TITLE) -fig_roc.suptitle(PLOT_TITLE) - axs = np.ravel(np.array(axs)) axs_log = np.ravel(np.array(axs_log)) axs_roc = np.ravel(np.array(axs_roc)) @@ -317,7 +310,6 @@ for clf_label, roc in rocs.items(): label=f"{clf_label}", ) plt.yscale("log") -plt.title(INPUT_PATH_FORMAT) plt.legend() plt.xlabel("Energy, GeV") plt.ylabel("Pion rejection at 95%") From a4bf3e64fdb602d81e48a7934b63b6fa986e6f7c Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 7 Dec 2024 01:30:30 -0500 Subject: [PATCH 05/11] backwards_ecal: adjust colors, don't connect points for resolution plot --- benchmarks/backwards_ecal/backwards_ecal.org | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/backwards_ecal/backwards_ecal.org b/benchmarks/backwards_ecal/backwards_ecal.org index 34be9b7b..7f49184e 100644 --- a/benchmarks/backwards_ecal/backwards_ecal.org +++ b/benchmarks/backwards_ecal/backwards_ecal.org @@ -209,12 +209,12 @@ for clf_label, sigma_rel_FWHM_cb in sigmas_rel_FWHM_cb.items(): energy_values, sigma_over_e, marker=".", + ls="none", label=f"{clf_label}" ) plt.plot( energy_values[cond], f(energy_values[cond], *par), - color="black", ls="--", lw=0.5, label=f"{clf_label}, ${np.ceil(stochastic * 10) / 10:.1f}\% / \sqrt{{E}} \oplus {np.ceil(constant * 10) / 10:.1f}\%$", From afe1d17b1037247cf5484d0c13e96a96e1c65e0f Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 7 Dec 2024 01:31:12 -0500 Subject: [PATCH 06/11] backwards_ecal: fix ylim top for resolution --- benchmarks/backwards_ecal/backwards_ecal.org | 1 + 1 file changed, 1 insertion(+) diff --git a/benchmarks/backwards_ecal/backwards_ecal.org b/benchmarks/backwards_ecal/backwards_ecal.org index 7f49184e..fb079216 100644 --- a/benchmarks/backwards_ecal/backwards_ecal.org +++ b/benchmarks/backwards_ecal/backwards_ecal.org @@ -225,6 +225,7 @@ plt.fill_between( np.sqrt((2 / np.sqrt(energy_values)) ** 2 + 3 ** 2), alpha=0.2, color="black", label=r"YR requirement $2\% / \sqrt{E} \oplus (1-3)\%$", ) +plt.ylim(top=10.) plt.legend() plt.xlabel("Energy, GeV", loc="right") plt.ylabel(r"$\sigma_{E} / E$ derived from FWHM, %", loc="top") From 7fdc56243485e78c8e8f009dc6bdf7d2731feb28 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 7 Dec 2024 01:53:46 -0500 Subject: [PATCH 07/11] backwards_ecal: drop full calorimeter hit sum plot --- benchmarks/backwards_ecal/backwards_ecal.org | 177 +++++++++---------- 1 file changed, 80 insertions(+), 97 deletions(-) diff --git a/benchmarks/backwards_ecal/backwards_ecal.org b/benchmarks/backwards_ecal/backwards_ecal.org index fb079216..ec359464 100644 --- a/benchmarks/backwards_ecal/backwards_ecal.org +++ b/benchmarks/backwards_ecal/backwards_ecal.org @@ -91,10 +91,8 @@ energies = [ "10GeV", "20GeV", ] - filter_name = [ "MCParticles.*", - "*EcalEndcapNRecHits*", "*EcalEndcapNClusters*", ] @@ -128,61 +126,53 @@ sigmas_rel_FWHM_cb = {} fractions_below = {} for ix, energy in enumerate(energies): - for use_clusters in [False, True]: - energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3")) - if use_clusters: - clf_label = "leading cluster" - else: - clf_label = "sum all hits" - def clf(events): - if use_clusters: - return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value - else: - return ak.sum(events["EcalEndcapNRecHits.energy"], axis=-1) / energy_value - e_pred = clf(e_eval[energy]) - - plt.sca(axs[ix]) - counts, bins, patches = plt.hist(e_pred, weights=np.full_like(e_pred, 1.0 / ak.num(e_pred, axis=0)), bins=np.linspace(0.01, 1.01, 101), label=rf"$e^-$ {clf_label}", hatch=None if use_clusters else r"xxx", alpha=0.8 if use_clusters else 1.) - plt.title(f"{energy}") + energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3")) + clf_label = "leading cluster" + def clf(events): + return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value + e_pred = clf(e_eval[energy]) + + plt.sca(axs[ix]) + counts, bins, patches = plt.hist(e_pred, weights=np.full_like(e_pred, 1.0 / ak.num(e_pred, axis=0)), bins=np.linspace(0.01, 1.01, 101), label=rf"$e^-$ {clf_label}") + plt.title(f"{energy}") + + e_over_p = (bins[1:] + bins[:-1]) / 2 + import scipy.stats + def f(x, n, beta, m, loc, scale): + return n * scipy.stats.crystalball.pdf(x, beta, m, loc, scale) + p0 = (np.sum(counts[10:]), 2., 3., 0.95, 0.05) - e_over_p = (bins[1:] + bins[:-1]) / 2 - import scipy.stats - def f(x, n, beta, m, loc, scale): - return n * scipy.stats.crystalball.pdf(x, beta, m, loc, scale) - p0 = (np.sum(counts[10:]), 2., 3., 0.95, 0.05) - - try: - import scipy.optimize - par, pcov = scipy.optimize.curve_fit(f, e_over_p[5:], counts[5:], p0=p0, maxfev=10000) - except RuntimeError: - par = None - plt.plot(e_over_p, f(e_over_p, *par), label=rf"Crystal Ball fit", color="tab:green" if use_clusters else "green", lw=0.8) - - def summarize_fit(par): - _, _, _, loc_cb, scale_cb = par - # Calculate FWHM - y_max = np.max(f(np.linspace(0., 1., 100), *par)) - f_prime = lambda x: f(x, *par) - y_max / 2 - x_plus, = scipy.optimize.root(f_prime, loc_cb + scale_cb).x - x_minus, = scipy.optimize.root(f_prime, loc_cb - scale_cb).x - color = "cyan" if use_clusters else "orange" - plt.axvline(x_minus, ls="--", lw=0.75, color=patches[0].get_facecolor(), label=r"$\mu - $FWHM") - plt.axvline(x_plus, ls=":", lw=0.75, color=patches[0].get_facecolor(), label=r"$\mu + $FWHM") - fwhm = (x_plus - x_minus) / loc_cb - sigma_rel_FWHM_cb = fwhm / 2 / np.sqrt(2 * np.log(2)) - - cutoff_x = loc_cb - fwhm - fraction_below = np.sum(counts[e_over_p < cutoff_x]) / ak.num(e_pred, axis=0) - - return sigma_rel_FWHM_cb, fraction_below - - sigma_rel_FWHM_cb, fraction_below = summarize_fit(par) - sigmas_rel_FWHM_cb.setdefault(clf_label, {})[energy] = sigma_rel_FWHM_cb - fractions_below.setdefault(clf_label, {})[energy] = fraction_below + try: + import scipy.optimize + par, pcov = scipy.optimize.curve_fit(f, e_over_p[5:], counts[5:], p0=p0, maxfev=10000) + except RuntimeError: + par = None + plt.plot(e_over_p, f(e_over_p, *par), label=rf"Crystal Ball fit", color="tab:green", lw=0.8) - plt.legend() - plt.xlabel("$E/p$", loc="right") - plt.ylabel("Event yield", loc="top") + def summarize_fit(par): + _, _, _, loc_cb, scale_cb = par + # Calculate FWHM + y_max = np.max(f(np.linspace(0., 1., 100), *par)) + f_prime = lambda x: f(x, *par) - y_max / 2 + x_plus, = scipy.optimize.root(f_prime, loc_cb + scale_cb).x + x_minus, = scipy.optimize.root(f_prime, loc_cb - scale_cb).x + plt.axvline(x_minus, ls="--", lw=0.75, color=patches[0].get_facecolor(), label=r"$\mu - $FWHM") + plt.axvline(x_plus, ls=":", lw=0.75, color=patches[0].get_facecolor(), label=r"$\mu + $FWHM") + fwhm = (x_plus - x_minus) / loc_cb + sigma_rel_FWHM_cb = fwhm / 2 / np.sqrt(2 * np.log(2)) + + cutoff_x = loc_cb - fwhm + fraction_below = np.sum(counts[e_over_p < cutoff_x]) / ak.num(e_pred, axis=0) + + return sigma_rel_FWHM_cb, fraction_below + + sigma_rel_FWHM_cb, fraction_below = summarize_fit(par) + sigmas_rel_FWHM_cb.setdefault(clf_label, {})[energy] = sigma_rel_FWHM_cb + fractions_below.setdefault(clf_label, {})[energy] = fraction_below + + plt.legend() + plt.xlabel("$E/p$", loc="right") + plt.ylabel("Event yield", loc="top") fig.savefig(output_dir / f"resolution_plots.pdf", bbox_inches="tight") fig.savefig(output_dir / f"resolution_plots.png", bbox_inches="tight") @@ -248,49 +238,42 @@ axs_roc = np.ravel(np.array(axs_roc)) rocs = {} for ix, energy in enumerate(energies): - for use_clusters in [False, True]: - energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3")) - if use_clusters: - clf_label = "leading cluster" - else: - clf_label = "sum all hits" - def clf(events): - if use_clusters: - return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value - else: - return ak.sum(events["EcalEndcapNRecHits.energy"], axis=-1) / energy_value - e_pred = clf(e_eval[energy]) - pi_pred = clf(pi_eval[energy]) - - for do_log, ax in [(False, axs[ix]), (True, axs_log[ix])]: - plt.sca(ax) - plt.hist(e_pred, weights=np.full_like(e_pred, 1.0 / ak.num(e_pred, axis=0)), bins=np.linspace(0., 1.01, 101), label=rf"$e^-$ {clf_label}", hatch=None if use_clusters else r"xxx", alpha=0.8 if use_clusters else 1.) - plt.hist(pi_pred, weights=np.full_like(pi_pred, 1.0 / ak.num(pi_pred, axis=0)), bins=np.linspace(0., 1.01, 101), label=rf"$\pi^-$ {clf_label}", histtype="step") - plt.title(f"{energy}") - plt.legend() - plt.xlabel("Classifier output") - plt.ylabel("Event yield") - if do_log: - plt.yscale("log") - - plt.sca(axs_roc[ix]) - fpr, tpr, _ = roc_curve( - np.concatenate([np.ones_like(e_pred), np.zeros_like(pi_pred)]), - np.concatenate([e_pred, pi_pred]), - ) - cond = fpr != 0 # avoid infinite rejection (region of large uncertainty) - cond &= tpr != 1 # avoid linear interpolation (region of large uncertainty) - def mk_interp(tpr, fpr): - def interp(eff): - return np.interp(eff, tpr, fpr) - return interp - rocs.setdefault(clf_label, {})[energy] = mk_interp(tpr, fpr) - plt.plot(tpr[cond] * 100, 1 / fpr[cond], label=f"{clf_label}") - plt.yscale("log") + energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3")) + clf_label = "leading cluster" + def clf(events): + return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value + e_pred = clf(e_eval[energy]) + pi_pred = clf(pi_eval[energy]) + + for do_log, ax in [(False, axs[ix]), (True, axs_log[ix])]: + plt.sca(ax) + plt.hist(e_pred, weights=np.full_like(e_pred, 1.0 / ak.num(e_pred, axis=0)), bins=np.linspace(0., 1.01, 101), label=rf"$e^-$ {clf_label}") + plt.hist(pi_pred, weights=np.full_like(pi_pred, 1.0 / ak.num(pi_pred, axis=0)), bins=np.linspace(0., 1.01, 101), label=rf"$\pi^-$ {clf_label}", histtype="step") plt.title(f"{energy}") - plt.legend(loc="lower left") - plt.xlabel("Electron efficiency, %") - plt.ylabel("Pion rejection factor") + plt.legend() + plt.xlabel("Classifier output") + plt.ylabel("Event yield") + if do_log: + plt.yscale("log") + + plt.sca(axs_roc[ix]) + fpr, tpr, _ = roc_curve( + np.concatenate([np.ones_like(e_pred), np.zeros_like(pi_pred)]), + np.concatenate([e_pred, pi_pred]), + ) + cond = fpr != 0 # avoid infinite rejection (region of large uncertainty) + cond &= tpr != 1 # avoid linear interpolation (region of large uncertainty) + def mk_interp(tpr, fpr): + def interp(eff): + return np.interp(eff, tpr, fpr) + return interp + rocs.setdefault(clf_label, {})[energy] = mk_interp(tpr, fpr) + plt.plot(tpr[cond] * 100, 1 / fpr[cond], label=f"{clf_label}") + plt.yscale("log") + plt.title(f"{energy}") + plt.legend(loc="lower left") + plt.xlabel("Electron efficiency, %") + plt.ylabel("Pion rejection factor") fig.savefig(output_dir / f"pred.pdf", bbox_inches="tight") fig.savefig(output_dir / f"pred.png", bbox_inches="tight") From 19021f03e52a4aacc487a8c361e6455819d0e62f Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 7 Dec 2024 01:54:09 -0500 Subject: [PATCH 08/11] backwards_ecal: smooth fit plot --- benchmarks/backwards_ecal/backwards_ecal.org | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/benchmarks/backwards_ecal/backwards_ecal.org b/benchmarks/backwards_ecal/backwards_ecal.org index ec359464..1d70e2c9 100644 --- a/benchmarks/backwards_ecal/backwards_ecal.org +++ b/benchmarks/backwards_ecal/backwards_ecal.org @@ -202,12 +202,15 @@ for clf_label, sigma_rel_FWHM_cb in sigmas_rel_FWHM_cb.items(): ls="none", label=f"{clf_label}" ) + xmin = np.min(energy_values[cond]) + xmax = np.max(energy_values[cond]) + xs = np.arange(xmin, xmax, 0.1) plt.plot( - energy_values[cond], - f(energy_values[cond], *par), + xs, + f(xs, *par), ls="--", lw=0.5, - label=f"{clf_label}, ${np.ceil(stochastic * 10) / 10:.1f}\% / \sqrt{{E}} \oplus {np.ceil(constant * 10) / 10:.1f}\%$", + label=f"Functional fit: ${np.ceil(stochastic * 10) / 10:.1f}\% / \sqrt{{E}} \oplus {np.ceil(constant * 10) / 10:.1f}\%$", ) plt.fill_between( energy_values, From 06349afd80a2e91cac49a52df15898030cf60418 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 7 Dec 2024 12:38:58 -0500 Subject: [PATCH 09/11] backwards_ecal: put PLOT_TITLE to legend --- benchmarks/backwards_ecal/Snakefile | 5 +++++ benchmarks/backwards_ecal/backwards_ecal.org | 4 ++-- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/benchmarks/backwards_ecal/Snakefile b/benchmarks/backwards_ecal/Snakefile index 9bcbb9d6..b9b19787 100644 --- a/benchmarks/backwards_ecal/Snakefile +++ b/benchmarks/backwards_ecal/Snakefile @@ -154,6 +154,11 @@ rule backwards_ecal: directory("results/backwards_ecal/{CAMPAIGN}/") shell: """ +if [[ "{wildcards.CAMPAIGN}" == "local" ]]; then +export PLOT_TITLE="Benchmark simulation" +else +export PLOT_TITLE="\\textbf{{ePIC}} Simulation {wildcards.CAMPAIGN}" +fi env \ MATPLOTLIBRC={input.matplotlibrc} \ DETECTOR_CONFIG=""" + DETECTOR_CONFIG + """ \ diff --git a/benchmarks/backwards_ecal/backwards_ecal.org b/benchmarks/backwards_ecal/backwards_ecal.org index 1d70e2c9..5ef59c82 100644 --- a/benchmarks/backwards_ecal/backwards_ecal.org +++ b/benchmarks/backwards_ecal/backwards_ecal.org @@ -127,7 +127,7 @@ fractions_below = {} for ix, energy in enumerate(energies): energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3")) - clf_label = "leading cluster" + clf_label = PLOT_TITLE def clf(events): return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value e_pred = clf(e_eval[energy]) @@ -242,7 +242,7 @@ rocs = {} for ix, energy in enumerate(energies): energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3")) - clf_label = "leading cluster" + clf_label = PLOT_TITLE def clf(events): return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value e_pred = clf(e_eval[energy]) From cfb458ca290283b3341280df2c8197e35f96d1d6 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 7 Dec 2024 12:39:42 -0500 Subject: [PATCH 10/11] backwards_ecal: smooth fil_between, adjust xlim --- benchmarks/backwards_ecal/backwards_ecal.org | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/benchmarks/backwards_ecal/backwards_ecal.org b/benchmarks/backwards_ecal/backwards_ecal.org index 5ef59c82..6ca1e916 100644 --- a/benchmarks/backwards_ecal/backwards_ecal.org +++ b/benchmarks/backwards_ecal/backwards_ecal.org @@ -212,13 +212,17 @@ for clf_label, sigma_rel_FWHM_cb in sigmas_rel_FWHM_cb.items(): lw=0.5, label=f"Functional fit: ${np.ceil(stochastic * 10) / 10:.1f}\% / \sqrt{{E}} \oplus {np.ceil(constant * 10) / 10:.1f}\%$", ) +xmin = np.min(energy_values) +xmax = np.max(energy_values) * 1.05 +xs = np.arange(xmin, xmax + 0.1, 0.1) plt.fill_between( - energy_values, - np.sqrt((2 / np.sqrt(energy_values)) ** 2 + 1 ** 2), - np.sqrt((2 / np.sqrt(energy_values)) ** 2 + 3 ** 2), - alpha=0.2, color="black", label=r"YR requirement $2\% / \sqrt{E} \oplus (1-3)\%$", + xs, + np.sqrt((2 / np.sqrt(xs)) ** 2 + 1 ** 2), + np.sqrt((2 / np.sqrt(xs)) ** 2 + 3 ** 2), + lw=0., alpha=0.2, color="black", label=r"YR requirement $2\% / \sqrt{E} \oplus (1-3)\%$", ) -plt.ylim(top=10.) +plt.xlim(0., xmax) +plt.ylim(top=6.) plt.legend() plt.xlabel("Energy, GeV", loc="right") plt.ylabel(r"$\sigma_{E} / E$ derived from FWHM, %", loc="top") @@ -296,6 +300,8 @@ for clf_label, roc in rocs.items(): marker=".", label=f"{clf_label}", ) +xmax = np.max(energy_values) * 1.05 +plt.xlim(0., xmax) plt.yscale("log") plt.legend() plt.xlabel("Energy, GeV") From 5dfe66983267e7e20b960bf71ba8ad43a6f5ecf4 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Tue, 10 Dec 2024 00:49:41 -0500 Subject: [PATCH 11/11] add bench:backwards_ecal_campaigns job, stream from xrootd by defaultwq --- benchmarks/backwards_ecal/Snakefile | 2 +- benchmarks/backwards_ecal/config.yml | 10 ++++++++++ 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/benchmarks/backwards_ecal/Snakefile b/benchmarks/backwards_ecal/Snakefile index b9b19787..5efaf844 100644 --- a/benchmarks/backwards_ecal/Snakefile +++ b/benchmarks/backwards_ecal/Snakefile @@ -77,7 +77,7 @@ rule backwards_ecal_local_sim_list: fp.write("\n".join(input)) -if False: +if config.get("stream_from_xrootd", True) not in [False, "", "0", "false"]: rule backwards_ecal_campaign_sim_list: output: "listing/backwards_ecal/{CAMPAIGN}/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}.lst", diff --git a/benchmarks/backwards_ecal/config.yml b/benchmarks/backwards_ecal/config.yml index 75d8f464..68a183e5 100644 --- a/benchmarks/backwards_ecal/config.yml +++ b/benchmarks/backwards_ecal/config.yml @@ -28,6 +28,16 @@ bench:backwards_ecal: - pip install -r benchmarks/backwards_ecal/requirements.txt - snakemake $SNAKEMAKE_FLAGS --cores 1 results/backwards_ecal/local +bench:backwards_ecal_campaigns: + extends: .det_benchmark + stage: benchmarks + when: manual + timeout: 4 hours + script: + - export PYTHONUSERBASE=$LOCAL_DATA_PATH/deps + - pip install -r benchmarks/backwards_ecal/requirements.txt + - snakemake $SNAKEMAKE_FLAGS --cores 1 results/backwards_ecal/24.10.1 + collect_results:backwards_ecal: extends: .det_benchmark stage: collect