Improve stability of fits (alisw#1019)

Author: qgp

machine_learning_hep/analysis/analyzer_jets.py

@@ -126,7 +126,7 @@ def __init__(self, datap, case, typean, period):
                     ("#it{#mu}", "#it{#sigma}", "significance", "#it{#chi}^{2}"),
                 )
             }
-            for level in ("data", "mc", "mcsig", "mcrefl")
+            for level in self.fit_levels
         }
         self.n_events = {}
         self.n_colls_read = {}
@@ -395,7 +395,7 @@ def _roofit_mass(self, level, hist, ipt, pdfnames, param_names, fitcfg, roows=No
         chi2 = frame.chiSquare()
         if "ptjet" not in filename:
             self.h_fit_results[level]["chi2"].SetBinContent(ipt + 1, chi2)
-        if chi2 > 5.0:
+        if chi2 > 5.0 and level != "predata":
             self.logger.error(
                 "Roofit fit is too bad: %s, ipt: %d, pthf: %g-%g, Chi2 = %g",
                 level,
@@ -590,7 +590,7 @@ def fit(self):
                                 level,
                                 ipt,
                                 self.bins_candpt[ipt],
-                                self.bins_candpt[ipt + 1]
+                                self.bins_candpt[ipt + 1],
                             )
                         # if level == 'mc':
                         #     roo_ws.Print()

machine_learning_hep/data/data_run3/database_ml_parameters_LcJet_pp.yml

@@ -410,206 +410,35 @@ LcJet_pp:
         double_gauss_sigma: "wide"
         fraction_refl: "frac_refl"
 
+      fit_levels: ["mc", "predata", "data"]
       mass_roofit:
         - level: mc
-          ptrange: [0., 2.]
           range: [2.18, 2.39]
           components:
-            sig:
-              fn: "Gaussian::peak(m[1., 10], mean[2.286, 2.283,2.289], sigma_g1[.007,.007,.012])"
-            wide:
-              fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))'
             model:
-              fn: "SUM::sig(f_peak[0.,1.]*peak, wide)"
-        - level: mc
-          ptrange: [2., 3.]
-          range: [2.18, 2.39]
-          components:
-            sig:
-              fn: "Gaussian::peak(m[1., 10], mean[2.286, 2.283,2.289], sigma_g1[.007,.007,.012])"
-            wide:
-              fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))'
-            model:
-              fn: "SUM::sig(f_peak[0.,1.]*peak, wide)"
-        - level: mc
-          ptrange: [3., 4.]
-          range: [2.18, 2.39]
-          components:
-            sig:
-              fn: "Gaussian::peak(m[1., 10], mean[2.286, 2.283,2.289], sigma_g1[.007,.007,.012])"
-            wide:
-              fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))'
-            model:
-              fn: "SUM::sig(f_peak[0.,1.]*peak, wide)"
-        - level: mc
-          ptrange: [4., 5.]
-          range: [2.18, 2.39]
-          components:
-            sig:
-              fn: "Gaussian::peak(m[1., 10], mean[2.286, 2.283,2.289], sigma_g1[.007,.007,.012])"
-            wide:
-              fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))'
-            model:
-              fn: "SUM::sig(f_peak[0.,1.]*peak, wide)"
-        - level: mc
-          ptrange: [5., 6.]
-          range: [2.18, 2.39]
-          components:
-            sig:
-              fn: "Gaussian::peak(m[1., 10], mean[2.286, 2.283,2.289], sigma_g1[.012,.01,.016])"
-            wide:
-              fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))'
-            model:
-              fn: "SUM::sig(f_peak[0.,1.]*peak, wide)"
-        - level: mc
-          ptrange: [6., 7.]
-          range: [2.17, 2.40]
-          components:
-            sig:
-              fn: "Gaussian::peak(m[1., 10], mean[2.286, 2.283,2.289], sigma_g1[.012,.01,.016])"
-            wide:
-              fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))'
-            model:
-              fn: "SUM::sig(f_peak[0.,1.]*peak, wide)"
-        - level: mc
-          ptrange: [7., 8.]
-          range: [2.10, 2.45]
-          components:
-            sig:
-              fn: "Gaussian::peak(m[1., 10], mean[2.2865, 2.283,2.289], sigma_g1[.013,.013,.020])"
-            wide:
-              fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))'
-            model:
-              fn: "SUM::sig(f_peak[0.,1.]*peak, wide)"
-        - level: mc
-          ptrange: [8., 10.]
-          range: [2.10, 2.45]
-          components:
-            sig:
-              fn: "Gaussian::peak(m[1., 10], mean[2.2865, 2.283,2.289], sigma_g1[.013,.013,.020])"
-            wide:
-              fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))'
-            model:
-              fn: "SUM::sig(f_peak[0.,1.]*peak, wide)"
-        - level: mc
-          ptrange: [10., 12.]
-          range: [2.10, 2.45]
-          components:
-            sig:
-              fn: "Gaussian::peak(m[1., 10], mean[2.2865, 2.283,2.289], sigma_g1[.013,.013,.020])"
-            wide:
-              fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))'
-            model:
-              fn: "SUM::sig(f_peak[0.,1.]*peak, wide)"
-        - level: mc
-          ptrange: [12., 16.]
-          range: [2.10, 2.45]
-          components:
-            sig:
-              fn: "Gaussian::peak(m[1., 10], mean[2.2865, 2.283,2.289], sigma_g1[.020,.020,.029])"
-            wide:
-              fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))'
-            model:
-              fn: "SUM::sig(f_peak[0.,1.]*peak, wide)"
-        - level: mc
-          ptrange: [16., 30.]
-          range: [2.10, 2.47]
-          components:
-            sig:
-              fn: "Gaussian::peak(m[1., 10], mean[2.2865, 2.283,2.289], sigma_g1[.020,.020,.029])"
-            wide:
-              fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))'
-            model:
-              fn: "SUM::sig(f_peak[0.,1.]*peak, wide)"
-        - ptrange: [0., 2.]
-          range: [2.21, 2.36] #2.216, 2.36
-          fix_params: ['n', 'f_peak']
-          components:
-            #sig:
-            #fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.005,.015])'
-            bkg:
-              fn: 'Polynomial::bkg(m, {a1[-0.15, -3, 3], a2[0.1, -3., 3]})'
-            model:
-              fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)'
-        - ptrange: [2., 3.]
-          range: [2.20, 2.37]
-          fix_params: ['n', 'f_peak']
-          components:
-            #sig:
-            #   fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.005,.015])'
-            bkg:
-              fn: 'Polynomial::bkg(m, {a1[-0.15, -3, 3], a2[0.1, -3., 3]})'
-            model:
-              fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)'
-        - ptrange: [3., 4.]
-          range: [2.18, 2.38]
-          fix_params: ['n', 'f_peak']
-          components:
-            #sig:
-            #fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.005,.015])'
-            bkg:
-              fn: 'Polynomial::bkg(m, {a1[-0.15, -3, 3], a2[0.1, -3., 3]})'
-            model:
-              fn: 'SUM::sum(f_sig[0., 1.]*sig, bkg)'
-        - ptrange: [4., 5.]
-          range: [2.19, 2.38]
-          fix_params: ['n', 'f_peak']
-          components:
-            # sig:
-            #   fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.005,.015])'
-            bkg:
-              fn: 'Polynomial::bkg(m, {a1[-0.15, -3, 3], a2[0.1, -3., 3]})'
-            model:
-              fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)'
-        - ptrange: [5., 6.]
-          range: [2.18, 2.39]
-          fix_params: ['n', 'f_peak']
-          components:
-            # sig:
-            #   fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.005,.015])'
-            bkg:
-              fn: 'Polynomial::bkg(m, {a1[-0.2 , -3, 3], a2[0., -3, 3.]})'
-            model:
-              fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)'
-        - ptrange: [6., 7.]
-          range: [2.16, 2.41]
-          fix_params: ['n', 'f_peak']
-          components:
-            # sig:
-            #   fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.03])'
-            bkg:
-              fn: 'Polynomial::bkg(m, {a1[-0.2 , -3, 3], a2[0., -3, 3.]})'
-            model:
-              fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)'
-        - ptrange: [7., 8.]
-          range: [2.16, 2.41]
-          fix_params: ['n', 'f_peak']
-          components:
-            # sig:
-            #   fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.03])'
-            bkg:
-              fn: 'Polynomial::bkg(m, {a1[-0.2 , -3, 3], a2[0., -3, 3.]})'
-            model:
-              fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)'
-        - ptrange: [8., 10.]
-          range: [2.1, 2.46]
-          fix_params: ['n', 'f_peak']
+              fn: "RooCrystalBall::sig(m[1., 4.], mean[2.286, 2.27, 2.30], sigma_g1[0.006, 0.035], alpha1[0., 2.], n1[0., 100.], alpha1, n1)"
+        - level: predata
+          range: [2.16, 2.39]
+          fix_params: ["alpha1", "n1"]
           components:
-            # sig:
-            #   fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.03])'
             bkg:
-              fn: 'Polynomial::bkg(m, {a1[-0.2 , -3, 3], a2[-0.2, -3, 3.]})'
+              fn: "Polynomial::prebkg(m, {a0[1.], a1[0., -.5, 0.]}, lowestOrder = 0)"
             model:
-              fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)'
-        - range: [2.1, 2.47]
-          fix_params: ['n', 'f_peak']
+              fn: "SUM::presum(f_sig[0.,1.]*sig, prebkg)"
+        - level: data
+          range: [2.16, 2.39]
+          fix_params: ["alpha1", "n1", "a1"]
           components:
-            # sig:
-            #   fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.03])'
             bkg:
-              fn: 'Polynomial::bkg(m, {a1[-0.2 , -3, 3], a2[-0.2, -3, 3.]})'
+              # fn: "Polynomial::bkg(m, {a1[0., -.5, .5], a2[0., -5., 5.]})"
+              # fn: "Polynomial::bkg(m, {expr('a1 - 2.*a2*mean', a1[0., -.5, .5], a2[0., -.1, .1], mean), a2})"
+              # fn: "Polynomial::bkg(m, {expr('a1 - 2.*a2*offset', a1[0., -.5, -.2], a2[0., 0., 5.], offset[3., 0., 100.]), a2})"
+              # fn: "Polynomial::bkg(m, {expr('a1 - 2.*a2*offset', a1, a2[10., 0., 5000.], offset[0., -100., 100.]), a2})"
+              fn: "Polynomial::bkg(m, {expr('a0 + a2*offset*offset', offset[0., 3., 100.], a0, a2[0., -100., 100.]), expr('a1 - 2.*a2*offset', a1, a2, offset), a2}, lowestOrder = 0)"
+              # fn: "Polynomial::bkg(m, {expr('c + a2*offset*offset', c[0., 0., 1.], offset[0., -1000., 1000.], a2[0., -5000., 5000.]), expr('- 2.*a2*offset', a2, offset), a2}, lowestOrder = 0)"
+              # fn: "Polynomial::bkg(expr('m-mean', m, mean), {a1[-0.15, -10000.0, 10000.], a2[0.1, -10000., 10000.]})"
             model:
-              fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)'
+              fn: "SUM::sum(f_sig[0.,1.]*sig, bkg)"
 
       #sidesub_per_ptjet: true
       sidesub:

machine_learning_hep/fitting/roofitter.py

@@ -17,6 +17,7 @@
 import ROOT
 from ROOT import RooAddPdf, RooArgList, RooArgSet, RooFit, RooRealVar, TPaveText
 
+USE_EXTMODEL = True
 
 # pylint: disable=too-few-public-methods, too-many-statements
 # (temporary until we add more functionality)
@@ -42,7 +43,7 @@ def fit_mass_new(self, hist, pdfnames, fit_spec, level, roows=None, plot=False):
                 model = fn
         m = ws.var(var_m)
 
-        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")
@@ -60,13 +61,13 @@ def fit_mass_new(self, hist, pdfnames, fit_spec, level, roows=None, plot=False):
             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':
+            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=1, MaxCalls=5000)
-            if level == 'data':
+            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)
@@ -108,7 +109,7 @@ def fit_mass_new(self, hist, pdfnames, fit_spec, level, roows=None, plot=False):
             # c.Modified()
             # c.Update()
 
-        if level == "data":
+        if level == "data" and USE_EXTMODEL:
             residuals = frame.residHist("data", "pdf_bkg")
             residual_frame = m.frame()
             residual_frame.addPlotable(residuals, "P")
@@ -153,6 +154,8 @@ def fit_mass(self, hist, fit_spec, plot=False):
 
 
 def calc_signif(roows, res, pdfnames, param_names, mean_sgn, sigma_sgn):
+    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()
@@ -239,6 +242,14 @@ def add_text_info_fit(text_info, frame, roows, param_names):
         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():.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"):
+        text_info.AddText(f"a1 = {a1.getVal():.3f} #pm {a1.getError():.3f}")
+    if a2 := roows.var("a2"):
+        text_info.AddText(f"a2 = {a2.getVal():.3f} #pm {a2.getError():.3f}")
+    if offset := roows.var("offset"):
+        text_info.AddText(f"offset = {offset.getVal():.3f} #pm {offset.getError():.3f}")
 
 
 def add_text_info_perf(text_info, sig, sig_err, bkg, bkg_err, s_over_b, s_over_b_err, signif, signif_err):