Goal: Correct for singles rate with $W$ calculation

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}$

This requires calculating $S_i$ and $S_j$ from the cced files. I need to rewrite my analysis from the beginning, or write another function that parses the cced file.

In this file, I will import the singles and bicorr data and calculate all $D_{i,j}$, $S_i$, $S_j$, and $W_{i,j}$.


In [1]:
import os
import sys
import matplotlib.pyplot as plt
import numpy as np
import imageio

In [2]:
import pandas as pd

In [3]:
import seaborn as sns
sns.set(style='ticks')

In [4]:
sys.path.append('../scripts/')

In [5]:
import bicorr as bicorr
import bicorr_e as bicorr_e
import bicorr_plot as bicorr_plot
import bicorr_sums as bicorr_sums
import bicorr_math as bicorr_math

In [6]:
%load_ext autoreload
%autoreload 2

Load some data

I'm going to work with the data from the combined data sets. The analysis for this data set is in analysis\Cf072115_to_Cf072215b.

The one limitation here is that this data has already cut out the fission chamber neighbors.

det_df without fission chamber neighbors


In [7]:
det_df = bicorr.load_det_df('../meas_info/det_df_pairs_angles.csv')
pair_is = bicorr.generate_pair_is(det_df,ignore_fc_neighbors_flag=True)
det_df = det_df.loc[pair_is].reset_index().rename(columns={'index':'index_og'}).copy()
det_df.head()


Out[7]:
index_og d1 d2 d1d2 angle
0 44 2 3 203 15.0
1 45 2 4 204 30.0
2 46 2 5 205 45.0
3 47 2 6 206 60.0
4 48 2 7 207 75.0

I am going to add in a new optional input parameter in bicorr.load_det_df that will let you provide this det_df without fission chamber neighbors directly. Try it out.


In [7]:
det_df = bicorr.load_det_df('../meas_info/det_df_pairs_angles.csv', remove_fc_neighbors=True)
det_df.head()


Out[7]:
index_og d1 d2 d1d2 angle
0 44 2 3 203 15.0
1 45 2 4 204 30.0
2 46 2 5 205 45.0
3 47 2 6 206 60.0
4 48 2 7 207 75.0

In [8]:
chList, fcList, detList, num_dets, num_det_pairs = bicorr.build_ch_lists()
dict_pair_to_index, dict_index_to_pair, dict_pair_to_angle = bicorr.build_dict_det_pair(det_df)

In [9]:
num_fissions = 2194651200.00

singles_hist.npz


In [10]:
singles_hist, dt_bin_edges_sh, dict_det_to_index, dict_index_to_det = bicorr.load_singles_hist(filepath='../analysis/Cf072115_to_Cf072215b/datap',plot_flag=True,show_flag=True)


<Figure size 576x396 with 0 Axes>

Load bhp_nn for all pairs

I'm going to skip a few steps in order to save memory. This data was produced in analysis_build_bhp_nn_by_pair_1_ns.ipynb and is stored in datap\bhp_nn_by_pair_1ns.npz. Load it now, as explained in the notebook.


In [11]:
npzfile = np.load('../analysis/Cf072115_to_Cf072215b/datap/bhp_nn_by_pair_1ns.npz')
pair_is = npzfile['pair_is']
bhp_nn_pos = npzfile['bhp_nn_pos']
bhp_nn_neg = npzfile['bhp_nn_neg']
dt_bin_edges = npzfile['dt_bin_edges']

The fission chamber neighbors have already been removed


In [12]:
len(pair_is)


Out[12]:
861

Specify energy range


In [12]:
emin = 0.62
emax = 12

Calculate sums

Singles- set up singles_df

I will store this in a pandas dataframe.

Columns:

  • Channel number
  • Sp - Singles counts, positive
  • Sn - Singles counts, negative
  • Sd - Singles counts, br-subtracted
  • Sd_err - Singles counts, br-subtracted, err

In [15]:
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 [16]:
singles_df['Sp']= 0.0
singles_df['Sn']= 0.0
singles_df['Sd']= 0.0
singles_df['Sd_err'] = 0.0

Singles- calculate sums


In [17]:
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 [18]:
singles_df.head()


Out[18]:
ch Sp Sn Sd Sd_err
1 2 5150169.0 76920.0 5073249.0 2286.282791
2 3 5024494.0 73426.0 4951068.0 2257.857391
3 4 5391714.0 76594.0 5315120.0 2338.441361
4 5 5717925.0 81136.0 5636789.0 2408.123959
5 6 5448282.0 71482.0 5376800.0 2349.417800

In [19]:
bicorr_plot.Sd_vs_angle_all(singles_df)


<matplotlib.figure.Figure at 0xccde6a0>

Doubles- set up det_df


In [20]:
det_df.head()


Out[20]:
index_og d1 d2 d1d2 angle
0 44 2 3 203 15.0
1 45 2 4 204 30.0
2 46 2 5 205 45.0
3 47 2 6 206 60.0
4 48 2 7 207 75.0

In [21]:
det_df['Cp'] = 0.0
det_df['Cn'] = 0.0
det_df['Cd'] = 0.0
det_df['Cd_err'] = 0.0
det_df['Np'] = 0.0
det_df['Nn'] = 0.0
det_df['Nd'] = 0.0
det_df['Nd_err'] = 0.0

In [22]:
det_df.head()


Out[22]:
index_og d1 d2 d1d2 angle Cp Cn Cd Cd_err Np Nn Nd Nd_err
0 44 2 3 203 15.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
1 45 2 4 204 30.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2 46 2 5 205 45.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
3 47 2 6 206 60.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
4 48 2 7 207 75.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

Doubles- Calculate sums


In [23]:
for index in det_df.index.values:
    Cp, Cn, Cd, err_Cd = bicorr.calc_nn_sum_br(bhp_nn_pos[index,:,:],
                                               bhp_nn_neg[index,:,:],
                                               dt_bin_edges,
                                               emin=emin, emax=emax)
    det_df.loc[index,'Cp'] = Cp
    det_df.loc[index,'Cn'] = Cn
    det_df.loc[index,'Cd'] = Cd
    det_df.loc[index,'Cd_err'] = err_Cd
    Np, Nn, Nd, err_Nd = bicorr.calc_nn_sum_br(bhp_nn_pos[index,:,:],
                                               bhp_nn_neg[index,:,:],
                                               dt_bin_edges,
                                               emin=emin, emax=emax,
                                               norm_factor = num_fissions)
    det_df.loc[index,'Np'] = Np
    det_df.loc[index,'Nn'] = Nn
    det_df.loc[index,'Nd'] = Nd
    det_df.loc[index,'Nd_err'] = err_Nd

In [24]:
det_df.head()


Out[24]:
index_og d1 d2 d1d2 angle Cp Cn Cd Cd_err Np Nn Nd Nd_err
0 44 2 3 203 15.0 22338.0 229.0 22109.0 150.223167 22338.0 229.0 22109.0 0.003207
1 45 2 4 204 30.0 14779.0 129.0 14650.0 122.098321 14779.0 129.0 14650.0 0.002606
2 46 2 5 205 45.0 14042.0 116.0 13926.0 118.987394 14042.0 116.0 13926.0 0.002540
3 47 2 6 206 60.0 11993.0 121.0 11872.0 110.063618 11993.0 121.0 11872.0 0.002349
4 48 2 7 207 75.0 11707.0 104.0 11603.0 108.678425 11707.0 104.0 11603.0 0.002320

In [25]:
bicorr_plot.counts_vs_angle_all(det_df, normalized=True)


<matplotlib.figure.Figure at 0xc8a2278>
<matplotlib.figure.Figure at 0xbfbd630>
<matplotlib.figure.Figure at 0xcbb1940>
yes
<matplotlib.figure.Figure at 0xcb78fd0>

Perform the correction

Now I am going to loop through all pairs and calculate $W$.

  • Loop through each pair
  • Identify $i$, $j$
  • Fetch $S_i$, $S_j$
  • Calculate $W$
  • Propagate error for $W_{err}$
  • Store in det_df

Add W, W_err columns to det_df


In [26]:
det_df['Sd1']     = 0.0
det_df['Sd1_err'] = 0.0
det_df['Sd2']     = 0.0
det_df['Sd2_err'] = 0.0

det_df['W']     = 0.0
det_df['W_err'] = 0.0

In [27]:
det_df.head()


Out[27]:
index_og d1 d2 d1d2 angle Cp Cn Cd Cd_err Np Nn Nd Nd_err Sd1 Sd1_err Sd2 Sd2_err W W_err
0 44 2 3 203 15.0 22338.0 229.0 22109.0 150.223167 22338.0 229.0 22109.0 0.003207 0.0 0.0 0.0 0.0 0.0 0.0
1 45 2 4 204 30.0 14779.0 129.0 14650.0 122.098321 14779.0 129.0 14650.0 0.002606 0.0 0.0 0.0 0.0 0.0 0.0
2 46 2 5 205 45.0 14042.0 116.0 13926.0 118.987394 14042.0 116.0 13926.0 0.002540 0.0 0.0 0.0 0.0 0.0 0.0
3 47 2 6 206 60.0 11993.0 121.0 11872.0 110.063618 11993.0 121.0 11872.0 0.002349 0.0 0.0 0.0 0.0 0.0 0.0
4 48 2 7 207 75.0 11707.0 104.0 11603.0 108.678425 11707.0 104.0 11603.0 0.002320 0.0 0.0 0.0 0.0 0.0 0.0

In [28]:
singles_df.head()


Out[28]:
ch Sp Sn Sd Sd_err
1 2 5150169.0 76920.0 5073249.0 2286.282791
2 3 5024494.0 73426.0 4951068.0 2257.857391
3 4 5391714.0 76594.0 5315120.0 2338.441361
4 5 5717925.0 81136.0 5636789.0 2408.123959
5 6 5448282.0 71482.0 5376800.0 2349.417800

Loop through det_df, store singles rates

Fill the S and S_err values for each channel in each detector pair.


In [30]:
# Fill S columns in det_df
for index in singles_df.index.values:
    ch = singles_df.loc[index,'ch']
    
    d1_indices = (det_df[det_df['d1'] == ch]).index.tolist()
    d2_indices = (det_df[det_df['d2'] == ch]).index.tolist()
    
    det_df.loc[d1_indices,'Sd1']     = singles_df.loc[index,'Sd']
    det_df.loc[d1_indices,'Sd1_err'] = singles_df.loc[index,'Sd_err']
    det_df.loc[d2_indices,'Sd2']     = singles_df.loc[index,'Sd']
    det_df.loc[d2_indices,'Sd2_err'] = singles_df.loc[index,'Sd_err']

In [38]:
# Calculate W, W_err from S columns
det_df['W'] = det_df['Cd']/(det_df['Sd1']*det_df['Sd2'])

In [39]:
det_df['W_err'] = det_df['W'] * np.sqrt((det_df['Cd_err']/det_df['Cd'])**2 + 
                                        (det_df['Sd1_err']/det_df['Sd1'])**2 + 
                                        (det_df['Sd2_err']/det_df['Sd2'])**2)

In [40]:
det_df.head()


Out[40]:
index_og d1 d2 d1d2 angle Cp Cn Cd Cd_err Np Nn Nd Nd_err Sd1 Sd1_err Sd2 Sd2_err W W_err
0 44 2 3 203 15.0 22338.0 229.0 22109.0 150.223167 22338.0 229.0 22109.0 0.003207 5073249.0 2286.282791 4951068.0 2257.857391 8.802054e-10 6.007264e-12
1 45 2 4 204 30.0 14779.0 129.0 14650.0 122.098321 14779.0 129.0 14650.0 0.002606 5073249.0 2286.282791 5315120.0 2338.441361 5.432983e-10 4.540952e-12
2 46 2 5 205 45.0 14042.0 116.0 13926.0 118.987394 14042.0 116.0 13926.0 0.002540 5073249.0 2286.282791 5636789.0 2408.123959 4.869770e-10 4.171833e-12
3 47 2 6 206 60.0 11993.0 121.0 11872.0 110.063618 11993.0 121.0 11872.0 0.002349 5073249.0 2286.282791 5376800.0 2349.417800 4.352250e-10 4.044147e-12
4 48 2 7 207 75.0 11707.0 104.0 11603.0 108.678425 11707.0 104.0 11603.0 0.002320 5073249.0 2286.282791 5564633.0 2394.397419 4.110055e-10 3.858153e-12

In [44]:
bicorr_plot.W_vs_angle_all(det_df)


<matplotlib.figure.Figure at 0xbfbe1d0>

This is much "tighter" than the raw counts.

Functionalize

Write functions to perform all of these calculations. Demo them here. The functions are in a new script called bicorr_sums.py.

You have to specify emin, emax.


In [17]:
emin = 0.62
emax = 12

Data you have to have loaded:

  • det_df
  • dict_index_to_det
  • singles_hist
  • dt_bin_edges_sh
  • bhp_nn_pos
  • bhp_nn_neg
  • dt_bin_edges
  • emin
  • emax
  • num_fissions
  • angle_bin_edges

Produce and fill singles_df:


In [15]:
singles_df = bicorr_sums.init_singles_df(dict_index_to_det)
singles_df.head()


Out[15]:
ch Sp Sn Sd Sd_err
1 2 NaN NaN NaN NaN
2 3 NaN NaN NaN NaN
3 4 NaN NaN NaN NaN
4 5 NaN NaN NaN NaN
5 6 NaN NaN NaN NaN

In [18]:
singles_df = bicorr_sums.fill_singles_df(dict_index_to_det, singles_hist, dt_bin_edges_sh, emin, emax)
singles_df.head()


Out[18]:
ch Sp Sn Sd Sd_err
1 2 5214283.0 80614.0 5133669.0 2301.064319
2 3 5081297.0 77135.0 5004162.0 2271.218175
3 4 5453342.0 80374.0 5372968.0 2352.385173
4 5 5774814.0 85302.0 5689512.0 2420.767647
5 6 5497475.0 74890.0 5422585.0 2360.585732

In [19]:
bicorr_plot.Sd_vs_angle_all(singles_df)


<Figure size 576x396 with 0 Axes>

Expand, fill det_df


In [20]:
det_df.head()


Out[20]:
index_og d1 d2 d1d2 angle
0 44 2 3 203 15.0
1 45 2 4 204 30.0
2 46 2 5 205 45.0
3 47 2 6 206 60.0
4 48 2 7 207 75.0

In [23]:
det_df = bicorr_sums.init_det_df_sums(det_df, t_flag = True)
det_df = bicorr_sums.fill_det_df_singles_sums(det_df, singles_df)
det_df = bicorr_sums.fill_det_df_doubles_t_sums(det_df, bhp_nn_pos, bhp_nn_neg, dt_bin_edges, emin, emax)
det_df = bicorr_sums.calc_det_df_W(det_df)
det_df.head()


Out[23]:
index_og d1 d2 d1d2 angle Cp Cn Cd Cd_err Sd1 Sd1_err Sd2 Sd2_err W W_err
0 44 2 3 203 15.0 22602.0 233.0 22369.0 151.112541 5133669.0 2301.064319 5004162.0 2271.218175 8.707377e-10 5.908386e-12
1 45 2 4 204 30.0 15049.0 138.0 14911.0 123.235547 5133669.0 2301.064319 5372968.0 2352.385173 5.405858e-10 4.480623e-12
2 46 2 5 205 45.0 14381.0 128.0 14253.0 120.453310 5133669.0 2301.064319 5689512.0 2420.767647 4.879816e-10 4.134986e-12
3 47 2 6 206 60.0 12284.0 129.0 12155.0 111.413644 5133669.0 2301.064319 5422585.0 2360.585732 4.366372e-10 4.011538e-12
4 48 2 7 207 75.0 12019.0 113.0 11906.0 110.145359 5133669.0 2301.064319 5631831.0 2410.150410 4.118020e-10 3.818220e-12

Condense into angle bins


In [24]:
angle_bin_edges = np.arange(8,190,10)

In [26]:
by_angle_df = bicorr_sums.condense_det_df_by_angle(det_df,angle_bin_edges)
by_angle_df.head()


Out[26]:
angle_bin_centers angle_bin_max angle_bin_min len pair_is std_angle W W_err std W
0 13.0 18 8 40.0 0.439220 8.543146e-10 8.983238e-13 2.815420e-11
1 23.0 28 18 17.0 0.705868 6.372004e-10 1.202704e-12 1.330271e-11
2 33.0 38 28 84.0 2.934340 5.486453e-10 4.920738e-13 2.081072e-11
3 43.0 48 38 57.0 1.684873 4.917936e-10 5.634607e-13 1.257825e-11
4 53.0 58 48 53.0 2.414525 4.648758e-10 5.762423e-13 2.569397e-11

In [27]:
bicorr_plot.W_vs_angle(det_df, by_angle_df, save_flag=False)


<Figure size 576x396 with 0 Axes>

Put all of this into one function

Returns: singles_df, det_df, by_angle_df


In [14]:
angle_bin_edges = np.arange(8,190,10)

In [15]:
singles_df, det_df, by_angle_df = bicorr_sums.perform_W_calcs(det_df,
                    dict_index_to_det, singles_hist, dt_bin_edges_sh,
                    bhp_nn_pos, bhp_nn_neg, dt_bin_edges,
                    num_fissions, emin, emax, angle_bin_edges)

In [16]:
det_df.head()


Out[16]:
index_og d1 d2 d1d2 angle Cp Cn Cd Cd_err Sd1 Sd1_err Sd2 Sd2_err W W_err
0 44 2 3 203 15.0 22602.0 233.0 22369.0 151.112541 5133669.0 2301.064319 5004162.0 2271.218175 8.707377e-10 5.908386e-12
1 45 2 4 204 30.0 15049.0 138.0 14911.0 123.235547 5133669.0 2301.064319 5372968.0 2352.385173 5.405858e-10 4.480623e-12
2 46 2 5 205 45.0 14381.0 128.0 14253.0 120.453310 5133669.0 2301.064319 5689512.0 2420.767647 4.879816e-10 4.134986e-12
3 47 2 6 206 60.0 12284.0 129.0 12155.0 111.413644 5133669.0 2301.064319 5422585.0 2360.585732 4.366372e-10 4.011538e-12
4 48 2 7 207 75.0 12019.0 113.0 11906.0 110.145359 5133669.0 2301.064319 5631831.0 2410.150410 4.118020e-10 3.818220e-12

In [17]:
bicorr_plot.W_vs_angle(det_df, by_angle_df, save_flag = False)


<Figure size 576x396 with 0 Axes>

In [ ]: