diff --git a/machine_learning_hep/analysis/analyzerdhadrons.py b/machine_learning_hep/analysis/analyzerdhadrons.py index 24672b0303..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, @@ -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,19 @@ 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 - res, ws, frame, residual_frame = self.fitter.fit_mass_new(hist, pdfnames, fitcfg, level, roows, True) + 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=roows, plot=True) + except ValueError: + 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() @@ -180,7 +193,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 +206,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 +214,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 +223,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,16 +289,16 @@ 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)) 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) @@ -304,6 +319,7 @@ def fit(self): self.lpt_finbinmax[ipt], lpt_probcutfin[ipt], ) + rfile.cd() h_invmass = rfile.Get("hmass" + suffix) # Rebin h_invmass.Rebin(self.p_rebin[ipt]) @@ -341,19 +357,21 @@ 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", ) 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 +389,27 @@ 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"]) + for pdf, outlabel in zip((bkg_pdf, sig_pdf, model), ("bkg", "sgn", "total"), + strict=False): + 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 yieldshistos.SetBinContent(ipt + 1, sig) yieldshistos.SetBinError(ipt + 1, sig_err) @@ -385,13 +421,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() - fileout.Close() + chihistos.Write() def yield_syst(self): # Enable ROOT batch mode and reset in the end @@ -418,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) @@ -470,6 +509,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: diff --git a/machine_learning_hep/analysis/analyzerdhadrons_mult.py b/machine_learning_hep/analysis/analyzerdhadrons_mult.py index cd62e163eb..fcdf26c466 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, @@ -105,6 +107,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,10 +201,19 @@ 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 - res, ws, frame, residual_frame = self.fitter.fit_mass_new(hist, pdfnames, fitcfg, level, roows, True) + 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=roows, plot=True) + except ValueError: + 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() @@ -208,7 +221,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 +234,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 +251,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: @@ -301,8 +316,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) @@ -315,6 +329,7 @@ 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)): @@ -342,6 +357,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]) @@ -391,13 +407,15 @@ 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, 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/mult_{multrange[0]}-{multrange[1]}/" f"h_mass_fitted_pthf-{ptrange[0]}-{ptrange[1]}" @@ -407,7 +425,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 +443,27 @@ 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"]) + for pdf, outlabel in zip((bkg_pdf, sig_pdf, model), ("bkg", "sgn", "total"), + strict=False): + 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 yieldshistos.SetBinContent(ipt + 1, sig) yieldshistos.SetBinError(ipt + 1, sig_err) @@ -439,13 +475,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() - fileout.Close() + chihistos.Write() def get_efficiency(self, ibin1, ibin2): fileouteff = TFile.Open(f"{self.d_resultsallpmc}/efficiencies{self.case}{self.typean}.root", "read") @@ -499,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) @@ -520,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) @@ -575,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) @@ -596,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) @@ -650,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) @@ -687,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]) @@ -860,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) diff --git a/machine_learning_hep/fitting/roofitter.py b/machine_learning_hep/fitting/roofitter.py index 41ef0f660c..09f3720759 100644 --- a/machine_learning_hep/fitting/roofitter.py +++ b/machine_learning_hep/fitting/roofitter.py @@ -32,8 +32,10 @@ def __init__(self): ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.WARNING) ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.ERROR) + # 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: @@ -53,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"]) @@ -70,7 +77,7 @@ def fit_mass_new( 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: + if level == "data" and USE_EXTMODEL: for v in ws.allVars(): v.setConstant(True) res = extmodel.fitTo( @@ -78,7 +85,7 @@ def fit_mass_new( ) else: res = model.fitTo(dh, Save=True, PrintLevel=-1, Strategy=1, MaxCalls=5000) - if level == 'data' and USE_EXTMODEL: + 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) @@ -138,7 +145,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""" @@ -200,7 +207,10 @@ 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) + 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 +221,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,