Measured_Calibration_tests.py

import h5py
import matplotlib.pyplot as plt
import numpy as np
import sys
import re
import ast
from scipy.optimize import curve_fit
import os
from scipy import interpolate

# Solo levanto algunos experimentos
Calib_Files = """000002231-LaserPowerCalibration.h5""" 

Meas_Files = """000002239-LaserPowerCalibration.h5
000002241-LaserPowerCalibration.h5"""

Test_Files = """000002231-LaserPowerCalibration.h5
000002232-LaserPowerCalibration.h5
000002233-LaserPowerCalibration.h5
000002237-LaserPowerCalibration.h5""" 


#carpeta pc nico labo escritorio:
#C:\Users\Usuario\Documents\artiq\artiq_experiments\analisis\plots\20210716_MeasuredCalibrationsTests\Data

def SeeKeys(files):
    for i, fname in enumerate(files.split()):
        data = h5py.File(fname, 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
        print(fname)
        print(list(data['datasets'].keys()))


def GetXToInterpolate(vec):
    return np.arange(min(vec), max(vec), 1e-2*np.abs((vec[1]-vec[0])))

def Test_Calibrate_amplitudes(PD_UV_counts, Calibration_freqs, Calibration_amps, Calibration_PD_Value, Experiment_frequencies):
    
    Interpolation_functions = []
    Interpolation_functions_2 = []

    Calibration_freqs = Calibration_freqs[0:int(len(PD_UV_counts)/len(Calibration_amps))]

    n = len(Calibration_freqs)
    PD_split = np.matrix([PD_UV_counts[i * n:(i + 1) * n] for i in range((len(PD_UV_counts) + n - 1) // n )])
    PD_T = np.transpose(PD_split)

    amps = Calibration_amps

    for i in range(len(PD_T)):
        #PD_values = np.array(PD_T[i][0])[0]
        PD_values = np.array(PD_T[i])[0]
        f = interpolate.interp1d(PD_values, amps, kind='cubic')
        g = interpolate.interp1d(amps, PD_values, kind='cubic')

        Interpolation_functions.append(f)
        Interpolation_functions_2.append(g)
        

    Calibrated_Amplitudes = []
    for fi in Interpolation_functions:
        try:
            Calibrated_Amplitudes.append(float(fi(Calibration_PD_Value)))
        except:
            Calibrated_Amplitudes.append(0)
            print('no, rey')


    Function_Calibs_vs_freq = interpolate.interp1d(Calibration_freqs, Calibrated_Amplitudes, kind='quadratic')    
    
    Experiment_amps = []
    
    for freq in Experiment_frequencies:
        calibamp = Function_Calibs_vs_freq(freq)
        Experiment_amps.append(calibamp)
    
    fs = GetXToInterpolate(Calibration_freqs)
    plt.figure()
    plt.plot(Calibration_freqs, Calibrated_Amplitudes, 'o')
    plt.plot(fs, Function_Calibs_vs_freq(fs))
    
    ExpectedFluo = []
    i = 0
    while i < len(Experiment_amps):
        if Experiment_amps[i] == 0:
            ExpectedFluo.append(0)
        else:
            ExpectedFluo.append(Interpolation_functions_2[i](Experiment_amps[i]))
           
        i = i + 1
        
    
    return Experiment_amps, ExpectedFluo
        

#%%    


Calibration_amps_vec = []
Calibration_freqs_vec = []

Experiment_amps_vec = []
Experiment_freqs_vec = []

PD_UV_counts_vec = []

for i, fname in enumerate(Calib_Files.split()):
    print(SeeKeys(Calib_Files))
    print(i)
    print(fname)
    data = h5py.File(fname, 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
    print(list(data['datasets'].keys()))
    Calibration_amps_vec.append(np.array(data['datasets']['Calibration_amps']))
    Calibration_freqs_vec.append(np.array(data['datasets']['Calibration_freqs']))
    Experiment_amps_vec.append(np.array(data['datasets']['Experiment_amps']))
    Experiment_freqs_vec.append(np.array(data['datasets']['Experiment_freqs']))
    PD_UV_counts_vec.append(np.array(data['datasets']['PD_UV_counts']))
    


Final_Experiment_amps_vec = []
Final_Experiment_freqs_vec = []
Measured_PD_UV_counts_vec = []

for i, fname in enumerate(Meas_Files.split()):
    print(i)
    print(fname)
    data = h5py.File(fname, 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
    print(list(data['datasets'].keys()))
    Final_Experiment_amps_vec.append(np.array(data['datasets']['Experiment_amps']))
    Final_Experiment_freqs_vec.append(np.array(data['datasets']['Experiment_freqs']))
    Measured_PD_UV_counts_vec.append(np.array(data['datasets']['Measured_PD_UV_counts']))
    

#%%

i = 0

Calibration_amps= Calibration_amps_vec[i]
Calibration_freqs = Calibration_freqs_vec[i]
PD_UV_counts = PD_UV_counts_vec[i]
Measured_PD_UV_counts = Measured_PD_UV_counts_vec[i]


PDvalue = 0.04

Experiment_frequencies = Final_Experiment_freqs_vec[0]

Test_Experiment_frequencies = Calibration_freqs[0:len(Calibration_amps)]
#Test_Experiment_frequencies = np.arange(0.9e8, 1.3e8, 0.05e8)

calibamps, expectedfluos = Test_Calibrate_amplitudes(PD_UV_counts, Calibration_freqs, Calibration_amps, PDvalue, Experiment_frequencies)

plt.figure()
plt.plot(Experiment_frequencies, calibamps, 'ro', label='amps')
plt.plot(Experiment_frequencies, expectedfluos, 'gx', label='Fluos')
#plt.plot(Experiment_frequencies, Measured_PD_UV_counts, 'b.', label='Measures')
plt.axhline(PDvalue)
plt.legend()

#%%
len_amps = len(Calibration_amps)
len_freqs = int(len(Calibration_freqs)/len_amps)

freqs = Calibration_freqs[0:len_freqs]

n = len(freqs)
PD_split = np.matrix([PD_UV_counts[i * n:(i + 1) * n] for i in range((len(PD_UV_counts) + n - 1) // n )])
PD_T = np.transpose(PD_split)

plt.figure()
for i in range(len(PD_split)):
    plt.plot([f*1e-6 for f in freqs], np.array(PD_split[i])[0], 'o')
plt.xlabel('Frecuencia (MHz)')
plt.ylabel('Voltaje fotodiodo (V)')
plt.title('Láser UV')


i = 10
plt.figure()
plt.plot(np.array(PD_T[i][0])[0], Calibration_amps, 'o', label=f"Freq: {round(freqs[i]*1e-6, 3)} MHz")
plt.ylabel('Amplitud Artiq')
plt.xlabel("Voltaje fotodiodo (V)")
plt.legend()

#%%

def RemoveZeros(calib, freqs):
    i, j = 0, 1
    new_calib = calib
    while i < len(calib):
        if new_calib[0] == 0:
            new_calib = new_calib[1:]
            i = i + 1
        else:            
            init_len = len(new_calib)
            final_calib = new_calib
            while j < init_len:
                if final_calib[-1] == 0:
                    final_calib = final_calib[:-1]
                    j = j + 1
                else:
                    if j == 1:
                        return (i, j), freqs[i:], final_calib
                    else:
                        return (i, j), freqs[i:-j+1], final_calib
        
lista = [0, 1, 0, 1, 0, 0, 1, 1, 0, 0]
freqs = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
tupla, new_freqs, calib = RemoveZeros(lista, freqs)

print(new_freqs)
print(calib)

#%%
i = 15

interpolations = []

interpolations_g = []

for i in range(len(PD_T)):
    PD = np.array(PD_T[i][0])[0]
    amps = Calibration_amps
    f = interpolate.interp1d(PD, amps, kind='quadratic')
    g = interpolate.interp1d(amps, PD, kind='quadratic')
    xs = np.arange(min(PD), max(PD), 1e-4)
    
    interpolations.append(f)
    interpolations_g.append(g)
    

Calibration_PD_value = 0.2

Calibs = []
for fi in interpolations:
    try:
        Calibs.append(float(fi(Calibration_PD_value)))
    except:
        Calibs.append(0)


fig, ax1 = plt.subplots()
ax2 = ax1.twinx()

ax1.plot([f*1e-6 for f in freqs], Calibs, 'o', label='Amplitudes a setear')


EffectivePD = []
for i in range(len(Calibs)):
    try:
        EffectivePD.append(interpolations_g[i](Calibs[i]))
    except:
        EffectivePD.append(0)

ax2.plot([f*1e-6 for f in freqs], EffectivePD, 'ro', label='Potencia')

ax1.set_xlabel('Frecuencia (MHz)')
ax1.set_ylabel('Amplitudes')
ax2.set_ylabel('Potencia PD (V)')
fig.legend()


h = interpolate.interp1d(freqs, Calibs, kind='quadratic')

xs2 = np.arange(min(freqs), max(freqs), 1e5)

plt.figure()
plt.plot(freqs, Calibs, 'o')
plt.plot(xs2, h(xs2))