CPT_plotter_20230321.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/20230321_heatingrate/Data/')

CPT_FILES = """000010420-IR_Scan_withcal_optimized
000010436-IR_Scan_withcal_optimized
000010437-IR_Scan_withcal_optimized
000010438-IR_Scan_withcal_optimized
000010439-IR_Scan_withcal_optimized
000010440-IR_Scan_withcal_optimized
000010441-IR_Scan_withcal_optimized
000010425-IR_Scan_withcal_optimized
000010426-IR_Scan_withcal_optimized
000010427-IR_Scan_withcal_optimized
000010428-IR_Scan_withcal_optimized
000010429-IR_Scan_withcal_optimized
000010430-IR_Scan_withcal_optimized
""" 

HEATING_FILES = """000010422-HeatingRate"""


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))
print(SeeKeys(HEATING_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']))


Counts_heating = []
Times = []

AmpTisa_heating = []
UVCPTAmp_heating = []
No_measures_heating = []

for i, fname in enumerate(HEATING_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.
    Times.append(np.array(data['datasets']['heating_times']))
    Counts_heating.append(np.array(data['datasets']['counts_spectrum']))
    #AmpTisa.append(np.array(data['datasets']['TISA_CPT_amp']))
    UVCPTAmp_heating.append(np.array(data['datasets']['UV_CPT_amp']))
    No_measures_heating.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()

#%%
from EITfit.threeLevel_2repumps_AnalysisFunctions import CalculoTeoricoDarkResonances_8levels, GetMinimaInfo, GetPlotsofFluovsAngle_8levels, PerformExperiment_8levels, FindDRFrequencies, FindRelativeFluorescencesOfDR, GenerateNoisyCPT, SmoothNoisyCPT, GetFinalMaps, GenerateNoisyCPT_fixedRabi, GenerateNoisyCPT_fit
from EITfit.threeLevel_2repumps_AnalysisFunctions import MeasureRelativeFluorescenceFromCPT, IdentifyPolarizationCoincidences, RetrieveAbsoluteCoincidencesBetweenMaps, GetClosestIndex
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

c=3e8
B = (u/(2*np.pi))/c

correccion = 6 #con 3 fitea bien

offsetxpi = 440+1+correccion
DetDoppler = -10.5-correccion

FreqsDRpi = [2*f*1e-6-offsetxpi+14 for f in Freqs[0]]
CountsDRpi = Counts[0]

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, scale, offset, temp):
    #temp = 2.9e-4
    MeasuredFreq, MeasuredFluo = GenerateNoisyCPT_fit(SG, sr, SP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, temp, alpha, phidoppler, titadoppler, phiprobe, [titaprobe], phirepump, titarepump, freqs, plot=False, solvemode=1, detpvec=None, noiseamplitude=noiseamplitude)
    #scale = 2.85164232e+04
    #offset = 2.16392191e+03
    FinalFluo = [f*scale + offset for f in MeasuredFluo]
    return FinalFluo


#con doppler broadening solo en uv da: SG, SP, scale, offset, temp
#[4.57165783e-01 8.05907940e+00 6.80017204e+04 3.61150824e-17, 6.77295995e-03] 
#con broadening en ambos laseres da
#[4.65604861e-01 8.06781613e+00 6.65914036e+04 4.12611199e-14, 1.80590400e-03]

#esos valores anteriores dan mal la ordenada al origen, cambie correccion y ahora da
#[4.71134671e-01, 7.63142299e+00, 7.30866544e+04, 1.80899906e+02, 1.20863371e-03]

popt_fullcpt, pcov_fullcpt = curve_fit(FitEITpi, FreqsDRpi, CountsDRpi, p0=[0.5, 4.5, 1e4, 1e3, 1e-3], bounds=((0, 0, 0, 0, 0), (2, 10, 1e5, 1e5, 10e-3)))

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


print(popt_fullcpt)

Sat = popt_fullcpt[0]
Det = popt_fullcpt[1]

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
"""


i_DR = 955

plt.figure()
plt.errorbar(FreqsDRpi, CountsDRpi, yerr=2*np.sqrt(CountsDRpi), fmt='o', capsize=2, markersize=2)
plt.plot(freqslongpi, FittedEITpi)
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()