"""
This script is used for averaging raw data with several different methods
which are described in https://doi.org/10.1063/1.4921479 and the document
we provided regarding the averaging of raw data (or master thesis).
**The experiment which was used:**
* Broadband Pump-Probe
* 1 x Chopper
* No Wobbler
* 1 pump interferometer arm blocked
"""
# Import tools -------------------------
import numpy as np
from numpy import ndarray
from matplotlib import pyplot as plt
import logging
logger = logging.getLogger(__name__)
import os
# Functions used in all algorithms --------------
[docs]def digitize_chopper(chopper_voltage: ndarray) -> ndarray:
# find the voltage that sits in the middle of the high and the low
# state of the chopper
separation_level = np.array(
[750]
) # (chopper_voltage.max() - chopper_voltage.min() ) / 2
logger.debug("separation level: {}".format(separation_level))
# Digitize chopper voltages generate array with 0 (chopper low) and
# 1 (chopper high) seperating the values in chopper_voltage at
# separation_level
chopper_states = np.digitize(chopper_voltage, separation_level)
return chopper_states
# --------------------------------------------------------
# Method 1 (m1) ------------------------
# Corresponds to equation 21 in https://doi.org/10.1063/1.4921479
# Currently corresponds to method 3.1 (Absorption averaged data) of our
# document
[docs]def evaluate_data_m1(data: ndarray) -> tuple:
"""
Absorption averaged data: First calculate shot to
shot difference spectrum and then take average and
variance of this.
Calculate average of raw data as described in
equation 21 of https://doi.org/10.1063/1.4921479.
This currently corresponds to method described in
section 3.1 of our document.
Args:
data (ndarray): raw *linearised* adc data. shape: (144, samples)
Returns:
ndarray: mean of shot to shot difference spectra. shape: (64)
ndarray: 1/variance of shot to shot difference spectra. shape: (64)
"""
chopper_states = digitize_chopper(
data[138]
) # this row contains to the chopper voltage
# calculate shot to shot (sts) transmission dividing the first 64
# pixels intensity by the last 64 pixels intensities
sts_transmission = data[:64] / data[64:128]
# calculate shot to shot (sts) difference (diff) spectra
sts_diff_spectra = -np.log10(sts_transmission[:, ::2]) + np.log10(
sts_transmission[:, 1::2]
)
# Invert sign if first chopper state was high
if chopper_states[0] == 1:
sts_diff_spectra *= -1
# calculate mean and variance of shot to shot difference spectra
mean_difference_spectrum = np.average(sts_diff_spectra, axis=1)
inv_variance = 1 / np.var(sts_diff_spectra, axis=1) # inv = inverse
return mean_difference_spectrum, inv_variance
# --------------------------------------------------------
# Method 2 (m2) ------------------------
# Corresponds to equation 22 in https://doi.org/10.1063/1.4921479
# Currently corresponds to method 3.2 (Transmission averaged data) of
# our document
[docs]def evaluate_data_m2(data: ndarray) -> tuple:
"""
Transmission averaged data: First calculate shot to
shot transmission and then take variance and average
of this. Afterwards calculate the difference spectrum
from this mean transmission.
Calculate average of raw data as described in
equation 22 of https://doi.org/10.1063/1.4921479.
This currently corresponds to method described in
section 3.2 of our document.
Args:
data (ndarray): raw *linearised* adc data. shape: (144, samples)
Returns:
ndarray: mean of shot to shot transmission. shape: (64, 2)
2 because of two different chopper states
ndarray: 1/variance of shot to shot difference spectra. shape: (64, 2)
2 because of two different chopper states
ndarray: difference spectrum calculated from mean transmission for
each chopper state. shape: (64)
"""
chopper_states = digitize_chopper(
data[138]
) # this row contains to the chopper voltage
# calculate shot to shot (sts) transmission dividing the first 64
# pixels intensity by the last 64 pixels intensities
sts_transmission = data[:64] / data[64:128]
# calculate average and variance of transmission for each chopper
# state
mean_transmission, inv_variance, _, _ = dp.sort_data(
sts_transmission, chopper_states, np.array(2)
)
# Calculate difference spectrum
difference_spectrum = -np.log10(mean_transmission[:, 0]) + np.log10(
mean_transmission[:, 1]
)
return mean_transmission, inv_variance, difference_spectrum
# --------------------------------------------------------
# Method 3 (m3) ------------------------ Corresponds to equation 23 in
# https://doi.org/10.1063/1.4921479 Currently corresponds to method 3.3
# (Intensity averaged data) of our document
[docs]def evaluate_data_m3(data: ndarray) -> tuple:
"""
Intensity averaged data: Calculate average and variance directly of
intensities. Afterwards calculate the difference spectrum from this
mean intensity for each each pixel and chopper state.
Calculate average of raw data as described in equation 23 of
https://doi.org/10.1063/1.4921479. This currently corresponds to
method described in section 3.3 of our document.
Args:
data (ndarray): raw *linearised* adc data. shape: (144, samples)
Returns:
ndarray: mean of intensity. shape: (128, 2)
2 because of two different chopper states
ndarray: 1/variance of intensity. shape: (128, 2)
2 because of two different chopper states
ndarray: difference spectrum calculated from mean intensity for
each chopper state and pixel. shape: (64)
"""
chopper_states = digitize_chopper(
data[138]
) # this row contains to the chopper voltage
# calculate average and variance of intensites for each chopper
# state
mean_intensity, inv_variance, _, _ = dp.sort_data(
data[:128], chopper_states, np.array(2)
)
# Calculate difference spectrum
transmission = mean_intensity[:64] / mean_intensity[64:128]
difference_spectrum = -np.log10(transmission[:, 0]) + np.log10(transmission[:, 1])
return mean_intensity, inv_variance, difference_spectrum
[docs]def process_data(name, path):
curr_dir = os.getcwd()
os.chdir(path)
# Yes I hardcoded it... find your own path
prl = dp.PixelResponseLinearization(
"/Users/arthun/Documents/Uni/Masterarbeit/Messsoftware/akb_software/hardware_config_files/h-lab_averaged_pixel_linearization_fit_parameters.json"
)
# Make directories to save processed data into
for i in range(1, 4):
path_name = "m{}".format(i)
if not os.path.exists(path_name):
os.mkdir(path_name)
# Make figure directories
if not os.path.exists("./figures"):
os.mkdir("./figures")
os.mkdir("./figures/histograms")
os.mkdir("./figures/difference_spectra")
for file in os.listdir():
if name in file:
raw_data = np.loadtxt(file, skiprows=1)
# Linearise pixel data
raw_data[:128] = prl.linearize(raw_data[:128])
# Visualize data now (histograms)
nrows = 43
ncols = 3
pixel_histograms, hist_ax = plt.subplots(
nrows=nrows, ncols=ncols, figsize=(8.3, 11.7 * 11)
)
pixel_histograms.suptitle(
"Histograms of linearized intensity for {}".format(file), y=0.995
)
for idx, data in enumerate(raw_data[:128]):
# Convert 1d index to 2d index
ax_idx = np.unravel_index(idx, shape=(nrows, ncols))
hist_ax[ax_idx].hist(data, bins=50)
hist_ax[ax_idx].set_xlabel("intensity [a.u.]")
hist_ax[ax_idx].set_ylabel("occurences")
hist_ax[ax_idx].grid()
hist_ax[ax_idx].set_title("Pixel {}".format(idx))
plt.tight_layout()
plt.subplots_adjust(top=0.99)
plt.savefig("./figures/histograms/{}.pdf".format(file))
plt.close(pixel_histograms)
# Calculate means by different methods
mean_m1, weight_m1 = evaluate_data_m1(raw_data)
mean_m2, weight_m2, diff_spec_m2 = evaluate_data_m2(raw_data)
mean_m3, weight_m3, diff_spec_m3 = evaluate_data_m3(raw_data)
# Save data - this would correspond to the resulting temp files
np.save("./m1/" + file + "_mean_m1", mean_m1)
np.save("./m1/" + file + "_weight_m1", weight_m1)
np.save("./m2/" + file + "_mean_m2", mean_m2)
np.save("./m2/" + file + "_weight_m2", weight_m2)
np.save("./m2/" + file + "_difference_spec_m2", diff_spec_m2)
np.save("./m3/" + file + "_mean_m3", mean_m3)
np.save("./m3/" + file + "_weight_m3", weight_m3)
np.save("./m3/" + file + "_difference_spec_m3", diff_spec_m3)
# Plot difference spectra and difference of difference spectra
spectra_fig, spec_ax = plt.subplots(nrows=4, figsize=(8.3, 11.7))
spectra_fig.suptitle(
"Difference spectra for {} calcluated with different methods".format(
file
)
)
spec_ax[0].plot(mean_m1, label="Method 1 / Equation 21")
spec_ax[0].plot(diff_spec_m2, label="Method 2 / Equation 22")
spec_ax[0].plot(diff_spec_m3, label="Method 3 / Equation 23")
spec_ax[1].plot(mean_m1 - mean_m1, label="Method 1 - Method 1")
spec_ax[1].plot(diff_spec_m2 - mean_m1, label="Method 2 - Method 1")
spec_ax[1].plot(diff_spec_m3 - mean_m1, label="Method 3 - Method 1")
spec_ax[2].plot(mean_m1 - diff_spec_m2, label="Method 1 - Method 2")
spec_ax[2].plot(diff_spec_m2 - diff_spec_m2, label="Method 2 - Method 2")
spec_ax[2].plot(diff_spec_m3 - diff_spec_m2, label="Method 3 - Method 2")
spec_ax[3].plot(mean_m1 - diff_spec_m3, label="Method 1 - Method 3")
spec_ax[3].plot(diff_spec_m2 - diff_spec_m3, label="Method 2 - Method 3")
spec_ax[3].plot(diff_spec_m3 - diff_spec_m3, label="Method 3 - Method 3")
[ax.grid() for ax in spec_ax]
[ax.legend() for ax in spec_ax]
[ax.set_xlabel("pixel") for ax in spec_ax]
[ax.set_ylabel("intensity [a.u.]") for ax in spec_ax]
plt.tight_layout()
plt.subplots_adjust(top=0.95)
plt.savefig("./figures/difference_spectra/{}.pdf".format(file))
plt.close(spectra_fig)
os.chdir(curr_dir)
if __name__ == "__main__":
# Change path to import self made data processing module
curr_dir = os.getcwd()
os.chdir("/Users/arthun/Documents/Uni/Masterarbeit/Messsoftware/akb_software")
import data_processing as dp
os.chdir(curr_dir)
# run it
process_data(
"200526_Probe_Big",
"/Users/arthun/Documents/Uni/Masterarbeit/Messsoftware/raw_data_20200527/probe_big",
)
process_data(
"200526_Probe_Only",
"/Users/arthun/Documents/Uni/Masterarbeit/Messsoftware/raw_data_20200527/probe_only",
)
process_data(
"200526_Probe_Small",
"/Users/arthun/Documents/Uni/Masterarbeit/Messsoftware/raw_data_20200527/probe_small",
)