Patricia Schuster
University of Michigan
January 11, 2018
I built the singles histogram in the other methods notebook, singles_histogram.ipynb. Now I am going to calculate the sum on that histogram.
Update: July 17, 2018. I am going to redo this method using the energy histograms.
The reason we care about the singles rates:
In order to correct for differences in detection efficiencies and solid angles, we will divide all of the doubles rates by the singles rates of the two detectors as follows:
$ W_{i,j} = \frac{D_{i,j}}{S_i*S_j}$
I can calculate $S_i$ and $S_j$ from the singles histogram. I will write a general function to calculate the histogram, and also demonstrate how to store those sums in a pandas dataframe.
In [1]:
import os
import sys
import matplotlib.pyplot as plt
import matplotlib.colors
import numpy as np
import os
import scipy.io as sio
import sys
import pandas as pd
from tqdm import *
# Plot entire array
np.set_printoptions(threshold=np.nan)
In [3]:
import seaborn as sns
sns.set_style(style='white')
In [89]:
sys.path.append('../scripts/')
import bicorr as bicorr
import bicorr_plot as bicorr_plot
import bicorr_e as bicorr_e
import bicorr_sums as bicorr_sums
In [90]:
%load_ext autoreload
%autoreload 2
In [7]:
det_df = bicorr.load_det_df()
chList, fcList, detList, num_dets, num_det_pairs = bicorr.build_ch_lists()
In [22]:
bicorr_e.load_singles_hist_both?
In [26]:
singles_hist_e_n, e_bin_edges, dict_det_to_index, dict_index_to_det = bicorr_e.load_singles_hist_both(filepath = r'../analysis/Cf072115_to_Cf072215b/datap', plot_flag=True, show_flag = True)
In [28]:
singles_hist_e_n.shape
Out[28]:
In [86]:
e_min = 0.62
e_max = 15
In [70]:
i_min = np.digitize(e_min,e_bin_edges)-1
i_max = np.digitize(e_max,e_bin_edges)-1
print(i_min, i_max)
In [65]:
e_bin_edges[i_max]
Out[65]:
In [49]:
e_bin_edges[np.digitize(e_max,e_bin_edges)-1]
Out[49]:
In [59]:
singles_hist_e_n[0,i_min:i_max]
Out[59]:
In [60]:
singles_hist_e_n[0,0:10]
Out[60]:
In [71]:
bicorr_sums.calc_n_sum_e(singles_hist_e_n, e_bin_edges, 0, e_min, e_max)
Out[71]:
In [75]:
singles_e_df = bicorr_sums.init_singles_e_df(dict_index_to_det)
singles_e_df.head()
Out[75]:
In [87]:
singles_e_df = bicorr_sums.fill_singles_e_df(dict_index_to_det, singles_hist_e_n, e_bin_edges, e_min, e_max)
singles_e_df.head()
Out[87]:
In [88]:
plt.errorbar(singles_e_df['ch'],
singles_e_df['Se'],
yerr=singles_e_df['Se_err'], fmt='.')
plt.xlabel('Detector channel')
plt.ylabel('Singles count rate')
plt.title('Singles count rate for each detector (errorbars shown)')
plt.show()
In [19]:
singles_hist, dt_bin_edges_sh, dict_det_to_index, dict_index_to_det = bicorr.load_singles_hist(filepath=r'../analysis/Cf072115_to_Cf072215b/datap/',
plot_flag = True, show_flag =True)
In [24]:
emin = 0.62
emax = 12
tmin = bicorr.convert_energy_to_time(emax)
tmax = bicorr.convert_energy_to_time(emin)
print(tmin,tmax)
In [15]:
dt_bin_edges_sh[-10:]
Out[15]:
In [17]:
i_min = np.min(np.argwhere(tmin<dt_bin_edges_sh))
i_max = np.min(np.argwhere(tmax<dt_bin_edges_sh))
print(i_min,i_max)
I also need to find the time ranges for the negative sum.
In [35]:
i_min_neg = np.min(np.argwhere(-tmax<dt_bin_edges_sh))
i_max_neg = np.min(np.argwhere(-tmin<dt_bin_edges_sh))
print(i_min_neg,i_max_neg)
In [40]:
dt_bin_centers = (dt_bin_edges_sh[:-1]+dt_bin_edges_sh[1:])/2
plt.plot(dt_bin_centers,np.sum(singles_hist[0,:,:],axis=(0)))
plt.plot(dt_bin_centers,np.sum(singles_hist[1,:,:],axis=(0)))
plt.axvline(dt_bin_edges_sh[i_min],color='r',linestyle='--')
plt.axvline(dt_bin_edges_sh[i_max],color='r',linestyle='--')
plt.axvline(dt_bin_edges_sh[i_min_neg],color='r',linestyle='--')
plt.axvline(dt_bin_edges_sh[i_max_neg],color='r',linestyle='--')
plt.xlabel('Time (ns)')
plt.ylabel('Number of events')
plt.title('Singles TOF distribution, all channels')
plt.legend(['N','G'])
plt.yscale('log')
plt.show()
For this calculation, I am only working with neutron events, so selected 0 in the first dimension of singles_hist.
For now, I will use all detector pairs. Start with the positive sum.
In [21]:
np.sum(singles_hist[[0],:,i_min:i_max])
Out[21]:
Now the negative sum
In [41]:
np.sum(singles_hist[[0],:,i_min_neg:i_max_neg])
Out[41]:
For now I am going to limit this to a single detector pair, since at the moment I only envision occasions where I have to calculate the sum for a specific pair and therefore the sums for each detector.
In [43]:
dict_det_to_index[4]
Out[43]:
In [44]:
ch_1 = 10
ch_2 = 11
print(np.sum(singles_hist[[0],dict_det_to_index[ch_1],i_min:i_max]))
print(np.sum(singles_hist[[0],dict_det_to_index[ch_2],i_min:i_max]))
In [45]:
help(bicorr.calc_n_sum_br)
In [46]:
bicorr.calc_n_sum_br(singles_hist, dt_bin_edges_sh, 1)
Out[46]:
Loop through all of the detectors and calculate it. Populate an array with these values
In [55]:
singles_rates = np.zeros((num_dets,2))
# Dimension 0: Detector index
# Dimension 1: 0- n_sum_br
# 1- n_sum_br_err
In [56]:
for det in detList:
det_i = dict_det_to_index[det]
singles_rates[det_i,:] = bicorr.calc_n_sum_br(singles_hist, dt_bin_edges_sh, det_i)
In [61]:
plt.errorbar(detList, singles_rates[:,0], yerr=singles_rates[:,1], fmt='.')
plt.xlabel('Detector channel')
plt.ylabel('Singles count rate')
plt.title('Singles count rate for each detector (errorbars shown)')
plt.show()
Columns:
Sp - Singles counts, positiveSn - Singles counts, negativeSd - Singles counts, br-subtractedSd_err - Singles counts, br-subtracted, err
In [20]:
singles_df = pd.DataFrame.from_dict(dict_index_to_det,orient='index',dtype=np.int8).rename(columns={0:'ch'})
chIgnore = [1,17,33]
singles_df = singles_df[~singles_df['ch'].isin(chIgnore)].copy()
In [21]:
singles_df['Sp']= 0.0
singles_df['Sn']= 0.0
singles_df['Sd']= 0.0
singles_df['Sd_err'] = 0.0
In [22]:
singles_df.head()
Out[22]:
In [25]:
for index in singles_df.index.values:
Sp, Sn, Sd, Sd_err = bicorr.calc_n_sum_br(singles_hist, dt_bin_edges_sh, index, emin=emin, emax=emax)
singles_df.loc[index,'Sp'] = Sp
singles_df.loc[index,'Sn'] = Sn
singles_df.loc[index,'Sd'] = Sd
singles_df.loc[index,'Sd_err'] = Sd_err
In [28]:
singles_df.head()
Out[28]:
Now use the plotting functions I developed.
In [26]:
bicorr_plot.Sd_vs_angle_all(singles_df)
I was going to save this to disk, but actually, the singles count rate depends on the timing windows so I need to provide the timing windows. Instead, I'll call the previous loop with the energy windows I desire.