Merge pull request #67 from nbiederbeck/fixes

Stack Counts, use Underscores

scripts/plot_counts_after_cuts.py

@@ -1,75 +1,39 @@
 from argparse import ArgumentParser
 
-import numpy as np
-from astropy import units as u
 from astropy.table import Table
 from matplotlib import pyplot as plt
 
 parser = ArgumentParser()
-parser.add_argument("-i", "--input-paths", required=True, nargs="+")
+parser.add_argument("-i", "--input-path", required=True)
 parser.add_argument("-o", "--output-path", required=True)
-parser.add_argument("--norm", default="none")
 args = parser.parse_args()
 
 
-def main(input_paths, output_path, norm):
- cuts_after_trigger = []
- cuts_after_gh = []
- cuts_after_gh_theta = []
- t_effective = []
- t_elapsed = []
-
- for path in input_paths:
- cuts = Table.read(path, path="cuts")
- cuts_after_trigger.append(cuts["after_trigger"])
- cuts_after_gh.append(cuts["after_gh"])
- cuts_after_gh_theta.append(cuts["after_gh_theta"])
- t_effective.append(cuts.meta["t_effective"])
- t_elapsed.append(cuts.meta["t_elapsed"])
-
- fig, ax = plt.subplots()
+def main(input_path, output_path):
+ cuts = Table.read(input_path, path="cuts")
 
  x = cuts["center"]
  xerr = (cuts["high"] - x, x - cuts["low"])
 
- t_eff = u.Quantity(t_effective).reshape(-1, 1)
- if norm == "none":
- norm = u.Quantity(1)
- elif norm == "eff":
- norm = np.sum(t_eff)
- else:
- raise NotImplementedError(f"Unsupported norm {norm}")
-
- cuts_after_trigger = np.sum(
- np.array(cuts_after_trigger) / norm,
- axis=0,
- )
- cuts_after_gh = np.sum(
- np.array(cuts_after_gh) / norm,
- axis=0,
- )
- cuts_after_gh_theta = np.sum(
- np.array(cuts_after_gh_theta) / norm,
- axis=0,
- )
+ fig, ax = plt.subplots()
 
  ax.errorbar(
  x,
- cuts_after_trigger,
+ cuts["after_trigger"],
  xerr=xerr,
  ls="",
  label="Trigger",
  )
  ax.errorbar(
  x,
- cuts_after_gh,
+ cuts["after_gh"],
  xerr=xerr,
  ls="",
  label="GH Cut",
  )
  ax.errorbar(
  x,
- cuts_after_gh_theta,
+ cuts["after_gh_theta"],
  xerr=xerr,
  ls="",
  label="Theta Cut",
@@ -78,7 +42,7 @@ def main(input_paths, output_path, norm):
  ax.set_yscale("log")
 
  ax.set_xlabel(rf"E_{{\text{{reco}}}} \:/\: {x.unit}")
- ax.set_ylabel(cuts_after_trigger.unit)
+ ax.set_ylabel(cuts["after_trigger"].unit)
 
  ax.legend(title="Counts after")
 

scripts/stack_counts_after_cuts.py

@@ -0,0 +1,64 @@
+from argparse import ArgumentParser
+
+import numpy as np
+from astropy import units as u
+from astropy.table import Table
+
+parser = ArgumentParser()
+parser.add_argument("-i", "--input-paths", required=True, nargs="+")
+parser.add_argument("-o", "--output-path", required=True)
+parser.add_argument("--norm", default="none")
+args = parser.parse_args()
+
+
+def main(input_paths, output_path, norm):
+ cuts_after_trigger = []
+ cuts_after_gh = []
+ cuts_after_gh_theta = []
+ t_effective = []
+ t_elapsed = []
+
+ for path in input_paths:
+ cuts = Table.read(path, path="cuts")
+ cuts_after_trigger.append(cuts["after_trigger"])
+ cuts_after_gh.append(cuts["after_gh"])
+ cuts_after_gh_theta.append(cuts["after_gh_theta"])
+ t_effective.append(cuts.meta["t_effective"])
+ t_elapsed.append(cuts.meta["t_elapsed"])
+
+ t_eff = u.Quantity(t_effective).reshape(-1, 1)
+ if norm == "none":
+ norm = u.Quantity(1)
+ elif norm == "eff":
+ norm = np.sum(t_eff)
+ else:
+ raise NotImplementedError(f"Unsupported norm {norm}")
+
+ cuts_after_trigger = np.sum(
+ np.array(cuts_after_trigger) / norm,
+ axis=0,
+ )
+ cuts_after_gh = np.sum(
+ np.array(cuts_after_gh) / norm,
+ axis=0,
+ )
+ cuts_after_gh_theta = np.sum(
+ np.array(cuts_after_gh_theta) / norm,
+ axis=0,
+ )
+
+ table = Table(
+ {
+ "after_trigger": cuts_after_trigger,
+ "after_gh": cuts_after_gh,
+ "after_gh_theta": cuts_after_gh_theta,
+ "center": cuts["center"],
+ "high": cuts["high"],
+ "low": cuts["low"],
+ },
+ )
+ table.write(output_path, path="cuts", overwrite=True, serialize_meta=True)
+
+
+if __name__ == "__main__":
+ main(**vars(args))

workflow/Snakefile

@@ -33,7 +33,7 @@ analyses = [
 dl3_plots = [
  build_dir / f"plots/{plot}/{run_id}.pdf"
  for run_id in RUN_IDS + ["stacked"]
- for plot in ["theta2", "skymap", "skymap-dl3"]
+ for plot in ["theta2", "skymap", "skymap_dl3"]
 ]
 
 

workflow/rules/create_dl3.smk

@@ -111,7 +111,7 @@ rule cuts_dl2_dl3:
  conda:
  lstchain_env
  output:
- build_dir / "dl3/counts_after_gh_theta_cut_Run{run_id}.h5",
+ build_dir / "dl3/counts/after_gh_theta_cut_{run_id}.h5",
  input:
  dl2=build_dir / "dl2/dl2_LST-1.Run{run_id}.h5",
  irf=build_dir / "irf/calculated/irf_Run{run_id}.fits.gz",
@@ -121,22 +121,35 @@ rule cuts_dl2_dl3:
  "python {input.script} --input-dl2 {input.dl2} --input-irf {input.irf} -c {input.config} -o {output}"
 
 
-rule plot_cuts_dl2_dl3:
+rule stack_cuts_dl2_dl3:
  conda:
  lstchain_env
  output:
- build_dir / "plots/counts_after_gh_theta_cut_{norm}.pdf",
+ build_dir / "dl3/counts/after_gh_theta_cut_{norm}_stacked.h5",
  input:
  data=expand(
- build_dir / "dl3/counts_after_gh_theta_cut_Run{run_id}.h5",
+ build_dir / "dl3/counts/after_gh_theta_cut_{run_id}.h5",
  run_id=RUN_IDS,
  ),
- script="scripts/plot_counts_after_cuts.py",
+ script="scripts/stack_counts_after_cuts.py",
  rc=os.environ.get("MATPLOTLIBRC", config_dir / "matplotlibrc"),
  shell:
  "MATPLOTLIBRC={input.rc} python {input.script} -i {input.data} -o {output} --norm {wildcards.norm}"
 
 
+rule plot_cuts_dl2_dl3:
+ conda:
+ lstchain_env
+ output:
+ build_dir / "plots/counts_after_gh_theta_cut_{norm}.pdf",
+ input:
+ data=build_dir / "dl3/counts/after_gh_theta_cut_{norm}.h5",
+ script="scripts/plot_counts_after_cuts.py",
+ rc=os.environ.get("MATPLOTLIBRC", config_dir / "matplotlibrc"),
+ shell:
+ "MATPLOTLIBRC={input.rc} python {input.script} -i {input.data} -o {output}"
+
+
 rule calc_skymap:
  resources:
  mem_mb="64G",
@@ -185,10 +198,10 @@ rule stack_skymaps_dl3:
  conda:
  lstchain_env
  output:
- build_dir / "dl3/skymap-dl3/stacked.fits",
+ build_dir / "dl3/skymap_dl3/stacked.fits",
  input:
  data=expand(
- build_dir / "dl3/skymap-dl3/{run_id}.fits",
+ build_dir / "dl3/skymap_dl3/{run_id}.fits",
  run_id=RUN_IDS,
  ),
  script="scripts/stack_skymap.py",

workflow/rules/post_dl3.smk

@@ -219,7 +219,7 @@ rule model_best_fit:
 
 rule calc_skymap_per_obs:
  output:
- build_dir / "dl3/skymap-dl3/{run_id}.fits",
+ build_dir / "dl3/skymap_dl3/{run_id}.fits",
  input:
  data=build_dir / "dl3/dl3_LST-1.Run{run_id}.fits.gz",
  script="scripts/calc_skymap_gammas.py",
@@ -238,9 +238,9 @@ rule calc_skymap_per_obs:
 
 rule plot_skymap_dl3:
  output:
- build_dir / "plots/skymap-dl3/{runid}.pdf",
+ build_dir / "plots/skymap_dl3/{runid}.pdf",
  input:
- data=build_dir / "dl3/skymap-dl3/{runid}.fits",
+ data=build_dir / "dl3/skymap_dl3/{runid}.fits",
  script="scripts/plot_skymap_dl3.py",
  rc=os.environ.get("MATPLOTLIBRC", config_dir / "matplotlibrc"),
  conda: