pruebas para comp microm

artiq_master/last_rid.pyon

-20217
+20231

artiq_master/repository/Calibrations/Micromotion compensation/Transitorios_osci.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
+
+#os.chdir('/home/nico/Documents/artiq_experiments/analisis/plots/20221006_transitoriosv2')
+
+os.chdir('/home/nico/Documents/artiq_experiments/analisis/plots/20241210_posibleosc/Data')
+
+def expo(T, tau, N0, C):
+    
+    global T0
+    return N0*np.exp(-(T-T0)/tau) + C
+
+def pow_from_amp(amp):
+    """Paso de amplitud urukul a potencia medida por Nico"""
+    # Forma altamente ineficiente de hacer esto, pero me salio asi
+    amplitudes_UV = np.flip(np.array([0.08, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.22, 0.24, 0.26, 0.28, 0.30]))
+    assert amp in amplitudes_UV
+    potencias_UV = np.flip(np.array([4, 10, 19, 32, 49, 71, 96, 125, 155, 183, 208, 229]))
+    return potencias_UV[np.where(amplitudes_UV == amp)][0]
+
+
+def SP_Bkgr_builder(amp_in, amp_fin, derivadainicio, derivadafin, longbins):
+    CalibCurve = []
+    j=0
+    while j<longbins:
+        if j<=derivadainicio:
+            CalibCurve.append(amp_in)
+        elif j>=derivadainicio and j<=derivadafin:
+            pendiente=(amp_fin-amp_in)/(derivadafin-derivadainicio)
+            CalibCurve.append(amp_in+pendiente*(j-derivadainicio))
+        else:
+            CalibCurve.append(amp_fin)
+        j=j+1
+    return CalibCurve
+
+    
+
+"""
+plt.plot(amplitudes_UV, potencias_UV, 'ko-', lw=0.2)
+plt.xlabel("Amplitud Urukul")
+plt.ylabel("Potencia /uW")
+plt.grid()
+"""
+#%%
+import scipy.fftpack
+
+BINW = 1e-2
+T0 = -0.4e-6
+
+files = [19613] #el 19613 es el que muestra el pico en 48.5 Hz
+
+SP_Heigths = []
+SP_Bins = []
+
+for i, fname in enumerate(files):
+    #print(i)
+    #print(fname)
+    data = h5py.File('Data/0000'+str(fname)+'-CRB.h5', 'r')
+    counts = np.array(data['datasets']['counts'])
+    bines = np.arange(counts.min(), counts.max()+BINW, BINW)
+
+    heigs, binsf  = np.histogram(counts, bines[bines>T0])
+    
+    SP_Heigths.append(heigs)
+    SP_Bins.append(binsf)
+    
+
+
+
+
+from scipy.optimize import curve_fit
+
+refe = len(SP_Bins[0])
+prop = 1
+
+RefBins = [t for t in SP_Bins[0][:int(prop*refe)]]
+
+# plt.figure()
+# plt.plot(RefBins[:-1], SP_Heigths[0][:int(prop*refe)],'-o')
+    
+x = RefBins[:-1]
+y = SP_Heigths[0][:int(prop*refe)]
+
+
+cc = np.ones(len(counts))
+fluo = np.cumsum(cc) - np.polyval(np.polyfit(counts,np.cumsum(cc),1),counts)
+plt.figure()
+plt.plot(counts,fluo)
+plt.xlabel('Tiempo (s)')
+plt.ylabel('Cumsum photons')
+
+
+y = np.array([f for f in fluo])
+dt = np.mean(np.diff(np.array(counts)))
+print(dt)
+
+plt.figure()
+# plt.plot(x, y, label='Señal original')
+plt.magnitude_spectrum(y*np.hanning(len(y)),Fs=1/dt)
+plt.xlabel('Frecuencia (Hz)')
+plt.ylabel('Amplitud')
+# plt.xlim(-10,200)
+plt.grid(True)
+# plt.legend()
+
+#%%
+
+"""
+Esto binea fotones asumienod que la cantidad de fotones que llegan en un intervalo tau
+esta modulada por una funcion periodica. Entonces a los tiempos de llegada entre
+0 y tau no les hace nada, a los tiempos entre tau y 2tau les resta tau,
+a los tiempos entre 2tau y 3tau les resta tau, y asi. Y despues
+binea y ajusta con una sinusoidal y devuelve la frecuencia del ajuste.
+
+Asi anda barbaro para ver la oscilacion a casi 50 Hz
+"""
+
+def bin_time_arrivals(arrival_times, tau):
+    """
+    Binea los tiempos de llegada de los fotones según el periodo tau.
+    
+    Parameters:
+        arrival_times (numpy array): Vector con los tiempos de llegada de los fotones.
+        tau (float): Periodo de bineado.
+        
+    Returns:
+        numpy array: Vector con los tiempos bienados.
+    """
+    return arrival_times - tau * (arrival_times // tau)
+
+taurf = 1/(22.135e6)
+
+taurf = 1/49.9
+
+tiemposarreglados = bin_time_arrivals(counts, taurf)
+
+b = taurf/100  # Ancho de bineo
+bins = np.arange(0, taurf, b)
+hist, bin_edges = np.histogram(tiemposarreglados, bins=bins)
+
+colors = ['b', 'g', 'r']  # Colores para cada repetición
+
+x_extended = np.concatenate([bin_edges[:-1] + i * taurf for i in range(3)])
+y_extended = np.tile(hist, 3)
+
+def sinusoidal(x, A, B, C, D):
+    return A * np.sin(B * x + C) + D
+
+x_extended_dense = np.arange(np.min(x_extended),np.max(x_extended),0.1*(x_extended[1]-x_extended[0]))
+
+params, _ = curve_fit(sinusoidal, x_extended, y_extended, p0=[max(y_extended), 2*np.pi/taurf, 0, np.mean(y_extended)])
+
+# Graficar los datos y el ajuste
+plt.figure(figsize=(8, 4))
+plt.plot(x_extended, y_extended, marker='o', linestyle='-', label="Datos", color='gray')
+plt.plot(x_extended_dense, sinusoidal(x_extended_dense, *params), linestyle='--', label="Ajuste sinusoidal", color='red')
+for i in range(3):
+    plt.plot(bin_edges[:-1] + i * taurf, hist, marker='o', linestyle='-', color=colors[i])
+
+plt.xlabel("Tiempo bineado (repetido 3 veces)")
+plt.ylabel("Frecuencia")
+plt.title("Distribución de tiempos bineados con ajuste sinusoidal")
+plt.legend()
+plt.show()
+
+print(f'Frecuencia del ajuste: {round(params[1]/(2*np.pi),2)} Hz')
+
+
+#%%
+
+"""
+Esto deberia hacerse asi para ver la oscilacion a la RF (micromocion)
+"""
+
+BINW = 1e-2
+T0 = -0.4e-6
+
+files = [19607] #el 19613 es el que muestra el pico en 48.5 Hz
+
+SP_Heigths = []
+SP_Bins = []
+
+for i, fname in enumerate(files):
+    #print(i)
+    #print(fname)
+    data = h5py.File('Data/0000'+str(fname)+'-CRB.h5', 'r')
+    counts = np.array(data['datasets']['counts'])
+    bines = np.arange(counts.min(), counts.max()+BINW, BINW)
+
+    heigs, binsf  = np.histogram(counts, bines[bines>T0])
+    
+    SP_Heigths.append(heigs)
+    SP_Bins.append(binsf)
+    
+    
+
+def bin_time_arrivals(arrival_times, tau):
+    """
+    Binea los tiempos de llegada de los fotones según el periodo tau.
+    
+    Parameters:
+        arrival_times (numpy array): Vector con los tiempos de llegada de los fotones.
+        tau (float): Periodo de bineado.
+        
+    Returns:
+        numpy array: Vector con los tiempos bienados.
+    """
+    return arrival_times - tau * (arrival_times // tau)
+
+freqrf = 22.135e6
+
+taurf = 1/freqrf
+
+tiemposarreglados = bin_time_arrivals(counts, taurf)
+
+b = taurf/50  # Ancho de bineo
+bins = np.arange(0, taurf, b)
+hist, bin_edges = np.histogram(tiemposarreglados, bins=bins)
+
+colors = ['b', 'g', 'r']  # Colores para cada repetición
+
+x_extended = np.concatenate([bin_edges[:-1] + i * taurf for i in range(3)])
+y_extended = np.tile(hist, 3)
+
+def sinusoidal(x, A, C, D):
+    return A * np.sin(2*np.pi*freqrf * x + C) + D
+
+x_extended_dense = np.arange(np.min(x_extended),np.max(x_extended),0.1*(x_extended[1]-x_extended[0]))
+
+params, _ = curve_fit(sinusoidal, x_extended, y_extended, p0=[max(y_extended),  0, np.mean(y_extended)])
+
+# Graficar los datos y el ajuste
+plt.figure(figsize=(8, 4))
+# plt.plot(x_extended, y_extended, marker='o', linestyle='-', label="Datos", color='gray')
+plt.plot(x_extended_dense, sinusoidal(x_extended_dense, *params), linestyle='--', label="Ajuste sinusoidal", color='black',linewidth=3,zorder=2)
+for i in range(3):
+    plt.plot(bin_edges[:-1] + i * taurf, hist, marker='o', linestyle='-', color=colors[i],markersize=3,zorder=1)
+
+plt.xlabel("Tiempo bineado (repetido 3 veces)")
+plt.ylabel("Frecuencia")
+plt.title("Distribución de tiempos bineados con ajuste sinusoidal")
+plt.ylim(0,np.max(hist)*1.2)
+plt.legend()
+plt.show()
+
+print(f'Amplitud normalizada: {np.abs(round(100*params[0]/params[2],3))}%')
+
+#%%
+# Parámetros
+n = len(x)  # Número de muestras
+dt = x[1] - x[0]  # Paso temporal
+
+# Transformada de Fourier
+fft_y = np.fft.fft(y)
+frequencies = np.fft.fftfreq(n, dt)  # Frecuencias asociadas
+
+# Solo nos quedamos con la mitad positiva (frecuencias reales)
+positive_frequencies = frequencies[:n // 2]
+positive_fft = np.abs(fft_y[:n // 2])
+
+# Graficar la señal original y su espectro de frecuencias
+plt.figure()
+# plt.plot(x, y, label='Señal original')
+plt.magnitude_spectrum((y-y.mean())*np.hanning(len(y)),Fs=1/dt)
+plt.title('Señal en el Dominio del Tiempo')
+plt.xlabel('Tiempo (s)')
+plt.ylabel('Amplitud')
+plt.grid(True)
+plt.legend()
+
+# Espectro de frecuencias
+plt.figure()
+plt.plot(positive_frequencies[1:], positive_fft[1:], '-o',label='Espectro de frecuencias')
+plt.title('Transformada de Fourier')
+plt.xlabel('Frecuencia (Hz)')
+plt.ylabel('Amplitud')
+plt.grid(True)
+# plt.xlim(0,160)
+#plt.axvline(50)
+# plt.ylim(-10,300)
+plt.legend()
+
+plt.tight_layout()
+plt.show()