In [1]:
%matplotlib notebook
In [2]:
import os
import sys
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import (MultipleLocator, FormatStrFormatter, AutoMinorLocator)
import pandas as pd
import scipy.io as sio
In [3]:
os.getcwd()
Out[3]:
In [4]:
os.chdir('..')
os.getcwd()
Out[4]:
In [5]:
sys.path.append('../scripts')
In [6]:
import bicorr_sim as bicorr_sim
import bicorr_plot as bicorr_plot
import bicorr_math as bicorr_math
In [7]:
by_angle_e_meas = pd.read_csv(r'Cf072115_to_Cf072215b/datap/by_angle_e_df.csv',index_col=0)
by_angle_e_cgmf = pd.read_csv(r'cgmf/datap/by_angle_e_df.csv',index_col=0)
by_angle_e_freya = pd.read_csv(r'freya/datap/by_angle_e_df.csv',index_col=0)
by_angle_e_ipol = pd.read_csv(r'ipol/datap/by_angle_e_df.csv',index_col=0)
by_angle_e_ipol_noct = pd.read_csv(r'ipol_noct/datap/by_angle_e_df.csv',index_col=0)
Load nps
In [8]:
num_fission_meas = int(int(sio.loadmat('Cf072115_to_Cf072215b/datap/num_fissions.mat')['num_fissions'])*float(sio.loadmat('Cf072115_to_Cf072215b/datap/fc_efficiency.mat')['fc_efficiency']))
num_fission_cgmf = int(sio.loadmat('cgmf/datap/num_fissions.mat')['num_fissions'])
num_fission_freya= int(sio.loadmat('freya/datap/num_fissions.mat')['num_fissions'])
num_fission_ipol = int(sio.loadmat('ipol/datap/num_fissions.mat')['num_fissions'])
num_fission_ipol_noct = int(sio.loadmat('ipol_noct/datap/num_fissions.mat')['num_fissions'])
num_fissions = [num_fission_meas,
num_fission_cgmf,
num_fission_freya,
num_fission_ipol,
num_fission_ipol_noct]
print(num_fissions)
In [21]:
by_angle_es = [by_angle_e_meas,
by_angle_e_cgmf,
by_angle_e_freya,
by_angle_e_ipol,
by_angle_e_ipol_noct]
legends =['Experiment', 'CGMF', 'FREYA', 'PoliMi', 'PoliMi-No CT']
fmts = ['x', 's', 'D', 'o', '^']
colors = ['#5d269b', '#dd673b', '#80bc31', '#3cbfe0', '#4242f4']
to_plot = [0,1, 2, 3]
In [22]:
line_thickness = 1
ebar_width = 3
In [23]:
fig = plt.figure(figsize=(4,4))
ax = plt.gca()
for i in to_plot:
by_angle_df = by_angle_es[i]
nps = num_fissions[i]
x = by_angle_df['angle_bin_centers']
y = by_angle_df['W']
yerr = by_angle_df['std W']
norm_factor = np.sum(y[x>20])
print(norm_factor)
y = y/norm_factor
yerr = yerr/norm_factor
plt.errorbar(x, y, yerr=yerr,
fmt = fmts[i],
markeredgewidth=1,
markerfacecolor='none',
elinewidth = line_thickness,
capthick = line_thickness,
capsize = ebar_width,
c = colors[i])
leg = plt.legend([legends[i] for i in to_plot])
leg.get_frame().set_edgecolor('w')
ax.axvspan(0,20,facecolor='gray', alpha=0.2)
ax.set_xlabel(r'$\theta$ (degrees)')
ax.set_ylabel(r'$\overline{W}(\theta)$ (arb. units)')
ax.set_xlim([0,180])
# Set up ticks
ax.tick_params(axis='both',
which='major',
direction='inout',
length=6,
color='k',
bottom=True, right=True, top=True, left=True)
ax.tick_params(axis='both',
which='minor',
direction='in',
length=3,
bottom=True, right=True, top=True, left=True)
# Major
ax.xaxis.set_major_locator(MultipleLocator(45))
#ax.yaxis.set_major_locator(MultipleLocator(0.02))
# Minor
ax.xaxis.set_minor_locator(MultipleLocator(15))
#ax.yaxis.set_minor_locator(MultipleLocator(0.005))
ax.text(45,0.12,'(a)', size=15, backgroundcolor='white')
plt.tight_layout()
In [19]:
os.getcwd()
Out[19]:
In [24]:
bicorr_plot.save_fig_to_folder('W_normd_by_integral',r'compare\fig')
In [25]:
def plot_calcs(by_angle_df):
x = by_angle_df['angle_bin_centers']
W = by_angle_df['W']
stdW = by_angle_df['std W']
norm_factor = np.sum(W[x>20])
y = W/norm_factor
yerr = stdW/norm_factor
return x, y, yerr
In [27]:
by_angle_df_exp = by_angle_es[0]
to_plot = [1,2,3]
fig = plt.figure(figsize=(4,4))
ax = plt.gca()
x_exp, y_exp, yerr_exp = plot_calcs(by_angle_df_exp)
for i in to_plot:
by_angle_df = by_angle_es[i]
x_sim, y_sim, yerr_sim = plot_calcs(by_angle_df)
y, yerr = bicorr_math.prop_err_division(y_sim,yerr_sim,y_exp,yerr_exp)
plt.errorbar(x, y, yerr=yerr,
fmt = fmts[i],
markeredgewidth=1,
markerfacecolor='none',
elinewidth = line_thickness,
capthick = line_thickness,
capsize = ebar_width,
c = colors[i])
leg = plt.legend([legends[i] for i in to_plot])
leg.get_frame().set_edgecolor('w')
plt.axhline(1.0,color='gray', linewidth=1,linestyle='--')
ax.axvspan(0,20,facecolor='gray', alpha=0.2)
ax.set_xlabel(r'$\theta$ (degrees)')
ax.set_ylabel(r'$\left[\overline{W}(\theta)\right]_{SIM} / \left[\overline{W}(\theta)\right]_{EXP}$')
ax.set_xlim([0,180])
# Set up ticks
ax.tick_params(axis='both',
which='major',
direction='inout',
length=6,
color='k',
bottom=True, right=True, top=True, left=True)
ax.tick_params(axis='both',
which='minor',
direction='in',
length=3,
bottom=True, right=True, top=True, left=True)
# Major
ax.xaxis.set_major_locator(MultipleLocator(45))
#ax.yaxis.set_major_locator(MultipleLocator(0.02))
# Minor
ax.xaxis.set_minor_locator(MultipleLocator(15))
#ax.yaxis.set_minor_locator(MultipleLocator(0.005))
ax.text(45,1.43,'(b)', size=15, backgroundcolor='white')
plt.tight_layout()
bicorr_plot.save_fig_to_folder('W_normd_diff',r'compare\fig')
In [29]:
In [ ]:
In [ ]:
In [15]:
by_angle_e_meas.head()
Out[15]:
In [16]:
by_angle_e_cgmf.head()
Out[16]:
In [17]:
by_angle_e_ipol.head()
Out[17]:
In [27]:
plt.figure(figsize=(4,4))
for i in to_plot:
by_angle_df = by_angle_es[i]
nps = num_fissions[i]
print(nps)
x = by_angle_df['angle_bin_centers']
# y = by_angle_df['W']*nps
y = by_angle_df['W']*np.sqrt(nps)
yerr = 0
plt.errorbar(x, y, yerr=yerr,
fmt = fmts[i],
markeredgewidth=1,
markerfacecolor='none',
elinewidth = line_thickness,
capthick = line_thickness,
capsize = ebar_width,
c = colors[i])
plt.legend([legends[i] for i in to_plot])
plt.xlabel(r'$\theta$ (degrees)')
plt.ylabel(r'$\overline{W}(\theta)$ (arb. units)')
plt.tight_layout()
In [13]:
by_angle_df.head()
Out[13]:
In [39]:
plt.figure(figsize=(4,4))
for i in to_plot:
by_angle_df = by_angle_es[i]
nps = num_fissions[i]
x = by_angle_df['angle_bin_centers']
y = np.multiply(by_angle_df['Sd1'], by_angle_df['Sd2'])/nps
yerr = 0
plt.errorbar(x, y, yerr=yerr,
fmt = fmts[i],
markeredgewidth=1,
markerfacecolor='none',
elinewidth = line_thickness,
capthick = line_thickness,
capsize = ebar_width,
c = colors[i])
plt.legend([legends[i] for i in to_plot])
plt.xlabel(r'$\theta$ (degrees)')
# plt.ylabel(r'$\overline{W}(\theta)$ (arb. units)')
plt.tight_layout()
In [ ]: