From 89bb7b52160e9e7228819b7dc53b360d784881e9 Mon Sep 17 00:00:00 2001
From: Felix Metzner <felixmetzner@outlook.com>
Date: Mon, 15 Apr 2024 18:12:53 +0200
Subject: [PATCH] Adding normalization option to shape sys eval pull plots and
 making shape sys eval plots more costumizable.

---
 .../dedicated_fit_routine.py                  | 40 ++++++++++++++-----
 .../dedicated_fit_approach/plotting_tools.py  | 33 ++++++++++++---
 2 files changed, 57 insertions(+), 16 deletions(-)

diff --git a/rdstar/offline_analysis/fitting/dedicated_fit_approach/dedicated_fit_routine.py b/rdstar/offline_analysis/fitting/dedicated_fit_approach/dedicated_fit_routine.py
index c4ba31311..3dcbdd066 100644
--- a/rdstar/offline_analysis/fitting/dedicated_fit_approach/dedicated_fit_routine.py
+++ b/rdstar/offline_analysis/fitting/dedicated_fit_approach/dedicated_fit_routine.py
@@ -1121,6 +1121,9 @@ class RDStarFitEvaluator:
         self,
         fit_status_container: RDStarFitStatusContainer,
         add_statistical_uncertainties: bool,
+        normalize_pulls: bool,
+        plot_subsets: bool,
+        scale_factor: Optional[float] = None,
     ) -> Generator[SpecificShapePlotInfoContainer, None, None]:
         sys_shape_plotter: SystematicsShapePlotter = SystematicsShapePlotter(
             base_output_dir_path=self.base_output_dir_path,
@@ -1136,19 +1139,23 @@ class RDStarFitEvaluator:
             relative_shape_error=fit_instance.ff_sig_rel_shape_sys_uncert_matrix,
             pure_bin_counts=fit_instance.squared_template_stat_error,
             components_to_consider=("BpDztau", "BzDmtau", "BpDzStau", "BzDmStau"),
-            scale_factor=10,
-            plot_subsets=False,
+            scale_factor=scale_factor,
+            plot_subsets=plot_subsets,
         )
 
         yield from sys_shape_plotter.create_systematics_shape_plots(
             sys_shape_info=example_sys_shape_info,
             add_statistical_uncertainties=add_statistical_uncertainties,
+            normalize_pulls=normalize_pulls,
         )
 
     def plot_norm_ff_sys_shape_effects(
         self,
         fit_status_container: RDStarFitStatusContainer,
         add_statistical_uncertainties: bool,
+        normalize_pulls: bool,
+        plot_subsets: bool,
+        scale_factor: Optional[float] = None,
     ) -> Generator[SpecificShapePlotInfoContainer, None, None]:
         sys_shape_plotter: SystematicsShapePlotter = SystematicsShapePlotter(
             base_output_dir_path=self.base_output_dir_path,
@@ -1164,19 +1171,23 @@ class RDStarFitEvaluator:
             relative_shape_error=fit_instance.ff_norm_rel_shape_sys_uncert_matrix,
             pure_bin_counts=fit_instance.squared_template_stat_error,
             components_to_consider=("BpDzl", "BzDml", "BpDzSl", "BzDmSl"),
-            scale_factor=10,
-            plot_subsets=False,
+            scale_factor=scale_factor,
+            plot_subsets=plot_subsets,
         )
 
         yield from sys_shape_plotter.create_systematics_shape_plots(
             sys_shape_info=example_sys_shape_info,
             add_statistical_uncertainties=add_statistical_uncertainties,
+            normalize_pulls=normalize_pulls,
         )
 
     def plot_l_id_sys_shape_effects(
         self,
         fit_status_container: RDStarFitStatusContainer,
         add_statistical_uncertainties: bool,
+        normalize_pulls: bool,
+        plot_subsets: bool,
+        scale_factor: Optional[float] = None,
     ) -> Generator[SpecificShapePlotInfoContainer, None, None]:
         sys_shape_plotter: SystematicsShapePlotter = SystematicsShapePlotter(
             base_output_dir_path=self.base_output_dir_path,
@@ -1192,19 +1203,23 @@ class RDStarFitEvaluator:
             relative_shape_error=fit_instance.lepton_id_rel_shape_sys_uncert_matrix,
             pure_bin_counts=fit_instance.squared_template_stat_error,
             components_to_consider=("BpDzl", "BzDml", "BpDzSl", "BzDmSl", "BBbarBKG_in_cB", "BBbarBKG_in_nB"),
-            scale_factor=10,
-            plot_subsets=False,
+            scale_factor=scale_factor,
+            plot_subsets=plot_subsets,
         )
 
         yield from sys_shape_plotter.create_systematics_shape_plots(
             sys_shape_info=example_sys_shape_info,
             add_statistical_uncertainties=add_statistical_uncertainties,
+            normalize_pulls=normalize_pulls,
         )
 
     def plot_tracking_sys_shape_effects(
         self,
         fit_status_container: RDStarFitStatusContainer,
         add_statistical_uncertainties: bool,
+        normalize_pulls: bool,
+        plot_subsets: bool,
+        scale_factor: Optional[float] = None,
     ) -> Generator[SpecificShapePlotInfoContainer, None, None]:
         sys_shape_plotter: SystematicsShapePlotter = SystematicsShapePlotter(
             base_output_dir_path=self.base_output_dir_path,
@@ -1220,19 +1235,23 @@ class RDStarFitEvaluator:
             relative_shape_error=fit_instance.tracking_rel_shape_sys_effect,
             pure_bin_counts=fit_instance.squared_template_stat_error,
             components_to_consider=None,
-            scale_factor=10,
-            plot_subsets=False,
+            scale_factor=scale_factor,
+            plot_subsets=plot_subsets,
         )
 
         yield from sys_shape_plotter.create_systematics_shape_plots(
             sys_shape_info=example_sys_shape_info,
             add_statistical_uncertainties=add_statistical_uncertainties,
+            normalize_pulls=normalize_pulls,
         )
 
     def plot_k_short_sys_shape_effects(
         self,
         fit_status_container: RDStarFitStatusContainer,
         add_statistical_uncertainties: bool,
+        normalize_pulls: bool,
+        plot_subsets: bool,
+        scale_factor: Optional[float] = None,
     ) -> Generator[SpecificShapePlotInfoContainer, None, None]:
         sys_shape_plotter: SystematicsShapePlotter = SystematicsShapePlotter(
             base_output_dir_path=self.base_output_dir_path,
@@ -1248,13 +1267,14 @@ class RDStarFitEvaluator:
             relative_shape_error=fit_instance.k_short_rel_shape_sys_uncert_matrix,
             pure_bin_counts=fit_instance.squared_template_stat_error,
             components_to_consider=None,
-            scale_factor=10,
-            plot_subsets=False,
+            scale_factor=scale_factor,
+            plot_subsets=plot_subsets,
         )
 
         yield from sys_shape_plotter.create_systematics_shape_plots(
             sys_shape_info=example_sys_shape_info,
             add_statistical_uncertainties=add_statistical_uncertainties,
+            normalize_pulls=normalize_pulls,
         )
 
     def dump_data_for_external_fitter(
diff --git a/rdstar/offline_analysis/fitting/dedicated_fit_approach/plotting_tools.py b/rdstar/offline_analysis/fitting/dedicated_fit_approach/plotting_tools.py
index b2e3feab5..1ec0538e4 100644
--- a/rdstar/offline_analysis/fitting/dedicated_fit_approach/plotting_tools.py
+++ b/rdstar/offline_analysis/fitting/dedicated_fit_approach/plotting_tools.py
@@ -215,8 +215,9 @@ class SpecificShapePlotInfoContainer:
     observable: FitObservableInfo
     reco_ch_info: RecoChannelInfo
     component_info: ComponentInfo
+    add_statistical_uncertainty: bool
+    normalize_pulls: bool
     scale_factor: Optional[float] = None
-    add_statistical_uncertainty: bool = True
 
     def __post_init__(self) -> None:
         assert len(self.normed_base_shape.shape) == 1, (len(self.normed_base_shape.shape), self.normed_base_shape.shape)
@@ -276,6 +277,7 @@ class SystematicsShapePlotter:
         reco_ch_info: RecoChannelInfo,
         component_info: ComponentInfo,
         add_statistical_uncertainty: bool,
+        normalize_pulls: bool,
     ) -> Generator[SpecificShapePlotInfoContainer, None, None]:
 
         assert len(selected_normed_base_shape.shape) == 1, selected_normed_base_shape.shape
@@ -360,8 +362,9 @@ class SystematicsShapePlotter:
                         observable=obs_info,
                         reco_ch_info=reco_ch_info,
                         component_info=component_info,
-                        scale_factor=sys_shape_info.scale_factor,
                         add_statistical_uncertainty=add_statistical_uncertainty,
+                        normalize_pulls=normalize_pulls,
+                        scale_factor=sys_shape_info.scale_factor,
                     )
 
             yield SpecificShapePlotInfoContainer(
@@ -374,14 +377,16 @@ class SystematicsShapePlotter:
                 observable=obs_info,
                 reco_ch_info=reco_ch_info,
                 component_info=component_info,
-                scale_factor=sys_shape_info.scale_factor,
                 add_statistical_uncertainty=add_statistical_uncertainty,
+                normalize_pulls=normalize_pulls,
+                scale_factor=sys_shape_info.scale_factor,
             )
 
     def create_systematics_shape_plots(
         self,
         sys_shape_info: ShapePlotInfoContainer,
         add_statistical_uncertainties: bool = True,
+        normalize_pulls: bool = True,
     ) -> Generator[SpecificShapePlotInfoContainer, None, None]:
         target_dir_path: PathType = os.path.join(self.output_dir_path, sys_shape_info.name)
         os.makedirs(target_dir_path, exist_ok=True)
@@ -439,6 +444,7 @@ class SystematicsShapePlotter:
                     reco_ch_info=reco_ch_info,
                     component_info=comp_info,
                     add_statistical_uncertainty=add_statistical_uncertainties,
+                    normalize_pulls=normalize_pulls,
                 ):
                     self.plot_systematics_shape_effect_for(
                         spec_sys_shape_info=_this_shape_info,
@@ -497,6 +503,8 @@ class SystematicsShapePlotter:
 
         bin_stat_uncert: np.ndarray = spec_sys_shape_info.stat_error
 
+        _bin_norm_values: np.ndarray = np.where(norm_bin_count > 0.0, norm_bin_count, 1.0)
+
         assert np.all(u_bin_shape_uncert * d_bin_shape_uncert <= 0.0), (
             np.sum(u_bin_shape_uncert * d_bin_shape_uncert > 0.0),
             u_bin_shape_uncert * d_bin_shape_uncert,
@@ -560,21 +568,34 @@ class SystematicsShapePlotter:
 
         # Pull Plot
         if spec_sys_shape_info.add_statistical_uncertainty:
+            if spec_sys_shape_info.normalize_pulls:
+                _bin_stat_uncert_values: np.ndarray = bin_stat_uncert / _bin_norm_values
+            else:
+                _bin_stat_uncert_values = bin_stat_uncert
+
             ax2.bar(
                 x=bin_mids,
-                height=2 * bin_stat_uncert,
+                height=2 * _bin_stat_uncert_values,
                 width=bin_widths,
-                bottom=-1.0 * bin_stat_uncert,
+                bottom=-1.0 * _bin_stat_uncert_values,
                 color=KITColors.grey,
                 alpha=0.4,
                 fill=True,
                 lw=0,
                 label="stat.",
             )
+
+        if spec_sys_shape_info.normalize_pulls:
+            _u_bin_shape_uncert_values: np.ndarray = u_bin_shape_uncert / _bin_norm_values
+            _d_bin_shape_uncert_values: np.ndarray = d_bin_shape_uncert / _bin_norm_values
+        else:
+            _u_bin_shape_uncert_values = u_bin_shape_uncert
+            _d_bin_shape_uncert_values = d_bin_shape_uncert
+
         ax2.hist(
             x=[bin_mids, bin_mids],
             bins=bin_edges,
-            weights=[u_bin_shape_uncert, d_bin_shape_uncert],
+            weights=[_u_bin_shape_uncert_values, _d_bin_shape_uncert_values],
             stacked=False,
             color=[KITColors.kit_blue, KITColors.kit_orange],
             lw=1.3,
-- 
GitLab