In [1]:
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')
The goal of this work is to design tools to:
bicorr_hist_master.npz file that contains the bicorr_hist_master numpy array, dt_bin_edges, and dict_pair_to_index, I will be building these tools using the bicorr_hist_master.npz data in 2017_03_22_pfs_combined_arrays_analysis. This data has 1 ns bins, which is nice to work with because the data will load and plot much faster than the 0.25 ns bins.
At the time of this work, the data includes the measurements Cf072115 and Cf072215a, which is about ten times as much data as in Cf072115 alone.
In [2]:
import numpy as np
import scipy.io as sio
import os
import sys
import matplotlib.pyplot as plt
import matplotlib.colors
import inspect
In [3]:
# Load the bicorr.py functions I have already developed
sys.path.append('../scripts')
import bicorr as bicorr
In [4]:
%load_ext autoreload
%autoreload 2
In [5]:
help(bicorr.load_sparse_bhm)
I'm loading some combined data so I have enough to work with.
In [5]:
sparse_bhm, dt_bin_edges = bicorr.load_sparse_bhm('../datar/')
det_df = bicorr.load_det_df()
bhm = bicorr.revive_sparse_bhm(sparse_bhm, det_df, dt_bin_edges)
In [6]:
print(sparse_bhm.shape)
print(bhm.shape)
In [7]:
det_df.head()
Out[7]:
bicorr_hist_plot, bhpIn producing the bicorrelation plot, I need to first generate a two-dimensional array that contains only the subset of data I wish to plot. I am going to separate this from the function that plots that array. The two-dimensional array will be called bicorr_hist_plot.
The technique for constructing bicorr_hist_plot should be as general and versatile as possible. I will build a function that returns bicorr_hist_plot from bicorr_hist_master and several selection options.
In [15]:
np.sum(bhm,axis=(0,1)).shape
Out[15]:
In [17]:
bhp = np.sum(bhm,axis=(0,1))
In [18]:
# Dimensions of histogram for a single event type
bhm[:,[1],:,:].shape
Out[18]:
Note: Here I use this bracket technique,
bicorr_hist_master[:,[1],:,:],
rather than the following alternative without brackets,
bicorr_hist_master[:,1,:,:],
because the bracket maintains the four-dimensional shape. If I did not use the bracket, the shape would change to (990,250,250). Then I would have to change the sum technique to only sum over axis 0. By using the bracket technique, I can keep my summation more general, always summing over axes 0 and 1.
In [19]:
# Now perform the sum
np.sum(bhm[:,[1],:,:],axis=(0,1)).shape
Out[19]:
In [20]:
bhp = np.sum(bhm[:,[1],:,:],axis=(0,1))
I am going to produce many plots with different subsets of the data. In order to keep it organized, I will use a naming convention that encodes what event types are in that distribution.
I will have four characters to code the interaction types that are being plotted:
allt at the end of the image filename.For instance, a plot of type 0 events only will be encoded as 0xxx. A plot of type 0, 1, and 2 events will be encoded as 012x.
The subset of data specified in the previous section would be encoded with the image filename extension _x1xx to make it clear that only interaction type 1 is included in that distribution.
In [33]:
dict_pair_to_index, dict_index_to_pair = bicorr.build_dict_det_pair(det_df)
In [35]:
det1 = 30
det2 = 31
pair_is = dict_pair_to_index[100*det1+det2]
np.sum(bhm[[pair_is],:,:,:],axis=(0,1)).shape
Out[35]:
In [37]:
bhp = np.sum(bhm[[pair_is],:,:,:],axis=(0,1))
Use the function I built in the analysis here:
https://github.com/pfschus/fission_bicorrelation/blob/master/analysis/detector_pair_angles.ipynb
In [6]:
print(inspect.getsource(bicorr.generate_pair_is))
In order to use this function, I need to provide the dictionary dict_pair_to_angle, which is built by the function bicorr.build_dict_det_pair.
In [7]:
help(bicorr.build_dict_det_pair)
In [42]:
dict_pair_to_angle = bicorr.build_dict_det_pair(det_df)[2]
Remind myself what the contents of this dictionary look like.
In [43]:
lists = sorted(dict_pair_to_angle.items()) # sorted by key, return a list of tuples
x, y = zip(*lists) # unpack a list of pairs into two tuples
plt.plot(x, y,'.k')
plt.xlabel('Detector pair in (100*det1ch+det2ch) notation')
plt.ylabel('Angle between detectors (degrees)')
plt.title('Contents of dict_pair_to_angles')
plt.show()
In [12]:
th_min = 0
th_max = 20
pair_is = bicorr.generate_pair_is(det_df,th_min,th_max)
This returned the detector pair corresponding indices in bicorr_hist_master.
In [57]:
pair_is
Out[57]:
In [61]:
plt.plot(pair_is,det_df.iloc[pair_is]['angle'],'.k')
plt.xlabel('Detector pair index')
plt.ylabel('Detector pair angle')
plt.title('Detector pairs in range {} to {} degrees'.format(th_min,th_max))
plt.show()
In [63]:
plt.plot(det_df.iloc[pair_is]['d1d2'],det_df.iloc[pair_is]['angle'],'.k')
plt.xlabel('Detector pair (d1d2 notation)')
plt.ylabel('Detector pair angle')
plt.title('Detector pairs in range {} to {} degrees'.format(th_min,th_max))
plt.show()
Looks good. Now select bicorr_hist_plot from this subset of detector pair indices pair_is.
In [65]:
bhp = np.sum(bhm[[pair_is],:,:,:],axis=(0,1))
bicorr_hist_plotI am going to write a function in my bicorr module that will generate bicorr_hist_plot from bicorr_hist_master and the specified subset of data requested by the user. The user should be able to specify the detector pairs as pair_is, in terms of the index of that pair in bicorr_hist_master, and the interaction types in terms of (0=nn, 1=np, 2=pn, 3=pp).
In order to normalize, I will follow the technique in my other notebook here:
https://github.com/pfschus/fission_bicorrelation/blob/master/analysis/coarsen_time_and_normalization.ipynb. The user will provide the number of fission events, and if the normalization flag norm_flag is set to True, then the distribution will be normalized by the number of fission events, the time bin width, and the number of detector pairs.
In [67]:
help(bicorr.build_bhm)
Take it for a test run. First, using the print_flag, without normalization:
In [8]:
bhp = bicorr.build_bhp(bhm,dt_bin_edges)[0]
In [9]:
bhp.shape
Out[9]:
In [10]:
np.sum(bhp)
Out[10]:
In [13]:
bhp = bicorr.build_bhp(bhm,dt_bin_edges, pair_is = pair_is)[0]
print(bhp.shape)
print(np.sum(bhp))
Now try it with normalization:
In [14]:
fission_rate = 2.9498e+05 # fissions/s
meas_length = 5*120 + 33*120
num_fissions = (fission_rate*meas_length)
In [15]:
bhp = bicorr.build_bhp(bhm, dt_bin_edges, num_fissions = num_fissions, print_flag=True)[0]
It looks like the process for constructing the bicorr_hist_plot array is complete. Now we need to plot it!
pcolormesh with logarithmic scaleThe technique for producing the 2D histogram, bicorr_hist_plot was covered in the previous section. Now I need to determine the best method for visualizing that histogram.
I will use the matplotlib function pcolormesh, which appears to be faster than the alternative pcolor.
Start by producing the desired bicorr_hist_plot. For now, look at all detector pairs and all interaction types.
Plot and compare log vs. linear scales for the colorbar.
In [15]:
bhp.shape
Out[15]:
In [17]:
plt.pcolormesh(dt_bin_edges, dt_bin_edges, bhp) # pcolormesh faster than pcolor
plt.colorbar()
plt.xlabel('$\Delta t_1$ (ns)')
plt.ylabel('$\Delta t_2$ (ns)')
plt.title('Bicorr plot for all pairs and all types, lin-scale')
plt.show()
plt.pcolormesh(dt_bin_edges, dt_bin_edges, bhp, norm=matplotlib.colors.LogNorm()) # pcolormesh faster than pcolor
plt.colorbar()
plt.xlabel('$\Delta t_1$ (ns)')
plt.ylabel('$\Delta t_2$ (ns)')
plt.title('Bicorr plot for all pairs and all types, log-scale')
plt.show()
For this data, it is important to use the logarithmic scale so that the non-gg features are visible.
In [20]:
plt.pcolormesh(dt_bin_edges, dt_bin_edges, bhp, norm=matplotlib.colors.LogNorm()
,vmin=1e-11,vmax=1e-7) # pcolormesh faster than pcolor
plt.colorbar()
plt.xlabel('$\Delta t_1$ (ns)')
plt.ylabel('$\Delta t_2$ (ns)')
plt.title('Bicorr plot for all pairs and all types, log-scale')
plt.show()
This "flattens" the colorbar by making anything outside of the range between vmin and vmax dark blue or dark red.
Using the versatility I built into bicorr.build_bicorr_hist_plot(), look at different conditions and subsets of data in bicorr_hist_plot.
First, look only at nn events.
In [22]:
bhp = bicorr.build_bhp(bhm, dt_bin_edges, type_is = [0], print_flag=True)[0]
plt.pcolormesh(dt_bin_edges, dt_bin_edges, bhp, norm=matplotlib.colors.LogNorm()) # pcolormesh faster than pcolor
plt.colorbar()
plt.xlabel('$\Delta t_1$ (ns)')
plt.ylabel('$\Delta t_2$ (ns)')
plt.title('Bicorr plot for all pairs and type 0 nn = 0xxx')
plt.show()
Now everything for a given pair.
In [16]:
dict_pair_to_index, dict_index_to_pair, dict_pair_to_angle = bicorr.build_dict_det_pair(det_df)
In [17]:
i = dict_pair_to_index[102]
bhp = bicorr.build_bhp(bhm, dt_bin_edges, pair_is = [i], print_flag=True)[0]
plt.pcolormesh(dt_bin_edges, dt_bin_edges, bhp, norm=matplotlib.colors.LogNorm()) # pcolormesh faster than pcolor
plt.colorbar()
plt.xlabel('$\Delta t_1$ (ns)')
plt.ylabel('$\Delta t_2$ (ns)')
plt.title('Bicorr plot for pair 102 (at index {}) and type 0 nn = 0xxx'.format(i))
plt.show()
Make a general bicorr_plot function that takes bicorr_hist_plot as an input parameter. View the contents here:
In [18]:
print(inspect.getsource(bicorr.bicorr_plot))
Try it out with lots of different options.
In [19]:
bhp = bicorr.build_bhp(bhm,dt_bin_edges,print_flag=True)[0]
In [20]:
bicorr.bicorr_plot(bhp,dt_bin_edges, show_flag = True)
In [21]:
bicorr.bicorr_plot(bhp,dt_bin_edges, show_flag = True, title='All events, all pairs'
, vmin = 1e2, vmax = 1e6)
Looks like it is working great.
There is a summary of my method for organizing detector pairs by their angles in my notebook here:
https://github.com/pfschus/fission_bicorrelation/blob/master/analysis/detector_pair_angles.ipynb
The steps in the process will be
1) Load detector pair angle information
2) Sort detector pairs by angle
3) Generate a bicorrelation plot at each angle, save to file
4) Stitch together in a .gif animation
In [22]:
dict_pair_to_index, dict_index_to_pair, dict_pair_to_angle = bicorr.build_dict_det_pair(det_df)
I also need to create vectors of the detector numbers and the angle between them, so I built the following function to unpack the dictionary into those vectors.
In [23]:
chList, fcList, detList, num_dets, num_det_pairs = bicorr.build_ch_lists()
I can see the angle between detector pairs by plotting det1ch, det2ch, and angle in a square scatter plot as shown below:
In [24]:
plt.scatter(det_df['d1'],det_df['d2'],c=det_df['angle'],s=13,marker='s',edgecolor='none')
plt.colorbar()
plt.xlim([0,50])
plt.ylim([0,50])
plt.xlabel('Detector 1')
plt.ylabel('Detector 2')
plt.title('Angle between detector pairs (degrees)')
plt.savefig('../fig/angle_btw_pairs.png')
plt.show()
What does the distribution of these detector angles look like?
In [25]:
plt.hist(det_df['angle'],20)
plt.xlabel('Angle between detector pairs (degrees)')
plt.ylabel('Number detector pairs')
plt.title('Distribution of detector pair angles')
plt.savefig('../fig/dist_angles.png')
plt.show()
In [26]:
th_bin_edges = np.arange(0,181,20)
print(th_bin_edges)
Use the function I wrote, bicorr.generate_pair_is_th_range, for generating pair_is given th_min and th_max (explained in
https://github.com/pfschus/fission_bicorrelation/blob/master/analysis/detector_pair_angles.ipynb).
In [27]:
help(bicorr.generate_pair_is)
In [41]:
th_min = th_bin_edges[0]
th_max = th_bin_edges[1]
pair_is = bicorr.generate_pair_is(det_df,th_min=th_min,th_max=th_max)
print(pair_is)
In [29]:
plt.plot(det_df.iloc[pair_is]['d1d2'],pair_is,'.k')
plt.xlabel('Pair number (100*det1ch+det2ch)')
plt.ylabel('Pair index')
plt.title('Pair indices within the given angular range')
plt.show()
In [30]:
plt.scatter(det_df.iloc[pair_is]['d1'],det_df.iloc[pair_is]['d2'],c=list(det_df.iloc[pair_is]['angle']),s=13,marker='s',edgecolor='none')
plt.colorbar()
plt.xlabel('Detector 1')
plt.ylabel('Detector 2')
plt.title('Angle between detector pairs (degrees)')
plt.show()
Loop through the angle ranges defined by th_bin_edges and produce a bicorrelation plot at each angle. For now, include all event types, and print the output.
The colorbar range will vary across all of the plots. I should generate bicorr_hist_plot for all angles and then calculate the overall max and min, then plot all bicorrelation plots with the same vmax and vmix.
Instead of making bicorr_hist_plot only two dimensions, I will add a third dimension for the angle range.
In [31]:
bhp = np.zeros([len(th_bin_edges)-1,len(dt_bin_edges)-1,len(dt_bin_edges)-1])
In [32]:
bhp.shape
Out[32]:
In [44]:
fission_rate = 2.9498e+05 # fissions/s
meas_length = 5*120 + 33*120
num_fissions = (fission_rate*meas_length)
for th_i in np.arange(0,len(th_bin_edges)-1):
th_min = th_bin_edges[th_i]
th_max = th_bin_edges[th_i+1]
pair_is = bicorr.generate_pair_is(det_df,th_min=th_min,th_max=th_max)
bhp[th_i,:,:] = bicorr.build_bhp(bhm, dt_bin_edges,num_fissions=num_fissions,pair_is = pair_is, print_flag=True)[0]
Now calculate vmin and vmax from all of the data in bicorr_hist_plot. There are many empty bins still with magnitude zero, so I need to calculate vmin as the minimum of the nonzero data.
In [45]:
vmin = np.min(bhp[np.nonzero(bhp)])
vmax = np.max(bhp)
print(vmin,vmax)
I am going to save the plots to .png files in a local folder, and then stitch them together in a .gif.
In [47]:
filenames = []
for th_i in np.arange(0,len(th_bin_edges)-1):
th_min = th_bin_edges[th_i]
th_max = th_bin_edges[th_i+1]
save_filename = r'bicorr_plot_{:03d}_{:03d}'.format(th_min,th_max)
filenames.append(save_filename)
bicorr.bicorr_plot(bhp[th_i,:], dt_bin_edges, title="All events for detector pairs from {} to {} degrees".format(th_min,th_max), vmin = vmin, vmax = vmax, save_flag = True, save_filename = save_filename, save_folder = '../fig/animate/', show_flag = True)
.gif animationTry the technique shown on this StackOverflow:
http://stackoverflow.com/questions/753190/programmatically-generate-video-or-animated-gif-in-python
I am essentially going to load each of the images I just saved and tack it on to the end of a .gif that I am building.
-Stitch-together-in-a-.gif-animation){kind=link}
In [48]:
print(filenames)
In [49]:
import imageio
help(imageio.mimsave)
In [50]:
images = []
for filename in filenames:
images.append(imageio.imread('../fig/animate/'+filename+'.png'))
imageio.mimsave('../fig/animate/all.gif', images, fps=1)
The .gif is now saved to file. Can I view it here in the jupyter notebook? Use a markdown cell and the following code:
<img src="fig/all.gif">
In [52]:
# Produce important variables for detector indices, pairs
det_df = bicorr.load_det_df()
dict_pair_to_index, dict_index_to_pair, dict_pair_to_angle = bicorr.build_dict_det_pair(det_df)
chList, fcList, detList, num_dets, num_det_pairs = bicorr.build_ch_lists()
In [53]:
# Specify angle bins
th_bin_edges = np.arange(0,181,20)
# Set up bicorr_hist_plot
bhp = np.zeros([len(th_bin_edges)-1,len(dt_bin_edges)-1,len(dt_bin_edges)-1])
In [57]:
# Calculate number of fissions based on measurement
fission_rate = 2.9498e+05 # fissions/s
meas_length = 5*120 + 33*120
num_fissions = (fission_rate*meas_length)
for th_i in np.arange(0,len(th_bin_edges)-1):
th_min = th_bin_edges[th_i]
th_max = th_bin_edges[th_i+1]
pair_is = bicorr.generate_pair_is(det_df,th_min=th_min,th_max=th_max)
bhp[th_i,:,:] = bicorr.build_bhp(bhm, dt_bin_edges,num_fissions=num_fissions,pair_is = pair_is, print_flag=False)[0]
# Calculate colorbar range
vmin = np.min(bhp[np.nonzero(bhp)])
vmax = np.max(bhp)
In [58]:
filenames = []
for th_i in np.arange(0,len(th_bin_edges)-1):
th_min = th_bin_edges[th_i]
th_max = th_bin_edges[th_i+1]
save_filename = r'bicorr_plot_{:03d}_{:03d}'.format(th_min,th_max)
filenames.append(save_filename)
bicorr.bicorr_plot(bhp[th_i,:], dt_bin_edges, title="All events for detector pairs from {} to {} degrees".format(th_min,th_max), vmin = vmin, vmax = vmax, save_flag = True, save_filename = save_filename, save_folder = '../fig/animate/', show_flag = True)
In [60]:
images = []
for filename in filenames:
images.append(imageio.imread('../fig/animate/'+filename+'.png'))
imageio.mimsave('../fig/animate/all.gif', images, fps=1)
<img src="fig/all.gif">
In [ ]: