agrego lo que falta del analisis

analisis/plots/20230510_MotionalSpectrum/MotionalElectric_new.py

@@ -14,7 +14,7 @@ from scipy import interpolate
 #os.chdir('/home/nico/Documents/artiq_experiments/analisis/plots/20230510_MotionalSpectrum/Data/')
 
 Data = []
-for ii in range(1,5):
+for ii in range(1,10):
     Data.append(f"ss0{ii}.dat")
 #for ii in range(10,20):
 #    Data.append(f"ss{ii}.dat")
@@ -35,8 +35,8 @@ for dd in Data:
     FluoVec.append([data[i][2] for i in range(len(data))])
 
 #%%
-
-ivec = [1,2,3]
+#los voltajes son 2, 2.5, 3, 5, 7 y 9
+ivec = [1,2,3,4, 5, 6, 7,8]
 
 plt.figure()
 for i in ivec:
@@ -49,6 +49,35 @@ for i in ivec:
     plt.plot([f*1e-6 for f in FrequencyVec[i]], FluoVec[i],'o')
 plt.xlim(0.87,0.89)    
 
+#%%
+#veo la menor radial en funcion del voltaje pp del driving
+
+ivec = [1,2,3,4, 5, 6, 7, 8]
+
+Vpp = [0.5, 2, 2.5, 3, 5, 2.5, 7, 9, 1]
+
+plt.figure()
+for i in ivec:
+    plt.plot([f*1e-6 for f in FrequencyVec[i]], FluoVec[i],'o', label=f'V={Vpp[i]} v')
+plt.xlim(0.78,0.8)    
+plt.legend()
+
+#%%
+#lo arreglo un poco
+plt.figure()
+plt.plot([f*1e-6+0.001700 for f in FrequencyVec[6]], FluoVec[6],'-o', label=f'V={Vpp[6]} v')
+plt.plot([f*1e-6+0.0025 for f in FrequencyVec[7]], FluoVec[7],'-o', label=f'V={Vpp[7]} v')
+plt.plot([f*1e-6+0.0013 for f in FrequencyVec[8]], FluoVec[8],'-o', label=f'V={Vpp[8]} v')
+plt.plot([f*1e-6 for f in FrequencyVec[3]], FluoVec[3],'-o', label=f'V={Vpp[3]} v')
+plt.plot([f*1e-6 for f in FrequencyVec[4]], FluoVec[4],'-o', label=f'V={Vpp[4]} v')
+plt.xlim(0.78,0.8)    
+plt.legend()
+
+
+
+
+
+
 #%%
 import scipy
 

analisis/plots/20230510_MotionalSpectrum/MotionalOptic_TransitoryEffect.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/20230510_MotionalSpectrum/DataOpticTransitory/')
+
+
+MOTIONAL_FILES = """
+000012033-AD9910RAM
+000012032-AD9910RAM
+000012031-AD9910RAM
+000012030-AD9910RAM
+000012025-AD9910RAM
+000012026-AD9910RAM
+000012027-AD9910RAM
+000012028-AD9910RAM
+000012029-AD9910RAM
+000012034-AD9910RAM
+"""
+
+
+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(MOTIONAL_FILES))
+
+#%%
+#carpeta pc nico labo escritorio:
+#C:\Users\Usuario\Documents\artiq\artiq_experiments\analisis\plots\20211101_CPT_DosLaseres_v03\Data
+
+Counts = []
+RealFreqs = []
+UV_amp_vec = []
+
+for i, fname in enumerate(MOTIONAL_FILES.split()):
+    print(str(i) + ' - ' + fname)
+    data = h5py.File(fname+'.h5', 'r')
+    RealFreqs.append(np.array(data['datasets']['real_freq']))
+    Counts.append(np.array(data['datasets']['counts']))
+    UV_amp_vec.append(np.array(data['datasets']['UV_amp']))
+
+
+Potencias = [7.4, 11.4, 16.8, 23.2, 31.1, 40.7, 51.5, 64, 78.2, 92.7, 108, 124, 134, 154, 167]
+
+PotenciasUsadas = [7.4, 11.4, 16.8, 23.2, 31.1, 40.7, 51.5, 64, 78.2]
+#%%
+
+"""
+Ploteo una curva para buscar su minimo
+"""
+
+jvec = [0,1,2]
+
+plt.figure()
+i = 0
+
+kmin = 106
+
+for j in jvec:
+    plt.errorbar([1*f*1e-3 for f in RealFreqs[j]], Counts[j], yerr=np.sqrt(Counts[j]), fmt='o', capsize=2, markersize=2, label=f"Pot: {Potencias[j]} uW")
+    plt.plot([1*f*1e-3 for f in RealFreqs[j]][kmin], Counts[j][kmin], 'o', markersize=15)
+    i = i + 1
+plt.xlabel('Frecuencia (kHz)')
+plt.ylabel('counts')
+#plt.xlim(782,787)
+#plt.ylim(2000,5000)
+plt.grid()
+plt.legend(loc='upper left')
+
+
+
+
+
+#%%
+
+"""
+Ploteo las curvas de referencia 
+"""
+
+jvec = [0, 1, 2, 3, 8]
+
+
+
+plt.figure()
+i = 0
+for j in jvec:
+    plt.errorbar([1*f*1e-3 for f in RealFreqs[j]], Counts[j], yerr=np.sqrt(Counts[j]), fmt='o', capsize=2, markersize=2, label=f"Pot: {Potencias[j]} uW")
+    i = i + 1
+plt.xlabel('Frecuencia (kHz)')
+plt.ylabel('counts')
+#plt.xlim(782,787)
+plt.ylim(2000,5000)
+plt.grid()
+plt.legend(loc='upper left')
+
+
+#%%
+"""
+Busco el cociente entre el minimo y el maximo
+"""
+kmins = [106, 106, 106, 106, 110, 110, 110, 112, 108]
+
+
+minabs = np.min(Counts[2])
+
+MinimosFluos = []
+MaximosFluos = []
+CocientesFluos = []
+ErrorCocientesFluos = []
+
+for m in range(9):
+    mini = Counts[m][kmins[m]]-minabs
+    maxi = Counts[m][0]-minabs
+    MinimosFluos.append(mini)
+    MaximosFluos.append(maxi)
+    CocientesFluos.append(mini/maxi)
+    ErrorCocientesFluos.append((mini/maxi)*(np.sqrt(mini)/mini + np.sqrt(maxi)/maxi))
+    
+    
+plt.figure()
+plt.errorbar(PotenciasUsadas, CocientesFluos, yerr=ErrorCocientesFluos, fmt='o', capsize=5, markersize=15)
+plt.axhline(1)
+plt.xlabel('Potencia IR (uW)')
+
+
+
+
+
+
+
+

analisis/plots/20230523_transitorioshighpower/Simulaciones/OBE3levels_chtime_fitting.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+import numpy as np
+import matplotlib.pyplot as plt
+from OBE3levels_functions import *
+from scipy.optimize import curve_fit
+from scipy.stats import poisson, norm
+import os
+import sys
+plt.ion()
+"""
+Este código simula la dinámica temporal de un ion de calcio sometido a un láser UV y dos láseres IR.
+
+"""
+def expo(T, tau, N0, C):
+    #global T0
+    return N0*np.exp(-(T-0e-6)/tau) + C
+
+ADD_NOISE = True
+NOISE_VAR = 30 # Varianza del ruido
+NOISE_ADJ = 78 # Escaleo del ruido versus el maximo
+NOISE_MEAN = 0 # Media del ruido
+NRO_REPS = 100 # Veces a promediar
+
+species = 'Calcium'
+#Parámetros de saturacion de los láseres rebombeo (IR), que son proporcionales a la intensidad de los mismos
+#Para reproducir SP, hay que poner 0 a los dos
+SatParRep1, SatParRep2 = 0, 0
+
+sat_intensity = 43.3e-5 # uW um-2
+g21 = 135591138.92893547
+
+#El número de los láseres IR aplicados en el experimento. Para el experimento de la transicion SP, es 0.
+Number_IR_lasers = 0
+
+#Detunings de los láseres en MHz. Si es 0, están en resonancia. En principio
+#podría ser un parámetro libre a ajustar porque influye mucho
+DetDoppler = -int(sys.argv[1])
+DetRep2 = 0
+DetRep1 = 0
+
+#Anchos de línea de los láseres, en MHz. No son tan cruciales, están entre 0.1 y 0.5 MHz
+DopplerLaserLinewidth, RepumpLaserLinewidth = 1e3, 1
+
+# Para trabajar con las curvas teóricas:
+# range_inicio = np.arange(sat_intensity/10, sat_intensity, sat_intensity/50)
+# range_fin = np.arange(sat_intensity, 30*sat_intensity, sat_intensity/5)
+# laser_intensities = np.append(range_inicio, range_fin)
+
+# Para trabajar con nuestros datos, saco las intensidades, a un radio de 100um:
+laser_intensities = np.array([0.00662085, 0.00582507, 0.0049338, 0.00397887, 0.00305577, 0.00226, 0.00155972, 0.00101859, 0.00060479, 0.00031831])
+
+
+OmegaRabi2 = (g21**2) * (laser_intensities / sat_intensity)
+SatParDopplerVec = OmegaRabi2 / ((DetDoppler*2*np.pi*1e6)**2 + g21**2)
+
+#Parámetros de la simulación, en us
+tfinOBE = 1
+pasoOBE = 1e-5
+
+#Estado electrónico inicial de la simulacion. Posibilidades: |S> = 1, |P> = 2, |D> = 3
+initial_state = 1
+
+#Corre y simula la dinámica. Devuelve elementos de matriz densidad
+# fig1, ax1 = plt.subplots()
+
+alltaus = list()
+allN0 = list()
+allmeans = list()
+allstds = list()
+
+temporal_taus = list()
+temporal_N0 = list()
+
+for i, SatParDoppler in enumerate(SatParDopplerVec):
+    t, RealRho11, RealRho22, RealRho33, AbsRho12, AbsRho13, AbsRho23 = Solve3LevelBloch(EquationSystemModel, species, SatParDoppler, SatParRep1, SatParRep2, DetDoppler, DetRep1, DetRep2, DopplerLaserLinewidth, RepumpLaserLinewidth, no_IR_lasers=Number_IR_lasers, initcond=initial_state, tfin=tfinOBE, paso=pasoOBE, printprogress=False)
+
+    k = int(0.1*len(t)) #parametro para que fitee solo la parte final
+    k = np.argmin(np.abs(t - 0.2e-6))
+
+    temporal_taus = list()
+    incert_estad = 0
+
+    for nro_rep in range(NRO_REPS):
+        poblacion_fit = RealRho22[k:].copy()
+        if ADD_NOISE:
+            # Genero ruido similar al medido
+            ruido_poisson = poisson.rvs(NOISE_VAR, size=len(poblacion_fit))
+            # Reescaleo el ruido para que sea compatible con las simulaciones
+            ruido_poisson *= np.max(poblacion_fit)/NOISE_ADJ
+
+            poblacion_fit += ruido_poisson
+            # poblacion_fit += norm.rvs(0, NOISE_VAR, size=len(poblacion_fit))*np.max(poblacion_fit)/NOISE_ADJ
+            poblacion_fit[np.where(poblacion_fit < 0 )] = 0 # Limpio posibles valores negativos
+
+        t_fit = t[k+1:]
+
+        popt, pcov = curve_fit(expo, t_fit, poblacion_fit, p0=(1e-6, 1e-2, 1e-2))
+
+        temporal_taus.append(popt[0])
+        incert_estad += np.sqrt(pcov[0,0])
+        # temporal_N0.append(popt[1])
+        if not ADD_NOISE: # si no quiere que agregue ruido, entonces solo hace el for una vez
+            break
+
+        print(f" {nro_rep}  ", end='\r')
+
+    if ADD_NOISE:
+        allmeans.append(np.mean(temporal_taus))
+        allstds.append(np.std(temporal_taus) + incert_estad/NRO_REPS)
+    else:
+        allmeans.append(np.mean(temporal_taus))
+
+
+    #ploteo en funcion de 100*Rho22 para que me devuelva en porcentaje la población del excitado que es proporcional a la fluorescencia detectada
+    # lbl='Detuning Rep1:' + str(DetRep1) + ' MHz'
+    # ax1.plot(t[:-1]*1e6, RealRho22, '.', label=lbl, markersize=1)
+    # ax1.plot(t[:-1]*1e6, 100*RealRho22, '.', label=lbl, markersize=1)
+    # ax1.plot(t_fit*1e6, 100*expo(t_fit, *popt))
+    print(f'({i:02d}/{len(SatParDopplerVec)}) Done {SatParDoppler}'.ljust(50, ' '))
+
+# ax1.legend(SatParDopplerVec)
+# ax1.set_xlabel('Time (us)')
+# ax1.set_ylabel('100*Rho22')
+
+#%% #####################################################################################
+### Guardo valores al CSV para analizar luego
+# nDet=f"Det{np.abs(DetDoppler)}"
+# CSV_FNAME = f"error_teorico_tau_{nDet}.csv"
+# try:
+#     import pandas as pd
+#     data = pd.read_csv(CSV_FNAME)
+#     data["meanval"] = allmeans
+#     data["stdval"] = allstds
+#     data.to_csv(CSV_FNAME, index=False)
+#     print(f"SAVED {nDet}".ljust(50, ' '))
+# except FileNotFoundError:
+#     data = pd.DataFrame({'I': laser_intensities, "meanval": allmeans, "stdval": allstds})
+#     data.to_csv(CSV_FNAME, index=False)
+#     print(f"CREATED CSV FILE".ljust(50, ' '))
+#     print(f"SAVED DETUNING {-DetDoppler}".ljust(50, ' '))
+
+#%% ####################################################################################
+### Plots de los valores de intensidades medidas para distintos parametros
+### Cada axes arma un eje horizontal distinto, es para comparar las posibilidades
+
+# fig4, ax4 = plt.subplots()
+# allmeans = [t*1e6 for t in allmeans]
+
+# ax4b = ax4.twinx()
+# ax4.plot(SatParDopplerVec, allmeans, '-o', lw=0.4)
+# ax4b.plot(SatParDopplerVec, allN0, 'k-^', lw=0.4)
+# ax4.set_xlabel("Parametro de saturación")
+# ax4.set_ylabel("Tau (circulo)")
+# ax4b.set_ylabel("Alturas (triang)")
+# ax4.set_title(f"SP; IR_Lasers=0; DetuningDoppler={DetDoppler}; LineWidth={DopplerLaserLinewidth}")
+
+# ax4by = ax4.twiny()
+# ax4by.plot(laser_intensities, allmeans, '-', lw=0)
+# ax4by.xaxis.set_label_position('bottom')
+# ax4by.xaxis.tick_bottom()
+# ax4by.spines['bottom'].set_position(("axes", -0.27))
+# ax4by.set_xlabel(r"Intensidad [ $\mu$W/$\mu$m$^2$ ]")
+
+# ax4cy = ax4.twiny()
+# ax4cy.plot(laser_intensities*np.pi*(35**2), allmeans, '-', lw=0)
+# ax4cy.xaxis.set_label_position('bottom')
+# ax4cy.xaxis.tick_bottom()
+# ax4cy.spines['bottom'].set_position(("axes", -0.58))
+# ax4cy.set_xlabel(r"Potencia (r=35$\mu$m) [ $\mu$W ]")
+
+# ax4dy = ax4.twiny()
+# ax4dy.plot(laser_intensities*np.pi*(50**2), allmeans, '-', lw=0)
+# ax4dy.xaxis.set_label_position('bottom')
+# ax4dy.xaxis.tick_bottom()
+# ax4dy.spines['bottom'].set_position(("axes", -0.9))
+# ax4dy.set_xlabel(r"Potencia (r=50$\mu$m) [ $\mu$W ]")

analisis/plots/20230523_transitorioshighpower/Simulaciones/error_teorico_tau_det20.csv

+I,meanval,stdval
+0.00662085,2.63563456542463e-07,3.62795752e-08
+0.00582507,2.674966392578803e-07,3.99478821e-08
+0.0049338,2.7426723419290064e-07,3.66249293e-08
+0.00397887,2.8263714050526124e-07,3.83554920e-08
+0.00305577,2.9770905441813384e-07,3.90291265e-08
+0.00226,3.1786383950851434e-07,4.27893803e-08
+0.00155972,3.562791436115684e-07,4.67339215e-08
+0.00101859,4.198979464288562e-07,5.28900825e-08
+0.00060479,5.457916989936301e-07,6.95843326e-08
+0.00031831,8.275611779916097e-07,1.30369576e-07

analisis/plots/20230523_transitorioshighpower/Simulaciones/medidos_powers_taus.csv

+Pow,Tau,N0
+208.0,4.107136373675053e-07,242.3949275128121
+183.0,3.994370315889412e-07,255.700355179914
+155.0,4.4882562779198994e-07,226.84461233128172
+125.0,4.727871415748306e-07,217.84774878119222
+96.0,5.254285404397582e-07,222.57732196795308
+71.0,6.079817946589069e-07,189.5426628192835
+49.0,7.274218641328625e-07,168.69185328968453
+32.0,8.823248220314772e-07,140.40570782757854
+19.0,1.280491601444429e-06,90.5188484395718
+10.0,1.96934320559994e-06,60.092743064494314

analisis/plots/20230523_transitorioshighpower/Simulaciones/plot_decaytimes_vs_intensity_detunings.py

+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+#plt.ion()
+
+#dataREAL = pd.read_csv("medidos_powers_taus.csv")
+dataSIM  = pd.read_csv("sim_intensidad_taus_full.csv")
+errsSIM  = pd.read_csv("error_teorico_tau_det20.csv")
+I_sat = 43.3*1e-5 #uW/um2
+
+#%% #####################################################################################
+### Plot de lineas teoricas tau-intensidad variando detuning + datos medidos
+
+
+UVpotVec=[4,6,12,17,26,36, 48, 77, 112, 151, 190, 226, 255, 279]
+
+Taus = [2.6474524026975392, 1.651982852985559, 1.0793431020294597, 0.8171909874107741, 0.6329192159654833, 0.5246427838474966, 0.468752585310023, 0.3943320198588099, 0.34233658233849174, 0.314035400291993, 0.30338246040910105,0.2891087574893047, 0.2861944000648665,0.2804274973493638]
+Taus_v2 = [1.6655917442406167, 1.1754160164890701, 1.0129299348502132, 0.6747311668331314, 0.5395244011650271, 0.4699144306437099, 0.41805942894076137, 0.3948601586682529, 0.3849453824721735, 0.3752507710597908]
+
+fig2, ax2 = plt.subplots(figsize=(4,3))
+detuning_lines = ['Det10', 'Det20', 'Det30', 'Det40', 'Det50', 'Det60', 'Det70', 'Det80', 'Det90']
+
+# Curvas de las simulaciones con tau en unidades de microseg
+for det in detuning_lines:
+    ax2.plot(dataSIM.I**1, [d for d in dataSIM[det]*1e6], 'k-', ms=1, lw=0.5)
+
+# Curvas de las mediciones cambiando el radio posible con tau en microseg
+
+for cambio in [0]: # esto es para shiftear la potencia medida
+    potencias = dataREAL.Pow - cambio
+    for rad in [128]: # esto evalua con distintos radios
+        ax2.errorbar([(p/(np.pi*(rad**2)))**1 for p in UVpotVec], [t for t in Taus], yerr=np.mean(1e6*errsSIM.stdval.values),
+                    fmt='.--', ms=6, lw=.5, \
+                    color="#3949ab",
+                    capsize=2)
+        ax2.errorbar([(p/(np.pi*(rad**2)))**1 for p in UVpotVec][4:], [t for t in Taus_v2], yerr=np.mean(1e6*errsSIM.stdval.values),
+                    fmt='.--', ms=6, lw=.5, \
+                    color="#3949ab",
+                    capsize=2)
+        
+fs = 9
+
+ax2.set_xlim([0.47e-4, 0.8e-2])
+#ax2.set_xlim([0., 0.000005])
+#ax2.set_ylim([2e-1, 4e1])
+ax2.set_ylabel(r"Characteristic time ($\mu$s)", fontname='STIXGeneral')
+#ax2.set_xlabel(r"$I = P / (\pi\,r^2)$ [$\mu$W/$\mu$m$^2$]")
+ax2.set_xlabel(r'UV intensity ($\mu$W/$\mu$m$^2$)', fontname='STIXGeneral')
+# ax2.set_title(r"datos(punteadas) $r \in [30; 190]\ \mu m$   |   simulacion(llenas) $\Delta \in -[10; 90]$ MHz")
+#ax2.legend(markerscale=2)
+#ax2.set_xticks([1e-4, 1e-3])
+ax2.set_xticklabels([1e-4, 1e-3], fontsize=fs, fontname='STIXGeneral')
+#ax2.set_yticks([1e0, 1e1])
+ax2.set_yticklabels([1e0, 1e1], fontsize=fs, fontname='STIXGeneral')
+ax2.grid(which='minor', alpha=0.2)
+ax2.set_xscale("log")
+ax2.set_yscale("log")
+
+plt.savefig('Fig03_varyingUVpower.pdf')
+plt.savefig('Fig03_varyingUVpower.svg')
+

analisis/plots/20230523_transitorioshighpower/Transitorios_03.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/20230523_transitorioshighpower/')
+# Solo levanto algunos experimentos
+SP_files = [12221, 12222, 12223, 12224]
+
+Random_files = [8749]
+
+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()
+"""
+#%%
+
+BINW = 1e-9
+T0 = -0.4e-6
+
+SP_Heigths = []
+SP_Bins = []
+
+for i, fname in enumerate(SP_files):
+    #print(i)
+    #print(fname)
+    data = h5py.File('Data/SP/0000'+str(fname)+'-SingleLine.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)
+    
+
+#%%
+"""
+Vectores de amplitudes y potencias
+"""
+
+UVampVec = [0.06, 0.07, 0.08, 0.09, 0.10, 0.11, 0.12, 0.14, 0.16, 0.18, 0.2, 0.22, 0.24, 0.26]
+UVpotVec = [4, 6, 12, 17, 26, 36, 48, 77, 112, 151, 190, 226, 255, 279]
+
+IRampVec = [0.10, 0.12, 0.14, 0.16, 0.18, 0.2, 0.22, 0.24, 0.26, 0.2, 0.08]
+
+
+#%%
+
+from scipy.optimize import curve_fit
+
+RefBins = [t*1e6 for t in SP_Bins[0][:-1]]
+
+plt.figure()
+for Height in [SP_Heigths[3]]:
+    plt.plot(RefBins, Height)
+    
+plt.xlim(0.15, 0.5)
+
+#%%
+
+"""
+Esto mira las curvas SP y les hace ajustes exponenciales a la ultima parte para plotear
+el tiempo caracteristico
+"""
+
+
+"""
+Primer detuning
+"""
+
+Taus = []
+Amps = []
+Offsets = []
+
+ErrorTaus = []
+
+
+popt_vec = []
+pcov_vec = []
+
+plt.figure()
+for j in range(len(SP_Heigths)):
+#for j in [12]:
+    print(j)
+    #BackgroundVector = SP_Bkgr_builder(np.mean(SP_Heigths[j][0:50]),np.mean(SP_Heigths[j][-50:]),59,67,965)
+    
+    #CorrectedSP_Height = [SP_Heigths[j][k]-BackgroundVector[k] for k in range(len(BackgroundVector))]
+    
+    CorrectedSP_Height = SP_Heigths[j]
+    
+    if j==3:
+        popt, pcov = curve_fit(expo, RefBins[90:], CorrectedSP_Height[90:], p0=(5, 1000, 100))
+        popt_vec.append(popt)
+        pcov_vec.append(pcov)
+        plt.plot(RefBins, CorrectedSP_Height)
+        plt.plot(RefBins[90:], [expo(r, *popt) for r in RefBins][90:])
+        Taus.append(popt[0])
+        Amps.append(popt[1])
+        Offsets.append(popt[2])
+        ErrorTaus.append(np.sqrt(pcov)[0][0])
+        print(Taus)        
+#%%
+plt.figure()
+plt.plot(UVpotVec[:-5], Taus[:-6],'o', markersize=10, color='purple')
+#plt.errorbar(UVpotVec[:-5], Taus[:-6], yerr=1e1*np.array(ErrorTaus[:-6]), fmt='.', capsize=2, markersize=2)
+plt.xlabel('UV power (mW)', fontsize=15)
+plt.ylabel('Characteristic time (us)', fontsize=15)
+plt.ylim(-0.1,3)
+plt.grid()
+
+
+#%%
+plt.figure()
+plt.plot(UVpotVec, Amps,'o')
+plt.xlabel('UV power (mW)')
+plt.ylabel('Characteristic time (us)')
+
+plt.figure()
+plt.plot(UVpotVec, Offsets,'o')
+
+
+#%%
+"""
+FIGURA PAPER SP CON AJUSTES
+"""
+
+import matplotlib
+import seaborn as sns
+
+matplotlib.rcParams['mathtext.fontset'] = 'stix'
+matplotlib.rcParams['font.family'] = 'STIXGeneral'
+plt.style.use('seaborn-bright')
+plt.rcParams.update({
+    "text.usetex": False,
+})
+
+plt.figure()
+plt.plot(Stat_Bins[0][:-1], Stat_Heigths[0])
+
+#plot figuras papers
+
+colors1=sns.color_palette("rocket", 10)
+colors2=sns.color_palette("mako", 10)
+color2 = colors2[1]
+color3 = colors2[4]
+color1 = colors2[8]
+
+
+#plt.figure(figsize = (3.8,2.8))
+#plt.plot(RefBins, )
+
+
+
+
+#%%
+
+"""
+Esto mira las curvas SP y les hace ajustes exponenciales a la ultima parte para plotear
+el tiempo caracteristico
+"""
+
+
+"""
+Segundo detuning
+"""
+
+RefBins = [t*1e6 for t in SP_Bins_v2[0][:-1]]
+
+bkgrvec = Calib_Bins[3][:-1]
+
+Taus_v2 = []
+Amps_v2 = []
+Offsets_v2 = []
+
+ErrorTaus_v2 = []
+
+
+popt_vec_v2 = []
+pcov_vec_v2 = []
+
+
+ki = 80 
+
+selectedj = [2,5,10]
+t0 = -0.2
+
+plt.figure()
+for j in range(len(SP_Heigths_v2)):
+#for j in [1]:
+    if j in selectedj:
+        #BackgroundVector = SP_Bkgr_builder(np.mean(SP_Heigths_v2[j][0:50]),np.mean(SP_Heigths_v2[j][-50:]),59,67,965)
+        
+        #CorrectedSP_Height_v2 = [SP_Heigths_v2[j][k]-BackgroundVector[k] for k in range(len(BackgroundVector))]
+        
+        CorrectedSP_Height_v2 = SP_Heigths_v2[j]
+        
+        popt, pcov = curve_fit(expo, RefBins[ki:], CorrectedSP_Height_v2[ki:], p0=(5, 1000, 100), bounds=((0, 0, 0), (50, 1e5, 1e5)))
+        popt_vec_v2.append(popt)
+        pcov_vec_v2.append(pcov)
+        plt.plot([r+t0 for r in RefBins], CorrectedSP_Height_v2)
+        plt.plot([r+t0 for r in RefBins[ki:]], [expo(r, *popt) for r in RefBins][ki:])
+        
+        print(popt[0])
+        
+        print(pcov[0][0])
+        
+        Taus_v2.append(popt[0])
+        Amps_v2.append(popt[1])
+        Offsets_v2.append(popt[2])
+        ErrorTaus_v2.append(np.sqrt(pcov)[0][0])
+    
+plt.xlim(-0.2,2)
+
+#%%
+
+plt.figure()
+plt.plot(UVpotVec, Taus_v2,'o')
+plt.xlabel('UV power (mW)')
+plt.ylabel('Characteristic time (us)')
+
+plt.figure()
+plt.plot(UVpotVec, Amps_v2,'o')
+plt.xlabel('UV power (mW)')
+plt.ylabel('Characteristic time (us)')
+
+plt.figure()
+plt.plot(UVpotVec, Offsets_v2,'o')
+
+#%%
+"""
+FIGURA PAPER
+"""
+import matplotlib
+import seaborn as sns
+
+matplotlib.rcParams['mathtext.fontset'] = 'stix'
+matplotlib.rcParams['font.family'] = 'STIXGeneral'
+plt.style.use('seaborn-bright')
+plt.rcParams.update({
+    "text.usetex": False,
+})
+
+jselected = [1, 5, 12]
+
+t0=0.18
+
+
+s1 = 0
+s2 = 700
+s3 = 2*s2
+
+kini = 75
+
+colors=sns.color_palette("rocket", 10)
+
+color3 = colors[0]
+color2 = colors[3]
+color1 = colors[8]
+
+
+plt.figure(figsize=(3.5,3))
+plt.plot([f*1e6-t0 for f in SP_Bins[jselected[0]][:-1]], [s+s1 for s in SP_Heigths[jselected[0]]],color=color1,alpha=0.7)
+plt.plot([r-t0 for r in RefBins[kini:]], [expo(r, *popt_vec[jselected[0]])+s1 for r in RefBins[kini:]],color=color1)
+plt.plot([f*1e6-t0 for f in SP_Bins[jselected[0]][:-1]], [s+s2 for s in SP_Heigths[jselected[1]]],color=color2,alpha=0.7)
+plt.plot([r-t0 for r in RefBins[kini:]], [expo(r, *popt_vec[jselected[1]])+s2 for r in RefBins[kini:]],color=color2)
+plt.plot([f*1e6-t0 for f in SP_Bins[jselected[0]][:-1]], [s+s3 for s in SP_Heigths[jselected[2]]],color=color3,alpha=0.7)
+plt.plot([r-t0 for r in RefBins[kini:]], [expo(r, *popt_vec[jselected[2]])+s3 for r in RefBins[kini:]],color=color3)
+plt.xlim(-0.3,2)
+plt.xlabel('Time (us)', fontname='STIXGeneral', fontsize=14)
+plt.ylabel('Counts', fontname='STIXGeneral', fontsize=14)
+plt.xticks([0, 0.5, 1, 1.5, 2], fontsize=12)
+plt.yticks([0, 1000, 2000, 3000, 4000], fontsize=12)
+plt.grid()
+plt.savefig('fig3_01.pdf')
+plt.savefig('fig3_01.svg')
+
+
+        #%%
+
+"""
+VEMOS UNA DP CON POTENCIA ALTA A VER SI DA LO QUE TIENE QUE DAR EL ANCHO DE LINEA
+"""
+
+BinsIRhigh = [t*1e6 for t in DP_Bins[-1][:-1]]
+HeightsIRhigh = DP_Heigths[-1]
+
+
+from scipy.optimize import curve_fit
+
+def expo(T, tau, N0, C, T0):
+    return N0*np.exp(-(T-T0)/tau) + C
+
+popt_IR, pcov_IR = curve_fit(expo, BinsIRhigh[int(0.1*len(BinsIRhigh)):], HeightsIRhigh[int(0.1*len(BinsIRhigh)):], p0=(5, 100, 200, 1))
+
+plt.figure()
+plt.plot(BinsIRhigh, HeightsIRhigh)
+plt.plot(BinsIRhigh, [expo(r, *popt_IR) for r in BinsIRhigh])
+
+tauIR = popt_IR[0]
+print(tauIR)
+
+
+#%%
+
+"""
+VEMOS UNA SP CON POTENCIA ALTA A VER SI DA LO QUE TIENE QUE DAR EL ANCHO DE LINEA
+"""
+
+# k=5
+# BinsUVhigh = [t*1e6 for t in Long_Bins[k][:-1]]
+# HeightsUVhigh = Long_Heigths[k]
+
+
+k=-3
+
+BinsUVhigh = [t*1e6 for t in SP_Bins[k][:-1]]
+HeightsUVhigh = SP_Heigths[k]
+
+
+from scipy.optimize import curve_fit
+
+def expo(T, tau, N0, C, T0):
+    return N0*np.exp(-(T-T0)/tau) + C
+
+initi=0.10
+
+popt_UV, pcov_UV = curve_fit(expo, BinsUVhigh[int(initi*len(BinsUVhigh)):], HeightsUVhigh[int(initi*len(BinsUVhigh)):], p0=(1, 1000, 1800, 1), bounds=((0,0,0,0),(20, 1e7, 1e5, 1e3)))
+
+plt.figure()
+plt.plot(BinsUVhigh, HeightsUVhigh,linewidth=4)
+plt.plot(BinsUVhigh[80:], [expo(r, *popt_UV) for r in BinsUVhigh][80:], color='fuchsia', linewidth=3, label='Exponential fit')
+#plt.plot(BinsUVhigh[int(initi*len(BinsUVhigh))], HeightsUVhigh[int(initi*len(BinsUVhigh))],'o',markersize=5)
+#plt.ylim(-1000,20000)
+plt.grid()
+plt.xlabel(r'Time ($\mu$s)', fontsize=15)
+plt.ylabel('Counts', fontsize=15)
+plt.xlim(-0.5, 5)
+plt.legend(fontsize=15)
+
+tauUV = popt_UV[0]
+print(tauUV)
+
+
+
+