Lets say that we have a two fits that we are trying to implement. I will assume that you have this version of double gaussian and exp class.
You really should use this one instead. It includes more sanity checks and allows you to pass iminuit limits as optional.
This means that you don't have to always put (-np.inf, np.inf) for each limit anymore.
class DoubleGaussian_plus_Exp:
def __init__(self, bins, nC, minuit_limits=None):
self.xMin, self.xMax, self.bin_width, self.x_vals, self.y_vals, self.y_errs = self.Fit_Setup(nC, bins)
self.minuit_limits = minuit_limits or {}
def Fit_Setup(self, nC, bins):
xMin = bins[0]
xMax = bins[-1]
#print('xMin, xMax = ',xMin,xMax)
bin_width = bins[1]-bins[0]
#print('bin_width = ',bin_width)
#x_vals = bins[0:len(bins)-1] + 0.5*bin_width
x_vals = 0.5 * (bins[:-1] + bins[1:])
#print("x_vals[0:10] = ",x_vals[0:10])
#print("x_vals shape = ",x_vals.shape)
y_vals = nC
#print("y_vals[0:10] = ",y_vals[0:10])
#print("y_vals shape = ",y_vals.shape)
y_errs = np.sqrt(nC)
return xMin, xMax, bin_width, x_vals, y_vals, y_errs
def gaussian(self, x, mu, sigma):
return self.bin_width * np.exp(-(x-mu)**2 / (2 * sigma**2)) / np.sqrt(2 * np.pi * sigma**2)
def expA(self, x, A, b):
integral = (A/b) * (1.0 - np.exp(-b*(self.xMax - self.xMin)))
norm = 1./integral
return self.bin_width * norm * A * np.exp(-b * (x - self.xMin))
def DoubleGaussian_plus_ExpA(self, x_vals, n_s, f, n_b, mu1, mu2, sigma1, sigma2, A, b):
n_s1 = n_s*f
n_s2 = n_s*(1-f)
return n_s1 * self.gaussian(x_vals, mu1, sigma1) + n_s2 * self.gaussian(x_vals, mu2, sigma2) + n_b * self.expA(x_vals, A, b)
def chi_squared(self, n_s, f, n_b, mu1, mu2, sigma1, sigma2, A, b):
mask = (0 != self.y_errs)
#print("mask shape = ",mask.shape)
prediction = self.DoubleGaussian_plus_ExpA(self.x_vals[mask],n_s, f, n_b, mu1, mu2, sigma1, sigma2, A, b)
ressq = (self.y_vals[mask] - prediction)**2 / np.square(self.y_errs[mask])
return ressq.sum()
def fit(self, init_pars):
m = Minuit(self.chi_squared, n_s=init_pars[0], f=init_pars[1], n_b=init_pars[2], mu1=init_pars[3], mu2=init_pars[4], sigma1=init_pars[5], sigma2=init_pars[6], A=init_pars[7], b=init_pars[8])
m.limits["n_s"] = self.minuit_limits.get("n_s", None)
m.limits["f"] = self.minuit_limits.get("f", None)
m.limits["n_b"] = self.minuit_limits.get("n_b", None)
m.limits["mu1"] = self.minuit_limits.get("mu1", None)
m.limits["mu2"] = self.minuit_limits.get("mu2", None)
m.limits["sigma1"] = self.minuit_limits.get("sigma1", None)
m.limits["sigma2"] = self.minuit_limits.get("sigma2", None)
m.limits["A"] = self.minuit_limits.get("A", None)
m.limits["b"] = self.minuit_limits.get("b", None)
m.migrad()
return m
def plot(self, m, bins, nC, title='Plot', xlabel='X', ylabel='Y', vlines=None, show_plot=True):
fig, ax = plt.subplots()
ax.set_title(title)
ax.hist(bins[:-1], bins=bins, weights=nC, label='Data')
mask = (0 != self.y_errs)
predictions = self.DoubleGaussian_plus_ExpA(self.x_vals[mask], *m.values)
ax.plot(self.x_vals[mask], predictions, 'r-', label='Fit')
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.axis([self.xMin, self.xMax, 0, 1.15*max(nC)])
ax.grid(True)
if vlines:
for vl in vlines:
ax.vlines(vl, 0., 0.6*max(nC), colors='yellow')
ax.legend()
if show_plot:
plt.show()
return fig, axNow lets try to fit both Lambda mass with both goodLambdaLL and goodLambdaDD cuts and then try to plot them together as subplot (just as an example)
The first step is that we are going to do the fit (and plot) like we used to do. We can do something like this
plt.ioff()
nC, bins, patches = plt.hist(L0_MM[good_LambdaLL], bins=100, range=(1100,1130))
plt.ion()
plt.close()
init_pars = [9120.817313121128, 0.929509388491892, 6548.478148791116, 1115.7823039513255, 1499.4428273473834, 1.6580107105881836, 4.957241526156779, -0.030793060590727006, -0.050084638862231065]
minuit_limits = {
"f": (0.0001, 0.9999),
}
fit_model_1 = DoubleGaussian_plus_Exp(bins, nC, minuit_limits)
fit_result_1 = fit_model_1.fit(init_pars)
fig_1, ax_1 = fit_model_1.plot(fit_result_1, bins, nC, title='L0_MM with goodLambdaLL ', xlabel='L0 Mass', ylabel='Entries per 0.5 MeV', vlines=[1114., 1118.], show_plot=False)
plt.close()
plt.ioff()
nC_2, bins_2, patches_2 = plt.hist(L0_MM[good_LambdaDD], bins=100, range=(1100,1130))
plt.ion()
plt.close()
init_pars_2 = [9120.817313121128, 0.929509388491892, 6548.478148791116, 1115.7823039513255, 1499.4428273473834, 1.6580107105881836, 4.957241526156779, -0.030793060590727006, -0.050084638862231065]
fit_model_2 = DoubleGaussian_plus_Exp(bins_2, nC_2, minuit_limits)
fit_result_2 = fit_model_2.fit(init_pars_2)
fig_2, ax_2 = fit_model_2.plot(fit_result_2, bins_2, nC_2, title='L0_MM with goodLambdaDD', xlabel='L0 Mass' , ylabel='Entries per 0.5 MeV', vlines=[1114., 1118.], show_plot=False)
plt.close()You might be seeing some new stuff here so let me explain them.
plt.ioff(): This line turns off the interactive mode in matplotlib.
In the interactive mode, plots are displayed as soon as they are created, while in the non-interactive mode,
they are displayed when plt.show() is called. And you can guess that plt.ion() turns the interactive mode back on.
There is a deep reason why do we have to use both plt.close, plt.ioff()` and plt.ion()` combination like this one.
I actually don't know but got it throught trial and error.
and for this part
fig_1, ax_1 = fit_model_1.plot(fit_result_1, bins, nC, title='L0_MM with goodLambdaLL ', xlabel='L0 Mass', ylabel='Entries per 0.5 MeV', vlines=[1114., 1118.], show_plot=False)Here, the plot method of the model is called. This plots the data, the fit results, and returns a figure object and an axes object. The plot is not shown because of show_plot=False.
we do the samething for `goodLambdaDD) also but we choose different variable names.
Now lets define a function that will take our two plots and combine them as a subplot.
import matplotlib.patches as mpatches
import matplotlib.lines as mlines
def subplot(from_ax, to_ax):
to_ax.set_xlim(from_ax.get_xlim())
to_ax.set_ylim(from_ax.get_ylim())
to_ax.set_title(from_ax.get_title())
to_ax.set_xlabel(from_ax.get_xlabel())
to_ax.set_ylabel(from_ax.get_ylabel())
for line in from_ax.get_lines():
new_line = mlines.Line2D(line.get_xdata(), line.get_ydata(), color=line.get_color())
to_ax.add_line(new_line)
for patch in from_ax.patches:
new_patch = mpatches.Rectangle((patch.get_x(), patch.get_y()), patch.get_width(), patch.get_height(), fill=patch.get_fill(), color=patch.get_facecolor())
to_ax.add_patch(new_patch)
handles, labels = from_ax.get_legend_handles_labels()
to_ax.legend(handles, labels)
to_ax.label_outer() Instead of trying to re-plot and the complications we will introduce. we can define a wrapper function that captures the current state of the plt module, uses that state to create a subplot, and then restores the original state. Or we can do better by creating a function that copies the lines, patches, and other attributes from one axes object to another. We create a new figure with two subplots using plt.subplots, and then use copy_axes to copy the data from the original plots to the subplots.
Here is a complete code for this function.
import matplotlib.patches as mpatches
import matplotlib.lines as mlines
def subplot(from_ax, to_ax):
to_ax.set_xlim(from_ax.get_xlim())
to_ax.set_ylim(from_ax.get_ylim())
to_ax.set_title(from_ax.get_title())
to_ax.set_xlabel(from_ax.get_xlabel())
to_ax.set_ylabel(from_ax.get_ylabel())
for line in from_ax.get_lines():
new_line = mlines.Line2D(line.get_xdata(), line.get_ydata(), color=line.get_color())
to_ax.add_line(new_line)
for patch in from_ax.patches:
new_patch = mpatches.Rectangle((patch.get_x(), patch.get_y()), patch.get_width(), patch.get_height(), fill=patch.get_fill(), color=patch.get_facecolor())
to_ax.add_patch(new_patch)
handles, labels = from_ax.get_legend_handles_labels()
to_ax.legend(handles, labels)
to_ax.label_outer() Lets go through each line
import matplotlib.patches as mpatches
import matplotlib.lines as mlinesThese lines import the patches and lines modules from matplotlib. Patches in matplotlib are used for 2D graphics like rectangles and circles, and lines for creating line segments.
def subplot(from_ax, to_ax):This line defines a function named subplot which takes two arguments: from_ax, the Axes object from which we want to copy elements, and to_ax, the Axes object to which we want to copy those elements.
to_ax.set_xlim(from_ax.get_xlim())
to_ax.set_ylim(from_ax.get_ylim())These lines copy the x-axis and y-axis limits from from_ax to to_ax.
to_ax.set_title(from_ax.get_title())
to_ax.set_xlabel(from_ax.get_xlabel())
to_ax.set_ylabel(from_ax.get_ylabel())These lines copy the title, x-label, and y-label from from_ax to to_ax.
for line in from_ax.get_lines()This loop goes through all line objects in from_ax, creates a new line with the same data and color, and adds it to to_ax.
for patch in from_ax.patchesThis loop goes through all patch objects in from_ax, creates a new patch (rectangle in this case) with the same properties, and adds it to to_ax.
handles, labels = from_ax.get_legend_handles_labels()
to_ax.legend(handles, labels)These lines extract the legend handles and labels from from_ax and use them to create a legend in to_ax.
to_ax.label_outer()This method is called to ensure labels for overlapping axes in the figure are displayed on the outermost subplots only.
So as we can deduce here. the purpose of this function is to create a duplicate plot by copying all relevant elements
(e.g., lines, patches, labels) from one Axes object (from_ax) to another (to_ax). Which provides a good way to solve our problem and make fit comparisons easier to visualize.