Commit 43cfb842 authored by feichtinger's avatar feichtinger
Browse files

Committing changes I did over a month ago

Doing this as I am preparing for the meeting tomorrow. I want to
restructure and clarify a bit, so that I directly can use the
notebook during my presentation.
parent ff5cc58b
......@@ -32,6 +32,7 @@ import pandas as pd
from pandas.plotting import scatter_matrix
import matplotlib.pyplot as plt
import seaborn as sns
import scipy
import scipy.stats
#matplotlib inline
# -
......@@ -167,7 +168,6 @@ def makeDataSetInterpolated(directory, excelFn, doInterpolate=True, dropBadPulse
def makeDataSetInterpolated2(directory, excelFn, doInterpolate=True, dropBadPulses=True, verbose=False,
CALCTthreshold=-50, CALCSthreshold=-50):
first = True
data = pd.DataFrame()
for filename in sorted(glob.glob(os.path.join(directory,'dp*-nomeans.csv'))):
mat = re.match('^dp(\d+).*', os.path.basename(filename))
......@@ -285,7 +285,7 @@ tmpdata = data.query(" | ".join(flist)) \
for col in colname:
sns.lmplot(x="index", y=col, hue="rawDataFile",
data=tmpdata,
fit_reg=False, size=4, aspect=3)
fit_reg=False, height=4, aspect=3)
# -
......@@ -324,7 +324,7 @@ for col in colname:
sns.lmplot(x="index", y=col, hue="rawDataFile",
data=tmpdata,
palette=sns.color_palette("Set1", len(flist)),
fit_reg=False, size=4, aspect=3)
fit_reg=False, height=4, aspect=3)
# -
......@@ -335,14 +335,14 @@ for col in colname:
sns.lmplot(x="PHOTON-ENERGY-PER-PULSE", y='CALCT', hue="rawDataFile",
data=data,
palette=sns.color_palette("Set1", len(flist)),
fit_reg=False, size=4, aspect=3, legend=False)
fit_reg=False, height=4, aspect=3, legend=False)
sns.lmplot(x="PEPP", y='CALCT', hue="rawDataFile",
data=data.rename(columns = {'PHOTON-ENERGY-PER-PULSE': 'PEPP'}) \
.query('180 < PEPP < 250'),
palette=sns.color_palette("Set1", len(flist)),
fit_reg=False, size=4, aspect=3, legend=True)
fit_reg=False, height=4, aspect=3, legend=True)
# So, the points are all part of file dp05
......@@ -477,6 +477,14 @@ physdf = pd.concat([data_mean.reset_index().set_index('SPECTRUM_CENTER'), tmpdf]
physdf.head()
# -
fig,ax = plt.subplots(figsize=(12,6))
physdf.plot(x='SPECTRUM_CENTER', y='cross-section',
linestyle='', marker='o', ax=ax, legend=False)
#ax.set_xscale('log')
#ax.set_yscale('log')
ax.set_ylabel('cross section')
ax.set_title("Xe cross section across data's photon energy range")
# The Electromultiplier signals like CALCS should be proportional to the number of photoions hitting the detector (multiplier was operated within the linear voltage range). We do not know the number of photons hitting the EM, but we know that it must be proportional to the number of photons $N_{ph}$ in the beam, which we obtain from the ion current detector.
#
# \begin{aligned} N_{ph} = \frac{\Phi}{\hbar \omega} = \frac{k}{e} \dot{} \frac{I_{\text {ion}} T}{\gamma \sigma z p} \end{aligned}
......@@ -488,7 +496,6 @@ physdf.head()
# \end{aligned}
physdf['Nphot_avg'] = physdf['PHOTON-ENERGY-PER-PULSE'] / physdf['SPECTRUM_CENTER']
physdf['CALCSUM'] = physdf.CALCS + physdf.CALCT
physdf.describe()
......@@ -521,9 +528,9 @@ agrid = sns.lmplot(x='Nphot_avg', y='CALCS', hue="rawDataFile",
for row in range(0, tmpdf.shape[0]):
#agrid.axes[0][0].text(tmpdf.Nphot_avg.iloc[0]+0.002, tmpdf.CALCS.iloc[0], "test")
agrid.ax.text(tmpdf.Nphot_avg.iloc[row]+0.0005, tmpdf.CALCS.iloc[row],
np.round(tmpdf.XeMultVoltag.iloc[row],1))
f'V:{np.round(tmpdf.XeMultVoltag.iloc[row],1)}')
agrid.ax.text(tmpdf.Nphot_avg.iloc[row]+0.0005, tmpdf.CALCS.iloc[row]-440,
np.round(tmpdf.SPECTRUM_CENTER.iloc[row],1))
f'E:{np.round(tmpdf.SPECTRUM_CENTER.iloc[row],1)}')
# +
......@@ -535,11 +542,15 @@ agrid = sns.lmplot(x='Nphot_avg', y='CALCS', hue="rawDataFile",
fit_reg=False, height=4, aspect=3)
for row in range(0, tmpdf.shape[0]):
#agrid.axes[0][0].text(tmpdf.Nphot_avg.iloc[0]+0.002, tmpdf.CALCS.iloc[0], "test")
agrid.ax.text(tmpdf.Nphot_avg.iloc[row]+0.0005, tmpdf.CALCS.iloc[row],
np.round(tmpdf.XeMultVoltag.iloc[row],1))
f'V:{np.round(tmpdf.XeMultVoltag.iloc[row],1)}')
agrid.ax.text(tmpdf.Nphot_avg.iloc[row]+0.0005, tmpdf.CALCS.iloc[row]-50,
np.round(tmpdf.SPECTRUM_CENTER.iloc[row],1))
f'E:{np.round(tmpdf.SPECTRUM_CENTER.iloc[row],1)}')
# -
# In the above graph, the dp01 file is a marked outlier. It has almost the same Voltage as dp32 (926). It is taken a lower photon energy (6101).
# +
tmpdf = physdf.query('1400 < XeMultVoltag < 1530')
......@@ -557,12 +568,9 @@ for row in range(0, tmpdf.shape[0]):
# -
import scipy
# ## fitting the gain equation
# \begin{aligned} {\text CALC} = N_{\text ph} \cdot A \cdot V^{B \cdot n} \end{aligned}
# \begin{aligned} {\text CALC} = N_{\text ph} \cdot A \cdot V^{B} \end{aligned}
#
# x input is composed of 2 arrays: [ nphot, voltage]
......@@ -606,32 +614,91 @@ popt, pcov = scipy.optimize.curve_fit(calcs_fn,
popt
# Let's plot the gain function based on these parameters. Note, that it is not the true gain function of the electromultiplier, because we have it in function of $N_{ph}$ and not of $N_{ion}$.
fig,ax = plt.subplots()
v_ax = np.linspace(physdf.XeMultVoltag.min(),physdf.XeMultVoltag.max(), 100)
# we get the gain by setting Nphot=1.0
ax.plot(v_ax, -calcs_fn([np.full(v_ax.shape, 1.0 ), v_ax], popt[0], popt[1]))
ax.set_yscale('log')
ax.set_xlabel('XeMultVoltag')
ax.set_ylabel('Gain')
# the predictions for CALCS given the experiments Nphot [from PEP / h \omega)]
# predictions for PEP based on the measured CALCS
physdf = physdf.assign(CALCS_pred = lambda x: calcs_fn([x.Nphot_avg, x.XeMultVoltag], popt[0], popt[1]),
PEP_pred = lambda x: x.SPECTRUM_CENTER * x.CALCS / (popt[0] * np.power(x.XeMultVoltag, popt[1])))
fig, ax = plt.subplots()
ax.plot(physdf.Nphot_avg, physdf.CALCS, label='data',linestyle='', marker='o')
ax.plot(physdf.Nphot_avg, calcs_fn([physdf.Nphot_avg, physdf.XeMultVoltag], popt[0], popt[1]),
ax.plot(physdf.Nphot_avg, physdf.CALCS_pred,
label='fit', linestyle='', marker='o')
ax.set_xlabel('Nphot')
ax.set_ylabel('CALCS')
ax.legend()
physdf.columns
# measured vs. predicted
fig,ax = plt.subplots(figsize=(8,4))
ax.plot(physdf.CALCS,
calcs_fn([physdf.Nphot_avg, physdf.XeMultVoltag], popt[0], popt[1]),
ax.plot(physdf.CALCS, physdf.CALCS_pred,
linestyle='', marker='o', alpha=0.6)
ax.set_xlabel('CALCS measured')
ax.set_ylabel('CALCS predicted')
# Let's plot the gain function based on these parameters. Note, that it is not the true gain function of the electromultiplier, because we have it in function of $N_{ph}$ and not of $N_{ion}$.
# measured vs. predicted
fig,ax = plt.subplots(figsize=(8,4))
ax.plot(physdf['PHOTON-ENERGY-PER-PULSE'], physdf.PEP_pred,
linestyle='', marker='o', alpha=0.6)
ax.set_xlabel('PEP measured')
ax.set_ylabel('PEP predicted')
# ## using cutoff for CALCS < 4000 uJ
physdf_cut = physdf.query('CALCS > -4000').copy()
popt_cut, pcov_cut = scipy.optimize.curve_fit(calcs_fn,
[ physdf_cut.Nphot_avg, physdf_cut.XeMultVoltag ],
physdf_cut.CALCS,
[-5e-16, 6.0],
maxfev=10000
)
# **Results for the fitted constants A and B of the gain equation**
popt_cut
# the predictions for CALCS given the experiments Nphot [from PEP / h \omega)]
# predictions for PEP based on the measured CALCS
physdf_cut = physdf_cut.assign(CALCS_pred = lambda x: calcs_fn([x.Nphot_avg, x.XeMultVoltag], popt[0], popt[1]),
PEP_pred = lambda x: x.SPECTRUM_CENTER * x.CALCS / (popt[0] * np.power(x.XeMultVoltag, popt[1])))
fig, ax = plt.subplots()
ax.plot(physdf_cut.Nphot_avg, physdf_cut.CALCS, label='data',linestyle='', marker='o')
ax.plot(physdf_cut.Nphot_avg, physdf_cut.CALCS_pred,
label='fit', linestyle='', marker='o')
ax.set_xlabel('Nphot')
ax.set_ylabel('CALCS')
ax.legend()
fig,ax = plt.subplots()
v_ax = np.linspace(physdf.XeMultVoltag.min(),physdf.XeMultVoltag.max(), 100)
# we get the gain by setting Nphot=1.0
ax.plot(v_ax, -calcs_fn([np.full(v_ax.shape, 1.0 ), v_ax], popt[0], popt[1]))
ax.set_yscale('log')
ax.set_xlabel('XeMultVoltag')
ax.set_ylabel('Gain')
# measured vs. predicted
fig,ax = plt.subplots(figsize=(8,4))
ax.plot(physdf_cut.CALCS,
physdf_cut.CALCS_pred,
linestyle='', marker='o', alpha=0.6)
ax.set_xlabel('CALCS measured')
ax.set_ylabel('CALCS predicted')
# measured vs. predicted
fig,ax = plt.subplots(figsize=(8,4))
ax.plot(physdf_cut['PHOTON-ENERGY-PER-PULSE'], physdf_cut.PEP_pred,
linestyle='', marker='o', alpha=0.6)
ax.set_xlabel('PEP measured')
ax.set_ylabel('PEP predicted')
physdf_cut
# # deduce information from CALCS distribution as function of voltage
# As described above, we can try to obtain the constant B from the gain formula by looking at the variation of the CALCS values in function of the supply voltage
......@@ -696,4 +763,52 @@ lreg.slope, lreg.intercept
# So, this method gave 6.12 for the parameter B, while the regression of the simple gain equation (earlier section) gave 6.5.
# # Looking at relationship of CALCS and CALCT
# Looking at the response from both electromultipliers. What is the expected relation if the gain equation is valid?
#
# \begin{aligned}
# \frac{\text CALCS}{\text CALCT} = \frac{N_{\text ph} \cdot A_S \cdot V^{B_S}}{N_{\text ph} \cdot A_T \cdot V^{B_T}}
# = \frac{A_S}{A_T} \cdot V^{B_S - B_T}
# \end{aligned}
#
# If both electromultipliers are very similar, then $B_S \approx B_T$ and there should be a linear relationship between CALCS and CALCT when operated at the same Voltage. In the data, we see a linear relationship over a wide range of the lower response function. Towards the higher values (abs(CALCS) > 15000) there is a "bend" in the relationship.
fig,ax = plt.subplots(figsize=(12,6))
ax.plot(data.CALCS, data.CALCT, linestyle='', marker='o',alpha=0.05)
ax.set_xlabel('CALCS')
ax.set_ylabel('CALCT')
fig,ax = plt.subplots(2,1,figsize=(12,6))
datatmp = data.query('CALCS > -15000') \
.assign(calcs_d_t = lambda x: x.CALCS/x.CALCT)
ax[0].plot(datatmp.CALCS, datatmp.CALCT, linestyle='', marker='o',alpha=0.05)
ax[0].set_xlabel('CALCS')
ax[0].set_ylabel('CALCT')
ax[1].plot(datatmp.XeMultVoltag,datatmp.calcs_d_t, linestyle='', marker='o', alpha=0.05)
ax[1].set_ylabel('CALCS/CALCT')
ax[1].set_xlabel('supply Voltage')
# +
#fig,ax = plt.subplots(2,1,figsize=(12,6))
#datatmp = data.query('CALCS > -15000') \
# .assign(calcs_d_t = lambda x: x.CALCS/x.CALCT)
#ax[0].plot(datatmp.CALCS, datatmp.CALCT, linestyle='', marker='o',alpha=0.05)
#ax[0].set_xlabel('CALCS')
#ax[0].set_ylabel('CALCT')
#ax[1].plot(datatmp.XeMultVoltag,datatmp.calcs_d_t, c=datatmp.SPECTRUM_CENTER, linestyle='', marker='o', alpha=0.05)
#ax[1].set_ylabel('CALCS/CALCT')
#ax[1].set_xlabel('supply Voltage')
# -
fig,ax = plt.subplots(figsize=(12,6))
ax.plot(data.CALCS/data.CALCT, data.XeMultVoltag, linestyle='', marker='o', alpha=0.05)
fig,ax = plt.subplots(figsize=(12,6))
datatmp = data.query('CALCS > -15000')
ax.plot(datatmp.CALCS/datatmp.CALCT, datatmp.SPECTRUM_CENTER, linestyle='', marker='o', alpha=0.05)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment