CPT_plotter_20230421.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

#Mediciones barriendo angulo del TISA y viendo kicking de resonancias oscuras

#C:\Users\Usuario\Documents\artiq\artiq_experiments\analisis\plots\20220106_CPT_DosLaseres_v08_TISA_DR\Data

os.chdir('/home/nico/Documents/artiq_experiments/analisis/plots/20230421_CPTconmicromocion/Data/')

CPT_FILES = """000011071-IR_Scan_withcal_optimized
000011072-IR_Scan_withcal_optimized
000011073-IR_Scan_withcal_optimized
""" 


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

print(SeeKeys(CPT_FILES))

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

Counts = []
Freqs = []

AmpTisa = []
UVCPTAmp = []
No_measures = []

for i, fname in enumerate(CPT_FILES.split()):
    print(str(i) + ' - ' + fname)
    #print(fname)
    data = h5py.File(fname+'.h5', 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'

    # Aca hago algo repugnante para poder levantar los strings que dejamos
    # que además tenian un error de tipeo al final. Esto no deberá ser necesario
    # cuando se solucione el error este del guardado.
    Freqs.append(np.array(data['datasets']['IR1_Frequencies']))
    Counts.append(np.array(data['datasets']['counts_spectrum']))
    #AmpTisa.append(np.array(data['datasets']['TISA_CPT_amp']))
    UVCPTAmp.append(np.array(data['datasets']['UV_CPT_amp']))
    No_measures.append(np.array(data['datasets']['no_measures']))


#%%

"""
Ploteo la cpt de referencia / plotting the reference CPT
"""

jvec = [0]

plt.figure()
i = 0
for j in jvec:
    plt.errorbar([2*f*1e-6 for f in Freqs[j]], Counts[j], yerr=np.sqrt(Counts[j]), fmt='o', capsize=2, markersize=2)
    i = i + 1
plt.xlabel('Frecuencia (MHz)')
plt.ylabel('counts')
plt.grid()
plt.legend()

#%%

"""
Intento mergear 0 1 y 2
"""

jvec = [0]

#f1 = 8.8
#f2 = 17.5
f1 = 8
f2 = 18

Freqs0 = Freqs[0]
Freqs1 = [f+f1*1e6 for f in Freqs[1]]
Freqs2 = [f+f2*1e6 for f in Freqs[2]]

plt.figure()

plt.errorbar([2*f*1e-6 for f in Freqs0], Counts[0], yerr=np.sqrt(Counts[0]), fmt='o', capsize=2, markersize=2)
plt.errorbar([2*f*1e-6 for f in Freqs1], Counts[1], yerr=np.sqrt(Counts[1]), fmt='o', capsize=2, markersize=2)
plt.errorbar([2*f*1e-6 for f in Freqs2], Counts[2], yerr=np.sqrt(Counts[2]), fmt='o', capsize=2, markersize=2)

plt.xlabel('Frecuencia (MHz)')
plt.ylabel('counts')
plt.grid()
plt.legend()


Freqs_merged = list(Freqs0) + Freqs1[125:] + Freqs2[126:]
Counts_merged = list(Counts[0]) + list(Counts[1][125:]) + list(Counts[2][126:])

plt.figure()

plt.errorbar([2*f*1e-6 for f in Freqs_merged], Counts_merged, yerr=np.sqrt(np.array(Counts_merged)), fmt='o', capsize=2, markersize=2)



#%%
#from EITfit.MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels_MM, GenerateNoisyCPT_MM_fit
from scipy.optimize import curve_fit

"""
Ajusto un cpt para obtener todos los parámetros relevantes primero.
I fit a cpt curve to retrieve all the relevant parameters first.
"""

phidoppler, titadoppler = 0, 90
phirepump, titarepump = 0,  0
phiprobe = 0
titaprobe = 90

gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 
alpha = 0
noiseamplitude = 0

T = 0.6e-3

sg = 0.544
sp = 4.5
sr = 0
DetRepump = 0


lw = 0.1
DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth = lw, lw, lw #ancho de linea de los laseres


u = 33.5e6

beta = 0

drivefreq = 2*np.pi*22.135e6

correccion = 12

offsetxpi = 458+correccion
DetDoppler = -2.5   -correccion

FreqsDRpi = [2*f*1e-6-offsetxpi+14 for f in Freqs_merged]
CountsDRpi = Counts_merged

freqslongpi = np.arange(min(FreqsDRpi), max(FreqsDRpi)+FreqsDRpi[1]-FreqsDRpi[0], 0.1*(FreqsDRpi[1]-FreqsDRpi[0]))

#[1.71811842e+04 3.34325038e-17]

def FitEITpi(freqs, SG, SP, BETA, scale, offset):
    temp = 1e-3
    MeasuredFreq, MeasuredFluo = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, BETA, drivefreq, min(freqs), max(freqs), freqs[1]-freqs[0])
    #MeasuredFreq, MeasuredFluo = GenerateNoisyCPT_MM_fit(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, BETA, drivefreq, freqs, plot=False, solvemode=1, detpvec=None, noiseamplitude=noiseamplitude)
    #scale = 0.3*2e+04
    #offset = 0.3*2e+03
    FinalFluo = [f*scale + offset for f in MeasuredFluo]
    return FinalFluo

popt_fullcpt, pcov_fullcpt = curve_fit(FitEITpi, FreqsDRpi[1:], CountsDRpi, p0=[0.5, 4.5, 1, 1e4, 1e3], bounds=((0, 0, 0, 1e1, 0), (2, 10, 1000, 1e5, 1e5)))

#print(f'Temperatura: ({round(1e3*popt_fullcpt[-1],2)} +- {round(1e3*np.sqrt(pcov_fullcpt[-1][-1]),2)}) mK')


print(popt_fullcpt)
#%%

FittedEITpi = FitEITpi(freqslongpi, popt_fullcpt[0], popt_fullcpt[1], 4, popt_fullcpt[3], popt_fullcpt[4])
#FittedEITpi = FitEITpi(freqslongpi, *popt_fullcpt)

"""
Ploteo la CPT de referencia junto al ajuste y a la resonancia oscura de interes
I plot the reference CPT along with the fit to the model and the dark resonance of interest
"""

DRs = [-26.5, -18, -11.5, -3]

plt.figure()
plt.errorbar(FreqsDRpi, Counts_merged, yerr=np.sqrt(np.array(Counts_merged)), fmt='o', capsize=2, markersize=2)
for dr in DRs:
    plt.axvline(dr, color='black', alpha=0.5)
    plt.axvline(dr-22.1, color='red', alpha=0.5)
    plt.axvline(dr+22.1, color='blue', alpha=0.2)

plt.plot(freqslongpi, FittedEITpi[1:])


#%%
#ignorar de aca para abajo por ahora

#%%
#i_DR = 955

DRs = [-31, -22, -16.5, -8]

plt.figure()
plt.errorbar(FreqsDRpi, CountsDRpi, yerr=2*np.sqrt(CountsDRpi), fmt='o', capsize=2, markersize=2)
plt.plot(freqslongpi[:-1], FittedEITpi)
for dr in DRs:
    dr = dr+4.7
    plt.axvline(dr, color='red',alpha=0.2)
    plt.axvline(dr+22.1, color='green',alpha=0.2)
    plt.axvline(dr-22.1, color='black',alpha=0.2)
#plt.axvline(DetDoppler-22.1)
#plt.axvline(DetDoppler+22.1)
#plt.plot(freqslongpi[i_DR], FittedEITpi[i_DR],'o', color='red', markersize=12)
plt.xlabel('Detuning (MHz)')
plt.ylabel('Counts')

#plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')



#%%
"""
Simulo CPTs con todos esos parámetros para distintas temperaturas
I simulate CPT curves with all the previous parameters but with varying temperatures
"""
TempVecTeorico = list(np.arange(0.3,1,0.1))+list(np.arange(1, 31, 1))
CurvasTeoricas = []

for tempi in TempVecTeorico:
    CurvasTeoricas.append(FitEITpi(freqslongpi, *popt_fullcpt[:-1], tempi*1e-3))



#%%
"""
Acá agarro la primera y busco el valor i_DR que corresponde a la resonancia oscura de interés
With the first one, I look for the value i_DR which corresponds to the dark resonance of interest
"""
curva_ref = CurvasTeoricas[0]

i_DR = 955

plt.figure()

plt.plot(freqslongpi, curva_ref)
plt.plot(freqslongpi[i_DR], curva_ref[i_DR],'o')

#%%
"""
ploteo algunos CPTs teoricos para algunas temperaturas
Plotting some theory cpt curves for some temperatures
"""

plt.plot(freqslongpi, CurvasTeoricas[0])
plt.plot(freqslongpi[i_DR], CurvasTeoricas[0][i_DR],'o',markersize=10)
plt.plot(freqslongpi, CurvasTeoricas[10])
plt.plot(freqslongpi[i_DR], CurvasTeoricas[10][i_DR],'o',markersize=10)
plt.plot(freqslongpi, CurvasTeoricas[20])
plt.plot(freqslongpi[i_DR], CurvasTeoricas[20][i_DR],'o',markersize=10)
plt.xlabel('Detuning (MHz)')
plt.ylabel('Fluorescence')
plt.grid()



#%%
"""
Ahora interpolo los valores teóricos de las profundidades de esas resonancias
y aplico la interpolación a las mediciones para obtener temperaturas.
Luego, grafico las temperaturas en función de los tiempos de calentamiento.

Now I interpolate the theoretical values of the depths of those resonances
and apply the interpolation to the measurements to obtain temperatures.
After that, I plot the temperatures with respect to the heating times
"""
from scipy.interpolate import interp1d

FluosDRTeo = [CurvasTeoricas[k][i_DR] for k in range(len(CurvasTeoricas))]

interpolado = interp1d(FluosDRTeo, TempVecTeorico) #creo funcion que interpola


meas = 0
maxi = 9 #valor maximo, dsp el ion se calento
 
Heating_tim = Times[meas][:maxi]
Heating_tim_ms = [t*1e3 for t in Heating_tim]
Heating_med = [2*c for c in Counts_heating[meas][:maxi]]

ErrorHeating_med = [2*np.sqrt(c) for c in Counts_heating[meas][:maxi]]

Temperaturas_interpoladas = [float(interpolado(h)) for h in Heating_med]
Error_Temperaturas_interpoladas = [float(interpolado(Heating_med[k]+0.5*ErrorHeating_med[k]))-Temperaturas_interpoladas[k] for k in range(len(Heating_med))]

plt.figure()
plt.plot(FluosDRTeo, [1*t for t in TempVecTeorico], 'o', color='orange')
plt.plot(np.linspace(FluosDRTeo[0],FluosDRTeo[-1],1000), interpolado(np.linspace(FluosDRTeo[0],FluosDRTeo[-1],1000)))
plt.xlabel('Cuentas de DR teoricas')
plt.ylabel('Temperatura (mK)')

plt.figure()
#plt.plot(Heating_med, Heating_tim, 'o', color='blue')
plt.errorbar([t*1e3 for t in Heating_tim], Heating_med, yerr=ErrorHeating_med, fmt='o', capsize=2, markersize=5)

plt.ylabel('Cuentas de DR medidas')
plt.xlabel('Heating time (s)')

def lineal(x,a,b):
    return a*x+b

#p1,p2 = curve_fit(lineal, Heating_tim_ms, Temperaturas_interpoladas, sigma=Error_Temperaturas_interpoladas)
p1,p2 = curve_fit(lineal, Heating_tim_ms, Temperaturas_interpoladas)


#%%
"""
Grafico finalmente el plot del heating rate de la trampa
Finally I plot the heating rate of the trap
"""


plt.figure()
#plt.plot(Heating_tim_ms,Temperaturas_interpoladas,'o')
plt.errorbar(Heating_tim_ms,Temperaturas_interpoladas, yerr=np.array(Error_Temperaturas_interpoladas), fmt='o', capsize=2, markersize=7, color='black')
plt.plot(Heating_tim_ms, lineal(np.array(Heating_tim_ms), *p1), color='red')
plt.xlabel('Heating time (ms)', fontname='STIXGeneral', fontsize=15)
plt.ylabel('Temperature (mK)', fontname='STIXGeneral', fontsize=15)
plt.grid()
plt.xticks([0, 5, 10, 15, 20 ,25, 30, 35], fontname='STIXGeneral', fontsize=15)
plt.yticks([0, 5, 10, 15], fontname='STIXGeneral', fontsize=15)
plt.title(f'Heating rate: ({round(p1[0],2)} +- {round(np.sqrt(p2[0][0]),2)}) mK/ms', fontname='STIXGeneral', fontsize=15)
plt.tight_layout()
plt.savefig('Fig_heatingrate.svg')
print(f'Heating rate: ({round(p1[0],2)} +- {round(np.sqrt(p2[0][0]),2)}) mK/ms')


#%%
"""
Ahora voy a ver CPT enteras con tiempos de calentamiento distintos.
Now I see whole CPT curves with different heating times
"""

jvec = [3, 4]

plt.figure()
i = 0
for j in jvec:
    if j==4:
        plt.errorbar([2*f*1e-6 for f in Freqs[j]], Counts[j], yerr=np.sqrt(Counts[j]), fmt='o', capsize=2, markersize=2, label='Without heating')
    elif j==3:
        plt.errorbar([2*f*1e-6 for f in Freqs[j]], Counts[j], yerr=np.sqrt(Counts[j]), fmt='o', capsize=2, markersize=2, label='5 ms heating')
    i = i + 1
plt.xlabel('Frecuencia (MHz)')
plt.ylabel('counts')
plt.ylim(1000,2900)
plt.grid()
plt.legend()


jvec = [1, 2]

plt.figure()
i = 0
for j in jvec:
    if j==2:
        plt.errorbar([2*f*1e-6 for f in Freqs[j]], Counts[j], yerr=np.sqrt(Counts[j]), fmt='o', capsize=2, markersize=2, label='Without heating')
    elif j==1:
        plt.errorbar([2*f*1e-6 for f in Freqs[j]], Counts[j], yerr=np.sqrt(Counts[j]), fmt='o', capsize=2, markersize=2, label='1 ms heating')
    i = i + 1
plt.xlabel('Frecuencia (MHz)')
plt.ylabel('counts')
plt.ylim(1000,2900)
plt.grid()
plt.legend()

#%%

"""
La siguiente curva probablemente no este bien medida ya que inmediatamente
despues, los laseres se deslockearon. La dejo por las dudas.

This curve is probably not well measured...

"""

jvec = [5, 6]

plt.figure()
i = 0
for j in jvec:
    if j==6:
        plt.errorbar([2*f*1e-6 for f in Freqs[j]], Counts[j], yerr=np.sqrt(Counts[j]), fmt='o', capsize=2, markersize=2, label='Without heating')
    elif j==5:
        plt.errorbar([2*f*1e-6 for f in Freqs[j]], Counts[j], yerr=np.sqrt(Counts[j]), fmt='o', capsize=2, markersize=2, label='10 ms heating')
    i = i + 1
plt.xlabel('Frecuencia (MHz)')
plt.ylabel('counts')
plt.ylim(1000,3900)
plt.grid()
plt.legend()
plt.title('Ojo: medicion condicionada por derivas')


#%%
"""
Ahora ploteo 6 curvas cpt para distintos valores de potencia del UV

This is a plot of 6 different cpt curves for 6 different UV powers. I should fit them
to obtain saturation parameters
"""


jvec = [7,8,9,10,11,12]

plt.figure()
for j in jvec:
    plt.errorbar([2*f*1e-6 for f in Freqs[j]], Counts[j], yerr=np.sqrt(Counts[j]), fmt='o', capsize=2, markersize=2, label='Without heating')
plt.xlabel('Frecuencia (MHz)')
plt.ylabel('counts')
#plt.ylim(1000,2900)
plt.grid()
#plt.legend()