Guide on flux and fluence measurements¶

Peak flux for an FRB¶

To calibrate the peak flux of an FRB, you need the gains file for the corresponding day. You also need to have installed the beam model library (here is a tutorial for the primary beam model) and be on a branch of baseband-analysis that contains the flux.py script. Below is an example of analysis steps to get the peak flux for event 172924963.

Python

import numpy as np
import matplotlib.pyplot as plt

from baseband_analysis.dev.flux import *
from baseband_analysis.core import *

# Load the data
data_path = "/data/chime/baseband/processed/2021/06/03/astro_172924963/singlebeam_172924963.h5"
data = BBData.from_file(data_path)
gains = read_gains("/data/chime/daily_gain_solutions/hdf5_files/gain_20210603T130623.102388Z_cyga_timing.h5")

Find the spectrum at the peak.

Python

spect_peak = get_frb_peak_spectrum(data, gains, thres=2.5, diagnostic_plots=True)

Load the primary beam model.

Python

import beam_model as bm

config = {
    "datapath_xpol" : "/usr/local/lib/python3.6/site-packages/beam_model/bm_data/beam_v1_SVD_XX.h5",
    "datapath_ypol" : "/usr/local/lib/python3.6/site-packages/beam_model/bm_data/beam_v1_SVD_YY.h5"
}

model = bm.primary.DataDrivenSVDPrimaryBeamModel(config)

Calibrate the data using the beam model.

Python

# Find the sensitivity at the FRB's position and frequencies
position = get_frb_position(data)
sens = model.get_sensitivity(position, spect_peak["f"])

# Calibrate
calib_x = spect_peak["spect_x"]/sens["sens_x"].squeeze()
calib_y = spect_peak["spect_y"]/sens["sens_y"].squeeze()

calibrated = (calib_x + calib_y)

Finally, the peak flux is just the average of the spectrum at the peak.

Python

# Calculate the peak flux
flux = np.nanmean(calibrated)
print(flux)

Fluence of an FRB¶

The fluence is the flux of the burst integrated over time. This calculation happens in three parts: first, we calibrate the spectra at every point in the main peak. Then, we get the average flux for each of these spectra (this gives a calibrated main peak in time). Finally, we take the integral of this peak over time. Using the same imports as above:

Python

spect = get_frb_spectra(data, gains)

Calibrate using the beam model. Note that we take the transpose of the sensitivities because the beam model returns sensitivities of shape (n positions, n frequencies), in this case (1, 1024), whereas the spectra are given in shape (n frequencies, n time points).

Python

# Find the sensitivity at the FRB's position and frequencies
position = get_frb_position(data)
sens = model.get_sensitivity(position, spect["f"])

# Calibrate
calib_x = spect["spect_x"]/sens["sens_x"].T
calib_y = spect["spect_y"]/sens["sens_y"].T

calibrated = (calib_x + calib_y)

Calculate the fluence by integrating under the pulse.

Python

# Calculate fluxes at every point in the burst
fluxes = np.nanmean(calibrated, axis=0)

# Calculate the fluence (Riemann sum)
downsampling_factor = spect["downsampling_factor"]
native_res = 2.56e-3 # in ms
t_res = native_res*downsampling_factor # time resolution of downsampled data

fluence = np.sum(fluxes*t_res)

print(fluence, "Jy ms")

Errors on the flux and fluence¶

The two main sources of error are from the data itself (noise) and from the primary beam model. From analyzing five observations of Cygnus A, it seems the errors that we find from the propagation are severely underestimated, but will likely improve as the new version of the primary beam model (the point source model) becomes available. Below is a process for getting the error on the peak flux and the fluence of an event.

Errors on the flux¶

Loading the data, same as above (this section also uses the same imports).

Python

# Load the data
data_path = "/data/chime/baseband/processed/2021/06/03/astro_172924963/singlebeam_172924963.h5"
data = BBData.from_file(data_path)
gains = read_gains("/data/chime/daily_gain_solutions/hdf5_files/gain_20210603T130623.102388Z_cyga_timing.h5")

# Find the clean spectrum at the peak
spect_peak = get_frb_peak_spectrum(data, gains, thres=2.5, diagnostic_plots=True)

Import the primary beam model with the error data.

Python

# Load the primary beam model
import beam_model as bm

config = {
    "datapath_xpol" : "/home/yuliya/CHIMEFRB/beam-model/beam_model/bm_data/beam_v1_SVD_XX.h5",
    "datapath_ypol" : "/home/yuliya/CHIMEFRB/beam-model/beam_model/bm_data/beam_v1_SVD_YY.h5",
    "datapath_xpol_resid" : "/home/yuliya/CHIMEFRB/beam-model/beam_model/bm_data/beam_v1_SVD_diff_XX.h5",
    "datapath_ypol_resid" : "/home/yuliya/CHIMEFRB/beam-model/beam_model/bm_data/beam_v1_SVD_diff_YY.h5"
}

model = bm.primary.DataDrivenSVDPrimaryBeamModel(config)

Then calibrate the spectrum, same as above.

Python

# Find the sensitivity at the FRB's position and frequencies
position = get_frb_position(data)
sens = model.get_sensitivity(position, spect_peak["f"])

# Calibrate
calib_x = spect_peak["spect_x"]/sens["sens_x"].squeeze()
calib_y = spect_peak["spect_y"]/sens["sens_y"].squeeze()

calibrated = (calib_x + calib_y)

# Calculate the peak flux
flux = np.nanmean(calibrated)
print(flux)

Find the error on the spectrum calibrated in flux using the propagation. From the expression of the spectrum calibration, where s_X and s_Y are the uncalibrated spectra in each polarization and m_X and m_Y are the beam model responses (sensitivities): \ spect

We can use basic error propagation to find the expression of the error on each point of the spectrum: \ spect_error

This formula is scripted below.

Python

# Find the error on the spectrum due to the data (first 2 terms of above formula)
err_data = (spect_peak["err_x"]/sens["sens_x"].squeeze())**2 + (spect_peak["err_y"]/sens["sens_y"].squeeze())**2

# Find the error on the spectrum due to the beam model
err_sens_x = model.get_error_on_sensitivity(position, spect_peak["f"])["err_x"].squeeze() # delta m_X
err_sens_y = model.get_error_on_sensitivity(position, spect_peak["f"])["err_y"].squeeze() # delta m_Y

err_model = (err_sens_x*spect_peak["spect_x"]/sens["sens_x"].squeeze()**2)**2 + (err_sens_y*spect_peak["spect_y"]/sens["sens_y"].squeeze()**2)**2

# Combine errors to find error on every spectrum point
d_f = np.sqrt(err_data + err_model)

# Note: Set either one of err_data or err_model to 0 to see how each influences the errors.
#    It is suspected that the errors on the beam model are underestimated or not used
#    properly, so more investigation is required.

We can now visualize the errors with a plot:

Python

plt.figure()
plt.title("Spectrum of FRB {} with errors from propagation".format(str(data.attrs["event_id"])))
plt.errorbar(x=spect_peak["f"], y=calibrated, yerr=d_f, fmt="k.", capsize=2, ecolor="xkcd:grey")
plt.xlabel("Frequency [MHz]")
plt.ylabel("Flux [Jy]")
plt.show()

FRB spectrum with (underestimated) error bars

Finally, find the error on the peak flux, again using the propagation. The peak flux is just the average of fluxes in the peak spectrum:\ peak_flux

Hence, the error on the peak flux can be found using the propagation. All the errors on the fluxes at individual frequencies were found in the previous step. N is the number of frequency channels.\ Peak flux error

Python

# Error on the peak flux
e = np.sqrt(np.nansum(np.square(d_f)))/len(d_f)
print(e)

We expect the output of the above to be 0.31711864919244825.

Errors on the fluence¶

Find the fluence (same as above).

Python

spect = get_frb_spectra(data, gains)

# Find the sensitivity at the FRB's position and frequencies
position = get_frb_position(data)
sens = model.get_sensitivity(position, spect["f"])

# Calibrate
calib_x = spect["spect_x"]/sens["sens_x"].T
calib_y = spect["spect_y"]/sens["sens_y"].T

calibrated = (calib_x + calib_y)

# Calculate fluxes at every point in the burst
fluxes = np.nanmean(calibrated, axis=0)

# Calculate the fluence (Riemann sum)
downsampling_factor = spect["downsampling_factor"]
native_res = 2.56e-3 # in ms

## Time resolution
t_res = native_res*downsampling_factor

fluence = np.sum(fluxes*t_res)

print(fluence, "Jy ms")

The following script allows you to find the error using error propagation. As above, this error is underestimated and should not be used to report results.

Python

## 1. Calculate the error on the spectrum, f(v, t)
# Find the error on the spectrum due to the data
err_data = (spect["err_x"]/sens["sens_x"])**2 + (spect["err_y"]/sens["sens_y"])**2

# Find the error on the spectrum due to the beam model
err_sens_x = model.get_error_on_sensitivity(position, spect["f"])["err_x"].T # delta m_X
err_sens_y = model.get_error_on_sensitivity(position, spect["f"])["err_y"].T # delta m_Y

err_model = (err_sens_x*spect["spect_x"]/sens["sens_x"].T**2)**2 + (err_sens_y*spect["spect_y"]/sens["sens_y"].T**2)**2

# Combine errors to get errors on each spectrum
d_spectra = np.sqrt(err_data.T + err_model)
print(d_spectra.shape) # should be: (n frequencies, n time points in the burst)

## 2. Calculate the error on the flux at each time point
d_f = np.sqrt(np.nansum(np.square(d_spectra), axis=0))/d_spectra.shape[0]
print(d_f.shape) # should be: (n time points)

## 3. Calculate the error on the fluence
d_F = t_res*np.sqrt(np.nansum(np.square(d_f)))
print(d_F) # Should be one value, the error on the fluence

Steady source spectrum calibration¶

Here is a worked example of calibrating five observations of Cygnus A from start to finish. First, load all files and imports.

Python

from baseband_analysis.core import BBData
from baseband_analysis.core.calibration import read_gains
from baseband_analysis.dev.flux import *

import numpy as np
import matplotlib.pyplot as plt

# files
data_paths = ["/data/chime/baseband/processed/2020/11/20/CygA_20201120000547/singlebeam_20201120000547_3beams.h5",
              "/data/chime/baseband/processed/2020/11/20/CygA_20201120235951/singlebeam_20201120235951_3beams.h5",
              "/data/chime/baseband/processed/2020/11/20/CygA_20201120000047/singlebeam_20201120000047_3beams.h5",
              "/data/chime/baseband/processed/2020/11/20/CygA_20201120235351/singlebeam_20201120235351_3beams.h5",
              "/data/chime/baseband/processed/2020/11/20/CygA_20201120235151/singlebeam_20201120235151_3beams.h5"
             ]

gains_paths = [
    "/data/chime/daily_gain_solutions/hdf5_files/gain_20201120T015901.218349Z_cyga_timing.h5",
]*len(data_paths)

data_files = []
gains_files = []
for n in range(len(data_paths)):
    data_files.append(BBData.from_file(data_paths[n]))
    gains_files.append(read_gains(gains_paths[n]))

Then, get the spectra. Set the plot argument to True to see all the diagnostic plots.

Python

spectra = []
for n in range(len(data_paths)):
    spectra.append(get_steady_source_spectrum(data_files[n], gains_files[n], remove_final=2, sep_pol=True, thres=2))

Find the positions of each observation of Cygnus A. As a sanity check, we expect all their y-positions to be the same. Note: get_frb_position can be used for steady sources as well as one-off events.

Python

positions = []
for n in range(len(data_files)):
    positions.append(get_frb_position(data_files[n]))

Load the primary beam model for calibration.

Python

import beam_model as bm

config = {
    "datapath_xpol" : "/home/yuliya/CHIMEFRB/beam-model/beam_model/bm_data/beam_v1_SVD_XX.h5",
    "datapath_ypol" : "/home/yuliya/CHIMEFRB/beam-model/beam_model/bm_data/beam_v1_SVD_YY.h5",
    "datapath_xpol_resid" : "/home/yuliya/CHIMEFRB/beam-model/beam_model/bm_data/beam_v1_SVD_diff_XX.h5",
    "datapath_ypol_resid" : "/home/yuliya/CHIMEFRB/beam-model/beam_model/bm_data/beam_v1_SVD_diff_YY.h5"
}

model = bm.primary.DataDrivenSVDPrimaryBeamModel(config)

sens = []
for n in range(len(data_paths)):
    sens.append(model.get_sensitivity(positions[n], spectra[n]["f"]))

Finally, calibrate the spectra.

Python

calib_svd = []
for n in range(len(data_paths)):
    # Can use find_source_from_spectrum() to find which beam to use, in this case 1
    x = spectra[n]["spect_x"][:,1] / sens[n]["sens_x"].squeeze()
    y = spectra[n]["spect_y"][:,1] / sens[n]["sens_y"].squeeze()
    calib_svd.append(x+y)

Use the cyg_a_polynomial function from flux.py to get the three best-fit parameters for the expected spectrum of Cygnus A (which follows a second-degree polynomial in frequency). This will be useful for plotting.

Python

[a,b,c] = cyg_a_polynomial()

You can now generate a nice plot!

Python

fig, (ax1, ax2) = plt.subplots(2,1, gridspec_kw={"height_ratios": [4,1]}, figsize=(10,7), sharex=True)

ax1.set_title("Cygnus A Spectrum") #
ax1.set_ylabel("Flux (Jy)")
ax1.set_xlabel("Frequency (MHz)")

f = np.linspace(400, 800, 1024)
ax1.plot(f, a*f**2 + b*f + c, "k--", label="Expected")
for n in range(len(data_files)):
    ax1.plot(spect[n]["f"], calib_svd[n], ".", label="x = " + '{:.3f}'.format(positions[n][0]))
ax1.legend(bbox_to_anchor=(1.02, 1), loc='upper left')
ax1.set_ylim(0, 7000)

ax2.set_ylabel("Residuals")
ax2.axhline(y=0, color='k', linestyle='--')
for n in range(len(data_files)):
    ax2.plot(spect[n]["f"], calib_svd[n] - (a*spect[n]["f"]**2 + b*spect[n]["f"] +c), ".")

plt.show()

Five observations of Cygnus A, calibrated using the primary beam model