A quick look to run 19¶

PSAU

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks, peak_widths
from lmfit.models import GaussianModel, ConstantModel
from scipy.optimize import curve_fit

from scipy.special import factorial

%matplotlib widget
In [2]:
def myLine(x, m, q):
    return m * x + q
In [3]:
data = np.loadtxt("run19-spettro.csv")
data.shape
Out[3]:
(8191,)
In [4]:
xVect = np.arange(0, 8191, 1)
xVect
Out[4]:
array([   0,    1,    2, ..., 8188, 8189, 8190], shape=(8191,))
In [5]:
minX = 2900
maxX = 3800

cond = (xVect > minX) & (xVect < maxX)
xSlice = xVect[cond]
ySlice = data[cond]

cond2 = (ySlice > 0)
xFit = xSlice[cond2]
yFit = ySlice[cond2]
In [6]:
# Peaks
peaks, _ = find_peaks(ySlice, height = 50, distance = 10, width = 5)

# FWHM
results_half = peak_widths(ySlice, peaks, rel_height=0.5)[0]
In [7]:
# Plot
fig, ax = plt.subplots()
fig.set_size_inches(12,5)

ax.plot(xSlice, ySlice, ds = "steps-mid", c = "darkgreen", lw = 2, label = "Run 19")
ax.fill_between(xSlice, ySlice, step = "mid", color = "lime", alpha = 1)


ax.grid()
ax.set_xlabel("Charge [ADC]", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
fig.suptitle("Multi Photon Spectrum", fontsize = 16)
ax.legend()
# ax.set_yscale("log")

ax.set_xlim(minX, maxX)

for p in peaks:
    ax.axvline(xSlice[p], c = "grey", ls = "--")


plt.show()
Figure
No description has been provided for this image

Fit¶

In [8]:
model = None
params = None

for i, p in enumerate(peaks):
    prefix = f"g{i}_"
    gauss = GaussianModel(prefix=prefix)

    if model is None:
        model = gauss
        params = gauss.make_params()
    else:
        model += gauss
        params.update(gauss.make_params())

    params[prefix+"center"].set(value = xSlice[p])
    params[prefix+"sigma"].set(value = results_half[i]/2.35, min=0)
    params[prefix+"amplitude"].set(value = ySlice[p], min=0)
    

# Baseline Gauss
prefix = f"g_b_"
gauss = GaussianModel(prefix=prefix)
params.update(gauss.make_params())

params[prefix+"center"].set(value = 3400)
params[prefix+"sigma"].set(value = 200, min=0)
params[prefix+"amplitude"].set(value = 300, min=0)

# # Baseline Const
# ConstantModel(prefix='c_')
# params.update(model.make_params(c_c=0))

# Fit
result = model.fit(yFit, params, x=xFit, weights=1/np.sqrt(yFit))
In [9]:
DEBUG = True
# DEBUG = False

# Plot
fig, ax = plt.subplots()
fig.set_size_inches(12,5)

ax.plot(xVect, data, ds = "steps-mid", c = "darkgreen", lw = 2, label = "Run 19")
ax.fill_between(xVect, data, step = "mid", color = "lime", alpha = 1)

plt.plot(xFit, result.best_fit, ls = ":", c = "k" , label="Fit (N Gauss + gauss)")

ax.grid()
ax.set_xlabel("Charge [ADC]", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
fig.suptitle("Multi Photon Spectrum with Multi Gaussian Fit", fontsize = 16)
ax.legend()
# ax.set_yscale("log")

ax.set_xlim(minX, maxX)

if DEBUG:
    for p in peaks:
        ax.axvline(xSlice[p], c = "grey", ls = "--")


plt.show()
Figure
No description has been provided for this image

Analysis on fit¶

In [10]:
# Extract mu from fit
centers = [] # Cannot use errors, most of them are undefined :(
for i in range(len(peaks)):
    cname = f"g{i}_center"
    centers.append(result.params[cname].value)
    
centers = np.array(centers)
gaussIdx = np.arange(0, len(centers) , 1)



# Fit with line
popt, pcov = curve_fit(myLine, gaussIdx, centers)
residuals = (centers - myLine(gaussIdx, *popt) ) / centers * 100



fig, ax = plt.subplots(2, 1, sharex = True,
                      gridspec_kw={"height_ratios":[3,1]})
fig.set_size_inches(12,5)

ax[0].scatter(gaussIdx, centers, marker = "*", color = "darkorchid", label = "$\mu$ extracted from fit")
ax[0].plot(gaussIdx, myLine(gaussIdx, *popt), ls = ":", c = "k", label = "p1 fit")

ax[1].scatter(gaussIdx, residuals, marker = "*", color = "darkorchid", label = r"Residuals ($y - y_\text{fit})$")


for a in ax:
    a.grid()
    a.legend()

ax[0].set_title("Fitted $\mu$ vs gaussian index", fontsize = 16)
ax[1].set_xlabel("Gaussian Index", fontsize = 14)
ax[0].set_ylabel("$\mu$", fontsize = 14)
ax[1].set_ylabel("Residuals [%]", fontsize = 14)

plt.show()
Figure
No description has been provided for this image
In [11]:
# Extract sigma from fit
sigmas = []
for i in range(len(peaks)):
    cname = f"g{i}_sigma"
    sigmas.append(result.params[cname].value)
    
sigmas = np.array(sigmas)[:-1]                        # Exclude last pt as it is patologic
gaussIdx = np.arange(0, len(centers) , 1)[:-1]



# Fit with line
popt, pcov = curve_fit(myLine, gaussIdx, sigmas)
residuals = (sigmas - myLine(gaussIdx, *popt) ) / sigmas * 100



fig, ax = plt.subplots(2, 1, sharex = True,
                      gridspec_kw={"height_ratios":[3,1]})
fig.set_size_inches(12,5)

ax[0].scatter(gaussIdx, sigmas, marker = "*", color = "darkorchid", label = "$\sigma$ extracted from fit")
ax[0].plot(gaussIdx, myLine(gaussIdx, *popt), ls = ":", c = "k", label = "p1 fit")

ax[1].scatter(gaussIdx, residuals, marker = "*", color = "darkorchid", label = r"Residuals ($y - y_\text{fit})$")


for a in ax:
    a.grid()
    a.legend()

fig.suptitle("Fitted $\sigma$ vs gaussian index", fontsize = 16)
ax[1].set_xlabel("Gaussian Index", fontsize = 14)
ax[0].set_ylabel("$\sigma$", fontsize = 14)
ax[1].set_ylabel("Residuals [%]", fontsize = 14)


plt.show()
Figure
No description has been provided for this image
In [12]:
diffs = np.diff(centers)
gaussIdx = np.arange(0, len(centers) -1 , 1) # Diff has one less point


# Fit with line
popt, pcov = curve_fit(myLine, gaussIdx, diffs)
residuals = (diffs - myLine(gaussIdx, *popt) ) / diffs * 100



fig, ax = plt.subplots(2, 1, sharex = True,
                      gridspec_kw={"height_ratios":[3,1]})
fig.set_size_inches(12,5)

ax[0].scatter(gaussIdx, diffs, marker = "*", color = "darkorchid", label = "$\mu$ extracted from fit")
ax[0].plot(gaussIdx, myLine(gaussIdx, *popt), ls = ":", c = "k", label = "p1 fit")

ax[1].scatter(gaussIdx, residuals, marker = "*", color = "darkorchid", label = r"Residuals ($y - y_\text{fit})$")


for a in ax:
    a.grid()
    a.legend()

fig.suptitle("Delta Peak to Peak - This is not good :(", fontsize = 16)
ax[1].set_xlabel("Gaussian Index", fontsize = 14)
ax[0].set_ylabel("$\Delta_\mu$", fontsize = 14)
ax[1].set_ylabel("Residuals [%]", fontsize = 14)

plt.show()
Figure
No description has been provided for this image

Poissonianity study¶

I try to do a sort of manual fit by minimizing the reduced $\Chi^2$ as a funcion of the $\mu$ parameter of a poissonian distribution, which is given by \begin{equation} P_\mu\left(n\right) = e^{-\mu}\,\dfrac{\mu^n}{n!} \end{equation}

But, as I don't know what is the real $n$ of the first peak, I do a spectra of such fits, varying each time $n_0$ in \begin{equation} P_\mu\left(n\right) = e^{-\mu}\,\dfrac{\mu^{(n+n_0)}}{(n+n_0)!} \end{equation}

I did this way as it is not that trivial to do a two-parameters manual fit

Furthermore, as the probability is normalized, we must normalize the amplitudes by diving for their sums. I know an appropriate normalization should take into account all the peaks, not just the one I am fitting. This is part of the error I accept to do

In [13]:
# Extract mu from fit
amplitudes = [] # Cannot use errors, most of them are undefined :(
for i in range(len(peaks)):
    cname = f"g{i}_amplitude"
    amplitudes.append(result.params[cname].value)
    
amplitudes = np.array(amplitudes)
amplitudes /= amplitudes.sum()
gaussIdx = np.arange(0, len(centers) , 1)
In [14]:
# Probability of having such counts for the n-th peak being mu the mean (to be estimated)
def myPoiss(n, mu, n0 = 0):
    return np.exp(-mu) * (mu **(n+n0)) / factorial((n+n0))

def myChi2(_bins, _h, _mu, n0_bypass):
    return np.sqrt(np.sum( ((_h - np.sum(_h) * myPoiss(_bins, _mu, n0_bypass))/np.sqrt(_h))**2 )) / (len(_bins)-1)
In [15]:
valoriScan = np.arange(10, 30, .5)
lstchi = []
cond = amplitudes > 0

lst_lstchi = []
lst_optmu = []
max_j = 8

for j in range(0, max_j):
    lstchi = np.array([ myChi2(gaussIdx[cond], amplitudes[cond], i, n0_bypass = j)  for i in  valoriScan])
    lst_lstchi.append(lstchi)
    lst_optmu.append(valoriScan[np.argmin(lstchi)])
    
    print(f"n0: {j} -- Opt mu: {lst_optmu[-1]} -- Min Chi: {np.min(lstchi):.4f}")

n0: 0 -- Opt mu: 11.5 -- Min Chi: 0.0076
n0: 1 -- Opt mu: 12.5 -- Min Chi: 0.0057
n0: 2 -- Opt mu: 13.5 -- Min Chi: 0.0039
n0: 3 -- Opt mu: 14.5 -- Min Chi: 0.0024
n0: 4 -- Opt mu: 15.5 -- Min Chi: 0.0020
n0: 5 -- Opt mu: 16.5 -- Min Chi: 0.0031
n0: 6 -- Opt mu: 17.5 -- Min Chi: 0.0046
n0: 7 -- Opt mu: 18.5 -- Min Chi: 0.0063
In [16]:
fig, ax  = plt.subplots()
fig.set_size_inches(12,5)


for j in range(0, max_j):
    ax.plot(valoriScan, lst_lstchi[j], label = f"n0 = {j}")

ax.set_title("Reduced chi2 vs $\mu$ for different $n_0$-s", fontsize = 16)
ax.set_xlabel("$\mu$", fontsize = 14)
ax.set_ylabel("Red. chi2", fontsize = 14)
ax.legend()

ax.grid()
plt.show()
Figure
No description has been provided for this image
In [17]:
fig, ax  = plt.subplots()
fig.set_size_inches(12,5)

# ax.plot(valoriScan, lstchi, "*g")

aVect = []
bVect = []

for j in range(0, max_j):
    aVect.append(j)
    bVect.append(lst_lstchi[j].min())
    
ax.plot(aVect, bVect, ls = ":", marker = "s", color = "darkorchid")
OPT_N0 = aVect[np.argmin(bVect)]

ax.set_title(f"Optimal Red. chi2 for each $n_0$ - Opt  $n_0$ is {OPT_N0}", fontsize = 16)
ax.set_xlabel("$n_0$", fontsize = 14)
ax.set_ylabel("Optimal Red. chi2", fontsize = 14)

ax.grid()
plt.show()
Figure
No description has been provided for this image
In [18]:
fig, ax = plt.subplots(2, 1, sharex = True,
                      gridspec_kw={"height_ratios":[3,1]})

fig.set_size_inches(12,5)

# Plot amplitude (poissonian) histogram
ax[0].plot(gaussIdx, amplitudes, ds = "steps-mid", c = "darkgreen", lw = 2, label = "Run 19")
ax[0].fill_between(gaussIdx, amplitudes, step = "mid", color = "lime", alpha = 1)

for j in range(0, max_j):
    tmpy = myPoiss(gaussIdx, lst_optmu[j], n0 = j)
    ax[0].plot(gaussIdx, tmpy,  label = f"Opt mu {lst_optmu[j]} - n0 = {j}")
    ax[1].plot(gaussIdx, (amplitudes - tmpy)/amplitudes, marker = "*")





for a in ax:
    a.grid()
    a.legend()

fig.suptitle("Optimal fitted poissonian for each $n_0$", fontsize = 16)
ax[1].set_xlabel("Gaussian Index", fontsize = 14)
ax[0].set_ylabel("$\Delta_\mu$", fontsize = 14)
ax[1].set_ylabel("Residuals [%]", fontsize = 14)

plt.show()

/tmp/ipykernel_2968125/1524409686.py:21: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  a.legend()
Figure
No description has been provided for this image
In [19]:
fig, ax = plt.subplots(2, 1, sharex = True,
                      gridspec_kw={"height_ratios":[3,1]})

fig.set_size_inches(12,5)

# Plot amplitude (poissonian) histogram
ax[0].plot(gaussIdx, amplitudes, ds = "steps-mid", c = "darkgreen", lw = 2, label = "Run 19")
ax[0].fill_between(gaussIdx, amplitudes, step = "mid", color = "lime", alpha = 1)

j = OPT_N0
tmpy = myPoiss(gaussIdx, lst_optmu[j], n0 = j)
ax[0].plot(gaussIdx, tmpy,  label = f"Opt mu {lst_optmu[j]} - n0 = {j}", c = "k", lw = 2)
ax[1].scatter(gaussIdx, (amplitudes - tmpy)/amplitudes, marker = "*", c = "k")





for a in ax:
    a.grid()
    a.legend()

fig.suptitle(f"n0 = {j}", fontsize = 16)
ax[1].set_xlabel("Gaussian Index", fontsize = 14)
ax[0].set_ylabel("$\Delta_\mu$", fontsize = 14)
ax[1].set_ylabel("Residuals [%]", fontsize = 14)

plt.show()

/tmp/ipykernel_2968125/916319418.py:21: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  a.legend()
Figure
No description has been provided for this image

After this section of analysis I can conclude that the first fitted peak seems to be $n=4$ even if by a close look of the spectrum I was tented to assume it to be $n=2$ But by looking closer at this histogram (and I did at the end of this analysis, to investigate what I obtained) i can confirm the first fitted peak is $n=4$

DATA ANALYSIS IS AWESOME

In [20]:
# Plot
fig, ax = plt.subplots()
fig.set_size_inches(12,5)

ax.plot(xVect, data, ds = "steps-mid", c = "darkgreen", lw = 2, label = "Run 19")
ax.fill_between(xVect, data, step = "mid", color = "lime", alpha = 1)
plt.plot(xFit, result.best_fit, ls = ":", c = "k" , label="Fit (N Gauss + gauss)")


ax.grid()
ax.set_xlabel("Charge [ADC]", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
fig.suptitle("Multi Photon Spectrum - Manual zoom until second fitted peak", fontsize = 16)
ax.legend()
ax.set_yscale("log")

ax.set_xlim(2500, 3000)

for p in peaks:
    ax.axvline(xSlice[p], c = "grey", ls = "--")


plt.show()
Figure
No description has been provided for this image

As a closing remark, the $\Chi^2$ I computed is not correclty dimensioned, as the counting error was adopted for a normalized distribution. Anyway, it is just rescaled for some constant I am to lazy to compute right now