reweighting_function.FluxReweighter Class Reference

Collaboration diagram for reweighting_function.FluxReweighter:

Public Member Functions
	__init__ (self, config_path="config.yaml")
	read_files (self, path)
	base_cut (self, df)
	compute_ratio (self)
	fit (self, degree)
	run (self)

Static Public Member Functions
	fit_poly_weighted (x, y, sigma, degree=1)

Public Attributes
	cfg
	bins_E
	E
	H_DPM
	s_DPM
	H_LC
	s_LC
	w
	sig

Static Public Attributes
dict	PID = {"nue": 12, "numu": 14, "nutau": 16}
list	COLS

Protected Member Functions
	_load_config (self, yaml_path)
	_report_selection (self, name, df, pid)
	_hist (self, df, mask, use_weight=True)
	_fit_mask (self, verbose=False)
	_make_plot (self, deg, poly, chi2red, fit_mask)

Static Protected Member Functions
	_require (d, key, section)
	_adaptive_edges (fine, dsw, dsw2, lsw, lsw2, min_neff)
	_poly_equation_text (poly, deg)

Detailed Description

Derive the reweighting function w(E) = fast-sim / DPMJET from flux files.

The class histograms the four input samples in neutrino energy, forms the
ratio bin by bin with propagated errors, and fits a polynomial to it. All
settings come from a YAML config passed to the constructor.

Typical use
-----------
    rw = FluxReweighter("Config_epos.yaml")
    rw.run()                       # diagnostic plot per degree
    poly, chi2red = rw.fit(2)      # or just the deg-2 weight function

Attributes populated by :meth:`compute_ratio`
---------------------------------------------
E      : bin centres [GeV]
w      : the ratio fast-sim / DPMJET per bin (this is w(E) sampled)
sig    : uncertainty on the ratio per bin
H_DPM  : summed DPMJET histogram (light + charm), s_DPM its error
H_LC   : summed fast-sim histogram (light + charm), s_LC its error

Definition at line 104 of file reweighting_function.py.

Constructor & Destructor Documentation

__init__()

reweighting_function.FluxReweighter.__init__	(		self,
			config_path = "config.yaml"
	)

Load the YAML config and prepare the default (uniform) bin edges.

Adaptive edges, if requested, are computed later in
:meth:`compute_ratio` because they depend on the selected events.

Definition at line 141 of file reweighting_function.py.

    def __init__(self, config_path="config.yaml"):
        """Load the YAML config and prepare the default (uniform) bin edges.
 
        Adaptive edges, if requested, are computed later in
        :meth:`compute_ratio` because they depend on the selected events.
        """
        self.cfg = self._load_config(config_path)
        self.bins_E = np.linspace(
            self.cfg["energy_min"], self.cfg["energy_max"], num=self.cfg["bins"]
        )
        # filled in by compute_ratio()
        self.E = self.w = self.sig = None
        self.H_DPM = self.s_DPM = self.H_LC = self.s_LC = None
 

Member Function Documentation

_adaptive_edges()

reweighting_function.FluxReweighter._adaptive_edges (   fine,   dsw,   dsw2,   lsw,   lsw2,   min_neff  )

staticprotected

Build variable-width bin edges with a minimum of effective entries.

Starting from a fine uniform grid, fine bins are accumulated left to
right and a bin is closed as soon as *both* the DPMJET and the fast-sim
side reach

    N_eff = (sum w)^2 / sum w^2  >=  min_neff

N_eff is the number of unweighted events that would give the same
statistical power, so this guarantees no bin is dominated by a handful
of large-weight events -- the usual failure mode at high energy, where
a single charm event with a huge weight otherwise makes the ratio jump
and its error explode.

Parameters
----------
fine : ndarray
    Edges of the fine starting grid.
dsw, dsw2 : ndarray
    Per-fine-bin sum of weights and sum of weights squared, DPMJET side.
lsw, lsw2 : ndarray
    The same for the fast-sim side.
min_neff : float
    Required effective entries per output bin.

Returns
-------
ndarray
    The merged bin edges (always closed at ``fine[-1]``).

Definition at line 163 of file reweighting_function.py.

    def _adaptive_edges(fine, dsw, dsw2, lsw, lsw2, min_neff):
        """Build variable-width bin edges with a minimum of effective entries.
 
        Starting from a fine uniform grid, fine bins are accumulated left to
        right and a bin is closed as soon as *both* the DPMJET and the fast-sim
        side reach
 
            N_eff = (sum w)^2 / sum w^2  >=  min_neff
 
        N_eff is the number of unweighted events that would give the same
        statistical power, so this guarantees no bin is dominated by a handful
        of large-weight events -- the usual failure mode at high energy, where
        a single charm event with a huge weight otherwise makes the ratio jump
        and its error explode.
 
        Parameters
        ----------
        fine : ndarray
            Edges of the fine starting grid.
        dsw, dsw2 : ndarray
            Per-fine-bin sum of weights and sum of weights squared, DPMJET side.
        lsw, lsw2 : ndarray
            The same for the fast-sim side.
        min_neff : float
            Required effective entries per output bin.
 
        Returns
        -------
        ndarray
            The merged bin edges (always closed at ``fine[-1]``).
        """
        edges = [fine[0]]
        adsw = adsw2 = alsw = alsw2 = 0.0
        for i in range(len(fine) - 1):
            adsw += dsw[i]; adsw2 += dsw2[i]
            alsw += lsw[i]; alsw2 += lsw2[i]
            nd = (adsw * adsw / adsw2) if adsw2 > 0 else 0.0
            nl = (alsw * alsw / alsw2) if alsw2 > 0 else 0.0
            if min(nd, nl) >= min_neff:
                edges.append(fine[i + 1]); adsw = adsw2 = alsw = alsw2 = 0.0
        if edges[-1] != fine[-1]:
            edges.append(fine[-1])
        return np.array(edges)
 

_fit_mask()

reweighting_function.FluxReweighter._fit_mask (   self,   verbose = False  )

protected

Boolean mask selecting which ratio bins enter the polynomial fit.

A bin is kept when it lies inside ``fit_range``, has a non-zero ratio
and error, passes the relative-error cap, and falls inside the optional
ratio window. With ``verbose`` the number of bins removed by the ratio
window is printed.

Definition at line 454 of file reweighting_function.py.

    def _fit_mask(self, verbose=False):
        """Boolean mask selecting which ratio bins enter the polynomial fit.
 
        A bin is kept when it lies inside ``fit_range``, has a non-zero ratio
        and error, passes the relative-error cap, and falls inside the optional
        ratio window. With ``verbose`` the number of bins removed by the ratio
        window is printed.
        """
        c = self.cfg
        rel_err = np.zeros_like(self.w)
        nz = (self.w > 0) & (self.sig > 0)
        rel_err[nz] = self.sig[nz] / self.w[nz]
 
        base = (
            (self.E >= c["fit_min"]) & (self.E <= c["fit_max"]) &
            (self.w > 0) & (self.sig > 0) &
            (rel_err <= c["rel_err_max"])
        )
 
        rmin = -np.inf if c["ratio_min"] is None else c["ratio_min"]
        rmax =  np.inf if c["ratio_max"] is None else c["ratio_max"]
        within = (self.w >= rmin) & (self.w <= rmax)
 
        if verbose and (c["ratio_min"] is not None or c["ratio_max"] is not None):
            dropped = int((base & ~within).sum())
            print(f"  ratio window [{rmin:g}, {rmax:g}]: dropped {dropped} of "
                  f"{int(base.sum())} otherwise-eligible bin(s)")
        return base & within
 

_hist()

reweighting_function.FluxReweighter._hist (   self,   df,   mask,   use_weight = True  )

protected

Weighted energy histogram of the selected rows, with its error.

Each event contributes its own ``weight in pb`` (or 1.0 when
``use_weight`` is false). The error is the standard weighted-histogram
one, sqrt(sum w^2) per bin.

Returns ``(counts, sigma)``.

Definition at line 347 of file reweighting_function.py.

    def _hist(self, df, mask, use_weight=True):
        """Weighted energy histogram of the selected rows, with its error.
 
        Each event contributes its own ``weight in pb`` (or 1.0 when
        ``use_weight`` is false). The error is the standard weighted-histogram
        one, sqrt(sum w^2) per bin.
 
        Returns ``(counts, sigma)``.
        """
        x = df.loc[mask, "E in GeV"]
        w = df.loc[mask, "weight in pb"] if use_weight else pd.Series(1.0, index=x.index)
        counts, _ = np.histogram(x, bins=self.bins_E, weights=w)
        var, _    = np.histogram(x, bins=self.bins_E, weights=w ** 2)
        return counts.astype(float), np.sqrt(var.astype(float))
 

_load_config()

reweighting_function.FluxReweighter._load_config (   self,   yaml_path  )

protected

Parse the reweighter YAML into a flat dict of settings.

Only the four input paths are mandatory; every other key has a default
so a minimal config stays valid. See Config_epos.yaml for the meaning
of each field.

Definition at line 207 of file reweighting_function.py.

    def _load_config(self, yaml_path):
        """Parse the reweighter YAML into a flat dict of settings.
 
        Only the four input paths are mandatory; every other key has a default
        so a minimal config stays valid. See Config_epos.yaml for the meaning
        of each field.
        """
        with open(yaml_path) as f:
            raw = yaml.safe_load(f) or {}
 
        inp = raw.get("inputs", {})
        eb  = raw.get("energy_binning", {})
        fr  = raw.get("fit_range", {})
        cut = raw.get("cuts", {})
        out = raw.get("output", {})
        wt  = raw.get("weights", {})
 
        return {
            # --- input flux files (mandatory) ---
            "light_dpm":     self._require(inp, "light_dpm", "inputs"),
            "charm_dpm":     self._require(inp, "charm_dpm", "inputs"),
            "light_fastsim": self._require(inp, "light_fastsim", "inputs"),
            "charm_fastsim": self._require(inp, "charm_fastsim", "inputs"),
 
            # --- neutrino flavour entering the ratio ---
            "flavour": raw.get("flavour", "mu"),
 
            # --- histogram binning in energy ---
            "energy_min": float(eb.get("min", 0.0)),
            "energy_max": float(eb.get("max", 4000.0)),
            "bins":       int(eb.get("bins", 51)),        # uniform mode: n edges
            "bin_mode":   eb.get("mode", "uniform"),      # "uniform" | "adaptive"
            "fine_bins":  int(eb.get("fine_bins", 400)),  # adaptive: start grid
            "min_neff":   float(eb.get("min_neff", 10.0)),# adaptive: entries/bin
 
            # --- energy window the polynomial is fitted over ---
            "fit_min": float(fr.get("min", 0.0)),
            "fit_max": float(fr.get("max", 4000.0)),
 
            # --- one diagnostic plot per degree ---
            "degrees": [int(d) for d in raw.get("polynomial_degrees", [1, 2, 3])],
 
            # --- quality gates on which bins enter the fit ---
            # Relative-error cap: drops bins whose sigma/ratio is too large.
            # This is the unbiased gate -- prefer it over the ratio window.
            "rel_err_max": float(raw.get("rel_err_max", 3.0)),
            # Optional window on the ratio VALUE. Keep it symmetric about the
            # data: a one-sided cap removes up-fluctuations only and biases the
            # fit low. None disables that side.
            "ratio_min": (None if raw.get("ratio_min") is None
                          else float(raw.get("ratio_min"))),
            "ratio_max": (None if raw.get("ratio_max") is None
                          else float(raw.get("ratio_max"))),
 
            # --- Selection cut ---
            "cut_energy_min": float(cut.get("energy_min", 10.0)),
            "cut_x_col":      cut.get("x_column", "x in m"),
            "cut_y_col":      cut.get("y_column", "y in m"),
            "cut_x_min":      float(cut.get("x_min", -46.0)),
            "cut_x_max":      float(cut.get("x_max",  -7.0)),
            "cut_y_min":      float(cut.get("y_min",  15.0)),
            "cut_y_max":      float(cut.get("y_max",  54.0)),
            # x/y are multiplied by this before the bounds are applied, so cm
            # bounds can be used with metre columns (unit_scale: 100).
            "cut_scale":      float(cut.get("unit_scale", 1.0)),
            # Set false when the inputs already passed the cut-based analysis.
            "cut_enabled":    bool(cut.get("enabled", True)),
 
            # --- whether 'weight in pb' is applied when histogramming ---
            # Both default to true. Disabling one side only puts numerator and
            # denominator on different scales and breaks the ratio's meaning.
            "weight_dpm":     bool(wt.get("dpm", True)),
            "weight_fastsim": bool(wt.get("fastsim", True)),
 
            # --- output ---
            "outpath": out.get("path", "./output"),
            "label":   out.get("label", "EPOS+POWHEG"),
        }
 

_make_plot()

reweighting_function.FluxReweighter._make_plot (   self,   deg,   poly,   chi2red,   fit_mask  )

protected

Write the diagnostic ratio plot for one polynomial degree.

Shows every ratio bin faintly, the subset entering the fit with error
bars, the fitted curve inside ``fit_range`` and its extrapolation
outside, plus a statistics box. Comparing the two point sets shows at a
glance what the quality gates removed.

Returns the path of the saved PDF.

Definition at line 532 of file reweighting_function.py.

    def _make_plot(self, deg, poly, chi2red, fit_mask):
        """Write the diagnostic ratio plot for one polynomial degree.
 
        Shows every ratio bin faintly, the subset entering the fit with error
        bars, the fitted curve inside ``fit_range`` and its extrapolation
        outside, plus a statistics box. Comparing the two point sets shows at a
        glance what the quality gates removed.
 
        Returns the path of the saved PDF.
        """
        c = self.cfg
        E, w, sig = self.E, self.w, self.sig
        E_fit, w_fit, sig_fit = E[fit_mask], w[fit_mask], sig[fit_mask]
 
        E_plot_fit = np.linspace(c["fit_min"], c["fit_max"], 300)
        E_plot_all = np.linspace(E.min(), E.max(), 400)
 
        fig, ax = plt.subplots(figsize=(9, 5.8))
        ax.errorbar(E, w, yerr=sig, fmt="o", ms=4, alpha=0.25, label="all data")
        ax.errorbar(E_fit, w_fit, yerr=sig_fit, fmt="o", ms=5, capsize=3, label="used in fit")
        ax.plot(E_plot_fit, poly(E_plot_fit), "-", lw=2.5, label=f"poly deg={deg} (fit range)")
        ax.plot(E_plot_all, poly(E_plot_all), "--", lw=1.8, alpha=0.6, label="poly extrapolation")
 
        ax.axvline(c["fit_min"], color="grey", ls=":", lw=1)
        ax.axvline(c["fit_max"], color="grey", ls=":", lw=1)
        ax.axhline(1.0, color="grey", ls="--", lw=1)   # ratio = 1: models agree
 
        ax.set_xlabel("Energy [GeV]")
        ax.set_ylabel("ratio")
        ax.set_title(c["label"])
        ax.grid(True, which="both", alpha=0.3)
        ax.legend(fontsize="small", loc="lower left")
 
        stats_text = "\n".join([
            rf"deg = {deg}",
            rf"$\chi^2$/ndf = {chi2red:.2f}",
            rf"$N_\mathrm{{points}}$ = {len(E_fit)}",
            rf"$E_\mathrm{{fit}}$: {c['fit_min']:.0f}-{c['fit_max']:.0f} GeV",
            self._poly_equation_text(poly, deg),
        ])
        ax.text(0.98, 0.95, stats_text, transform=ax.transAxes,
                ha="right", va="top", fontsize=10,
                bbox=dict(boxstyle="round", facecolor="white", alpha=0.85))
 
        fig.tight_layout()
        outfile = os.path.join(c["outpath"],
                               f"reweight_poly_deg{deg}_{c['flavour']}.pdf")
        fig.savefig(outfile)
        plt.close(fig)
        return outfile
 

_poly_equation_text()

reweighting_function.FluxReweighter._poly_equation_text (   poly,   deg  )

staticprotected

Format w(E) as a LaTeX string for the plot's statistics box.

Definition at line 520 of file reweighting_function.py.

    def _poly_equation_text(poly, deg):
        """Format w(E) as a LaTeX string for the plot's statistics box."""
        terms = []
        for p, coef in zip(range(deg, -1, -1), poly.coefficients):
            if p == 0:
                terms.append(rf"({coef:.3e})")
            elif p == 1:
                terms.append(rf"({coef:.3e})\,E")
            else:
                terms.append(rf"({coef:.3e})\,E^{p}")
        return r"$w(E) = " + " + ".join(terms) + r"$"
 

_report_selection()

reweighting_function.FluxReweighter._report_selection (   self,   name,   df,   pid  )

protected

Print how many rows survive each selection stage for one file.

Printed before the ratio is built. If a stage removes everything the
actual x/y data ranges are shown next to the cut bounds, which makes a
unit mismatch (cut in cm, data in metres) immediately obvious.

Returns the number of surviving rows.

Definition at line 316 of file reweighting_function.py.

    def _report_selection(self, name, df, pid):
        """Print how many rows survive each selection stage for one file.
 
        Printed before the ratio is built. If a stage removes everything the
        actual x/y data ranges are shown next to the cut bounds, which makes a
        unit mismatch (cut in cm, data in metres) immediately obvious.
 
        Returns the number of surviving rows.
        """
        c = self.cfg
        fl = (df["Neutrino"].abs() == pid)
        if not c["cut_enabled"]:
            print(f"  [{name:13s}] N={len(df):>7} | flavour={int(fl.sum()):>7} "
                  f"| cuts disabled")
            return int(fl.sum())
        e  = fl & (df["E in GeV"] > c["cut_energy_min"])
        xs = df[c["cut_x_col"]] * c["cut_scale"]
        ys = df[c["cut_y_col"]] * c["cut_scale"]
        x  = e & xs.between(c["cut_x_min"], c["cut_x_max"])
        y  = x & ys.between(c["cut_y_min"], c["cut_y_max"])
        n = int(y.sum())
        print(f"  [{name:13s}] N={len(df):>7} | flavour={int(fl.sum()):>7} "
              f"| +E>{c['cut_energy_min']:g}={int(e.sum()):>7} "
              f"| +x={int(x.sum()):>7} | +y={n:>7}")
        if n == 0:
            print(f"      0 survive -> {c['cut_x_col']}*{c['cut_scale']:g} range "
                  f"[{xs.min():.3g}, {xs.max():.3g}] vs cut [{c['cut_x_min']}, {c['cut_x_max']}]; "
                  f"{c['cut_y_col']}*{c['cut_scale']:g} range "
                  f"[{ys.min():.3g}, {ys.max():.3g}] vs cut [{c['cut_y_min']}, {c['cut_y_max']}]")
        return n
 

_require()

reweighting_function.FluxReweighter._require (   d,   key,   section  )

staticprotected

Return a mandatory config value, or raise naming the missing field.

Definition at line 156 of file reweighting_function.py.

    def _require(d, key, section):
        """Return a mandatory config value, or raise naming the missing field."""
        if d is None or key not in d or d[key] is None:
            raise ValueError(f"Missing required config field: {section}.{key}")
        return d[key]
 

base_cut()

reweighting_function.FluxReweighter.base_cut	(		self,
			df
	)

Boolean mask for the fiducial (detector acceptance) selection.

Applies a minimum energy and an x/y box, with the x/y columns scaled by
``cuts.unit_scale`` first. Returns all-true when ``cuts.enabled`` is
false, i.e. when the inputs already went through the cut-based analysis.

Definition at line 298 of file reweighting_function.py.

    def base_cut(self, df):
        """Boolean mask for the fiducial (detector acceptance) selection.
 
        Applies a minimum energy and an x/y box, with the x/y columns scaled by
        ``cuts.unit_scale`` first. Returns all-true when ``cuts.enabled`` is
        false, i.e. when the inputs already went through the cut-based analysis.
        """
        c = self.cfg
        if not c["cut_enabled"]:
            return pd.Series(True, index=df.index)
        xs = df[c["cut_x_col"]] * c["cut_scale"]
        ys = df[c["cut_y_col"]] * c["cut_scale"]
        return (
            (df["E in GeV"] > c["cut_energy_min"]) &
            xs.between(c["cut_x_min"], c["cut_x_max"]) &
            ys.between(c["cut_y_min"], c["cut_y_max"])
        )
 

compute_ratio()

reweighting_function.FluxReweighter.compute_ratio

(

self

)

Build the fast-sim / DPMJET ratio versus energy.

Steps: read the four files, print the selection report, optionally
derive adaptive bin edges, histogram each sample, sum light + charm on
each side, and divide.

Light and charm are histogrammed separately and the histograms added,
which is equivalent to weighting every event individually and filling
one histogram -- each event always keeps its own ``weight in pb``. The
errors are combined in quadrature, sqrt(s_light^2 + s_charm^2).

Populates ``E``, ``w``, ``sig``, ``H_DPM``, ``s_DPM``, ``H_LC``,
``s_LC`` and returns ``self`` so calls can be chained.

Definition at line 365 of file reweighting_function.py.

    def compute_ratio(self):
        """Build the fast-sim / DPMJET ratio versus energy.
 
        Steps: read the four files, print the selection report, optionally
        derive adaptive bin edges, histogram each sample, sum light + charm on
        each side, and divide.
 
        Light and charm are histogrammed separately and the histograms added,
        which is equivalent to weighting every event individually and filling
        one histogram -- each event always keeps its own ``weight in pb``. The
        errors are combined in quadrature, sqrt(s_light^2 + s_charm^2).
 
        Populates ``E``, ``w``, ``sig``, ``H_DPM``, ``s_DPM``, ``H_LC``,
        ``s_LC`` and returns ``self`` so calls can be chained.
        """
        pid = self.PID[self.cfg["flavour"]]
 
        df_ld = self.read_files(self.cfg["light_dpm"])
        df_cd = self.read_files(self.cfg["charm_dpm"])
        df_lf = self.read_files(self.cfg["light_fastsim"])
        df_cf = self.read_files(self.cfg["charm_fastsim"])
 
        # Selection breakdown per file; also reveals unit mismatches.
        print(f"Selection report ({self.cfg['label']}, flavour={self.cfg['flavour']}):")
        n_ld = self._report_selection("light_dpm", df_ld, pid)
        n_cd = self._report_selection("charm_dpm", df_cd, pid)
        n_lf = self._report_selection("light_fastsim", df_lf, pid)
        n_cf = self._report_selection("charm_fastsim", df_cf, pid)
        if (n_ld + n_cd) == 0 or (n_lf + n_cf) == 0:
            raise ValueError(
                "Selection leaves no events on one side of the ratio "
                f"(DPMJET={n_ld + n_cd}, fast-sim={n_lf + n_cf}). See the report "
                "above: the most common cause is the x/y fiducial cut being in "
                "different units than the data (cut in cm, 'x in m' column in "
                "metres). Fix via cuts.unit_scale (e.g. 100 for m->cm) or by "
                "putting the bounds in the data's units.")
 
        # Optional adaptive binning. Uses the same flavour + fiducial selection
        # as the histograms below, and overwrites self.bins_E before filling.
        if self.cfg["bin_mode"] == "adaptive":
            fine = np.linspace(self.cfg["energy_min"], self.cfg["energy_max"],
                               self.cfg["fine_bins"] + 1)
 
            def _fine(df, use_w):
                """Sum of w and w^2 per fine bin for one selected sample."""
                m = (df["Neutrino"].abs() == pid) & self.base_cut(df)
                E = df.loc[m, "E in GeV"].to_numpy()
                w = (df.loc[m, "weight in pb"].to_numpy()
                     if use_w else np.ones(int(m.sum())))
                sw, _  = np.histogram(E, bins=fine, weights=w)
                sw2, _ = np.histogram(E, bins=fine, weights=w ** 2)
                return sw.astype(float), sw2.astype(float)
 
            wd, wf = self.cfg["weight_dpm"], self.cfg["weight_fastsim"]
            aLd = _fine(df_ld, wd); aCd = _fine(df_cd, wd)
            aLf = _fine(df_lf, wf); aCf = _fine(df_cf, wf)
            self.bins_E = self._adaptive_edges(
                fine, aLd[0] + aCd[0], aLd[1] + aCd[1],
                aLf[0] + aCf[0], aLf[1] + aCf[1], self.cfg["min_neff"])
            print(f"adaptive binning -> {len(self.bins_E) - 1} bins "
                  f"(min_neff={self.cfg['min_neff']:g})")
 
        # Denominator: DPMJET total = light + charm.
        H_L_dpm, s_L_dpm = self._hist(df_ld, (df_ld["Neutrino"].abs() == pid) & self.base_cut(df_ld), use_weight=self.cfg["weight_dpm"])
        H_C_dpm, s_C_dpm = self._hist(df_cd, (df_cd["Neutrino"].abs() == pid) & self.base_cut(df_cd), use_weight=self.cfg["weight_dpm"])
        self.H_DPM = H_L_dpm + H_C_dpm
        self.s_DPM = np.sqrt(s_L_dpm ** 2 + s_C_dpm ** 2)
 
        # Numerator: fast-sim total = light (EPOS) + charm (POWHEG).
        H_L_fs, s_L_fs = self._hist(df_lf, (df_lf["Neutrino"].abs() == pid) & self.base_cut(df_lf), use_weight=self.cfg["weight_fastsim"])
        H_C_fs, s_C_fs = self._hist(df_cf, (df_cf["Neutrino"].abs() == pid) & self.base_cut(df_cf), use_weight=self.cfg["weight_fastsim"])
        self.H_LC = H_L_fs + H_C_fs
        self.s_LC = np.sqrt(s_L_fs ** 2 + s_C_fs ** 2)
 
        # Ratio with relative errors added in quadrature. Bins empty on either
        # side are left at zero and are excluded from the fit by _fit_mask().
        centers = 0.5 * (self.bins_E[:-1] + self.bins_E[1:])
        R  = np.zeros_like(self.H_DPM)
        sR = np.zeros_like(self.H_DPM)
 
        m = (self.H_DPM > 0) & (self.H_LC > 0)
        R[m] = self.H_LC[m] / self.H_DPM[m]
        rel2_num = (self.s_LC[m]  / self.H_LC[m])  ** 2
        rel2_den = (self.s_DPM[m] / self.H_DPM[m]) ** 2
        sR[m] = R[m] * np.sqrt(rel2_num + rel2_den)
 
        self.E, self.w, self.sig = centers, R, sR
        return self
 

fit()

reweighting_function.FluxReweighter.fit	(		self,
			degree
	)

Return the reweighting polynomial w(E) for one degree.

Computes the ratio on first use. Evaluate the result per event as
``poly(E)``; it is only meaningful inside ``fit_range``.

Returns ``(poly, chi2red)``. Raises ``ValueError`` if the selection
leaves fewer points than the degree requires.

Definition at line 495 of file reweighting_function.py.

    def fit(self, degree):
        """Return the reweighting polynomial w(E) for one degree.
 
        Computes the ratio on first use. Evaluate the result per event as
        ``poly(E)``; it is only meaningful inside ``fit_range``.
 
        Returns ``(poly, chi2red)``. Raises ``ValueError`` if the selection
        leaves fewer points than the degree requires.
        """
        if self.E is None:
            self.compute_ratio()
        mask = self._fit_mask(verbose=True)
        n_fit = int(mask.sum())
        if n_fit < degree + 1:
            raise ValueError(
                f"Not enough fit points ({n_fit}) for degree {degree}")
        poly, _chi2, chi2red = self.fit_poly_weighted(
            self.E[mask], self.w[mask], self.sig[mask], degree=degree
        )
        return poly, chi2red
 

fit_poly_weighted()

reweighting_function.FluxReweighter.fit_poly_weighted (   x,   y,   sigma,   degree = 1  )

static

Least-squares polynomial fit weighted by 1/sigma.

Returns ``(poly, chi2, chi2/ndf)`` where ``poly`` is a ``np.poly1d``
with coefficients ordered highest power first.

Definition at line 484 of file reweighting_function.py.

    def fit_poly_weighted(x, y, sigma, degree=1):
        """Least-squares polynomial fit weighted by 1/sigma.
 
        Returns ``(poly, chi2, chi2/ndf)`` where ``poly`` is a ``np.poly1d``
        with coefficients ordered highest power first.
        """
        poly  = np.poly1d(np.polyfit(x, y, deg=degree, w=1.0 / sigma))
        chi2  = np.sum(((y - poly(x)) / sigma) ** 2)
        ndof  = max(len(x) - (degree + 1), 1)
        return poly, chi2, chi2 / ndof
 

read_files()

reweighting_function.FluxReweighter.read_files	(		self,
			path
	)

Read one flux text file into a DataFrame with named columns.

POWHEG samples carry an extra trailing ``iEvent`` column; this is
detected from the file name so the column names line up either way.

Definition at line 289 of file reweighting_function.py.

    def read_files(self, path):
        """Read one flux text file into a DataFrame with named columns.
 
        POWHEG samples carry an extra trailing ``iEvent`` column; this is
        detected from the file name so the column names line up either way.
        """
        cols = self.COLS + ["iEvent"] if "powheg" in path.lower() else self.COLS
        return pd.read_csv(path, sep=r"\s+", header=None, names=cols)
 

run()

reweighting_function.FluxReweighter.run

(

self

)

Fit every degree in ``polynomial_degrees`` and save one plot each.

This is the "reweight" mode entry point: it is the diagnostic pass used
to choose a degree and check the binning and quality gates. Degrees
with too few surviving points are skipped with a message. Returns
``self``.

Definition at line 586 of file reweighting_function.py.

    def run(self):
        """Fit every degree in ``polynomial_degrees`` and save one plot each.
 
        This is the "reweight" mode entry point: it is the diagnostic pass used
        to choose a degree and check the binning and quality gates. Degrees
        with too few surviving points are skipped with a message. Returns
        ``self``.
        """
        os.makedirs(self.cfg["outpath"], exist_ok=True)   # create BEFORE saving
        if self.E is None:
            self.compute_ratio()
 
        fit_mask = self._fit_mask(verbose=True)
        n_fit = int(fit_mask.sum())
 
        for deg in self.cfg["degrees"]:
            if n_fit < deg + 1:
                print(f"[skip] deg={deg}: only {n_fit} fit points (need {deg + 1})")
                continue
            poly, chi2, chi2red = self.fit_poly_weighted(
                self.E[fit_mask], self.w[fit_mask], self.sig[fit_mask], degree=deg
            )
            print("Saved:", self._make_plot(deg, poly, chi2red, fit_mask))
 
        print("Output dir:", self.cfg["outpath"])
        return self
 
 
# ==========================================================================
# Run-level configuration (Analysis.yaml)
# ==========================================================================
 
#: PDG codes used when selecting the incoming neutrino in the ROOT files.

Member Data Documentation

bins_E

reweighting_function.FluxReweighter.bins_E

Definition at line 148 of file reweighting_function.py.

cfg

reweighting_function.FluxReweighter.cfg

Definition at line 147 of file reweighting_function.py.

COLS

list reweighting_function.FluxReweighter.COLS

static

Initial value:

=  [
        "Neutrino", "Parent",
        "x in m", "y in m", "z in m",
        "x angle in rad", "y angle in rad",
        "E in GeV", "weight in pb",
    ]

Definition at line 131 of file reweighting_function.py.

E

reweighting_function.FluxReweighter.E

Definition at line 152 of file reweighting_function.py.

H_DPM

reweighting_function.FluxReweighter.H_DPM

Definition at line 153 of file reweighting_function.py.

H_LC

reweighting_function.FluxReweighter.H_LC

Definition at line 436 of file reweighting_function.py.

PID

dict reweighting_function.FluxReweighter.PID = {"nue": 12, "numu": 14, "nutau": 16}

static

Definition at line 127 of file reweighting_function.py.

s_DPM

reweighting_function.FluxReweighter.s_DPM

Definition at line 431 of file reweighting_function.py.

s_LC

reweighting_function.FluxReweighter.s_LC

Definition at line 437 of file reweighting_function.py.

sig

reweighting_function.FluxReweighter.sig

Definition at line 451 of file reweighting_function.py.

w

reweighting_function.FluxReweighter.w

Definition at line 451 of file reweighting_function.py.

---

The documentation for this class was generated from the following file:

• macro/reweighting_function.py