Propagate POWHEG weights to luminosity scaling

Author: qgp

machine_learning_hep/analysis/analyzer_jets.py

@@ -1126,7 +1126,10 @@ def estimate_feeddown(self):
                                 project_hist(h3_fd_gen_orig[var], [0, 2], {}),
                                 f"fd/h_ptjet-{var}_feeddown_genonly_noeffscaling.png",
                             )
-                powheg_xsection_scale_factor = 0.  # FIXME: retrieve cross section
+                    h_norm = rfile.Get("histonorm")
+                    powheg_xsection_scale_factor = h_norm.GetBinContent(5)
+                    self.logger.info("powheg_xsection_scale_factor = %f", powheg_xsection_scale_factor)
+                    self.logger.info("POWHEG luminosity (mb^{-1}): %g", 1.0 / powheg_xsection_scale_factor)
 
             case fd_input:
                 self.logger.critical("Invalid feeddown input %s", fd_input)

machine_learning_hep/data/data_run3/database_ml_parameters_D0Jet_pp.yml

@@ -70,8 +70,9 @@ D0Jet_pp:
           O2hfd0mccollbase: [fPosZ]
         filter: "abs(fPosZ) < 10."
 
-      # powhegcoll:
-          # O2d0cmcpjetmcco: ...
+      coll_powheg_fd:
+        trees:
+          O2d0cmcpjetmcco: [fEventWeight]
       collcnt_mc:
         level: mc
         trees:
@@ -234,6 +235,9 @@ D0Jet_pp:
         source: colldata
         file: AnalysisResultsEvt.parquet
 
+      coll_powheg_fd:
+        level: fd
+        file: AnalysisResultsPowheg.parquet
       evtorig_fd:
         level: fd
         source: collfd
@@ -311,6 +315,7 @@ D0Jet_pp:
     namefile_reco_applieddata: AnalysisResultsRecoAppliedData.parquet
     namefile_reco_appliedmc: AnalysisResultsRecoAppliedMC.parquet
     namefile_mcweights: mcweights.root
+    namefile_wgt: AnalysisResultsPowheg.parquet
     treeoutput: "D0tree"
     histofilename: "masshisto.root"
     efffilename: "effhisto.root"

machine_learning_hep/processer.py

@@ -143,6 +143,7 @@ def __init__(
         self.n_evt = datap["files_names"]["namefile_evt"]
         self.n_collcnt = datap["files_names"].get("namefile_collcnt")
         self.n_bccnt = datap["files_names"].get("namefile_bccnt")
+        self.n_wgt = datap["files_names"].get("namefile_wgt", None)
         self.n_evtorig = datap["files_names"].get("namefile_evtorig")
         self.n_evt_count_ml = datap["files_names"].get("namefile_evt_count", "evtcount.yaml")
         self.n_gen = datap["files_names"]["namefile_gen"]
@@ -195,6 +196,7 @@ def __init__(
         self.l_evtorig = createlist(self.d_pkl, self.l_path, self.n_evtorig)
         self.l_collcnt = createlist(self.d_pkl, self.l_path, self.n_collcnt) if self.datatype != "fd" else None
         self.l_bccnt = createlist(self.d_pkl, self.l_path, self.n_bccnt) if self.datatype != "fd" else None
+        self.l_wgt = createlist(self.d_pkl, self.l_path, self.n_wgt) if self.datatype == "fd" and self.n_wgt else None
         self.l_histomass = createlist(self.d_results, self.l_path, self.n_filemass)
         self.l_histoeff = createlist(self.d_results, self.l_path, self.n_fileeff)
         # self.l_historesp = createlist(self.d_results, self.l_path, self.n_fileresp)
@@ -384,6 +386,8 @@ def dfread(rdir, trees, cols, idx_name=None):
                     trees = [trees]
                     cols = [cols]
                 # if all(type(var) is str for var in vars): vars = [vars]
+                if not all((name in rdir for name in trees)):
+                    self.logger.critical("Missing trees: %s", trees)
                 df = None
                 for tree, col in zip([rdir[name] for name in trees], cols):
                     try:
@@ -536,10 +540,13 @@ def dfuse(df_spec):
             for df_name, df_spec in self.df_write.items():
                 if dfuse(df_spec):
                     src = df_spec.get("source", df_name)
-                    dfo = dfquery(dfs[src], df_spec.get("filter", None))
-                    path = os.path.join(self.d_pkl, self.l_path[file_index], df_spec["file"])
-                    self.logger.info("writing %s to %s with info %s", df_name, path, dfo.info())
-                    write_df(dfo, path)
+                    if src in dfs:
+                        dfo = dfquery(dfs[src], df_spec.get("filter", None))
+                        path = os.path.join(self.d_pkl, self.l_path[file_index], df_spec["file"])
+                        self.logger.info("writing %s to %s with info %s", df_name, path, dfo.info())
+                        write_df(dfo, path)
+                    else:
+                        self.logger.error("could not write tree, missing source %s", src)
 
     def skim(self, file_index):
         dfreco = read_df(self.l_reco[file_index]) if self.datatype != "fd" else None

machine_learning_hep/processer_jet.py

@@ -241,7 +241,7 @@ def process_histomass_single(self, index):
         print(f"Opening file {self.l_histomass[index]}", flush=True)
         with TFile.Open(self.l_histomass[index], "recreate") as _:
             dfevtorig = read_df(self.l_evtorig[index])
-            histonorm = TH1F("histonorm", "histonorm", 4, 0, 4)
+            histonorm = TH1F("histonorm", "histonorm", 5, 0, 5)
             histonorm.SetBinContent(1, len(dfquery(dfevtorig, self.s_evtsel)))
             if self.l_collcnt:
                 dfcollcnt = read_df(self.l_collcnt[index])
@@ -257,10 +257,18 @@ def process_histomass_single(self, index):
                 ser_bccnt = dfbccnt[self.cfg("counter_tvx")]
                 bccnt_tvx = functools.reduce(lambda x, y: float(x) + float(y), (ar[0] for ar in ser_bccnt))
                 histonorm.SetBinContent(4, bccnt_tvx)
+            if self.l_wgt:
+                self.logger.info("Filling event weights")
+                dfwgt = read_df(self.l_wgt[index])
+                print(dfwgt.info())
+                histonorm.SetBinContent(5, dfwgt["fEventWeight"].sum())
+            else:
+                self.logger.warning("No event weights found, empty list: %s", self.l_wgt)
             get_axis(histonorm, 0).SetBinLabel(1, "N_{evt}")
             get_axis(histonorm, 0).SetBinLabel(2, "N_{coll}")
             get_axis(histonorm, 0).SetBinLabel(3, "N_{coll}^{TVX}")
             get_axis(histonorm, 0).SetBinLabel(4, "N_{BC}^{TVX}")
+            get_axis(histonorm, 0).SetBinLabel(5, "w_{POWHEG}")
             histonorm.Write()
 
             if self.datatype != "fd":