From 7f7d96d7dbb41384799f6a730b9d6eb43d26f3bb Mon Sep 17 00:00:00 2001 From: saganatt <8majak8@gmail.com> Date: Wed, 24 Sep 2025 12:14:46 +0200 Subject: [PATCH 1/7] Pick roofitter.py changes from pass6 --- machine_learning_hep/fitting/roofitter.py | 178 ++++++++++++++++------ 1 file changed, 133 insertions(+), 45 deletions(-) diff --git a/machine_learning_hep/fitting/roofitter.py b/machine_learning_hep/fitting/roofitter.py index 41ef0f660c..a7a00b1b24 100644 --- a/machine_learning_hep/fitting/roofitter.py +++ b/machine_learning_hep/fitting/roofitter.py @@ -14,12 +14,12 @@ """Definition of the RooFitter class and helper functions""" -from math import sqrt +from math import sqrt, isnan import ROOT from ROOT import RooAddPdf, RooArgList, RooArgSet, RooFit, RooRealVar, TPaveText -USE_EXTMODEL = True +USE_EXTMODEL = False # pylint: disable=too-few-public-methods, too-many-statements # (temporary until we add more functionality) @@ -28,21 +28,66 @@ class RooFitter: def __init__(self): ROOT.gErrorIgnoreLevel = ROOT.kError - ROOT.RooMsgService.instance().setSilentMode(True) + #ROOT.RooMsgService.instance().setSilentMode(True) ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.WARNING) ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.ERROR) + def find_best_a0(self, ws, m, dh, model, old_res, range_m = None): + kwargs = {"Save": True, + "PrintLevel": -1, + "Strategy": 2} + #"MaxCalls": 5000} + #"Minimizer": "Minuit2"} + if range_m: + kwargs["Range"] = (range_m[0], range_m[1]) + tmp_frame = m.frame() + dh.plotOn(tmp_frame, ROOT.RooFit.Name("data")) + model.plotOn(tmp_frame) + chi2 = tmp_frame.chiSquare() + a0_var = ws.var("a0") # Get the a0 parameter + attempt = 0 + a0_values = [10, 20, 30, 40, 50, 80, 100, 120, 150, 200, 250, 300, 350, 400, 450, 500, + 700, 1000, 1500, 2000, 2500, 3000, 5000, 10000] + + res = old_res + chi_threshold = 6. + while (chi2 > chi_threshold or isnan(chi2)) and attempt < len(a0_values): + print(f"Attempt {attempt+1}: Setting a0 to {a0_values[attempt]}") + a0_var.setVal(a0_values[attempt]) # Change a0 value + attempt += 1 + + res = model.fitTo(dh, **kwargs) + tmp_frame = m.frame() + dh.plotOn(tmp_frame) + model.plotOn(tmp_frame) + chi2 = tmp_frame.chiSquare() + + if chi2 <= chi_threshold: + print(f"Fit improved: chi2 = {chi2}, stopping adjustments.") + + return res, model + + # pylint: disable=too-many-branches def fit_mass_new( - self, hist, pdfnames: dict, fit_spec: dict, level: str, roows: ROOT.RooWorkspace = None, plot: bool = False + self, hist, pdfnames: dict, param_names: dict, fit_spec: dict, level: str, + fixed_sigma: bool = False, fixed_sigma_val: float = 0., roows: ROOT.RooWorkspace = None, plot: bool = False ): """New fit method""" + print(f"fit spec:\n{fit_spec}") + print(f'fitting for level {level} pt: {fit_spec.get("ptrange", "")}') if hist.GetEntries() == 0: raise UserWarning("Cannot fit histogram with no entries") ws = roows or ROOT.RooWorkspace("ws") var_m = fit_spec.get("var", "m") - n_signal = RooRealVar("n_signal", "Number of signal events", 1e7, 0, 1e10) - n_background = RooRealVar("n_background", "Number of background events", 1e7, 0, 1e10) + hist_integral = hist.Integral(*(hist.FindBin(mmass) for mmass in fit_spec.get("range"))) + if "data" in level: + print(f"hist integral: {hist_integral}") + #n_signal = RooRealVar("n_signal", "Number of signal events", 0.3 * hist_integral, 0., 1.2 * hist_integral) + #n_background = RooRealVar("n_background", "Number of background events", + #0.3 * hist_integral, 0., 1.2 * hist_integral) + n_signal = RooRealVar("n_signal", "Number of signal events", 1000, 100, 100000000) + n_background = RooRealVar("n_background", "Number of background events", 1000, 100, 100000000) model = None for comp, spec in fit_spec.get("components", {}).items(): @@ -53,35 +98,48 @@ def fit_mass_new( raise ValueError("model not set") m = ws.var(var_m) + print(f"m var: {m}") - if level == "data" and USE_EXTMODEL: + if level == "mc": + sigma_sgn = ws.var(param_names["gauss_sigma"]) + if fixed_sigma: + sigma_sgn.setVal(fixed_sigma_val) + sigma_sgn.setConstant(True) + + if level == "data": signal_pdf = ws.pdf(pdfnames["pdf_sig"]) if not signal_pdf: raise ValueError("sig PDF not found") background_pdf = ws.pdf(pdfnames["pdf_bkg"]) if not background_pdf: raise ValueError("bkg pdf not found") - extmodel = RooAddPdf( + model = RooAddPdf( "model", "Total model", RooArgList(signal_pdf, background_pdf), RooArgList(n_signal, n_background) ) dh = ROOT.RooDataHist("dh", "dh", [m], Import=hist) if range_m := fit_spec.get("range"): m.setRange("fit", *range_m) - # print(f'using fit range: {range_m}, var range: {m.getRange("fit")}') - res = model.fitTo(dh, Range=(range_m[0], range_m[1]), Save=True, PrintLevel=-1, Strategy=1, MaxCalls=5000) - if level == 'data' and USE_EXTMODEL: - for v in ws.allVars(): - v.setConstant(True) - res = extmodel.fitTo( - dh, Range=(range_m[0], range_m[1]), Save=True, PrintLevel=-1, Strategy=1, MaxCalls=5000 - ) + print(f'using fit range: {range_m}, var range: {m.getRange("fit")}') + res = model.fitTo(dh, Range=(range_m[0], range_m[1]), Save=True, PrintLevel=-1, Strategy=2) + ret_model = model + #if level == "data" and USE_EXTMODEL: + # for v in ws.allVars(): + # v.setConstant(True) + # res, ret_model = self.find_best_a0(ws, m, dh, extmodel, res, range_m) + #elif + if level == "data": + res, ret_model = self.find_best_a0(ws, m, dh, model, res, range_m) else: - res = model.fitTo(dh, Save=True, PrintLevel=-1, Strategy=1, MaxCalls=5000) - if level == 'data' and USE_EXTMODEL: - for v in ws.allVars(): - v.setConstant(True) - res = extmodel.fitTo(dh, Save=True, PrintLevel=-1, Strategy=1, MaxCalls=5000) + res = model.fitTo(dh, Save=True, PrintLevel=-1, Strategy=2) + ret_model = model + #if level == "data" and USE_EXTMODEL: + # for v in ws.allVars(): + # v.setConstant(True) + # res, ret_model = self.find_best_a0(ws, m, dh, extmodel, res) + #elif + if level == "data": + res, ret_model = self.find_best_a0(ws, m, dh, model, res) frame = None residual_frame = None if plot: @@ -90,8 +148,8 @@ def fit_mass_new( c.cd() frame = m.frame() dh.plotOn(frame, ROOT.RooFit.Name("data")) - model.plotOn(frame) - model.paramOn(frame, Layout=(0.65, 1.0, 0.9)) + ret_model.plotOn(frame) + ret_model.paramOn(frame, Layout=(0.65, 1.0, 0.9)) frame.getAttText().SetTextFont(42) frame.getAttText().SetTextSize(0.001) if range_m: @@ -99,9 +157,9 @@ def fit_mass_new( frame.SetAxisRange(0.0, frame.GetMaximum() + (frame.GetMaximum() * 0.3), "Y") try: - for pdf in model.pdfList(): + for pdf in ret_model.pdfList(): pdf_name = pdf.GetName() - model.plotOn( + ret_model.plotOn( frame, ROOT.RooFit.Components(pdf), ROOT.RooFit.Name(f"pdf_{pdf_name}"), @@ -110,7 +168,7 @@ def fit_mass_new( ROOT.RooFit.LineWidth(1), ) # model.SetName("bkg") - model.plotOn(frame, ROOT.RooFit.Name("model")) + ret_model.plotOn(frame, ROOT.RooFit.Name("model")) except: # pylint: disable=bare-except # noqa: E722 pass # for comp in fit_spec.get('components', {}): @@ -120,7 +178,7 @@ def fit_mass_new( # c.Modified() # c.Update() - if level == "data" and USE_EXTMODEL and frame is not None: + if level == "data": #and USE_EXTMODEL and frame is not None: residuals = frame.residHist("data", "pdf_bkg") residual_frame = m.frame() residual_frame.addPlotable(residuals, "P") @@ -138,7 +196,7 @@ def fit_mass_new( residual_frame.SetAxisRange(range_m[0], range_m[1], "X") residual_frame.SetYTitle("Residuals") - return (res, ws, frame, residual_frame) + return (res, ws, frame, residual_frame, dh, ret_model) def fit_mass(self, hist, fit_spec, plot=False): """Old fit method""" @@ -171,10 +229,23 @@ def fit_mass(self, hist, fit_spec, plot=False): return (res, ws, frame) +def estimate_signal(hist, fit_spec, n_bkg, bkg_integral): + bin_min = hist.FindBin(fit_spec["mass_mean"] - 3 * fit_spec["sigma_signal"]) + bin_max = hist.FindBin(fit_spec["mass_mean"] + 3 * fit_spec["sigma_signal"]) + msum = 0. + for ind in range(bin_min, bin_max + 1): + msum += hist.GetBinContent(ind) + bkg = calculate_background(n_bkg, bkg_integral) + return msum - bkg + + +def calculate_background(n_bkg, bkg_integral): + return n_bkg.getVal() * bkg_integral + def calc_signif(roows, res, pdfnames, param_names, mean_sgn, sigma_sgn): """Calculate significance, signal, background, signal/background ratio.""" - if not USE_EXTMODEL: - return (0., 0., 0., 0., 0., 0, 0, 0.) + #if not USE_EXTMODEL: + # return (0., 0., 0., 0., 0., 0, 0, 0.) f_sig = roows.pdf(pdfnames["pdf_sig"]) n_signal = res.floatParsFinal().find("n_signal").getVal() sigma_n_signal = res.floatParsFinal().find("n_signal").getError() @@ -200,7 +271,12 @@ def calc_signif(roows, res, pdfnames, param_names, mean_sgn, sigma_sgn): n_signal_signal = signal_integral.getVal() * n_signal n_bkg_signal = bkg_integral.getVal() * n_bkg - significance = n_signal_signal / sqrt(n_signal_signal + n_bkg_signal) + print(f"n signal signal: {n_signal_signal} n bkg signal: {n_bkg_signal}") + + if n_signal_signal + n_bkg_signal == 0.: + significance = 0.0 + else: + significance = n_signal_signal / sqrt(n_signal_signal + n_bkg_signal) # Calculate the error on the signal and bkg integrals using the covariance matrix sigma_signal_integral = signal_integral.getPropagatedError(res) @@ -211,17 +287,29 @@ def calc_signif(roows, res, pdfnames, param_names, mean_sgn, sigma_sgn): ) sigma_n_bkg_signal = sqrt((bkg_integral.getVal() * sigma_n_bkg) ** 2 + (n_bkg * sigma_bkg_integral) ** 2) - dS_dS = 1 / sqrt(n_signal_signal + n_bkg_signal) - ( - n_signal_signal / (2 * (n_signal_signal + n_bkg_signal) ** (3 / 2)) - ) - dS_dB = -n_signal_signal / (2 * (n_signal_signal + n_bkg_signal) ** (3 / 2)) - significance_err = sqrt((dS_dS * sigma_n_signal_signal) ** 2 + (dS_dB * sigma_n_bkg_signal) ** 2) - - # Signal to bkg ratio - s_over_b = n_signal_signal / n_bkg_signal - s_over_b_err = s_over_b * sqrt( - (sigma_n_signal_signal / n_signal_signal) ** 2 + (sigma_n_bkg_signal / n_bkg_signal) ** 2 - ) + if n_signal_signal + n_bkg_signal == 0.: + dS_dS = 0.0 + dS_dB = 0.0 + else: + dS_dS = (1 / sqrt(n_signal_signal + n_bkg_signal) - + (n_signal_signal / (2 * (n_signal_signal + n_bkg_signal)**(3/2)))) + dS_dB = -n_signal_signal / (2 * (n_signal_signal + n_bkg_signal)**(3/2)) + significance_err = sqrt( + (dS_dS * sigma_n_signal_signal) ** 2 + + (dS_dB * sigma_n_bkg_signal) ** 2) + + #Signal to bkg ratio + if n_bkg_signal == 0.: + s_over_b = 0.0 + s_over_b_err = 0.0 # as S/B is ill-defined + elif n_signal_signal == 0.: + s_over_b = 0.0 + s_over_b_err = s_over_b * sqrt((sigma_n_bkg_signal / n_bkg_signal) ** 2) + else: + s_over_b = n_signal_signal / n_bkg_signal + s_over_b_err = ( + s_over_b * sqrt((sigma_n_signal_signal / n_signal_signal) ** 2 + + (sigma_n_bkg_signal / n_bkg_signal) ** 2)) return ( n_signal_signal, @@ -258,11 +346,11 @@ def add_text_info_fit(text_info, frame, roows, param_names): refl_frac = roows.var(param_names["fraction_refl"]) text_info.AddText(f"#chi^{{2}}/ndf = {chi2:.2f}") text_info.AddText(f"#mu = {mean_sgn.getVal():.3f} #pm {mean_sgn.getError():.3f}") - text_info.AddText(f"#sigma = {sigma_sgn.getVal():.3f} #pm {sigma_sgn.getError():.3f}") + text_info.AddText(f"#sigma = {sigma_sgn.getVal():.4f} #pm {sigma_sgn.getError():.4f}") if sigmawide_sgn: - text_info.AddText(f"#sigma wide = {sigmawide_sgn.getVal():.3f} #pm {sigmawide_sgn.getError():.3f}") + text_info.AddText(f"#sigma wide = {sigmawide_sgn.getVal():.4f} #pm {sigmawide_sgn.getError():.4f}") if refl_frac: - text_info.AddText(f"refl.frac. = {refl_frac.getVal():.3f} #pm {refl_frac.getError():.3f}") + text_info.AddText(f"refl.frac. = {refl_frac.getVal():.4f} #pm {refl_frac.getError():.4f}") if a0 := roows.var("a0"): text_info.AddText(f"a0 = {a0.getVal():.3f} #pm {a0.getError():.3f}") if a1 := roows.var("a1"): From b77d689e7b541076c0b3cb27783a95e7f40a2dd4 Mon Sep 17 00:00:00 2001 From: saganatt <8majak8@gmail.com> Date: Wed, 24 Sep 2025 12:34:01 +0200 Subject: [PATCH 2/7] Take straight the fitting changes from pass6 in analysis/ --- .../analysis/analyzerdhadrons.py | 84 +++++++++++++++---- .../analysis/analyzerdhadrons_mult.py | 33 +++++--- 2 files changed, 90 insertions(+), 27 deletions(-) diff --git a/machine_learning_hep/analysis/analyzerdhadrons.py b/machine_learning_hep/analysis/analyzerdhadrons.py index 24672b0303..7ebc357050 100644 --- a/machine_learning_hep/analysis/analyzerdhadrons.py +++ b/machine_learning_hep/analysis/analyzerdhadrons.py @@ -37,6 +37,8 @@ gStyle, kBlue, kCyan, + RooConstVar, + RooArgSet, ) from machine_learning_hep.analysis.analyzer import Analyzer @@ -100,6 +102,8 @@ def __init__(self, datap, case, typean, period): self.n_fileff = os.path.join(self.d_resultsallpmc, self.n_fileff) self.p_bin_width = datap["analysis"][self.typean]["bin_width"] self.p_rebin = datap["analysis"][self.typean]["n_rebin"] + self.p_fixed_sigma = datap["analysis"][self.typean]["fixed_sigma"] + self.p_fixed_sigma_val = datap["analysis"][self.typean]["fixed_sigma_val"] self.p_pdfnames = datap["analysis"][self.typean]["pdf_names"] self.p_param_names = datap["analysis"][self.typean]["param_names"] @@ -169,10 +173,18 @@ def _save_hist(self, hist, filename, option=""): self.rfigfile.WriteObject(hist, rfilename) # region fitting - def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=None, filename=None): + def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, fixed_sigma, fixed_sigma_val, # pylint: disable=too-many-arguments + roows=None, filename=None): if fitcfg is None: + return None, None, None, None, None + try: + res, ws, frame, residual_frame, data_hist, model = self.fitter.fit_mass_new(hist, pdfnames, param_names, + fitcfg, level, + fixed_sigma, fixed_sigma_val, + roows, True) + except ValueError: + self.logger.error("Could not do fitting on %s {level} for pt %d - %d", level, self.bins_candpt[ipt], self.bins_candpt[ipt+1]) return None, None - res, ws, frame, residual_frame = self.fitter.fit_mass_new(hist, pdfnames, fitcfg, level, roows, True) frame.SetTitle(f"inv. mass for p_{{T}} {self.bins_candpt[ipt]} - {self.bins_candpt[ipt + 1]} GeV/c") c = TCanvas() @@ -180,7 +192,7 @@ def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=No add_text_info_fit(textInfoRight, frame, ws, param_names) textInfoLeft = create_text_info(0.12, 0.68, 0.6, 0.89) - if level == "data": + if res and level == "data": mean_sgn = ws.var(self.p_param_names["gauss_mean"]) sigma_sgn = ws.var(self.p_param_names["gauss_sigma"]) (sig, sig_err, bkg, bkg_err, signif, signif_err, s_over_b, s_over_b_err) = calc_signif( @@ -193,7 +205,7 @@ def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=No textInfoRight.Draw() textInfoLeft.Draw() - if res.status() == 0: + if res and res.status() == 0: self._save_canvas(c, filename) else: self.logger.warning("Invalid fit result for %s", hist.GetName()) @@ -201,7 +213,7 @@ def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=No filename = filename.replace(".png", "_invalid.png") self._save_canvas(c, filename) - if level == "data": + if level == "data" and residual_frame: residual_frame.SetTitle( f"inv. mass for p_{{T}} {self.bins_candpt[ipt]} - {self.bins_candpt[ipt + 1]} GeV/c" ) @@ -210,7 +222,9 @@ def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=No filename = filename.replace(".png", "_residual.png") self._save_canvas(cres, filename) - return res, ws + chi = frame.chiSquare() + + return res, ws, chi, data_hist, model def _fit_mass(self, hist, filename=None): if hist.GetEntries() == 0: @@ -274,22 +288,23 @@ def fit(self): fileout_name = self.make_file_path( self.d_resultsallpdata, self.yields_filename, "root", None, [self.case, self.typean] ) - fileout = TFile(fileout_name, "RECREATE") + # fileout = TFile(fileout_name, "RECREATE") yieldshistos = TH1F("hyields0", "", len(self.lpt_finbinmin), array("d", self.bins_candpt)) meanhistos = TH1F("hmean0", "", len(self.lpt_finbinmin), array("d", self.bins_candpt)) sigmahistos = TH1F("hsigmas0", "", len(self.lpt_finbinmin), array("d", self.bins_candpt)) signifhistos = TH1F("hsignifs0", "", len(self.lpt_finbinmin), array("d", self.bins_candpt)) soverbhistos = TH1F("hSoverB0", "", len(self.lpt_finbinmin), array("d", self.bins_candpt)) + chihistos = TH1F("hchi0", "", len(self.lpt_finbinmin), array("d", self.bins_candpt)) lpt_probcutfin = [None] * self.nbins - with TFile(rfilename) as rfile: + with TFile(rfilename) as rfile, TFile(fileout_name, "RECREATE") as fileout: for ipt in range(len(self.lpt_finbinmin)): lpt_probcutfin[ipt] = self.lpt_probcutfin_tmp[self.bin_matching[ipt]] self.logger.debug("fitting %s - %i", level, ipt) roows = self.roows.get(ipt) if self.mltype == "MultiClassification": - suffix = "%s%d_%d_%.2f%.2f%.2f" % ( + suffix = "%s%d_%d_%.2f%.2f%.3f" % ( self.v_var_binning, self.lpt_finbinmin[ipt], self.lpt_finbinmax[ipt], @@ -304,12 +319,16 @@ def fit(self): self.lpt_finbinmax[ipt], lpt_probcutfin[ipt], ) + rfile.cd() h_invmass = rfile.Get("hmass" + suffix) # Rebin + self.logger.info("suffix: %s ipt %d", suffix, ipt) h_invmass.Rebin(self.p_rebin[ipt]) if h_invmass.GetEntries() < 100: # TODO: reconsider criterion self.logger.error("Not enough entries to fit for %s bin %d", level, ipt) continue + if "data" in level: + self.logger.info("hist entries: %d, class: %s", h_invmass.GetEntries(), h_invmass.ClassName()) ptrange = (self.bins_candpt[ipt], self.bins_candpt[ipt + 1]) if self.cfg("mass_fit"): @@ -341,19 +360,23 @@ def fit(self): roows.var(fixpar).setConstant(True) if h_invmass.GetEntries() == 0: continue - roo_res, roo_ws = self._roofit_mass( + roo_res, roo_ws, chi, dh, model = self._roofit_mass( level, h_invmass, ipt, self.p_pdfnames, self.p_param_names, fitcfg, + self.p_fixed_sigma[ipt], + self.p_fixed_sigma_val[ipt], roows, f"roofit/h_mass_fitted_pthf-{ptrange[0]}-{ptrange[1]}_{level}.png", ) + if roo_res: + roo_res.Print() self.roo_ws[level][ipt] = roo_ws self.roows[ipt] = roo_ws - if roo_res.status() == 0: + if roo_res and roo_res.status() == 0: if level in ("data", "mc_sig"): self.fit_mean[level][ipt] = roo_ws.var(self.p_param_names["gauss_mean"]).getValV() self.fit_sigma[level][ipt] = roo_ws.var(self.p_param_names["gauss_sigma"]).getValV() @@ -371,9 +394,34 @@ def fit(self): if level == "data": mean_sgn = roo_ws.var(self.p_param_names["gauss_mean"]) sigma_sgn = roo_ws.var(self.p_param_names["gauss_sigma"]) - (sig, sig_err, _, _, signif, signif_err, s_over_b, s_over_b_err) = calc_signif( - roo_ws, roo_res, self.p_pdfnames, self.p_param_names, mean_sgn, sigma_sgn - ) + if roo_res: + (sig, sig_err, _, _, + signif, signif_err, s_over_b, s_over_b_err + ) = calc_signif(roo_ws, roo_res, self.p_pdfnames, self.p_param_names, mean_sgn, sigma_sgn) + fileout.cd() + one = RooConstVar("one", "constant 1.0", 1.0) + bkg_pdf = roo_ws.pdf(self.p_pdfnames["pdf_bkg"]) + sig_pdf = roo_ws.pdf(self.p_pdfnames["pdf_sig"]) + if not model: + self.logger.info("Model is null") + for pdf, outlabel in zip((bkg_pdf, sig_pdf, model), ("bkg", "sgn", "total")): + if not pdf: + self.logger.info("Pdf null") + continue + self.logger.info("Pdf %s", pdf) + obs = pdf.getObservables(dh) + self.logger.info("Observables %s", obs) + params = pdf.getParameters(dh) + self.logger.info("Parameters %s", params) + fit_func = pdf.asTF(obs, params, RooArgSet(one)) + self.logger.info("TF fit func %s", fit_func) + fit_func.Write(f"{outlabel}TF_{ptrange[0]:.0f}_{ptrange[1]:.0f}") + self.logger.info("Wrote TF func") + + h_invmass.Write(f"hmass_{ipt}") + self.logger.info("Wrote hist mass") + else: + sig = sig_err = signif = signif_err = s_over_b = s_over_b_err = 0.0 yieldshistos.SetBinContent(ipt + 1, sig) yieldshistos.SetBinError(ipt + 1, sig_err) @@ -385,13 +433,16 @@ def fit(self): signifhistos.SetBinError(ipt + 1, signif_err) soverbhistos.SetBinContent(ipt + 1, s_over_b) soverbhistos.SetBinError(ipt + 1, s_over_b_err) + chihistos.SetBinContent(ipt + 1, chi) + chihistos.SetBinError(ipt + 1, 0) fileout.cd() yieldshistos.Write() meanhistos.Write() sigmahistos.Write() signifhistos.Write() soverbhistos.Write() - fileout.Close() + chihistos.Write() + #fileout.Close() def yield_syst(self): # Enable ROOT batch mode and reset in the end @@ -470,6 +521,7 @@ def efficiency(self): legeffFD.Draw() cEffFD.SaveAs(f"{self.d_resultsallpmc}/EffFD{self.case}{self.typean}.eps") + @staticmethod def calculate_norm(logger, hevents, hselevents): # TO BE FIXED WITH EV SEL if not hevents: @@ -516,8 +568,6 @@ def makenormyields(self): # pylint: disable=import-outside-toplevel, too-many-b histonorm.SetBinContent(1, selnorm) self.logger.warning("Number of events %d", norm) - self.logger.warning("Number of events after event selection %d", selnorm) - if self.p_dobkgfromsideband: fileoutbkg = TFile.Open(f"{self.d_resultsallpdata}/Background_fromsidebands_{self.case}_{self.typean}.root") hbkg = fileoutbkg.Get("hbkg_fromsidebands") diff --git a/machine_learning_hep/analysis/analyzerdhadrons_mult.py b/machine_learning_hep/analysis/analyzerdhadrons_mult.py index cd62e163eb..596e7b34b5 100644 --- a/machine_learning_hep/analysis/analyzerdhadrons_mult.py +++ b/machine_learning_hep/analysis/analyzerdhadrons_mult.py @@ -200,7 +200,11 @@ def _save_hist(self, hist, filename, option=""): def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=None, filename=None): if fitcfg is None: return None, None - res, ws, frame, residual_frame = self.fitter.fit_mass_new(hist, pdfnames, fitcfg, level, roows, True) + try: + res, ws, frame, residual_frame = self.fitter.fit_mass_new(hist, pdfnames, param_names, fitcfg, level, roows=roows, plot=True) + except ValueError: + self.logger.error(f"Could not do fitting on {level} for pt {self.bins_candpt[ipt]} - {self.bins_candpt[ipt+1]}") + return None, None frame.SetTitle(f"inv. mass for p_{{T}} {self.bins_candpt[ipt]} - {self.bins_candpt[ipt + 1]} GeV/c") c = TCanvas() @@ -208,7 +212,7 @@ def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=No add_text_info_fit(textInfoRight, frame, ws, param_names) textInfoLeft = create_text_info(0.12, 0.68, 0.6, 0.89) - if level == "data": + if res and level == "data": mean_sgn = ws.var(self.p_param_names["gauss_mean"]) sigma_sgn = ws.var(self.p_param_names["gauss_sigma"]) (sig, sig_err, bkg, bkg_err, signif, signif_err, s_over_b, s_over_b_err) = calc_signif( @@ -221,7 +225,7 @@ def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=No textInfoRight.Draw() textInfoLeft.Draw() - if res.status() == 0: + if res and res.status() == 0: self._save_canvas(c, filename) else: self.logger.warning("Invalid fit result for %s", hist.GetName()) @@ -238,7 +242,9 @@ def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=No filename = filename.replace(".png", "_residual.png") self._save_canvas(cres, filename) - return res, ws + chi = frame.chiSquare() + + return res, ws, chi def _fit_mass(self, hist, filename=None): if hist.GetEntries() == 0: @@ -315,8 +321,8 @@ def fit(self): soverbhistos = TH1F( "hSoverB%d" % (ibin2), "", len(self.lpt_finbinmin), array("d", self.bins_candpt) ) + chihistos = TH1F("hchi%d" % (ibin2), "", len(self.lpt_finbinmin), array("d", self.bins_candpt)) - lpt_probcutfin = [None] * self.nbins for ipt in range(len(self.lpt_finbinmin)): lpt_probcutfin[ipt] = self.lpt_probcutfin_tmp[self.bin_matching[ipt]] self.logger.debug("fitting %s - %i - %i", level, ipt, ibin2) @@ -391,7 +397,7 @@ def fit(self): # Create the directory if it doesn't exist directory_path.mkdir(parents=True, exist_ok=True) - roo_res, roo_ws = self._roofit_mass( + roo_res, roo_ws, chi = self._roofit_mass( level, h_invmass, ipt, @@ -407,7 +413,7 @@ def fit(self): # roo_ws.Print() self.roo_ws[level][ipt] = roo_ws self.roows[ipt] = roo_ws - if roo_res.status() == 0: + if roo_res and roo_res.status() == 0: if level in ("data", "mc_sig"): self.fit_mean[level][ipt] = roo_ws.var(self.p_param_names["gauss_mean"]).getValV() self.fit_sigma[level][ipt] = roo_ws.var(self.p_param_names["gauss_sigma"]).getValV() @@ -425,9 +431,13 @@ def fit(self): if level == "data": mean_sgn = roo_ws.var(self.p_param_names["gauss_mean"]) sigma_sgn = roo_ws.var(self.p_param_names["gauss_sigma"]) - (sig, sig_err, _, _, signif, signif_err, s_over_b, s_over_b_err) = calc_signif( - roo_ws, roo_res, self.p_pdfnames, self.p_param_names, mean_sgn, sigma_sgn - ) + if roo_res: + (sig, sig_err, _, _, + signif, signif_err, s_over_b, s_over_b_err + ) = calc_signif(roo_ws, roo_res, self.p_pdfnames, \ + self.p_param_names, mean_sgn, sigma_sgn) + else: + sig = sig_err = signif = signif_err = s_over_b = s_over_b_err = 0.0 yieldshistos.SetBinContent(ipt + 1, sig) yieldshistos.SetBinError(ipt + 1, sig_err) @@ -439,12 +449,15 @@ def fit(self): signifhistos.SetBinError(ipt + 1, signif_err) soverbhistos.SetBinContent(ipt + 1, s_over_b) soverbhistos.SetBinError(ipt + 1, s_over_b_err) + chihistos.SetBinContent(ipt + 1, chi) + chihistos.SetBinError(ipt + 1, 0) fileout.cd() yieldshistos.Write() meanhistos.Write() sigmahistos.Write() signifhistos.Write() soverbhistos.Write() + chihistos.Write() fileout.Close() def get_efficiency(self, ibin1, ibin2): From a6c19bdb066e72cd56d132833894126e093fc926 Mon Sep 17 00:00:00 2001 From: saganatt <8majak8@gmail.com> Date: Wed, 1 Oct 2025 15:57:50 +0200 Subject: [PATCH 3/7] Drop the extended changes --- .../analysis/analyzerdhadrons.py | 44 ++---- .../analysis/analyzerdhadrons_mult.py | 25 +++- machine_learning_hep/fitting/roofitter.py | 138 +++++------------- 3 files changed, 62 insertions(+), 145 deletions(-) diff --git a/machine_learning_hep/analysis/analyzerdhadrons.py b/machine_learning_hep/analysis/analyzerdhadrons.py index 7ebc357050..5bcac1bbb1 100644 --- a/machine_learning_hep/analysis/analyzerdhadrons.py +++ b/machine_learning_hep/analysis/analyzerdhadrons.py @@ -37,8 +37,6 @@ gStyle, kBlue, kCyan, - RooConstVar, - RooArgSet, ) from machine_learning_hep.analysis.analyzer import Analyzer @@ -102,8 +100,6 @@ def __init__(self, datap, case, typean, period): self.n_fileff = os.path.join(self.d_resultsallpmc, self.n_fileff) self.p_bin_width = datap["analysis"][self.typean]["bin_width"] self.p_rebin = datap["analysis"][self.typean]["n_rebin"] - self.p_fixed_sigma = datap["analysis"][self.typean]["fixed_sigma"] - self.p_fixed_sigma_val = datap["analysis"][self.typean]["fixed_sigma_val"] self.p_pdfnames = datap["analysis"][self.typean]["pdf_names"] self.p_param_names = datap["analysis"][self.typean]["param_names"] @@ -173,18 +169,14 @@ def _save_hist(self, hist, filename, option=""): self.rfigfile.WriteObject(hist, rfilename) # region fitting - def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, fixed_sigma, fixed_sigma_val, # pylint: disable=too-many-arguments - roows=None, filename=None): + def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=None, filename=None): if fitcfg is None: - return None, None, None, None, None + return None, None, None try: - res, ws, frame, residual_frame, data_hist, model = self.fitter.fit_mass_new(hist, pdfnames, param_names, - fitcfg, level, - fixed_sigma, fixed_sigma_val, - roows, True) + res, ws, frame, residual_frame = self.fitter.fit_mass_new(hist, pdfnames, fitcfg, level, roows=roows, plot=True) except ValueError: - self.logger.error("Could not do fitting on %s {level} for pt %d - %d", level, self.bins_candpt[ipt], self.bins_candpt[ipt+1]) - return None, None + self.logger.error(f"Could not do fitting on {level} for pt {self.bins_candpt[ipt]} - {self.bins_candpt[ipt+1]}") + return None, None, None frame.SetTitle(f"inv. mass for p_{{T}} {self.bins_candpt[ipt]} - {self.bins_candpt[ipt + 1]} GeV/c") c = TCanvas() @@ -224,7 +216,7 @@ def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, fixed_si chi = frame.chiSquare() - return res, ws, chi, data_hist, model + return res, ws, chi def _fit_mass(self, hist, filename=None): if hist.GetEntries() == 0: @@ -288,7 +280,6 @@ def fit(self): fileout_name = self.make_file_path( self.d_resultsallpdata, self.yields_filename, "root", None, [self.case, self.typean] ) - # fileout = TFile(fileout_name, "RECREATE") yieldshistos = TH1F("hyields0", "", len(self.lpt_finbinmin), array("d", self.bins_candpt)) meanhistos = TH1F("hmean0", "", len(self.lpt_finbinmin), array("d", self.bins_candpt)) @@ -304,7 +295,7 @@ def fit(self): self.logger.debug("fitting %s - %i", level, ipt) roows = self.roows.get(ipt) if self.mltype == "MultiClassification": - suffix = "%s%d_%d_%.2f%.2f%.3f" % ( + suffix = "%s%d_%d_%.2f%.2f%.2f" % ( self.v_var_binning, self.lpt_finbinmin[ipt], self.lpt_finbinmax[ipt], @@ -322,13 +313,10 @@ def fit(self): rfile.cd() h_invmass = rfile.Get("hmass" + suffix) # Rebin - self.logger.info("suffix: %s ipt %d", suffix, ipt) h_invmass.Rebin(self.p_rebin[ipt]) if h_invmass.GetEntries() < 100: # TODO: reconsider criterion self.logger.error("Not enough entries to fit for %s bin %d", level, ipt) continue - if "data" in level: - self.logger.info("hist entries: %d, class: %s", h_invmass.GetEntries(), h_invmass.ClassName()) ptrange = (self.bins_candpt[ipt], self.bins_candpt[ipt + 1]) if self.cfg("mass_fit"): @@ -360,20 +348,16 @@ def fit(self): roows.var(fixpar).setConstant(True) if h_invmass.GetEntries() == 0: continue - roo_res, roo_ws, chi, dh, model = self._roofit_mass( + roo_res, roo_ws, chi = self._roofit_mass( level, h_invmass, ipt, self.p_pdfnames, self.p_param_names, fitcfg, - self.p_fixed_sigma[ipt], - self.p_fixed_sigma_val[ipt], roows, f"roofit/h_mass_fitted_pthf-{ptrange[0]}-{ptrange[1]}_{level}.png", ) - if roo_res: - roo_res.Print() self.roo_ws[level][ipt] = roo_ws self.roows[ipt] = roo_ws if roo_res and roo_res.status() == 0: @@ -402,24 +386,15 @@ def fit(self): one = RooConstVar("one", "constant 1.0", 1.0) bkg_pdf = roo_ws.pdf(self.p_pdfnames["pdf_bkg"]) sig_pdf = roo_ws.pdf(self.p_pdfnames["pdf_sig"]) - if not model: - self.logger.info("Model is null") for pdf, outlabel in zip((bkg_pdf, sig_pdf, model), ("bkg", "sgn", "total")): if not pdf: - self.logger.info("Pdf null") continue - self.logger.info("Pdf %s", pdf) obs = pdf.getObservables(dh) - self.logger.info("Observables %s", obs) params = pdf.getParameters(dh) - self.logger.info("Parameters %s", params) fit_func = pdf.asTF(obs, params, RooArgSet(one)) - self.logger.info("TF fit func %s", fit_func) fit_func.Write(f"{outlabel}TF_{ptrange[0]:.0f}_{ptrange[1]:.0f}") - self.logger.info("Wrote TF func") h_invmass.Write(f"hmass_{ipt}") - self.logger.info("Wrote hist mass") else: sig = sig_err = signif = signif_err = s_over_b = s_over_b_err = 0.0 @@ -442,7 +417,6 @@ def fit(self): signifhistos.Write() soverbhistos.Write() chihistos.Write() - #fileout.Close() def yield_syst(self): # Enable ROOT batch mode and reset in the end @@ -568,6 +542,8 @@ def makenormyields(self): # pylint: disable=import-outside-toplevel, too-many-b histonorm.SetBinContent(1, selnorm) self.logger.warning("Number of events %d", norm) + self.logger.warning("Number of events after event selection %d", selnorm) + if self.p_dobkgfromsideband: fileoutbkg = TFile.Open(f"{self.d_resultsallpdata}/Background_fromsidebands_{self.case}_{self.typean}.root") hbkg = fileoutbkg.Get("hbkg_fromsidebands") diff --git a/machine_learning_hep/analysis/analyzerdhadrons_mult.py b/machine_learning_hep/analysis/analyzerdhadrons_mult.py index 596e7b34b5..1ead07ae5f 100644 --- a/machine_learning_hep/analysis/analyzerdhadrons_mult.py +++ b/machine_learning_hep/analysis/analyzerdhadrons_mult.py @@ -199,12 +199,12 @@ def _save_hist(self, hist, filename, option=""): # region fitting def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=None, filename=None): if fitcfg is None: - return None, None + return None, None, None try: - res, ws, frame, residual_frame = self.fitter.fit_mass_new(hist, pdfnames, param_names, fitcfg, level, roows=roows, plot=True) + res, ws, frame, residual_frame = self.fitter.fit_mass_new(hist, pdfnames, fitcfg, level, roows=roows, plot=True) except ValueError: self.logger.error(f"Could not do fitting on {level} for pt {self.bins_candpt[ipt]} - {self.bins_candpt[ipt+1]}") - return None, None + return None, None, None frame.SetTitle(f"inv. mass for p_{{T}} {self.bins_candpt[ipt]} - {self.bins_candpt[ipt + 1]} GeV/c") c = TCanvas() @@ -307,8 +307,7 @@ def fit(self): fileout_name = self.make_file_path( self.d_resultsallpdata, self.yields_filename, "root", None, [self.case, self.typean] ) - fileout = TFile(fileout_name, "RECREATE") - with TFile(rfilename) as rfile: + with TFile(rfilename) as rfile, TFile(fileout_name, "RECREATE") as fileout: for ibin2 in range(len(self.lvar2_binmin)): yieldshistos = TH1F( "hyields%d" % (ibin2), "", len(self.lpt_finbinmin), array("d", self.bins_candpt) @@ -323,6 +322,7 @@ def fit(self): ) chihistos = TH1F("hchi%d" % (ibin2), "", len(self.lpt_finbinmin), array("d", self.bins_candpt)) + lpt_probcutfin = [None] * self.nbins for ipt in range(len(self.lpt_finbinmin)): lpt_probcutfin[ipt] = self.lpt_probcutfin_tmp[self.bin_matching[ipt]] self.logger.debug("fitting %s - %i - %i", level, ipt, ibin2) @@ -348,6 +348,7 @@ def fit(self): self.lvar2_binmin[ibin2], self.lvar2_binmax[ibin2], ) + rfile.cd() h_invmass = rfile.Get("hmass" + suffix) # Rebin h_invmass.Rebin(self.p_rebin[ipt]) @@ -436,6 +437,19 @@ def fit(self): signif, signif_err, s_over_b, s_over_b_err ) = calc_signif(roo_ws, roo_res, self.p_pdfnames, \ self.p_param_names, mean_sgn, sigma_sgn) + fileout.cd() + one = RooConstVar("one", "constant 1.0", 1.0) + bkg_pdf = roo_ws.pdf(self.p_pdfnames["pdf_bkg"]) + sig_pdf = roo_ws.pdf(self.p_pdfnames["pdf_sig"]) + for pdf, outlabel in zip((bkg_pdf, sig_pdf, model), ("bkg", "sgn", "total")): + if not pdf: + continue + obs = pdf.getObservables(dh) + params = pdf.getParameters(dh) + fit_func = pdf.asTF(obs, params, RooArgSet(one)) + fit_func.Write(f"{outlabel}TF_{ptrange[0]:.0f}_{ptrange[1]:.0f}") + + h_invmass.Write(f"hmass_{ipt}") else: sig = sig_err = signif = signif_err = s_over_b = s_over_b_err = 0.0 @@ -458,7 +472,6 @@ def fit(self): signifhistos.Write() soverbhistos.Write() chihistos.Write() - fileout.Close() def get_efficiency(self, ibin1, ibin2): fileouteff = TFile.Open(f"{self.d_resultsallpmc}/efficiencies{self.case}{self.typean}.root", "read") diff --git a/machine_learning_hep/fitting/roofitter.py b/machine_learning_hep/fitting/roofitter.py index a7a00b1b24..2b231c8283 100644 --- a/machine_learning_hep/fitting/roofitter.py +++ b/machine_learning_hep/fitting/roofitter.py @@ -14,12 +14,12 @@ """Definition of the RooFitter class and helper functions""" -from math import sqrt, isnan +from math import sqrt import ROOT from ROOT import RooAddPdf, RooArgList, RooArgSet, RooFit, RooRealVar, TPaveText -USE_EXTMODEL = False +USE_EXTMODEL = True # pylint: disable=too-few-public-methods, too-many-statements # (temporary until we add more functionality) @@ -28,66 +28,22 @@ class RooFitter: def __init__(self): ROOT.gErrorIgnoreLevel = ROOT.kError - #ROOT.RooMsgService.instance().setSilentMode(True) + ROOT.RooMsgService.instance().setSilentMode(True) ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.WARNING) ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.ERROR) - def find_best_a0(self, ws, m, dh, model, old_res, range_m = None): - kwargs = {"Save": True, - "PrintLevel": -1, - "Strategy": 2} - #"MaxCalls": 5000} - #"Minimizer": "Minuit2"} - if range_m: - kwargs["Range"] = (range_m[0], range_m[1]) - tmp_frame = m.frame() - dh.plotOn(tmp_frame, ROOT.RooFit.Name("data")) - model.plotOn(tmp_frame) - chi2 = tmp_frame.chiSquare() - a0_var = ws.var("a0") # Get the a0 parameter - attempt = 0 - a0_values = [10, 20, 30, 40, 50, 80, 100, 120, 150, 200, 250, 300, 350, 400, 450, 500, - 700, 1000, 1500, 2000, 2500, 3000, 5000, 10000] - - res = old_res - chi_threshold = 6. - while (chi2 > chi_threshold or isnan(chi2)) and attempt < len(a0_values): - print(f"Attempt {attempt+1}: Setting a0 to {a0_values[attempt]}") - a0_var.setVal(a0_values[attempt]) # Change a0 value - attempt += 1 - - res = model.fitTo(dh, **kwargs) - tmp_frame = m.frame() - dh.plotOn(tmp_frame) - model.plotOn(tmp_frame) - chi2 = tmp_frame.chiSquare() - - if chi2 <= chi_threshold: - print(f"Fit improved: chi2 = {chi2}, stopping adjustments.") - - return res, model - # pylint: disable=too-many-branches def fit_mass_new( - self, hist, pdfnames: dict, param_names: dict, fit_spec: dict, level: str, - fixed_sigma: bool = False, fixed_sigma_val: float = 0., roows: ROOT.RooWorkspace = None, plot: bool = False + self, hist, pdfnames: dict, fit_spec: dict, level: str, roows: ROOT.RooWorkspace = None, plot: bool = False ): """New fit method""" - print(f"fit spec:\n{fit_spec}") - print(f'fitting for level {level} pt: {fit_spec.get("ptrange", "")}') if hist.GetEntries() == 0: raise UserWarning("Cannot fit histogram with no entries") ws = roows or ROOT.RooWorkspace("ws") var_m = fit_spec.get("var", "m") - hist_integral = hist.Integral(*(hist.FindBin(mmass) for mmass in fit_spec.get("range"))) - if "data" in level: - print(f"hist integral: {hist_integral}") - #n_signal = RooRealVar("n_signal", "Number of signal events", 0.3 * hist_integral, 0., 1.2 * hist_integral) - #n_background = RooRealVar("n_background", "Number of background events", - #0.3 * hist_integral, 0., 1.2 * hist_integral) - n_signal = RooRealVar("n_signal", "Number of signal events", 1000, 100, 100000000) - n_background = RooRealVar("n_background", "Number of background events", 1000, 100, 100000000) + n_signal = RooRealVar("n_signal", "Number of signal events", 1e7, 0, 1e10) + n_background = RooRealVar("n_background", "Number of background events", 1e7, 0, 1e10) model = None for comp, spec in fit_spec.get("components", {}).items(): @@ -98,48 +54,35 @@ def fit_mass_new( raise ValueError("model not set") m = ws.var(var_m) - print(f"m var: {m}") - - if level == "mc": - sigma_sgn = ws.var(param_names["gauss_sigma"]) - if fixed_sigma: - sigma_sgn.setVal(fixed_sigma_val) - sigma_sgn.setConstant(True) - if level == "data": + if level == "data" and USE_EXTMODEL: signal_pdf = ws.pdf(pdfnames["pdf_sig"]) if not signal_pdf: raise ValueError("sig PDF not found") background_pdf = ws.pdf(pdfnames["pdf_bkg"]) if not background_pdf: raise ValueError("bkg pdf not found") - model = RooAddPdf( + extmodel = RooAddPdf( "model", "Total model", RooArgList(signal_pdf, background_pdf), RooArgList(n_signal, n_background) ) dh = ROOT.RooDataHist("dh", "dh", [m], Import=hist) if range_m := fit_spec.get("range"): m.setRange("fit", *range_m) - print(f'using fit range: {range_m}, var range: {m.getRange("fit")}') - res = model.fitTo(dh, Range=(range_m[0], range_m[1]), Save=True, PrintLevel=-1, Strategy=2) - ret_model = model - #if level == "data" and USE_EXTMODEL: - # for v in ws.allVars(): - # v.setConstant(True) - # res, ret_model = self.find_best_a0(ws, m, dh, extmodel, res, range_m) - #elif - if level == "data": - res, ret_model = self.find_best_a0(ws, m, dh, model, res, range_m) + # print(f'using fit range: {range_m}, var range: {m.getRange("fit")}') + res = model.fitTo(dh, Range=(range_m[0], range_m[1]), Save=True, PrintLevel=-1, Strategy=1, MaxCalls=5000) + if level == "data" and USE_EXTMODEL: + for v in ws.allVars(): + v.setConstant(True) + res = extmodel.fitTo( + dh, Range=(range_m[0], range_m[1]), Save=True, PrintLevel=-1, Strategy=1, MaxCalls=5000 + ) else: - res = model.fitTo(dh, Save=True, PrintLevel=-1, Strategy=2) - ret_model = model - #if level == "data" and USE_EXTMODEL: - # for v in ws.allVars(): - # v.setConstant(True) - # res, ret_model = self.find_best_a0(ws, m, dh, extmodel, res) - #elif - if level == "data": - res, ret_model = self.find_best_a0(ws, m, dh, model, res) + res = model.fitTo(dh, Save=True, PrintLevel=-1, Strategy=1, MaxCalls=5000) + if level == "data" and USE_EXTMODEL: + for v in ws.allVars(): + v.setConstant(True) + res = extmodel.fitTo(dh, Save=True, PrintLevel=-1, Strategy=1, MaxCalls=5000) frame = None residual_frame = None if plot: @@ -148,8 +91,8 @@ def fit_mass_new( c.cd() frame = m.frame() dh.plotOn(frame, ROOT.RooFit.Name("data")) - ret_model.plotOn(frame) - ret_model.paramOn(frame, Layout=(0.65, 1.0, 0.9)) + model.plotOn(frame) + model.paramOn(frame, Layout=(0.65, 1.0, 0.9)) frame.getAttText().SetTextFont(42) frame.getAttText().SetTextSize(0.001) if range_m: @@ -157,9 +100,9 @@ def fit_mass_new( frame.SetAxisRange(0.0, frame.GetMaximum() + (frame.GetMaximum() * 0.3), "Y") try: - for pdf in ret_model.pdfList(): + for pdf in model.pdfList(): pdf_name = pdf.GetName() - ret_model.plotOn( + model.plotOn( frame, ROOT.RooFit.Components(pdf), ROOT.RooFit.Name(f"pdf_{pdf_name}"), @@ -168,7 +111,7 @@ def fit_mass_new( ROOT.RooFit.LineWidth(1), ) # model.SetName("bkg") - ret_model.plotOn(frame, ROOT.RooFit.Name("model")) + model.plotOn(frame, ROOT.RooFit.Name("model")) except: # pylint: disable=bare-except # noqa: E722 pass # for comp in fit_spec.get('components', {}): @@ -178,7 +121,7 @@ def fit_mass_new( # c.Modified() # c.Update() - if level == "data": #and USE_EXTMODEL and frame is not None: + if level == "data" and USE_EXTMODEL and frame is not None: residuals = frame.residHist("data", "pdf_bkg") residual_frame = m.frame() residual_frame.addPlotable(residuals, "P") @@ -196,7 +139,7 @@ def fit_mass_new( residual_frame.SetAxisRange(range_m[0], range_m[1], "X") residual_frame.SetYTitle("Residuals") - return (res, ws, frame, residual_frame, dh, ret_model) + return (res, ws, frame, residual_frame) def fit_mass(self, hist, fit_spec, plot=False): """Old fit method""" @@ -229,23 +172,10 @@ def fit_mass(self, hist, fit_spec, plot=False): return (res, ws, frame) -def estimate_signal(hist, fit_spec, n_bkg, bkg_integral): - bin_min = hist.FindBin(fit_spec["mass_mean"] - 3 * fit_spec["sigma_signal"]) - bin_max = hist.FindBin(fit_spec["mass_mean"] + 3 * fit_spec["sigma_signal"]) - msum = 0. - for ind in range(bin_min, bin_max + 1): - msum += hist.GetBinContent(ind) - bkg = calculate_background(n_bkg, bkg_integral) - return msum - bkg - - -def calculate_background(n_bkg, bkg_integral): - return n_bkg.getVal() * bkg_integral - def calc_signif(roows, res, pdfnames, param_names, mean_sgn, sigma_sgn): """Calculate significance, signal, background, signal/background ratio.""" - #if not USE_EXTMODEL: - # return (0., 0., 0., 0., 0., 0, 0, 0.) + if not USE_EXTMODEL: + return (0., 0., 0., 0., 0., 0, 0, 0.) f_sig = roows.pdf(pdfnames["pdf_sig"]) n_signal = res.floatParsFinal().find("n_signal").getVal() sigma_n_signal = res.floatParsFinal().find("n_signal").getError() @@ -271,8 +201,6 @@ def calc_signif(roows, res, pdfnames, param_names, mean_sgn, sigma_sgn): n_signal_signal = signal_integral.getVal() * n_signal n_bkg_signal = bkg_integral.getVal() * n_bkg - print(f"n signal signal: {n_signal_signal} n bkg signal: {n_bkg_signal}") - if n_signal_signal + n_bkg_signal == 0.: significance = 0.0 else: @@ -346,11 +274,11 @@ def add_text_info_fit(text_info, frame, roows, param_names): refl_frac = roows.var(param_names["fraction_refl"]) text_info.AddText(f"#chi^{{2}}/ndf = {chi2:.2f}") text_info.AddText(f"#mu = {mean_sgn.getVal():.3f} #pm {mean_sgn.getError():.3f}") - text_info.AddText(f"#sigma = {sigma_sgn.getVal():.4f} #pm {sigma_sgn.getError():.4f}") + text_info.AddText(f"#sigma = {sigma_sgn.getVal():.3f} #pm {sigma_sgn.getError():.3f}") if sigmawide_sgn: - text_info.AddText(f"#sigma wide = {sigmawide_sgn.getVal():.4f} #pm {sigmawide_sgn.getError():.4f}") + text_info.AddText(f"#sigma wide = {sigmawide_sgn.getVal():.3f} #pm {sigmawide_sgn.getError():.3f}") if refl_frac: - text_info.AddText(f"refl.frac. = {refl_frac.getVal():.4f} #pm {refl_frac.getError():.4f}") + text_info.AddText(f"refl.frac. = {refl_frac.getVal():.3f} #pm {refl_frac.getError():.3f}") if a0 := roows.var("a0"): text_info.AddText(f"a0 = {a0.getVal():.3f} #pm {a0.getError():.3f}") if a1 := roows.var("a1"): From 5943d7c80a72accbd22cbe96d86054152d08c4c7 Mon Sep 17 00:00:00 2001 From: saganatt <8majak8@gmail.com> Date: Fri, 3 Oct 2025 12:39:21 +0200 Subject: [PATCH 4/7] Return the model to save TF functions --- machine_learning_hep/analysis/analyzerdhadrons.py | 11 ++++++----- .../analysis/analyzerdhadrons_mult.py | 11 ++++++----- machine_learning_hep/fitting/roofitter.py | 2 +- 3 files changed, 13 insertions(+), 11 deletions(-) diff --git a/machine_learning_hep/analysis/analyzerdhadrons.py b/machine_learning_hep/analysis/analyzerdhadrons.py index 5bcac1bbb1..1e3b386dc6 100644 --- a/machine_learning_hep/analysis/analyzerdhadrons.py +++ b/machine_learning_hep/analysis/analyzerdhadrons.py @@ -171,12 +171,13 @@ def _save_hist(self, hist, filename, option=""): # region fitting def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=None, filename=None): if fitcfg is None: - return None, None, None + return None, None, None, None, None try: - res, ws, frame, residual_frame = self.fitter.fit_mass_new(hist, pdfnames, fitcfg, level, roows=roows, plot=True) + res, ws, frame, residual_frame, data_hist, model = self.fitter.fit_mass_new(hist, pdfnames, fitcfg, level, + roows=roows, plot=True) except ValueError: self.logger.error(f"Could not do fitting on {level} for pt {self.bins_candpt[ipt]} - {self.bins_candpt[ipt+1]}") - return None, None, None + return None, None, None, None, None frame.SetTitle(f"inv. mass for p_{{T}} {self.bins_candpt[ipt]} - {self.bins_candpt[ipt + 1]} GeV/c") c = TCanvas() @@ -216,7 +217,7 @@ def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=No chi = frame.chiSquare() - return res, ws, chi + return res, ws, chi, data_hist, model def _fit_mass(self, hist, filename=None): if hist.GetEntries() == 0: @@ -348,7 +349,7 @@ def fit(self): roows.var(fixpar).setConstant(True) if h_invmass.GetEntries() == 0: continue - roo_res, roo_ws, chi = self._roofit_mass( + roo_res, roo_ws, chi, dh, model = self._roofit_mass( level, h_invmass, ipt, diff --git a/machine_learning_hep/analysis/analyzerdhadrons_mult.py b/machine_learning_hep/analysis/analyzerdhadrons_mult.py index 1ead07ae5f..89aa2a9a2f 100644 --- a/machine_learning_hep/analysis/analyzerdhadrons_mult.py +++ b/machine_learning_hep/analysis/analyzerdhadrons_mult.py @@ -199,12 +199,13 @@ def _save_hist(self, hist, filename, option=""): # region fitting def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=None, filename=None): if fitcfg is None: - return None, None, None + return None, None, None, None, None try: - res, ws, frame, residual_frame = self.fitter.fit_mass_new(hist, pdfnames, fitcfg, level, roows=roows, plot=True) + res, ws, frame, residual_frame, data_hist, model = self.fitter.fit_mass_new(hist, pdfnames, fitcfg, level, + roows=roows, plot=True) except ValueError: self.logger.error(f"Could not do fitting on {level} for pt {self.bins_candpt[ipt]} - {self.bins_candpt[ipt+1]}") - return None, None, None + return None, None, None, None, None frame.SetTitle(f"inv. mass for p_{{T}} {self.bins_candpt[ipt]} - {self.bins_candpt[ipt + 1]} GeV/c") c = TCanvas() @@ -244,7 +245,7 @@ def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=No chi = frame.chiSquare() - return res, ws, chi + return res, ws, chi, data_hist, model def _fit_mass(self, hist, filename=None): if hist.GetEntries() == 0: @@ -398,7 +399,7 @@ def fit(self): # Create the directory if it doesn't exist directory_path.mkdir(parents=True, exist_ok=True) - roo_res, roo_ws, chi = self._roofit_mass( + roo_res, roo_ws, chi, dh, model = self._roofit_mass( level, h_invmass, ipt, diff --git a/machine_learning_hep/fitting/roofitter.py b/machine_learning_hep/fitting/roofitter.py index 2b231c8283..28cf8272e1 100644 --- a/machine_learning_hep/fitting/roofitter.py +++ b/machine_learning_hep/fitting/roofitter.py @@ -139,7 +139,7 @@ def fit_mass_new( residual_frame.SetAxisRange(range_m[0], range_m[1], "X") residual_frame.SetYTitle("Residuals") - return (res, ws, frame, residual_frame) + return (res, ws, frame, residual_frame, dh, model) def fit_mass(self, hist, fit_spec, plot=False): """Old fit method""" From 565a49f0cd4ae46aa3a4fb4c5083bd20744f9032 Mon Sep 17 00:00:00 2001 From: saganatt <8majak8@gmail.com> Date: Fri, 3 Oct 2025 12:52:14 +0200 Subject: [PATCH 5/7] Add the possibility of setting fixed sigma --- machine_learning_hep/analysis/analyzerdhadrons.py | 11 +++++++++-- .../analysis/analyzerdhadrons_mult.py | 11 +++++++++-- machine_learning_hep/fitting/roofitter.py | 8 +++++++- 3 files changed, 25 insertions(+), 5 deletions(-) diff --git a/machine_learning_hep/analysis/analyzerdhadrons.py b/machine_learning_hep/analysis/analyzerdhadrons.py index 1e3b386dc6..cb7e0dd24e 100644 --- a/machine_learning_hep/analysis/analyzerdhadrons.py +++ b/machine_learning_hep/analysis/analyzerdhadrons.py @@ -100,6 +100,8 @@ def __init__(self, datap, case, typean, period): self.n_fileff = os.path.join(self.d_resultsallpmc, self.n_fileff) self.p_bin_width = datap["analysis"][self.typean]["bin_width"] self.p_rebin = datap["analysis"][self.typean]["n_rebin"] + self.p_fixed_sigma = datap["analysis"][self.typean]["fixed_sigma"] + self.p_fixed_sigma_val = datap["analysis"][self.typean]["fixed_sigma_val"] self.p_pdfnames = datap["analysis"][self.typean]["pdf_names"] self.p_param_names = datap["analysis"][self.typean]["param_names"] @@ -169,11 +171,14 @@ def _save_hist(self, hist, filename, option=""): self.rfigfile.WriteObject(hist, rfilename) # region fitting - def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=None, filename=None): + def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, fixed_sigma, fixed_sigma_val, # pylint: disable=too-many-arguments + roows=None, filename=None): if fitcfg is None: return None, None, None, None, None try: - res, ws, frame, residual_frame, data_hist, model = self.fitter.fit_mass_new(hist, pdfnames, fitcfg, level, + res, ws, frame, residual_frame, data_hist, model = self.fitter.fit_mass_new(hist, pdfnames, param_names, + fitcfg, level, + fixed_sigma, fixed_sigma_val, roows=roows, plot=True) except ValueError: self.logger.error(f"Could not do fitting on {level} for pt {self.bins_candpt[ipt]} - {self.bins_candpt[ipt+1]}") @@ -356,6 +361,8 @@ def fit(self): self.p_pdfnames, self.p_param_names, fitcfg, + self.p_fixed_sigma[ipt], + self.p_fixed_sigma_val[ipt], roows, f"roofit/h_mass_fitted_pthf-{ptrange[0]}-{ptrange[1]}_{level}.png", ) diff --git a/machine_learning_hep/analysis/analyzerdhadrons_mult.py b/machine_learning_hep/analysis/analyzerdhadrons_mult.py index 89aa2a9a2f..b0ea3a758c 100644 --- a/machine_learning_hep/analysis/analyzerdhadrons_mult.py +++ b/machine_learning_hep/analysis/analyzerdhadrons_mult.py @@ -105,6 +105,8 @@ def __init__(self, datap, case, typean, period): self.p_bin_width = datap["analysis"][self.typean]["bin_width"] self.p_rebin = datap["analysis"][self.typean]["n_rebin"] + self.p_fixed_sigma = datap["analysis"][self.typean]["fixed_sigma"] + self.p_fixed_sigma_val = datap["analysis"][self.typean]["fixed_sigma_val"] self.p_pdfnames = datap["analysis"][self.typean]["pdf_names"] self.p_param_names = datap["analysis"][self.typean]["param_names"] @@ -197,11 +199,14 @@ def _save_hist(self, hist, filename, option=""): self.rfigfile.WriteObject(hist, rfilename) # region fitting - def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=None, filename=None): + def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, fixed_sigma, fixed_sigma_val, # pylint: disable=too-many-arguments + roows=None, filename=None): if fitcfg is None: return None, None, None, None, None try: - res, ws, frame, residual_frame, data_hist, model = self.fitter.fit_mass_new(hist, pdfnames, fitcfg, level, + res, ws, frame, residual_frame, data_hist, model = self.fitter.fit_mass_new(hist, pdfnames, param_names, + fitcfg, level, + fixed_sigma, fixed_sigma_val, roows=roows, plot=True) except ValueError: self.logger.error(f"Could not do fitting on {level} for pt {self.bins_candpt[ipt]} - {self.bins_candpt[ipt+1]}") @@ -406,6 +411,8 @@ def fit(self): self.p_pdfnames, self.p_param_names, fitcfg, + self.p_fixed_sigma[ipt], + self.p_fixed_sigma_val[ipt], roows, f"roofit/mult_{multrange[0]}-{multrange[1]}/" f"h_mass_fitted_pthf-{ptrange[0]}-{ptrange[1]}" diff --git a/machine_learning_hep/fitting/roofitter.py b/machine_learning_hep/fitting/roofitter.py index 28cf8272e1..09f3720759 100644 --- a/machine_learning_hep/fitting/roofitter.py +++ b/machine_learning_hep/fitting/roofitter.py @@ -34,7 +34,8 @@ def __init__(self): # pylint: disable=too-many-branches def fit_mass_new( - self, hist, pdfnames: dict, fit_spec: dict, level: str, roows: ROOT.RooWorkspace = None, plot: bool = False + self, hist, pdfnames: dict, param_names: dict, fit_spec: dict, level: str, + fixed_sigma: bool = False, fixed_sigma_val: float = 0., roows: ROOT.RooWorkspace = None, plot: bool = False ): """New fit method""" if hist.GetEntries() == 0: @@ -54,6 +55,11 @@ def fit_mass_new( raise ValueError("model not set") m = ws.var(var_m) + if level == "mc": + sigma_sgn = ws.var(param_names["gauss_sigma"]) + if fixed_sigma: + sigma_sgn.setVal(fixed_sigma_val) + sigma_sgn.setConstant(True) if level == "data" and USE_EXTMODEL: signal_pdf = ws.pdf(pdfnames["pdf_sig"]) From fe5884d5f46b8958ba0ddc6fea8c9daa9f3d326e Mon Sep 17 00:00:00 2001 From: saganatt <8majak8@gmail.com> Date: Fri, 3 Oct 2025 16:31:14 +0200 Subject: [PATCH 6/7] Manual fix of related pylint errors --- machine_learning_hep/analysis/analyzerdhadrons.py | 14 ++++++++++---- .../analysis/analyzerdhadrons_mult.py | 8 ++++++-- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/machine_learning_hep/analysis/analyzerdhadrons.py b/machine_learning_hep/analysis/analyzerdhadrons.py index cb7e0dd24e..ec010dae7b 100644 --- a/machine_learning_hep/analysis/analyzerdhadrons.py +++ b/machine_learning_hep/analysis/analyzerdhadrons.py @@ -29,6 +29,8 @@ TF1, TH1, TH1F, + RooArgSet, + RooConstVar, TCanvas, TFile, TLegend, @@ -181,7 +183,8 @@ def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, fixed_si fixed_sigma, fixed_sigma_val, roows=roows, plot=True) except ValueError: - self.logger.error(f"Could not do fitting on {level} for pt {self.bins_candpt[ipt]} - {self.bins_candpt[ipt+1]}") + self.logger.error("Could not do fitting on %d for pt %.1f - %.1f", + level, self.bins_candpt[ipt], self.bins_candpt[ipt+1]) return None, None, None, None, None frame.SetTitle(f"inv. mass for p_{{T}} {self.bins_candpt[ipt]} - {self.bins_candpt[ipt + 1]} GeV/c") c = TCanvas() @@ -389,12 +392,14 @@ def fit(self): if roo_res: (sig, sig_err, _, _, signif, signif_err, s_over_b, s_over_b_err - ) = calc_signif(roo_ws, roo_res, self.p_pdfnames, self.p_param_names, mean_sgn, sigma_sgn) + ) = calc_signif(roo_ws, roo_res, self.p_pdfnames, self.p_param_names, + mean_sgn, sigma_sgn) fileout.cd() one = RooConstVar("one", "constant 1.0", 1.0) bkg_pdf = roo_ws.pdf(self.p_pdfnames["pdf_bkg"]) sig_pdf = roo_ws.pdf(self.p_pdfnames["pdf_sig"]) - for pdf, outlabel in zip((bkg_pdf, sig_pdf, model), ("bkg", "sgn", "total")): + for pdf, outlabel in zip((bkg_pdf, sig_pdf, model), ("bkg", "sgn", "total"), + strict=False): if not pdf: continue obs = pdf.getObservables(dh) @@ -451,7 +456,8 @@ def efficiency(self): print(self.n_fileff) lfileeff = TFile.Open(self.n_fileff) lfileeff.ls() - fileouteff = TFile.Open(f"{self.d_resultsallpmc}/{self.efficiency_filename}{self.case}{self.typean}.root", "recreate") + fileouteff = TFile.Open(f"{self.d_resultsallpmc}/{self.efficiency_filename}{self.case}{self.typean}.root", + "recreate") cEff = TCanvas("cEff", "The Fit Canvas") cEff.SetCanvasSize(1900, 1500) cEff.SetWindowSize(500, 500) diff --git a/machine_learning_hep/analysis/analyzerdhadrons_mult.py b/machine_learning_hep/analysis/analyzerdhadrons_mult.py index b0ea3a758c..30ab47597f 100644 --- a/machine_learning_hep/analysis/analyzerdhadrons_mult.py +++ b/machine_learning_hep/analysis/analyzerdhadrons_mult.py @@ -29,6 +29,8 @@ TF1, TH1, TH1F, + RooArgSet, + RooConstVar, TCanvas, TFile, TLegend, @@ -209,7 +211,8 @@ def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, fixed_si fixed_sigma, fixed_sigma_val, roows=roows, plot=True) except ValueError: - self.logger.error(f"Could not do fitting on {level} for pt {self.bins_candpt[ipt]} - {self.bins_candpt[ipt+1]}") + self.logger.error("Could not do fitting on %d for pt %.1f - %.1f", + level, self.bins_candpt[ipt], self.bins_candpt[ipt+1]) return None, None, None, None, None frame.SetTitle(f"inv. mass for p_{{T}} {self.bins_candpt[ipt]} - {self.bins_candpt[ipt + 1]} GeV/c") c = TCanvas() @@ -449,7 +452,8 @@ def fit(self): one = RooConstVar("one", "constant 1.0", 1.0) bkg_pdf = roo_ws.pdf(self.p_pdfnames["pdf_bkg"]) sig_pdf = roo_ws.pdf(self.p_pdfnames["pdf_sig"]) - for pdf, outlabel in zip((bkg_pdf, sig_pdf, model), ("bkg", "sgn", "total")): + for pdf, outlabel in zip((bkg_pdf, sig_pdf, model), ("bkg", "sgn", "total"), + strict=False): if not pdf: continue obs = pdf.getObservables(dh) From fa23bc549638d029e03d79e65a33181ca096e898 Mon Sep 17 00:00:00 2001 From: saganatt <8majak8@gmail.com> Date: Fri, 3 Oct 2025 17:06:57 +0200 Subject: [PATCH 7/7] Ruff autofix --- .../analysis/analyzerdhadrons_mult.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/machine_learning_hep/analysis/analyzerdhadrons_mult.py b/machine_learning_hep/analysis/analyzerdhadrons_mult.py index 30ab47597f..fcdf26c466 100644 --- a/machine_learning_hep/analysis/analyzerdhadrons_mult.py +++ b/machine_learning_hep/analysis/analyzerdhadrons_mult.py @@ -537,7 +537,7 @@ def efficiency(self): legsl.AddEntry(h_sel_pr_sl, legeffstring, "LEP") h_sel_pr_sl.GetXaxis().SetTitle("#it{p}_{T} (GeV/#it{c})") - h_sel_pr_sl.GetYaxis().SetTitle("Signal loss (prompt) %s" % (self.p_latexnhadron)) + h_sel_pr_sl.GetYaxis().SetTitle(f"Signal loss (prompt) {self.p_latexnhadron}") h_sel_pr_sl.SetMinimum(0.7) h_sel_pr_sl.SetMaximum(1.0) @@ -558,7 +558,7 @@ def efficiency(self): h_sel_pr.Write() legeff.AddEntry(h_sel_pr, legeffstring, "LEP") h_sel_pr.GetXaxis().SetTitle("#it{p}_{T} (GeV/#it{c})") - h_sel_pr.GetYaxis().SetTitle("Acc x efficiency (prompt) %s" % (self.p_latexnhadron)) + h_sel_pr.GetYaxis().SetTitle(f"Acc x efficiency (prompt) {self.p_latexnhadron}") h_sel_pr.SetMinimum(0.0004) h_sel_pr.SetMaximum(0.4) @@ -613,7 +613,7 @@ def efficiency(self): legslFD.AddEntry(h_sel_fd_sl, legeffstring, "LEP") h_sel_fd_sl.GetXaxis().SetTitle("#it{p}_{T} (GeV/#it{c})") - h_sel_fd_sl.GetYaxis().SetTitle("Signal loss (feeddown) %s" % (self.p_latexnhadron)) + h_sel_fd_sl.GetYaxis().SetTitle(f"Signal loss (feeddown) {self.p_latexnhadron}") h_sel_fd_sl.SetMinimum(0.7) h_sel_fd_sl.SetMaximum(1.0) @@ -634,7 +634,7 @@ def efficiency(self): h_sel_fd.Write() legeffFD.AddEntry(h_sel_fd, legeffFDstring, "LEP") h_sel_fd.GetXaxis().SetTitle("#it{p}_{T} (GeV/#it{c})") - h_sel_fd.GetYaxis().SetTitle("Acc x efficiency feed-down %s" % (self.p_latexnhadron)) + h_sel_fd.GetYaxis().SetTitle(f"Acc x efficiency feed-down {self.p_latexnhadron}") h_sel_fd.SetMinimum(0.0004) h_sel_fd.SetMaximum(0.4) @@ -688,7 +688,7 @@ def plotter(self): norm = 2 * self.p_br * self.p_nevents / (self.p_sigmamb * 1e12) hcross.Scale(1.0 / norm) fileoutcross.cd() - hcross.GetXaxis().SetTitle("#it{p}_{T} %s (GeV/#it{c})" % self.p_latexnhadron) + hcross.GetXaxis().SetTitle(f"#it{{p}}_{{T}} {self.p_latexnhadron} (GeV/#it{{c}})") hcross.GetYaxis().SetTitle(f"d#sigma/d#it{{p}}_{{T}} ({self.p_latexnhadron}) {self.typean}") hcross.SetName("hcross%d" % imult) hcross.GetYaxis().SetRangeUser(1e1, 1e10) @@ -725,7 +725,7 @@ def plotter(self): print("pt", ipt) for imult in range(self.p_nbin2): hcrossvsvar2[ipt].SetLineColor(ipt + 1) - hcrossvsvar2[ipt].GetXaxis().SetTitle("%s" % self.p_latexbin2var) + hcrossvsvar2[ipt].GetXaxis().SetTitle(f"{self.p_latexbin2var}") hcrossvsvar2[ipt].GetYaxis().SetTitle(self.p_latexnhadron) hcrossvsvar2[ipt].SetBinContent(imult + 1, listvalues[imult][ipt]) hcrossvsvar2[ipt].SetBinError(imult + 1, listvalueserr[imult][ipt]) @@ -898,7 +898,7 @@ def plotternormyields(self): hcross.Scale(1.0 / (self.p_sigmamb * 1e12)) hcross.SetLineColor(imult + 1) hcross.SetMarkerColor(imult + 1) - hcross.GetXaxis().SetTitle("#it{p}_{T} %s (GeV/#it{c})" % self.p_latexnhadron) + hcross.GetXaxis().SetTitle(f"#it{{p}}_{{T}} {self.p_latexnhadron} (GeV/#it{{c}})") hcross.GetYaxis().SetTitleOffset(1.3) hcross.GetYaxis().SetTitle(f"Corrected yield/events ({self.p_latexnhadron}) {self.typean}") hcross.GetYaxis().SetRangeUser(1e-10, 1)