Beta
Table of Contents
The outline of your notebook will show up here. You can include headings in any text cell by starting a line with # , ## , ### , etc., depending on the desired title hierarchy.
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
from scipy import integrate
from scipy.optimize import fsolve
import math
import sympy as sp
from scipy.interpolate import interp1d
from numpy.polynomial.polynomial import Polynomial
#x1 = 800
#x2 = 1200
x1 = 900
x2 = 1050
#x1_ = 1150
x2_ = 1300
x1_ = 1200
#x2_ = 1300
c = 3*10**8 #m/s
xdat = np.linspace(x1, x2, 10000)
xdat_953 = np.linspace(940, 965, 10000)
xdat_964 = np.linspace(950, 975, 10000)
xdat_ = np.linspace(x1_, x2_, 10000)
xdat_1232 = np.linspace(1222, 1242, 10000)
xdat_1245 = np.linspace(1230, 1260, 10000)
t_data = np.array([10, 60, 120, 150, 290])*10**-6 # Updated time points
theta_ring1 = 157.6
#2+ -> 0+
# Setting name of the file that the data is to be extracted from in python
#file_JG1 = 'APPACOMM_EuBa_60Co.txt'
file_JG1_ring1_10um = 'Mg_Nb_10um_ring1.txt'
file_JG1_ring2_10um = 'Mg_Nb_10um_ring2.txt'
file_JG1_ring1_30um = 'Mg_Nb_30um_ring1.txt'
file_JG1_ring2_30um = 'Mg_Nb_30um_ring2.txt'
file_JG1_ring1_60um = 'Mg_Nb_60um_ring1.txt'
file_JG1_ring2_60um = 'Mg_Nb_60um_ring2.txt'
file_JG1_ring1_100um = 'Mg_Nb_100um_ring1.txt'
file_JG1_ring2_100um = 'Mg_Nb_100um_ring2.txt'
file_JG1_ring1_120um = 'Mg_Nb_120um_ring1.txt'
file_JG1_ring2_120um = 'Mg_Nb_120um_ring2.txt'
file_JG1_ring1_150um = 'Mg_Nb_150um_ring1.txt'
file_JG1_ring2_150um = 'Mg_Nb_150um_ring2.txt'
file_JG1_ring1_290um = 'Mg_Nb_290um_ring1.txt'
file_JG1_ring2_290um = 'Mg_Nb_290um_ring2.txt'
# Loading file data into numpy array and storing it in variable called data_collected
JG1_data_10um_ring1 = np.loadtxt(file_JG1_ring1_10um)
JG1_data_10um_ring2 = np.loadtxt(file_JG1_ring2_10um)
JG1_data_30um_ring1 = np.loadtxt(file_JG1_ring1_30um)
JG1_data_30um_ring2 = np.loadtxt(file_JG1_ring2_30um)
JG1_data_60um_ring1 = np.loadtxt(file_JG1_ring1_60um)
JG1_data_60um_ring2 = np.loadtxt(file_JG1_ring2_60um)
JG1_data_100um_ring1 = np.loadtxt(file_JG1_ring1_100um)
JG1_data_100um_ring2 = np.loadtxt(file_JG1_ring2_100um)
JG1_data_120um_ring1 = np.loadtxt(file_JG1_ring1_120um)
JG1_data_120um_ring2 = np.loadtxt(file_JG1_ring2_120um)
JG1_data_150um_ring1 = np.loadtxt(file_JG1_ring1_150um)
JG1_data_150um_ring2 = np.loadtxt(file_JG1_ring2_150um)
JG1_data_290um_ring1 = np.loadtxt(file_JG1_ring1_290um)
JG1_data_290um_ring2 = np.loadtxt(file_JG1_ring2_290um)
#regions for 2+ -> 0+ transition
x_10 = JG1_data_10um_ring1[x1:x2,0]
y_10 = JG1_data_10um_ring1[x1:x2,1]
x_30 = JG1_data_30um_ring1[x1:x2,0]
y_30 = JG1_data_30um_ring1[x1:x2,1]
x_60 = JG1_data_60um_ring1[x1:x2,0]
y_60 = JG1_data_60um_ring1[x1:x2,1]
x_100 = JG1_data_100um_ring1[x1:x2,0]
y_100 = JG1_data_100um_ring1[x1:x2,1]
x_120 = JG1_data_120um_ring1[x1:x2,0]
y_120 = JG1_data_120um_ring1[x1:x2,1]
x_150 = JG1_data_150um_ring1[x1:x2,0]
y_150 = JG1_data_150um_ring1[x1:x2,1]
x_290 = JG1_data_290um_ring1[x1:x2,0]
y_290 = JG1_data_290um_ring1[x1:x2,1]
#regions for 4+ -> 2+ transition
x_f_10 = JG1_data_10um_ring1[x1_:x2_,0]
y_f_10 = JG1_data_10um_ring1[x1_:x2_,1]
x_f_30 = JG1_data_30um_ring1[x1_:x2_,0]
y_f_30 = JG1_data_30um_ring1[x1_:x2_,1]
x_f_60 = JG1_data_60um_ring1[x1_:x2_,0]
y_f_60 = JG1_data_60um_ring1[x1_:x2_,1]
x_f_100 = JG1_data_100um_ring1[x1_:x2_,0]
y_f_100 = JG1_data_100um_ring1[x1_:x2_,1]
x_f_120 = JG1_data_120um_ring1[x1_:x2_,0]
y_f_120 = JG1_data_120um_ring1[x1_:x2_,1]
x_f_150 = JG1_data_150um_ring1[x1_:x2_,0]
y_f_150 = JG1_data_150um_ring1[x1_:x2_,1]
x_f_290 = JG1_data_290um_ring1[x1_:x2_,0]
y_f_290 = JG1_data_290um_ring1[x1_:x2_,1]
##FUNCTIONS
def bkg(x, a, b):
return a*x + b
def gaussian(x, amp, mu, sig, m, c):#, a, b):
fun_val = amp * np.exp(-np.power(x-mu, 2) / (2*np.power(sig, 2))) + m * x + c
return fun_val
def four_gaussian_with_background(x, A1, mu1, sigma1, A2, mu2, sigma2, A3, mu3, sigma3, A4, mu4, sigma4, a, b):
return (A1 * np.exp(-(x - mu1)**2 / (2 * sigma1**2)) +
A2 * np.exp(-(x - mu2)**2 / (2 * sigma2**2)) +
A3 * np.exp(-(x - mu3)**2 / (2 * sigma3**2)) +
A4 * np.exp(-(x - mu4)**2 / (2 * sigma4**2)) +
a*x + b)
def six_gaussian_with_background(x, A1, mu1, sigma1, A2, mu2, sigma2, A3, mu3, sigma3, A4, mu4, sigma4, A5, mu5, sigma5, A6, mu6, sigma6, a, b):
return (A1 * np.exp(-(x - mu1)**2 / (2 * sigma1**2)) +
A2 * np.exp(-(x - mu2)**2 / (2 * sigma2**2)) +
A3 * np.exp(-(x - mu3)**2 / (2 * sigma3**2)) +
A4 * np.exp(-(x - mu4)**2 / (2 * sigma4**2)) +
A5 * np.exp(-(x - mu5)**2 / (2 * sigma5**2)) +
A6 * np.exp(-(x - mu6)**2 / (2 * sigma6**2)) +
a*x + b)
#FIXED PARAMETERS (E[keV] & sig[keV] for shifted and degraded 2+ -> 0+ transition peak)
fixed_953 = [953.78506621, 2.22422379]
fixed_964 = [963.81743382, 1.90161346]
#params_fixed_290 = [1231.8887296913715, 2.6925183648441835]
#params_fixed_10 = [1245.2793648, 2.51969264]
#fixed_1232 = [1231.896458563245, 2.764334269653125] #28.3.2024
#fixed_1245 = [1245.3540131724, 2.7225153082197826] #28.3.2024
fixed_1245 = [1245.3122002818593, 2.5730777214538882]
fixed_1232 = [1231.8787334405167, 2.638551806792974]
#2+ -> 0+ transition
#10 um
initial_guess_four_f_10 = [50, 1222, 1, 50, 1232, 3, 150, 1245, 3, 5, 1252.5, 1.2, -0.145, 110]
bnd_four_f_10 = ((0, 1220, 0, 0, 1231, 1, 0, 1242, 1, 0, 1252, 1, -0.1475, 100), (1000, 1225, 2.5, 1000, 1233, 5, 2000, 1246, 4, 150, 1256, 1.5, -0.1395, 232))
#initial_guess_four_f_10 = [50, 1208, 1, 40, 1220, 3, 30, 1233, 3, 200, 1245, 2, -0.067, 50]
#bnd_four_f_10 = ((10, 1207, 0, 5, 1218, 0, 0, 1230, 0, 20, 1243, 1, -0.07, 0),
# (300, 1209, 5, 150, 1225, 5, 1000, 1235, 5, 300, 1248, 4, -0.065, 140))
params_four_f_10, covariance_f_10 = curve_fit(four_gaussian_with_background, x_f_10, y_f_10, p0=initial_guess_four_f_10, bounds = bnd_four_f_10)
yfit_four_f_10 = four_gaussian_with_background(xdat_, *params_four_f_10)
fixed_1232_ = [1233.2732224475262, 1.4523473313953956]
params_four_fixed_both_f_10 = np.concatenate((params_four_f_10[0:4], fixed_1232_, [params_four_f_10[6]], fixed_1245, params_four_f_10[9:14]))
comb_1232_10 = np.concatenate(([params_four_f_10[3]],fixed_1232_, params_four_f_10[12:14]))
comb_1245_10 = np.concatenate(([params_four_f_10[6]],fixed_1245, params_four_f_10[12:14]))
plt.title('1232 keV (4+ -> 2+); 10 um')
plt.step(x_f_10, y_f_10, label='data (4+ -> 2+)')
plt.plot(xdat_1232, gaussian(xdat_1232, *comb_1232_10), 'r-', label='1232 keV')
plt.plot(xdat_1245, gaussian(xdat_1245, *comb_1245_10), 'k-', label='1245 keV')
plt.plot(xdat_, four_gaussian_with_background(xdat_, *params_four_fixed_both_f_10), 'r--')
plt.plot(xdat_, bkg(xdat_, *params_four_fixed_both_f_10[12:14]), 'k--')
plt.xlabel("E[keV]")
plt.show()
#4+ -> 2+
#10 um
I_1232_10 = integrate.simpson(gaussian(xdat_1232, *comb_1232_10), xdat_1232)
I_bkg_1232_10 = integrate.simpson(bkg(xdat_1232, *params_four_f_10[12:14]), xdat_1232)
I_net_1232_10 = I_1232_10 - I_bkg_1232_10
I_1245_10 = integrate.simpson(gaussian(xdat_1245, *comb_1245_10), xdat_1245)
I_bkg_1245_10 = integrate.simpson(bkg(xdat_1245, *params_four_f_10[12:14]), xdat_1245)
I_net_1245_10 = I_1245_10 - I_bkg_1245_10
b_10_f = I_net_1245_10/(I_net_1245_10 + I_net_1232_10)
unc_I_1232_10 = np.sqrt(I_1232_10)
unc_I_1245_10 = np.sqrt(I_1245_10)
unc_I_bkg_1232_10 = np.sqrt(I_bkg_1232_10)
unc_I_bkg_1245_10 = np.sqrt(I_bkg_1245_10)
unc_I_net_1232_10 = np.sqrt(unc_I_1232_10**2 + unc_I_bkg_1232_10**2)
unc_I_net_1245_10 = np.sqrt(unc_I_1245_10**2 + unc_I_bkg_1245_10**2)
#unc_b_10_f = b_10_f * np.sqrt((unc_I_net_1245_10 / I_net_1245_10)**2 + (unc_I_net_1232_10 / I_net_1232_10)**2)
unc_b_10 = b_10 * np.sqrt((unc_I_1232_10 / I_1232_10)**2 + (unc_I_1245_10 / I_1245_10)**2 + (unc_I_bkg_1232_10 / I_bkg_1232_10)**2 + (unc_I_bkg_1245_10 / I_bkg_1245_10)**2)
print('b_10_f =', b_10_f)
print('unc_b_10', unc_b_10)