In [1]:
import os
import sys
root_folder = os.path.dirname(os.getcwd())
sys.path.append(root_folder)
from ResoFit._pulse_shape import NeutronPulse
from ResoFit._pulse_shape import ProtonPulse
from ResoFit.experiment import Experiment
import numpy as np
from ResoFit.simulation import Simulation
import scipy.signal as ss
import matplotlib.pyplot as plt
import scipy

In [2]:
%matplotlib notebook

In [3]:
overwrite_csv = False

folder = 'data/IPTS_20440'
spectra_file = 'spectra.txt'
data_file1 = 'Ta_80C_12pC.csv'
norm_to_file = 'OB_80C_12pC.csv'
baseline = False
baseline_deg = 3

norm_factor = 1

image_start = None  # Can be omitted or =None
image_end = None  # Can be omitted or =None

source_to_detector_m = 16.46
offset_us = 0

In [4]:
experiment1 = Experiment(data_file=data_file1,
                         spectra_file=spectra_file,
                         folder=folder,
                         baseline=baseline,
                         baseline_deg=baseline_deg,
                         source_to_detector_m=source_to_detector_m,
                         offset_us=offset_us
                        )
experiment1.norm_to(file=norm_to_file, norm_factor=norm_factor)
experiment1.slice(start=image_start, end=image_end)

simulation = Simulation(energy_min=7.5, energy_max=15, energy_step=0.01, database='ENDF_VIII')
simulation.add_layer(layer='Ta', thickness_mm=0.127)

In [5]:
simulation._convolve_beam_shapes(source_to_detector_m=source_to_detector_m, conv_proton=False,
#                                  proton_params={'sigma':500},
                                )
x_n = simulation.x_tof_us
y_n = simulation.y_att


/Users/Shawn/anaconda3/envs/py37/lib/python3.7/site-packages/lmfit/models.py:34: FutureWarning:


The current behaviour of 'Series.argmin' is deprecated, use 'idxmin'
instead.
The behavior of 'argmin' will be corrected to return the positional
minimum in the future. For now, use 'series.values.argmin' or
'np.argmin(np.array(values))' to get the position of the minimum
row.

Fitted params for raw proton pulse shape:
Name          Value      Min      Max   Stderr     Vary     Expr Brute_Step
amplitude  1.984e+04     -inf      inf    36.58     True     None     None
center         1909     -inf      inf   0.3126     True     None     None
fwhm          345.7     -inf      inf   0.7362    False 2.3548200*sigma     None
height         53.9     -inf      inf  0.09941    False 0.3989423*amplitude/max(2.220446049250313e-16, sigma)     None
sigma         146.8        0      inf   0.3126     True     None     None
✅ '/Users/Shawn/Documents/GitHub/ResoFit/ResoFit/result/neutron_pulse/ikeda_carpenter/Neutron_fitted_params_1eV_500eV_ikeda_carpenter.csv' exists...
Fitted parameters file loaded.
✅ '/Users/Shawn/Documents/GitHub/ResoFit/ResoFit/result/neutron_pulse/ikeda_carpenter/Loglog_linear_within_1eV_500eV_ikeda_carpenter.csv' exists...
Parameters linear fitted file loaded.
✅ '/Users/Shawn/Documents/GitHub/ResoFit/ResoFit/result/neutron_pulse/ikeda_carpenter/TOF_shape_eV_7.5_15.0_0.01_us_0.05_183.7_0.01_for_sum_16.46m_ikeda_carpenter.csv' exists...
TOF neutron beam shape file loaded.

In [6]:
simulation._convolve_beam_shapes(source_to_detector_m=source_to_detector_m, conv_proton=True)
x_np = simulation.x_tof_us
y_np = simulation.y_att


Fitted params for raw proton pulse shape:
Name          Value      Min      Max   Stderr     Vary     Expr Brute_Step
amplitude  1.984e+04     -inf      inf    36.58     True     None     None
center         1909     -inf      inf   0.3126     True     None     None
fwhm          345.7     -inf      inf   0.7362    False 2.3548200*sigma     None
height         53.9     -inf      inf  0.09941    False 0.3989423*amplitude/max(2.220446049250313e-16, sigma)     None
sigma         146.8        0      inf   0.3126     True     None     None
✅ '/Users/Shawn/Documents/GitHub/ResoFit/ResoFit/result/neutron_pulse/ikeda_carpenter/Neutron_fitted_params_1eV_500eV_ikeda_carpenter.csv' exists...
Fitted parameters file loaded.
✅ '/Users/Shawn/Documents/GitHub/ResoFit/ResoFit/result/neutron_pulse/ikeda_carpenter/Loglog_linear_within_1eV_500eV_ikeda_carpenter.csv' exists...
Parameters linear fitted file loaded.
✅ '/Users/Shawn/Documents/GitHub/ResoFit/ResoFit/result/neutron_pulse/ikeda_carpenter/TOF_shape_eV_7.5_15.0_0.01_us_0.05_183.7_0.01_proton_for_sum_16.46m_ikeda_carpenter.csv' exists...
TOF neutron beam shape file loaded.

In [7]:
simulation._convolve_beam_shapes(source_to_detector_m=source_to_detector_m, conv_proton=True,
                                 proton_params={'sigma':300},
                                )
x_np300 = simulation.x_tof_us
y_np300 = simulation.y_att


Fitted params for raw proton pulse shape:
Name          Value      Min      Max   Stderr     Vary     Expr Brute_Step
amplitude  1.984e+04     -inf      inf    36.58     True     None     None
center         1909     -inf      inf   0.3126     True     None     None
fwhm          345.7     -inf      inf   0.7362    False 2.3548200*sigma     None
height         53.9     -inf      inf  0.09941    False 0.3989423*amplitude/max(2.220446049250313e-16, sigma)     None
sigma         146.8        0      inf   0.3126     True     None     None
✅ '/Users/Shawn/Documents/GitHub/ResoFit/ResoFit/result/neutron_pulse/ikeda_carpenter/Neutron_fitted_params_1eV_500eV_ikeda_carpenter.csv' exists...
Fitted parameters file loaded.
✅ '/Users/Shawn/Documents/GitHub/ResoFit/ResoFit/result/neutron_pulse/ikeda_carpenter/Loglog_linear_within_1eV_500eV_ikeda_carpenter.csv' exists...
Parameters linear fitted file loaded.
✅ '/Users/Shawn/Documents/GitHub/ResoFit/ResoFit/result/neutron_pulse/ikeda_carpenter/TOF_shape_eV_7.5_15.0_0.01_us_0.05_183.7_0.01_proton_sigma_300_for_sum_16.46m_ikeda_carpenter.csv' exists...
TOF neutron beam shape file loaded.

In [8]:
simulation._convolve_beam_shapes(source_to_detector_m=source_to_detector_m, conv_proton=True,
                                 proton_params={'sigma':500},
                                )
x_np500 = simulation.x_tof_us
y_np500 = simulation.y_att


Fitted params for raw proton pulse shape:
Name          Value      Min      Max   Stderr     Vary     Expr Brute_Step
amplitude  1.984e+04     -inf      inf    36.58     True     None     None
center         1909     -inf      inf   0.3126     True     None     None
fwhm          345.7     -inf      inf   0.7362    False 2.3548200*sigma     None
height         53.9     -inf      inf  0.09941    False 0.3989423*amplitude/max(2.220446049250313e-16, sigma)     None
sigma         146.8        0      inf   0.3126     True     None     None
✅ '/Users/Shawn/Documents/GitHub/ResoFit/ResoFit/result/neutron_pulse/ikeda_carpenter/Neutron_fitted_params_1eV_500eV_ikeda_carpenter.csv' exists...
Fitted parameters file loaded.
✅ '/Users/Shawn/Documents/GitHub/ResoFit/ResoFit/result/neutron_pulse/ikeda_carpenter/Loglog_linear_within_1eV_500eV_ikeda_carpenter.csv' exists...
Parameters linear fitted file loaded.
✅ '/Users/Shawn/Documents/GitHub/ResoFit/ResoFit/result/neutron_pulse/ikeda_carpenter/TOF_shape_eV_7.5_15.0_0.01_us_0.05_183.7_0.01_proton_sigma_500_for_sum_16.46m_ikeda_carpenter.csv' exists...
TOF neutron beam shape file loaded.

In [15]:
import pandas as pd
from scipy.interpolate import interp1d

In [31]:
d = 5.2

In [14]:
df1 = pd.DataFrame()
# df1['x_n'] = x_n+1.95-d
df1['x_n'] = x_n
df1['y_n'] = y_n
df1.dropna(inplace=True)
df1.reset_index(inplace=True, drop=True)
df1


Out[14]:
x_n y_n
0 128.246541 0.018150
1 128.247443 0.018150
2 128.247697 0.018150
3 128.247710 0.018150
4 128.248458 0.018150
... ... ...
36953 163.971566 0.005759
36954 163.977361 0.005759
36955 163.980180 0.005796
36956 163.985568 0.005796
36957 163.988796 0.005796

36958 rows × 2 columns


In [15]:
df2 = pd.DataFrame()
df2['x_np'] = x_np+0.03-d
df2['y_np'] = y_np
df2.dropna(inplace=True)
df2.reset_index(inplace=True, drop=True)
df2


Out[15]:
x_np y_np
0 127.245751 0.018150
1 127.245848 0.018150
2 127.246532 0.018150
3 127.246550 0.018150
4 127.247066 0.018150
... ... ...
35055 170.285171 0.005687
35056 170.293539 0.005722
35057 170.301909 0.005759
35058 170.310281 0.005759
35059 170.318654 0.005796

35060 rows × 2 columns


In [16]:
df3 = pd.DataFrame()
df3['x_np300'] = x_np300-d
df3['y_np300'] = y_np300
df3.dropna(inplace=True)
df3.reset_index(inplace=True, drop=True)
df3


Out[16]:
x_np300 y_np300
0 120.376140 0.018150
1 120.376690 0.018150
2 120.377151 0.018150
3 120.377499 0.018150
4 120.377680 0.018150
... ... ...
38149 164.297029 0.005687
38150 164.305405 0.005722
38151 164.313783 0.005759
38152 164.322163 0.005759
38153 164.330544 0.005796

38154 rows × 2 columns


In [17]:
df4 = pd.DataFrame()
df4['x_np500'] = x_np500-d
df4['y_np500'] = y_np500
df4.dropna(inplace=True)
df4.reset_index(inplace=True, drop=True)
df4


Out[17]:
x_np500 y_np500
0 120.182229 0.018150
1 120.190239 0.017839
2 120.198250 0.018058
3 120.206263 0.017765
4 120.214277 0.017601
... ... ...
38347 164.389256 0.005687
38348 164.397650 0.005722
38349 164.406045 0.005722
38350 164.414442 0.005759
38351 164.422841 0.005796

38352 rows × 2 columns


In [18]:
df = pd.DataFrame()
df['x_n'] = df1['x_n']
df['y_n'] = df1['y_n']
df['x_np'] = df2['x_np']
df['y_np'] = df2['y_np']
df['x_np300'] = df3['x_np300']
df['y_np300'] = df3['y_np300']
df['x_np500'] = df4['x_np500']
df['y_np500'] = df4['y_np500']
df.dropna(inplace=True)
df


Out[18]:
x_n y_n x_np y_np x_np300 y_np300 x_np500 y_np500
0 128.246541 0.018150 127.245751 0.018150 120.376140 0.018150 120.182229 0.018150
1 128.247443 0.018150 127.245848 0.018150 120.376690 0.018150 120.190239 0.017839
2 128.247697 0.018150 127.246532 0.018150 120.377151 0.018150 120.198250 0.018058
3 128.247710 0.018150 127.246550 0.018150 120.377499 0.018150 120.206263 0.017765
4 128.248458 0.018150 127.247066 0.018150 120.377680 0.018150 120.214277 0.017601
... ... ... ... ... ... ... ... ...
35055 155.131795 0.040106 170.285171 0.005687 147.709181 0.040488 146.982212 0.040725
35056 155.134417 0.040105 170.293539 0.005722 147.710011 0.040488 146.985921 0.040724
35057 155.140131 0.040105 170.301909 0.005759 147.717303 0.040486 146.990510 0.040723
35058 155.142443 0.040105 170.310281 0.005759 147.718446 0.040486 146.994517 0.040722
35059 155.148468 0.040102 170.318654 0.005796 147.725427 0.040484 146.998810 0.040721

35060 rows × 8 columns


In [19]:
t_max = 131.87
t_min = 128.54
t_step = 0.01
nbr_point = int((t_max - t_min) / t_step + 1)
ls_x = ['x_n', 'x_np', 'x_np300', 'x_np500']
ls_y = ['y_n', 'y_np', 'y_np300', 'y_np500']
df_out = pd.DataFrame()
x_axis = np.linspace(t_min, t_max, nbr_point).round(5)
df_out['x_interp'] = x_axis
for i, each in enumerate(ls_y):
    y_axis_function = interp1d(x=list(df[ls_x[i]]), y=list(df[each]), kind='cubic')
    df_out[each] = y_axis_function(x_axis)

In [20]:
df_out


Out[20]:
x_interp y_n y_np y_np300 y_np500
0 128.54 0.021277 0.027801 0.045665 0.045717
1 128.55 0.021771 0.028229 0.045663 0.045714
2 128.56 0.022300 0.028665 0.045660 0.045711
3 128.57 0.022864 0.029109 0.045657 0.045709
4 128.58 0.023463 0.029561 0.045654 0.045706
... ... ... ... ... ...
329 131.83 0.021529 0.022256 0.044742 0.044796
330 131.84 0.021797 0.022400 0.044740 0.044793
331 131.85 0.022080 0.022557 0.044737 0.044791
332 131.86 0.022378 0.022727 0.044734 0.044788
333 131.87 0.022689 0.022911 0.044731 0.044785

334 rows × 5 columns


In [55]:
df_out.to_clipboard()

In [18]:


In [16]:
simulation.export(x_axis='time', y_axis='attenuation', source_to_detector_m=source_to_detector_m, offset_us=2.78)


'Time (us)' was obtained with the following:
source_to_detector_m=16.45
offset_us=2.78
Exporting to clipboard completed.

In [18]:
experiment1.export(x_type='time', y_type='attenuation', t_unit='us')


Out[18]:
x y
0 0.96 0.032947
1 1.12 0.028682
2 1.28 0.023784
3 1.44 0.020349
4 1.60 0.022276
5 1.76 0.022747
6 1.92 0.020167
7 2.08 0.018862
8 2.24 0.023378
9 2.40 0.028210
10 2.56 0.027037
11 2.72 0.024648
12 2.88 0.020323
13 3.04 0.019840
14 3.20 0.019051
15 3.36 0.014892
16 3.52 0.011814
17 3.68 0.007386
18 3.84 0.006987
19 4.00 0.013964
20 4.16 0.014205
21 4.32 0.016009
22 4.48 0.014003
23 4.64 0.015825
24 4.80 0.017421
25 4.96 0.015155
26 5.12 0.012540
27 5.28 0.012222
28 5.44 0.009548
29 5.60 0.010404
... ... ...
2776 445.12 0.005094
2777 445.28 0.003596
2778 445.44 0.004051
2779 445.60 0.009737
2780 445.76 0.007203
2781 445.92 0.013659
2782 446.08 0.010629
2783 446.24 0.011610
2784 446.40 0.007048
2785 446.56 0.009207
2786 446.72 0.008866
2787 446.88 0.003544
2788 447.04 0.010973
2789 447.20 0.007840
2790 447.36 0.009436
2791 447.52 0.003415
2792 447.68 0.008944
2793 447.84 0.009021
2794 448.00 0.013826
2795 448.16 0.007227
2796 448.32 0.012424
2797 448.48 0.011139
2798 448.64 0.009720
2799 448.80 0.005232
2800 448.96 0.011359
2801 449.12 0.010899
2802 449.28 0.008012
2803 449.44 0.010435
2804 449.60 0.009637
2805 449.76 0.024209

2806 rows × 2 columns


In [9]:
ax = simulation.plot(x_type='time', y_type='attenuation', source_to_detector_m=source_to_detector_m, offset_us=offset_us,
                          mixed=False, all_layers=True)
ax = experiment1.plot(x_type='time', y_type='attenuation', t_unit='us', ax_mpl=ax)
ax.plot(x_n, y_n, label='Conv. Neutron')
ax.plot(x_np, y_np, label='Conv. Neutron*Proton(\u03C3=146.8)')
ax.plot(x_np300, y_np300, label='Conv. Neutron*Proton(\u03C3=300)')
ax.plot(x_np500, y_np500, label='Conv. Neutron*Proton(\u03C3=500)')
# experiment1.plot_raw(x_type='time', y_type='transmission', time_unit='us', ax_mpl=fig.axes[0])
# plt.plot(simulation.x_tof_us-5.2, 1-simulation.y_att, label='Conv. N*P(\u03C3=146.8)')
# ax.set_xlim(left=1, right=6)
ax.legend()
ax.set_title('Gd resonance dip at ~4eV')
ax.grid()


'Time (us)' was obtained with the following:
source_to_detector_m=16.46
offset_us=0

In [11]:
ax = simulation.plot(x_type='time', y_type='attenuation', source_to_detector_m=source_to_detector_m, offset_us=offset_us,
                          mixed=False, all_layers=True)
ax = experiment1.plot(x_type='time', y_type='attenuation', t_unit='us', ax_mpl=ax)
ax.plot(x_n, y_n, label='Conv. Neutron')
ax.plot(x_np, y_np/2, label='Conv. Neutron*Proton(\u03C3=146.8)')
ax.plot(x_np300, y_np300, label='Conv. Neutron*Proton(\u03C3=300)')
ax.plot(x_np500, y_np500, label='Conv. Neutron*Proton(\u03C3=500)')
# experiment1.plot_raw(x_type='time', y_type='transmission', time_unit='us', ax_mpl=fig.axes[0])
# plt.plot(simulation.x_tof_us-5.2, 1-simulation.y_att, label='Conv. N*P(\u03C3=146.8)')
# ax.set_xlim(left=1, right=6)
ax.legend()
ax.set_title('Ta resonance dip from 7-15eV')
ax.grid()


'Time (us)' was obtained with the following:
source_to_detector_m=16.46
offset_us=0

In [10]:
simulation.neutron_pulse.proton_pulse.plot()


Out[10]:
<matplotlib.axes._subplots.AxesSubplot at 0x1107fdcc0>

In [9]:
total_fig = simulation.neutron_pulse.plot_shape_total()



In [ ]: