**plot_fitted_gaussian(data, ax, bin_num, data_label=None, spec_limit=None, rss=None)**
* data (1D numpy array): data to be fitted and plotted
* ax (matplotlib axis): axis object for the plotted data
* bin_num (integer): bin number for the histogram
* data_label (optional, string): title for the plot
* spec_limit (float tuple): upper/lower limit locations
* rss (float tuple): rss locations

def get_sigma_bin(bins, mu, sigma, sigma_num):
    left_bin_ind = np.argmin(abs(bins - (mu - sigma_num * sigma)))
    right_bin_ind = np.argmin(abs(bins - (mu + sigma_num * sigma)))
    return left_bin_ind, right_bin_ind

def plot_fitted_gaussian(data, ax, bin_num, data_label=None, spec_limit=None, rss=None):
    usl = None
    lsl = None
    l_rss = None
    u_rss = None
    if spec_limit:
        lsl, usl = spec_limit
    if rss:
        l_rss, u_rss = rss
    sample_size = data.size
    trials, bins, _ = ax.hist(data, bins=bin_num, edgecolor="black", fill=False, linewidth=0.5)
    if data_label:
        ax.set_xlabel(data_label)   

    mu, sigma = stats.norm.fit(data)
    p = stats.norm.pdf(bins, mu, sigma)
    p_norm =  p/p.sum() * sample_size    

    # gaussian sigma area  
    left_sigma_1_ind, right_sigma_1_ind = get_sigma_bin(bins, mu, sigma, 1)
    left_sigma_2_ind, right_sigma_2_ind = get_sigma_bin(bins, mu, sigma, 2)
    print(f'Fitted Gaussian Info: mu:{mu}, sigma:{sigma}, 1 sigma sum: {p_norm[left_sigma_1_ind:right_sigma_1_ind + 1].sum()}, 2 sigma sum: {p_norm[left_sigma_2_ind:right_sigma_2_ind + 1].sum()}')
    ax.plot(bins, p_norm ,'o--b', linewidth=1, markersize=1.5, label='Fitted Gaussian')   
    # ax.fill_between(bins[left_sigma_1_ind:right_sigma_1_ind + 1], p_norm[left_sigma_1_ind:right_sigma_1_ind + 1], y2=0, facecolor='orange', alpha=0.45, label=r'1$\sigma$')
    
    ax.fill_between(bins[left_sigma_1_ind:right_sigma_1_ind + 1], p_norm[left_sigma_1_ind:right_sigma_1_ind + 1], y2=0, facecolor='orange', alpha=0.45, label=f'$1\/\sigma\hspace{{1}}({sci_tex(sigma)})$')
    ax.fill_between(bins[left_sigma_2_ind:left_sigma_1_ind + 1], p_norm[left_sigma_2_ind:left_sigma_1_ind + 1], y2=0, facecolor='green', alpha=0.45, label=f'$2\/\sigma\hspace{{1}}({sci_tex(2 * sigma)})$')
    ax.fill_between(bins[right_sigma_1_ind:right_sigma_2_ind + 1], p_norm[right_sigma_1_ind:right_sigma_2_ind + 1], y2=0, facecolor='green', alpha=0.45)

    # LSL & USL
    ylim = ax.get_ylim()
    if lsl:
        ax.plot([lsl, lsl], ax.get_ylim(), c='r', linestyle='--', linewidth=1, label='LSL')
    if usl:
        ax.plot([usl, usl], ax.get_ylim(), c='r', linestyle=':', linewidth=1, label='USL')
    
    # RSS
    if l_rss:
        ax.plot([l_rss, l_rss], ax.get_ylim(), c='black', linestyle='-.', linewidth=1, label='RSS')
    if u_rss:
        ax.plot([u_rss, u_rss], ax.get_ylim(), c='black', linestyle='-.', linewidth=1)
    
    ax.set_ylim(ylim)
    ax.legend()
    return

**Example**

from helpers import plot_fitted_gaussian
import numpy.random
import numpy as np
import matplotlib.pyplot as plt

data = numpy.random.randn(3, 10000)
fig, ax =  plt.subplots(3, 1)
for i, r in enumerate(data[:, :]):
    plot_fitted_gaussian(r, bin_num=24, ax=ax[i])

**Output**