Extend jet analyzer to read feeddown info from processor

Author: qgp

machine_learning_hep/analysis/analyzer_jets.py

@@ -66,10 +66,12 @@ def __init__(self, datap, case, typean, period):
         self.d_resultsallpdata = (
             self.cfg(f"data.results.{period}") if period is not None else self.cfg("data.resultsallp")
         )
+        self.d_resultsallpfd = self.cfg(f"fd.results.{period}") if period is not None else self.cfg("fd.resultsallp")
 
         # input directories (processor output)
         self.d_resultsallpmc_proc = self.d_resultsallpmc
         self.d_resultsallpdata_proc = self.d_resultsallpdata
+        self.d_resultsallpfd_proc = self.d_resultsallpfd
         # use a different processor output
         if "data_proc" in datap["analysis"][typean]:
             self.d_resultsallpdata_proc = (
@@ -79,11 +81,16 @@ def __init__(self, datap, case, typean, period):
             self.d_resultsallpmc_proc = (
                 self.cfg(f"mc_proc.results.{period}") if period is not None else self.cfg("mc_proc.resultsallp")
             )
+        if "fd_proc" in datap["analysis"][typean]:
+            self.d_resultsallpfd_proc = (
+                self.cfg(f"fd_proc.results.{period}") if period is not None else self.cfg("fd_proc.resultsallp")
+            )
 
         # input files
         n_filemass_name = datap["files_names"]["histofilename"]
         self.n_filemass = os.path.join(self.d_resultsallpdata_proc, n_filemass_name)
         self.n_filemass_mc = os.path.join(self.d_resultsallpmc_proc, n_filemass_name)
+        self.n_filemass_fd = os.path.join(self.d_resultsallpfd_proc, n_filemass_name)
         self.n_fileeff = datap["files_names"]["efffilename"]
         self.n_fileeff = os.path.join(self.d_resultsallpmc_proc, self.n_fileeff)
         self.n_fileresp = datap["files_names"]["respfilename"]
@@ -1050,57 +1057,73 @@ def _extract_signal(self, hist, var, mcordata, ipt):
         # hres.Sumw2() # TODO: check if we should do this here
         return hres
 
-    # region feeddown
-    # pylint: disable=too-many-statements
-    def estimate_feeddown(self):
-        self.logger.info("Estimating feeddown")
-
-        with TFile(self.cfg("fd_root")) as rfile:
-            powheg_xsection = rfile.Get("fHistXsection")
-            powheg_xsection_scale_factor = powheg_xsection.GetBinContent(1) / powheg_xsection.GetEntries()
-        self.logger.info("POWHEG luminosity (mb^{-1}): %g", 1.0 / powheg_xsection_scale_factor)
 
-        df = pd.read_parquet(self.cfg("fd_parquet"))
-        col_mapping = {"dr": "delta_r_jet", "zpar": "z"}  # TODO: check mapping
+    def estimate_feeddown(self):
+        """Estimate feeddown from legacy Run 2 trees or gen-only simulation"""
+        match self.cfg("fd_input", "tree"):
+            case "tree":
+                with TFile(self.cfg("fd_root")) as rfile:
+                    powheg_xsection = rfile.Get("fHistXsection")
+                    powheg_xsection_scale_factor = powheg_xsection.GetBinContent(1) / powheg_xsection.GetEntries()
+                self.logger.info("POWHEG luminosity (mb^{-1}): %g", 1.0 / powheg_xsection_scale_factor)
+
+                df = pd.read_parquet(self.cfg("fd_parquet"))
+                col_mapping = {"dr": "delta_r_jet", "zpar": "z"}  # TODO: check mapping
+
+                # TODO: generalize to higher dimensions
+                h3_fd_gen_orig = {}
+                for var in self.observables["all"]:
+                    bins_ptjet = np.asarray(self.cfg("bins_ptjet"), "d")
+                    # TODO: generalize or derive from histogram?
+                    bins_obs = {}
+                    if binning := self.cfg(f"observables.{var}.bins_gen_var"):
+                        bins_tmp = np.asarray(binning, "d")
+                    elif binning := self.cfg(f"observables.{var}.bins_gen_fix"):
+                        bins_tmp = bin_array(*binning)
+                    elif binning := self.cfg(f"observables.{var}.bins_var"):
+                        bins_tmp = np.asarray(binning, "d")
+                    elif binning := self.cfg(f"observables.{var}.bins_fix"):
+                        bins_tmp = bin_array(*binning)
+                    else:
+                        self.logger.error("no binning specified for %s, using defaults", var)
+                        bins_tmp = bin_array(10, 0.0, 1.0)
+                    bins_obs[var] = bins_tmp
+
+                    colname = col_mapping.get(var, f"{var}_jet")
+                    if f"{colname}" not in df:
+                        if var is not None:
+                            self.logger.error("No feeddown information for %s (%s), cannot estimate feeddown", var, colname)
+                            # print(df.info(), flush=True)
+                        continue
+
+                    # TODO: derive histogram
+                    h3_fd_gen_orig[var] = create_hist(
+                        "h3_feeddown_gen",
+                        f";p_{{T}}^{{jet}} (GeV/#it{{c}});p_{{T}}^{{HF}} (GeV/#it{{c}});{var}",
+                        bins_ptjet,
+                        self.bins_candpt,
+                        bins_obs[var],
+                    )
+                    fill_hist_fast(h3_fd_gen_orig, df[["pt_jet", "pt_cand", f"{colname}"]])
+                    self._save_hist(project_hist(h3_fd_gen_orig, [0, 2], {}), f"fd/h_ptjet-{var}_feeddown_gen_noeffscaling.png")
+
+            case "sim":
+                # TODO: recover cross section
+                h3_fd_gen_orig = {}
+                with TFile(self.n_filemass_fd) as rfile:
+                    for var in self.observables["all"]:
+                        self.logger.info("Running feeddown analysis for obs. %s", var)
+                        label = f"-{var}" if var else ""
+                        if fh := rfile.Get(f"h_mass-ptjet-pthf{label}"):
+                            h3_fd_gen_orig[var] = project_hist(fh, list(range(1, get_dim(fh))), {})
+                            ensure_sumw2(h3_fd_gen_orig[var])
+
+            case fd_input:
+                self.logger.critical("Invalid feeddown input %s", fd_input)
 
-        # TODO: generalize to higher dimensions
         for var in self.observables["all"]:
-            bins_ptjet = np.asarray(self.cfg("bins_ptjet"), "d")
-            # TODO: generalize or derive from histogram?
-            bins_obs = {}
-            if binning := self.cfg(f"observables.{var}.bins_gen_var"):
-                bins_tmp = np.asarray(binning, "d")
-            elif binning := self.cfg(f"observables.{var}.bins_gen_fix"):
-                bins_tmp = bin_array(*binning)
-            elif binning := self.cfg(f"observables.{var}.bins_var"):
-                bins_tmp = np.asarray(binning, "d")
-            elif binning := self.cfg(f"observables.{var}.bins_fix"):
-                bins_tmp = bin_array(*binning)
-            else:
-                self.logger.error("no binning specified for %s, using defaults", var)
-                bins_tmp = bin_array(10, 0.0, 1.0)
-            bins_obs[var] = bins_tmp
-
-            colname = col_mapping.get(var, f"{var}_jet")
-            if f"{colname}" not in df:
-                if var is not None:
-                    self.logger.error("No feeddown information for %s (%s), cannot estimate feeddown", var, colname)
-                    # print(df.info(), flush=True)
-                continue
-
-            # TODO: derive histogram
-            h3_fd_gen_orig = create_hist(
-                "h3_feeddown_gen",
-                f";p_{{T}}^{{jet}} (GeV/#it{{c}});p_{{T}}^{{HF}} (GeV/#it{{c}});{var}",
-                bins_ptjet,
-                self.bins_candpt,
-                bins_obs[var],
-            )
-            fill_hist_fast(h3_fd_gen_orig, df[["pt_jet", "pt_cand", f"{colname}"]])
-            self._save_hist(project_hist(h3_fd_gen_orig, [0, 2], {}), f"fd/h_ptjet-{var}_feeddown_gen_noeffscaling.png")
-
             # new method
-            h3_fd_gen = h3_fd_gen_orig.Clone()
+            h3_fd_gen = h3_fd_gen_orig[var].Clone()
             ensure_sumw2(h3_fd_gen)
             self._save_hist(project_hist(h3_fd_gen, [0, 2], {}), f"fd/h_ptjet-{var}_fdnew_gen.png")
             # apply np efficiency

machine_learning_hep/data/data_run3/database_ml_parameters_D0Jet_pp.yml

@@ -308,7 +308,7 @@ D0Jet_pp:
   multi:
     data:
       nprocessesparallel: 80
-      maxfiles: [1] #list of periods
+      maxfiles: [-1] #list of periods
       chunksizeunp: [100] #list of periods
       chunksizeskim: [100] #list of periods
       fracmerge: [.1] #list of periods
@@ -325,7 +325,7 @@ D0Jet_pp:
       mcreweights: [../Analyses] #list of periods
     mc:
       nprocessesparallel: 80
-      maxfiles: [1] #list of periods
+      maxfiles: [-1] #list of periods
       chunksizeunp: [100] #list of periods
       chunksizeskim: [1000] #list of periods
       fracmerge: [1.] #list of periods
@@ -342,7 +342,7 @@ D0Jet_pp:
       mcreweights: [../Analyses] #list of periods
     fd:
       nprocessesparallel: 80
-      maxfiles: [1] #list of periods
+      maxfiles: [-1] #list of periods
       chunksizeunp: [100] #list of periods
       chunksizeskim: [1000] #list of periods
       fracmerge: [1.] #list of periods

machine_learning_hep/submission/processor.yml

@@ -63,7 +63,7 @@ analysis:
     histomass: true # processer: process_histomass
   mc:
     histomass: true # processer: process_histomass
-    efficiency: false # processer: process_efficiency
+    efficiency: true # processer: process_efficiency
   fd:
     histomass: true
   steps: