Fix feeddown folding

Author: qgp

machine_learning_hep/analysis/analyzer_jets.py

@@ -226,7 +226,7 @@ def calculate_efficiencies(self):
                 # gen-level efficiency for feeddown estimation
                 h_eff_gen = h_genmatch[cat].Clone()
                 h_eff_gen.Divide(h_gen[cat])
-                self._save_hist(h_eff_gen, f'eff/h_effgen_{cat}.png')
+                self._save_hist(h_eff_gen, f'eff/h_effgen_{cat}.png', 'text')
                 self.hcandeff_gen[cat] = h_eff_gen
 
                 # matching loss
@@ -947,18 +947,18 @@ def estimate_feeddown(self):
                                          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()
             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
             for ipt in range(get_nbins(h3_fd_gen, 1)):
-                eff_np = self.hcandeff_gen['np'].GetBinContent(ipt+1)
-                for iptjet, ishape in itertools.product(
-                        range(get_nbins(h3_fd_gen, 0)), range(get_nbins(h3_fd_gen, 2))):
-                    scale_bin(h3_fd_gen, eff_np, iptjet+1, ipt+1, ishape+1)
+                for iptjet in range(get_nbins(h3_fd_gen, 0)):
+                    eff_np = self.hcandeff_gen['np'].GetBinContent(iptjet + 1, ipt + 1)
+                    print(f'scaling with {eff_np=} for {iptjet=}, {ipt=}', flush=True)
+                    for ishape in range(get_nbins(h3_fd_gen, 2)):
+                        scale_bin(h3_fd_gen, eff_np, iptjet+1, ipt+1, ishape+1)
             self._save_hist(project_hist(h3_fd_gen, [0, 2], {}), f'fd/h_ptjet-{var}_fdnew_gen_geneff.png')
 
             # 3d folding incl. kinematic efficiencies
@@ -970,6 +970,8 @@ def estimate_feeddown(self):
                     rfile.Get(f'h_effkine_fd_det_nocuts_{var}'),
                     rfile.Get(f'h_effkine_fd_det_cut_{var}'))
                 h_response = rfile.Get(f'h_response_fd_{var}')
+                h_response.Print('v')
+                print(f'fd folding for {var=}')
                 h_response_norm = norm_response(h_response, 3)
                 h3_fd_gen.Multiply(h_effkine_gen)
                 self._save_hist(project_hist(h3_fd_gen, [0, 2], {}), f'fd/h_ptjet-{var}_fdnew_gen_genkine.png')
@@ -987,14 +989,15 @@ def estimate_feeddown(self):
                 for iptjet, ishape in itertools.product(
                         range(get_nbins(h3_fd_det, 0)), range(get_nbins(h3_fd_det, 2))):
                     scale_bin(h3_fd_det, 1./eff_pr, iptjet+1, ipt+1, ishape+1)
-            self._save_hist(project_hist(h3_fd_det, [0, 2], {}), f'fd/h_ptjet-{var}_fdnew_det_deteff.png')
 
             # project to 2d (ptjet-shape)
             h_fd_det = project_hist(h3_fd_det, [0, 2], {})
+            self._save_hist(h_fd_det, f'fd/h_ptjet-{var}_fdnew_det_deteff.png')
 
             # old method
             h3_fd_gen = h3_fd_gen_orig.Clone()
             ensure_sumw2(h3_fd_gen)
+            self._save_hist(project_hist(h3_fd_gen, [0, 2], {}), f'fd/h_ptjet-{var}_feeddown_gen_noeffscaling.png')
             for ipt in range(get_nbins(h3_fd_gen, 1)):
                 eff_pr = self.hcandeff['pr'].GetBinContent(ipt+1)
                 eff_np = self.hcandeff['np'].GetBinContent(ipt+1)

machine_learning_hep/processer_jet.py

@@ -429,7 +429,7 @@ def _prepare_response_fd(self, dfi, h_effkine, h_response, var):
                     (df[f'{var}_gen'] >= axis_var_gen.GetXmin()) & (df[f'{var}_gen'] < axis_var_gen.GetXmax())]
         fill_hist(h_effkine[('det', 'cut', var)], df[['fJetPt', 'fPt', var]])
 
-        fill_hist(h_response[var], df[['fJetPt', 'fPt', f'{var}', 'fJetPt_gen', 'fPt_gen', f'{var}_gen']])
+        fill_hist(h_response[var], df[['fJetPt_gen', 'fPt_gen', f'{var}_gen', 'fJetPt_gen', 'fPt_gen', f'{var}_gen']])
 
         df = dfi
         df = df.loc[(df.fJetPt_gen >= axis_ptjet_gen.GetXmin()) & (df.fJetPt_gen < axis_ptjet_gen.GetXmax()) &

machine_learning_hep/utils/hist.py

@@ -293,13 +293,19 @@ def norm_response(response, dim_out):
     for bin_in in itertools.product(*(range(1, get_nbins(response_norm, iaxis) + 1)
                                       for iaxis in range(dim_out, get_dim(response_norm)))):
         for iaxis, val in enumerate(bin_in, dim_out):
-            get_axis(response_norm, iaxis).SetRange(val, val)
-        norm = response_norm.Projection(0).Integral()
+            get_axis(response, iaxis).SetRange(val, val)
+        # norm = response.Projection(0).Integral()
+        norm = 0.
+        for bin_out in itertools.product(*(range(1, get_nbins(response, i) + 1) for i in range(dim_out))):
+            norm += get_bin_val(response, bin_out + bin_in)
         if np.isclose(norm, 0.):
             continue
-        for bin_out in itertools.product(*(range(1, get_nbins(response_norm, i)+1) for i in range(dim_out))):
-            set_bin_val(response_norm, bin_out + bin_in, get_bin_val(response_norm, bin_out + bin_in) / norm)
-            set_bin_err(response_norm, bin_out + bin_in, get_bin_err(response_norm, bin_out + bin_in) / norm)
+        total = 0.
+        for bin_out in itertools.product(*(range(1, get_nbins(response_norm, i) + 1) for i in range(dim_out))):
+            set_bin_val(response_norm, bin_out + bin_in, get_bin_val(response, bin_out + bin_in) / norm)
+            set_bin_err(response_norm, bin_out + bin_in, get_bin_err(response, bin_out + bin_in) / norm)
+            total += get_bin_val(response_norm, bin_out + bin_in)
+        print(f'distributed {bin_in=} to {total=} counts')
     return response_norm
 
 
@@ -309,10 +315,18 @@ def fold_hist(hist, response):
     dim_out = get_dim(response) - get_dim(hist)
     axes_spec = list(np.array(get_axis(response, i).GetXbins(), 'd') for i in range(dim_out))
     hfold = create_hist('test', 'test', *axes_spec)
-    for bin_out in itertools.product(*(range(1, get_nbins(hfold, i)+1) for i in range(get_dim(hfold)))):
+    # TODO: setup axes
+    for bin_in in itertools.product(*(range(1, get_nbins(hist, i) + 1) for i in range(get_dim(hist)))):
+        total = 0.
+        for bin_out in itertools.product(*(range(1, get_nbins(hfold, i) + 1) for i in range(get_dim(hfold)))):
+            total += get_bin_val(response, bin_out + bin_in)
+        print(f'redistributed {bin_in=} to {total=} counts')
+
+    for bin_out in itertools.product(*(range(1, get_nbins(hfold, i) + 1) for i in range(get_dim(hfold)))):
         val = 0.
         err = 0.
-        for bin_in in itertools.product(*(range(1, get_nbins(hist, i)+1) for i in range(get_dim(hist)))):
+        for bin_in in itertools.product(*(range(1, get_nbins(hist, i) + 1) for i in range(get_dim(hist)))):
+            print(f'{bin_out=} collecting {bin_in=} with weight {get_bin_val(response, bin_out + bin_in)}')
             val += get_bin_val(hist, bin_in) * get_bin_val(response, bin_out + bin_in)
             err += get_bin_err(hist, bin_in)**2 * get_bin_val(response, bin_out + bin_in)**2
         set_bin_val(hfold, bin_out, val)