Merge branch 'master' of https://code.df.uba.ar/nnunez/artiq_experiments

analisis/plots/20221006_transitoriosv2/Simulaciones/plot_decaytimes_vs_intensity_detunings.py

@@ -11,7 +11,7 @@ I_sat = 43.3*1e-5 #uW/um2
 #%% #####################################################################################
 ### Plot de lineas teoricas tau-intensidad variando detuning + datos medidos
 
-fig2, ax2 = plt.subplots()
+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
@@ -22,11 +22,15 @@ for det in detuning_lines:
 
 for cambio in [0]: # esto es para shiftear la potencia medida
     potencias = dataREAL.Pow - cambio
-    for rad in [85.7]: # esto evalua con distintos radios
+    for rad in [85.7*1.5]: # esto evalua con distintos radios
         ax2.errorbar([p/(np.pi*(rad**2)) for p in UVpotVec], Taus, yerr=np.mean(1e6*errsSIM.stdval.values),
                     fmt='.--', ms=6, lw=.5, \
                     color="#3949ab",
-                    capsize=3)
+                    capsize=2)
+        ax2.errorbar([p/(np.pi*(rad**2)) for p in UVpotVec][4:], Taus_v2[4:], yerr=np.mean(1e6*errsSIM.stdval.values),
+                    fmt='.--', ms=6, lw=.5, \
+                    color="#3949ab",
+                    capsize=2)
     
     
         #ax2.errorbar(potencias/(np.pi*(rad**2)), dataREAL.Tau*1e6, \
@@ -36,13 +40,19 @@ for cambio in [0]: # esto es para shiftear la potencia medida
         #             capsize=3,
         #             label=f'r={rad}; cambio={cambio}')
 
-ax2.set_xlim([1.2e-4, 0.2e-1])
-ax2.set_ylim([2e-1, 0.8e1])
-ax2.set_ylabel(r"$\tau$ [$\mu$s]")
+fs = 9
+
+ax2.set_xlim([0.47e-4, 0.8e-2])
+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$)')
+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")

analisis/plots/20221006_transitoriosv2/Transitorios_02.py

@@ -7,14 +7,17 @@ import ast
 from scipy.optimize import curve_fit
 import os
 
-
+os.chdir('/home/nico/Documents/artiq_experiments/analisis/plots/20221006_transitoriosv2')
 # Solo levanto algunos experimentos
 SP_files = [8445, 8457, 8458, 8459, 8460, 8461, 8462, 8463, 8464, 8465, 8466, 8467, 8468, 8469]
+SP_files_v2 = [8763, 8764, 8765, 8766, 8767, 8768, 8769, 8770, 8771, 8772, 8773, 8774, 8775, 8776]
 DP_files = [8472, 8473, 8474, 8475, 8476, 8477, 8478, 8479, 8480, 8481, 8482]
 
-Long_files = ['DP_long', 'SP_long', 'DP_long_v2', 'DP_long_v3']
+Long_files = ['DP_long', 'SP_long', 'DP_long_v2', 'DP_long_v3', 'DP_long_new', 'SP_long_new']
 Calib_files = ['Cal_DP_5M', 'Cal_SP_5M', 'Cal_DP_25M', 'Cal_SP_25M']
 
+Random_files = [8749]
+
 def expo(T, tau, N0, C):
     global T0
     return N0*np.exp(-(T-T0)/tau) + C
@@ -27,6 +30,23 @@ def pow_from_amp(amp):
     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")
@@ -36,7 +56,7 @@ plt.grid()
 #%%
 
 BINW = 10e-9
-T0 = 0.4e-6
+T0 = -0.4e-6
 
 SP_Heigths = []
 SP_Bins = []
@@ -53,6 +73,22 @@ for i, fname in enumerate(SP_files):
     SP_Heigths.append(heigs)
     SP_Bins.append(binsf)
     
+SP_Heigths_v2 = []
+SP_Bins_v2 = []
+
+for i, fname in enumerate(SP_files_v2):
+    #print(i)
+    #print(fname)
+    data = h5py.File('Data/SP/00000'+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_v2.append(heigs)
+    SP_Bins_v2.append(binsf)    
+
+    
 DP_Heigths = []
 DP_Bins = []
 
@@ -69,8 +105,9 @@ for i, fname in enumerate(DP_files):
     DP_Bins.append(binsf)    
 
 
+
 BINW_long = 10e-9
-T0_long = -0.5e-6
+T0_long = -0.4e-6
 
 Long_Heigths = []
 Long_Bins = []
@@ -87,6 +124,10 @@ for i, fname in enumerate(Long_files):
     Long_Heigths.append(heigs)
     Long_Bins.append(binsf) 
 
+
+BINW_calib = 10e-9
+T0_calib = -0.4e-6
+
 Calib_Heigths = []
 Calib_Bins = []
 
@@ -95,13 +136,31 @@ for i, fname in enumerate(Calib_files):
     #print(fname)
     data = h5py.File('Data/Calibrations/'+str(fname)+'.h5', 'r')
     counts = np.array(data['datasets']['counts'])
-    bines = np.arange(counts.min(), counts.max()+BINW_long, BINW_long)
+    bines = np.arange(counts.min(), counts.max()+BINW_calib, BINW_calib)
 
-    heigs, binsf  = np.histogram(counts, bines[bines>T0_long])
+    heigs, binsf  = np.histogram(counts, bines[bines>T0_calib])
     
     Calib_Heigths.append(heigs)
     Calib_Bins.append(binsf) 
 
+BINW_Random = 10e-9
+T0_Random = -0.5e-6
+
+
+Random_Heigths = []
+Random_Bins = []
+
+for i, fname in enumerate(Random_files):
+    #print(i)
+    #print(fname)
+    data = h5py.File('Data/Random/00000'+str(fname)+'-SingleLine.h5', 'r')
+    counts = np.array(data['datasets']['counts'])
+    bines = np.arange(counts.min(), counts.max()+BINW_Random, BINW_Random)
+
+    heigs, binsf  = np.histogram(counts, bines[bines>T0_Random])
+    
+    Random_Heigths.append(heigs)
+    Random_Bins.append(binsf) 
 
 
 #plt.figure()
@@ -137,6 +196,11 @@ Esto mira las curvas SP y les hace ajustes exponenciales a la ultima parte para
 el tiempo caracteristico
 """
 
+
+"""
+Primer detuning
+"""
+
 bkgrvec = Calib_Bins[3][:-1]
 
 Taus = []
@@ -145,19 +209,30 @@ Offsets = []
 
 ErrorTaus = []
 
+
+popt_vec = []
+pcov_vec = []
+
 plt.figure()
 for j in range(len(SP_Heigths)):
+#for j in [12]:
+    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]
     
-    popt, pcov = curve_fit(expo, RefBins, CorrectedSP_Height, p0=(5, 1000, 100))
+    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, [expo(r, *popt) for r in RefBins])
+    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])
     
+#%%
 plt.figure()
 plt.plot(UVpotVec, Taus,'o')
 plt.xlabel('UV power (mW)')
@@ -169,39 +244,176 @@ plt.xlabel('UV power (mW)')
 plt.ylabel('Characteristic time (us)')
 
 plt.figure()
-plt.plot(UVampVec, Offsets,'o')
+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
+
+plt.figure()
+for j in range(len(SP_Heigths_v2)):
+#for j in [1]:
+    #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(RefBins, CorrectedSP_Height_v2)
+    plt.plot(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.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
+"""
 
+jselected = [1, 5, 12]
+
+plt.figure()
+plt.plot([f*1e6 for f in SP_Bins[jselected[0]][:-1]], SP_Heigths[jselected[0]])
+plt.plot(RefBins[90:], [expo(r, *popt_vec[jselected[0]]) for r in RefBins[90:]])
+plt.plot([f*1e6 for f in SP_Bins[jselected[0]][:-1]], SP_Heigths[jselected[1]])
+plt.plot(RefBins[90:], [expo(r, *popt_vec[jselected[1]]) for r in RefBins[90:]])
+plt.plot([f*1e6 for f in SP_Bins[jselected[0]][:-1]], SP_Heigths[jselected[2]])
+plt.plot(RefBins[90:], [expo(r, *popt_vec[jselected[2]]) for r in RefBins[90:]])
+plt.xlim(-0.5,4)
 
 #%%
 """
 Defino los vectores de las mediciones largas para obtener branching fractions, 25 M de mediciones
 """
 
-Bins_long = [t*1e6 for t in Long_Bins[0][:-1]]
-DP_long_raw = Long_Heigths[3] #la 0 y la 2 salieron mal pero las dejo por las dudas. uso la 3
-SP_long_raw = Long_Heigths[1]
+#Bins_long = [t*1e6 for t in Long_Bins[0][:-1]]
+Bins_long = [t*1e6 for t in Long_Bins[-1][:-1]]
+#DP_long_raw = Long_Heigths[3] #la 0 y la 2 salieron mal pero las dejo por las dudas. uso la 3
+#SP_long_raw = Long_Heigths[1]
+
+DP_long_raw = Long_Heigths[-2] #la 0 y la 2 salieron mal pero las dejo por las dudas. uso la 3
+SP_long_raw = Long_Heigths[-1]
+
 
 DP_bkg_long_raw = Calib_Heigths[2] #calibracion del DP, esta al final ni la uso
 SP_bkg_long_raw = Calib_Heigths[3] #calibracion del SP
 
+
+#Genero el vector de background para el SP
+BackgroundVectorLong = SP_Bkgr_builder(np.mean(SP_long_raw[0:50]),np.mean(SP_long_raw[-50:]),59,67,965)
+
+SP_long_corrected = [SP_long_raw[k]-BackgroundVectorLong[k] for k in range(len(SP_long_raw))]
+DP_long_corrected = [DP_long_raw[k]-np.mean(DP_long_raw[-100:]) for k in range(len(DP_long_raw))]
+
+plt.figure()
+plt.plot(SP_long_raw)
+plt.plot(BackgroundVectorLong,'o')
+plt.xlim(0,100)
 #%%
 
 """
 Ploteo las mediciones largas
 """
 plt.figure()
-plt.plot(Bins_long, SP_long_raw)
-plt.plot(Bins_long[0:475], DP_long_raw)
-plt.xlim(-1, 2)
+plt.plot(Bins_long, SP_long_corrected)
+plt.plot([b+0.25 for b in Bins_long], DP_long_corrected)
+plt.grid()
+#plt.xlim(0, 1.5)
+#plt.ylim(-100,100)
 
+stotal = np.sum(SP_long_corrected)
+dtotal = np.sum(DP_long_corrected)
+print(stotal/(dtotal+stotal))
 
 #%%
 """
 aca se plotea la SP con su curva de calibracion sin ion para restar
 """
 plt.figure()
-plt.plot(Bins_long, SP_long_raw)
+#plt.plot(Bins_long, SP_long_raw)
 plt.plot(Bins_long, SP_bkg_long_raw)
 plt.xlim(-1, 2)
 #%%
@@ -211,16 +423,34 @@ apagado porque en ningun momento esta prendido el UV
 """
 #Pruebo restar background directamente
 #DP_long_subs = [DP_long_raw[j]-DP_bkg_long_raw[j] for j in range(len(DP_long_raw))]
-DP_long_subs = [d-np.mean(DP_long_raw[-100:]) for d in DP_long_raw]
+DP_long_subs = [d-int(np.mean(DP_long_raw[-100:])) for d in DP_long_raw]
 SP_long_subs = [SP_long_raw[j]-SP_bkg_long_raw[j] for j in range(len(SP_long_raw))]
 
+stot = np.sum(SP_long_subs)
+dtot = np.sum(DP_long_subs)
+
+prob = stot/(stot+dtot)
+
+eff_deteccion = dtot/25e6
+
+print(f'p={prob}')
+print(f'Eff de deteccion: {eff_deteccion*100}')
+
 #%%
 plt.figure()
+plt.plot(Bins_long, SP_long_subs)
 plt.plot([b+0.25 for b in Bins_long[0:475]], DP_long_subs)
-#plt.plot(Bins_long, SP_long_subs)
-#plt.xlim(-0.2, 2)
+plt.xlim(-0.2, 2)
 plt.grid()
 
+#%%
+
+plt.figure()
+plt.plot(Random_Bins[0][:-1], Random_Heigths[0],'o')
+
+
+
+
 
 
 

analisis/plots/20221017_IonStatistics/IonStatistics.py

@@ -9,11 +9,8 @@ import os
 import scipy.stats as sts
 import seaborn as sns
 
-colors1=sns.color_palette("rocket", 10)
-colors2=sns.color_palette("mako", 10)
-color1 = colors1[5]
-color2 = colors1[3]
-color3 = colors1[1]
+
+os.chdir('/home/nico/Documents/artiq_experiments/analisis/plots/20221017_IonStatistics')
 
 plt.rcParams.update({
     "font.family": "STIXGeneral"
@@ -22,7 +19,7 @@ plt.rcParams.update({
     "font.size": 14
 })
 # Solo levanto algunos experimentos
-Stat_files = [8731, 8738, 8745]
+Stat_files = [8731, 8738, 8745, 8819, 8820, 8822]
 
 def expo(T, tau, N0, C):
     global T0
@@ -50,7 +47,10 @@ T0 = 0.0e-6
 Stat_Heigths = []
 Stat_Bins = []
 
-for i, fname in enumerate(Stat_files):
+OnOff_Heights = []
+OnOff_Bins = []
+
+for i, fname in enumerate(Stat_files[0:3]):
     #print(i)
     #print(fname)
     data = h5py.File('Data/00000'+str(fname)+'-IonStatistics.h5', 'r')
@@ -62,42 +62,125 @@ for i, fname in enumerate(Stat_files):
     Stat_Heigths.append(heigs)
     Stat_Bins.append(binsf)
     
+for i, fname in enumerate(Stat_files[3:6]):
+    #print(i)
+    #print(fname)
+    try:
+        data = h5py.File('Data/00000'+str(fname)+'-SingleLine.h5', 'r')
+    except:
+        data = h5py.File('Data/00000'+str(fname)+'-StationaryFluo.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])
+    
+    OnOff_Heights.append(heigs)
+    OnOff_Bins.append(binsf)
+
+    
 #%%
+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
 
-bins1 = np.arange(100,350, 1)
-bins2 = np.arange(10,50,1)
+colors1=sns.color_palette("rocket", 10)
+colors2=sns.color_palette("mako", 10)
+color2 = colors2[1]
+color3 = colors2[4]
+color1 = colors2[8]
+
+bins1 = np.arange(150,300, 1)
+bins2 = np.arange(0,50,1)
 bins3 = np.arange(30,100,1)
 
-plt.figure(figsize = (4.5,3))
-plt.hist(Stat_Heigths[0], bins=bins1, histtype='step',density = True,color = color1)#,label = 'BG')
-#plt.hist(Stat_Heigths[1], bins=bins2, histtype='step',density = True)
-plt.hist(Stat_Heigths[1], bins=bins2, histtype='step',density = True,color = color2)#,label = 'UV laser')
-#plt.hist(Stat_Heigths[3], bins=bins2, histtype='step',density = True)
-#plt.hist(Stat_Heigths[4], bins=bins2, histtype='step',density = True)
-plt.hist(Stat_Heigths[2], bins=bins3, histtype='step',density = True,color = color3)#,label = 'Ion')
+
+
+
+"""
+FIGURA ESTADISTICA POISSON
+"""
+
+plt.figure(figsize = (3.8,2.8))
+plt.hist(Stat_Heigths[0], bins=bins1, histtype='step',density = True,color = color1, alpha = 0.6)#,label = 'BG')
+plt.hist(Stat_Heigths[1], bins=bins2, histtype='step',density = True,color = color2, alpha = 0.6)#,label = 'UV laser')
+plt.hist(Stat_Heigths[2], bins=bins3, histtype='step',density = True,color = color3, alpha = 0.6)#,label = 'Ion')
 
 
 poisson1= sts.poisson.pmf(bins1,np.mean(Stat_Heigths[0]))
 poisson2= sts.poisson.pmf(bins2,np.mean(Stat_Heigths[1]))
 poisson3= sts.poisson.pmf(bins3,np.mean(Stat_Heigths[2]))
 
-
-plt.plot(bins2+0.5,poisson2,color = color2,alpha = 0.9,label = 'BG')
-plt.plot(bins3+0.5,poisson3,color = color3,alpha = 0.9,label = 'UV laser')
-plt.plot(bins1+0.5,poisson1,color = color1,alpha = 0.9,label = 'Ion')
-plt.legend()
+plt.plot(bins2+0.5,poisson2,color = color2,label = 'Background')
+plt.plot(bins3+0.5,poisson3,color = color3,label = 'UV laser')
+plt.plot(bins1+0.5,poisson1,color = color1,label = 'Ion')
+plt.legend(loc=(0.15,0.65), prop={'size': 11})
 plt.grid()
 plt.tight_layout()
-plt.xlabel('Counts')
-plt.ylabel('Event frequency')
+plt.xticks([0, 100, 200, 300], fontname='STIXGeneral')
+plt.yticks([0,0.02, 0.04, 0.06, 0.08], fontname='STIXGeneral')
+plt.xlabel('Counts', fontname='STIXGeneral')
+plt.ylabel('Event frequency', fontname='STIXGeneral')
 
 #poissoneidad =  np.var(Stat_Heigths)/np.mean(Stat_Heigths)
 #plt.title('Varianza/media = {:.4f}'.format(poissoneidad))
-plt.savefig('bg_laser_ion_stats.pdf',dpi = 600 )
+#plt.savefig('bg_laser_ion_stats.pdf',dpi = 600 )
+
+name='fig01a'
+
+plt.savefig('/home/nico/Nextcloud/G_liaf/Publicaciones/Papers/2022 Transient Phenomena JOSA B/Figures/'+name+'.pdf')
+plt.savefig('/home/nico/Nextcloud/G_liaf/Publicaciones/Papers/2022 Transient Phenomena JOSA B/Figures/'+name+'.svg')
+
+
+
+#%%
+import matplotlib
+
+matplotlib.rcParams['mathtext.fontset'] = 'stix'
+matplotlib.rcParams['font.family'] = 'STIXGeneral'
+plt.style.use('seaborn-bright')
+plt.rcParams.update({
+    "text.usetex": False,
+})
+
+"""
+FIGURA BACKGROUND, ENCENDIDO AOM Y ENCENDIDO ION
+"""
+
+plt.figure(figsize=(4.5,3))
+plt.plot([s*1e6 for s in OnOff_Bins[1][:-1]],OnOff_Heights[1], color=color2)
+plt.plot([s*1e6 for s in OnOff_Bins[0][:-1]],[(OnOff_Heights[0][j]-OnOff_Heights[1][j])*3.13+OnOff_Heights[1][j] for j in range(len(OnOff_Heights[0]))], color=color3)
+plt.plot([s*1e6 for s in OnOff_Bins[2][:-1]],OnOff_Heights[2], color=color1)
+plt.xlim(0,1)
+plt.grid()
+plt.xlabel(r'Time ($\mu$s)', fontname='STIXGeneral')
+plt.ylabel(r'Counts /10$~\mu$s', fontname='STIXGeneral')
+plt.xticks([0,0.2, 0.4, 0.6, 0.8, 1], fontname='STIXGeneral')
+plt.yticks([0,5000, 10000], fontname='STIXGeneral')
+plt.tight_layout()
+
+name='fig01b'
+
+plt.savefig('/home/nico/Nextcloud/G_liaf/Publicaciones/Papers/2022 Transient Phenomena JOSA B/Figures/'+name+'.pdf')
+plt.savefig('/home/nico/Nextcloud/G_liaf/Publicaciones/Papers/2022 Transient Phenomena JOSA B/Figures/'+name+'.svg')
+
+
+
+
+
+
+
 
 

analisis/plots/20221024_BranchingFractionMeasurement/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/20221024_BranchingFractionMeasurement/Simulaciones/OBE3levels_functions.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri May 15 18:11:00 2020
+
+@author: oem
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from odeintw import odeintw
+import time
+#matplotlib.use('Qt5Agg')
+
+
+"""
+odeintw es un wrapper de la funcion odeint de scipy.integrate que sirve para integrar
+sistemas de ecuaciones con parámetros complejos.
+En este script se resuelven las ecuaciones
+de Bloch de un sistema de 3 niveles en configuración lambda (1 y 3 son ground, 2 es excitado)
+para ion de Bario o de Calcio.
+Este script considera que hay dos láseres repump prendidos con sus respectivas frecuencias de
+rabi y detuning.
+Además, se usan dos distintos sistemas de referencia rotantes, condensados en las funciones
+RepumpReferenceSystemModel 1 y 2. La 2 es más útil para prender de a poco el láser repump2, y
+cuando su frec. de rabi es cero, no aparece en las ecuaciones.
+https://github.com/WarrenWeckesser/odeintw
+
+
+"""
+
+
+#CHEQUEADOS LOS DAMPINGS
+
+
+def EquationSystemModel(z, t, rabi12, rabi23, g21, g23, glg, glr, detr, detg, no_IR_lasers):
+    """
+    Se plantean acá las ecuaciones dinámicas de la matriz densidad (OBE: Optical Bloch Equations)
+    para la interacción entre un átomo y un campo electromagnético.
+        -rabi12, rabi23: frecuencias de rabi, proporcionales a la intensidad de los láseres
+        -g21, g23: anchos de línea naturales de las transiciones
+        -glg, glr: anchos de línea de los láseres doppler (UV) y repump (IR) respectivamente
+        -detg, detr: detunings de transición doppler (UV) y repump (IR) respectivamente
+        -no_IR_lasers: número de láseres IR aplicados
+    """
+
+    rho33, rho22, rho11, rho12, rho13, rho23 = z
+
+    rho21 = np.conjugate(rho12)
+    rho31 = np.conjugate(rho13)
+    rho32 = np.conjugate(rho23)
+
+    drho11dt = (1j)*0.5*rabi12*(rho12 - rho21) + rho22*g21
+    drho22dt = (-1j)*0.5*rabi12*(rho12 - rho21) - 1j*0.5*(np.conjugate(rabi23)*rho32 - rabi23*rho23) - rho22*(g21 + g23)
+    drho33dt = (-1j)*0.5*(rabi23*rho23 - np.conjugate(rabi23)*rho32) + g23*rho22
+
+    drho12dt = (-1j)*0.5*rabi12*(rho22 - rho11) - (1j)*(detg*rho12 - 0.5*rabi23*rho13) - rho12*(0.5*g21 + 0.5*g23 + glg)
+    drho13dt = (-1j)*(detg*rho13 + 0.5*rabi12*rho23 - 0.5*np.conjugate(rabi23)*rho12 - detr*rho13) - rho13*(glg + no_IR_lasers*glr)
+    drho23dt = (-1j)*0.5*rabi12*rho13 - (1j)*0.5*np.conjugate(rabi23)*(rho33-rho22) + (1j)*detr*rho23 - rho23*(0.5*g21 + 0.5*g23 + no_IR_lasers*glr)
+
+    dzdt = [drho33dt, drho22dt, drho11dt, drho12dt, drho13dt, drho23dt]
+
+    return dzdt
+
+
+def RepumpReferenceSystemModel(rabi23repump1, rabi23repump2, Detuningrepump1, Detuningrepump2, t):
+    """
+    Cuando hay más de un láser en una transición se producen batidos entre los mismos, haciendo que
+    la frecuencia de rabi sea un parámetro que varía en el tiempo. Para resolver eso numéricamente se
+    implementa de esta manera, calculando la frecuencia de rabi efectiva antes de resolver las OBE.
+    Al resolver el sistema, se pasa a un marco de referencia rotante respecto a los láseres. Como los dos láseres
+    IR excitan la misma transición, hay que elegir una de ambas frecuencias para rotar. Esto elige al repump1
+    como tal referencia, pero es posible elegir a repump2, o incluso a un promedio entre ambas frecuencias (por eso se llama detuningmedio).
+    Es decir, esto considera a la matriz de rotación como:
+        U = exp(-i.wg.t)|1><1| + |2><2| + exp(-i.wr1.t)|3><3|
+    """
+
+    rabi23 = [rabi23repump1 + rabi23repump2*np.exp(-1*(Detuningrepump1-Detuningrepump2)*1j*tiemp) for tiemp in t]
+    DetuningMedioRepump = Detuningrepump1
+    return rabi23, DetuningMedioRepump
+
+
+def GetTransitionLinewidths(species):
+    """
+    Devuelve los anchos de línea naturales según la especie iónica
+    """
+
+
+    if species == 'Calcium':
+        g21, g23, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 #anchos de linea de las transiciones
+    elif species == 'Barium':
+        g21, g23, = 2*np.pi*15.1e6, 2*np.pi*5.3e6 #anchos de linea de las transiciones
+    else:
+        print('Species not implemented')
+        raise
+
+    return g21, g23
+
+
+def Solve3LevelBloch(EquationSystemModel, species, SatParDoppler, SatParRep1, SatParRep2, DetDoppler, DetRep1, DetRep2,
+                     glg, glr, no_IR_lasers, initcond=1, tfin=10, paso=1e-3, printprogress=False):
+    """
+    Resuelve las ecuaciones de bloch del modelo dado. Los inputs son:
+        glg, glr: anchos de línea de doppler y repump en MHz
+        SatPar: parámetro de saturación de cada láser
+        Det: detunings en MHz
+        rabi12, rabi23 = sr*g21, sp*g23 #frecuencias de rabi
+    """
+
+    t = np.arange(0, tfin*1e-6, paso*1e-6)
+
+    g21, g23 = GetTransitionLinewidths(species)
+
+
+    #convierto todo a MHz
+    DopplerLaserLinewidth = glg*2*np.pi*1e6
+    RepumpLaserLinewidth = glr*2*np.pi*1e6
+
+    DetuningDoppler = DetDoppler*2*np.pi*1e6
+    Detuningrepump1 = DetRep1*2*np.pi*1e6
+    Detuningrepump2 = DetRep2*2*np.pi*1e6
+
+    rabi12 = np.sqrt(SatParDoppler*(4*(DetuningDoppler**2) + g21**2)/2)
+    rabi23repump1 = np.sqrt(SatParRep1*(4*(Detuningrepump1**2) + g23**2)/2)
+    rabi23repump2 = np.sqrt(SatParRep2*(4*(Detuningrepump2**2) + g23**2)/2)
+
+    rabi23vst, _ = RepumpReferenceSystemModel(rabi23repump1, rabi23repump2, Detuningrepump1, Detuningrepump2, t)
+
+    rho33sol = []
+    rho22sol = []
+    rho11sol = []
+
+    rho12sol = []
+    rho13sol = []
+    rho23sol = []
+
+
+    #condicion inicial: poblacion en el nivel 1
+    if initcond==1:
+        z0 = [0+0j, 0+0j, 1+0j, 0+0j, 0+0j, 0+0j] #formato: [rho33, rho22, rho11, rho12, rho13, rho23]
+    elif initcond==2:
+        z0 = [0+0j, 1+0j, 0+0j, 0+0j, 0+0j, 0+0j]
+    elif initcond==3:
+        z0 = [1+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j]
+    else:
+        print('Mala condicion inicial')
+        raise
+
+    #z0 = [0.5+0j, 0+0j, 0.5+0j, -0.5+0j, 0+0j, 0+0j] DARK STATE
+
+    for i in range(1,len(t)):
+        # span for next time step
+        tspan = [t[i-1],t[i]]
+        # solve for next step
+        #z = odeintw(EquationSystemModel, z0, tspan, args=(rabi12, rabi23vst[i], g21, g23, DopplerLaserLinewidth, RepumpLaserLinewidth, DetuningMedioRepump, DetuningDoppler), full_output=True)
+        z = odeintw(EquationSystemModel, z0, tspan, args=(rabi12, rabi23repump1, g21, g23, DopplerLaserLinewidth, RepumpLaserLinewidth, Detuningrepump1, DetuningDoppler, no_IR_lasers), full_output=True)
+        # store solution for plotting
+        rho22sol.append(z[0][1][1])
+        rho33sol.append(z[0][1][0])
+        rho11sol.append(z[0][1][2])
+
+        rho12sol.append(z[0][1][3])
+        rho13sol.append(z[0][1][4])
+        rho23sol.append(z[0][1][5])
+
+        # next initial condition
+        z0 = z[0][1]
+        if printprogress:
+            print(i, '/', len(t), end="\r", flush=True)
+
+    #Tomo la parte real de los diagonales y el módulo de las coherencias
+    RealRho22 = [r.real for r in rho22sol]
+    RealRho11 = [r.real for r in rho11sol]
+    RealRho33 = [r.real for r in rho33sol]
+    AbsRho12 = [np.abs(r) for r in rho12sol]
+    AbsRho13 = [np.abs(r) for r in rho13sol]
+    AbsRho23 = [np.abs(r) for r in rho23sol]
+
+    return t, RealRho11, RealRho22, RealRho33, AbsRho12, AbsRho13, AbsRho23
+
+
+
+def RepumpSpectrum(EquationSystemModel, species, SatParDoppler, SatParRep1, SatParRep2, DetDoppler, DetRep2,
+                     glg, glr, tfinOBE=10, pasoOBE=1e-3, frecini=-100, frecfin=100, pasofrec=1, plotSpectrum=True, sweepingrepump=True):
+    """
+    Hace un barrido en la frecuencia del repump.
+    Considerando que puede haber oscilaciones de la población por haber dos repumps,
+    """
+
+    g21, g23 = GetTransitionLinewidths(species)
+
+    FluorescenceVector = []
+
+    t = np.arange(0, tfinOBE*1e-6, pasoOBE*1e-6)
+
+    DopplerLaserLinewidth = glg*2*np.pi*1e6
+    RepumpLaserLinewidth = glr*2*np.pi*1e6
+
+    DetuningDoppler = DetDoppler*2*np.pi*1e6
+    Detuningrepump2 = DetRep2*2*np.pi*1e6
+
+    rabi12 = np.sqrt(SatParDoppler*(4*(DetuningDoppler**2) + g21**2)/2)
+    rabi23repump1 = np.sqrt(SatParRep1*(4*(Detuningrepump1**2) + g23**2)/2)
+    rabi23repump2 = np.sqrt(SatParRep2*(4*(Detuningrepump2**2) + g23**2)/2)
+
+    Repump1DetuningVector = np.arange(frecini*2*np.pi*1e6, frecfin*2*np.pi*1e6, pasofrec*2*np.pi*1e6)
+
+    counter = 1
+
+    FluorescenceVector = []
+
+    for DetRep1 in Repump1DetuningVector:
+        tini = time.time()
+
+        rho33sol = []
+        rho22sol = []
+        rho11sol = []
+
+        rho12sol = []
+        rho13sol = []
+        rho23sol = []
+
+        rabi23vst, DetuningMedioRepump = RepumpReferenceSystemModel(rabi23repump1, rabi23repump2, DetRep1, Detuningrepump2, t)
+
+        z0 = [0+0j, 0+0j, 1+0j, 0+0j, 0+0j, 0+0j] #formato: [rho33, rho22, rho11, rho12, rho13, rho23]
+
+
+        for i in range(1,len(t)):
+            # span for next time step
+            tspan = [t[i-1],t[i]]
+            # solve for next step
+            z = odeintw(EquationSystemModel,z0,tspan, args=(rabi12, rabi23vst[i], g21, g23, DopplerLaserLinewidth, RepumpLaserLinewidth, DetuningMedioRepump, DetuningDoppler), full_output=True)
+            # store solution for plotting
+
+            rho22sol.append(z[0][1][1])
+            rho33sol.append(z[0][1][0])
+            rho11sol.append(z[0][1][2])
+
+            rho12sol.append(z[0][1][3])
+            rho13sol.append(z[0][1][4])
+            rho23sol.append(z[0][1][5])
+
+            # next initial condition
+            z0 = z[0][1]
+    #        print(i, '/', len(t))
+        print
+
+        FluorescenceVector.append(np.real(rho12sol[-1]))
+
+        tfin = time.time()
+        print(counter, '/', len(Repump1DetuningVector), ', elapsed time: ', round(tfin-tini, 2), ' s')
+        counter = counter + 1
+
+    if plotSpectrum:
+        plt.plot([1e-6*d/(2*np.pi) for d in Repump1DetuningVector], [100*f for f in FluorescenceVector], 'o')
+        plt.xlabel('Repump detuning (MHz')
+        plt.ylabel('Fluorescence (A.U.)')
+        plt.title('Doppler detuning: ' + str(DetDoppler) + ' MHz, Repump 2 detuning: ' + str(DetRep2) + ' MHz')
+
+    return Repump1DetuningVector, FluorescenceVector
+
+
+#Funciones auxiliares para el análisis
+def FindFluoAtTime(fluo, time, timetarget, pasoti):
+    """
+    Encuentra la fluorescencia para determinado tiempo (en us).
+    El pasoti es la tolerancia temporal para el cual busca esa fluorescencia
+    """
+    ti = timetarget*1e-6
+    pasoti = pasoti*1e-6
+    i = 0
+    fluoattivec = []
+    while i < len(fluo):
+        if abs(time[i]-ti) < 2*pasoti:
+            fluoattivec.append(fluo[i])
+        i = i + 1
+    return np.mean(fluoattivec)
+
+def FindTimeForFluo(fluo, time, fluotarget, tolerancefluo=1e-3):
+    i = 0
+    timesforfluo = []
+    while i < len(time)-1:
+        if np.abs(fluo[i]-fluotarget) < 2*tolerancefluo:
+            timesforfluo.append(time[i])
+        i = i + 1
+    if len(timesforfluo) == 0:
+        print('La fluorescencia no toma ese valor')
+        return 0
+    elif len(timesforfluo) == 1:
+        return timesforfluo[0]
+    else:
+#        return np.mean(timesforfluo[1:])
+        return timesforfluo[-1]
+
+def FindCharacteristicDecayTime(fluo, time):
+    pass
+
+
+def IntegrateWindow(fluo, time, ti, window):
+    ti = ti*1e-6
+    window = window*1e-6
+    i = 0
+    fluowindow = []
+    while i<len(fluo):
+        if time[i]>ti and time[i]<=ti+window:
+            fluowindow.append(fluo[i])
+        i = i + 1
+        if time[i]>ti+window:
+            return np.mean(fluowindow)
+
+
+def Lorentz(x, x0, gamma, ymax):
+    scale = ymax*np.pi*gamma*0.5
+    return (scale/np.pi)*(0.5*gamma)/((x-x0)**2 + (0.5*gamma)**2)
+
+def rho23model(rabip, gamma, gammas, detuningp, rho0, t):
+    rho23 = ((1j)*rho0*rabip/(2*gamma))*(np.exp(-gamma*t) - np.exp(-(gammas+(rabip**2)/(4*gamma))*t -(1j)*detuningp*t))
+    return [np.real(rho23), np.abs(np.imag(rho23))]
+
+
+def dopplerBroadening(T, wldopp=0.397e-6, wlrep=0.866e-6, alpha=0, mcalcio = 6.655e-23*1e-3):
+    """
+    Calcula el broadening extra semiclásico por temperatura considerando que el ion atrapado se mueve.
+    wlr es la longitud de onda doppler, wlp la longitud de onda repump, T la temperatura del ion en kelvin, y alpha (en rads) el ángulo
+    que forman ambos láseres.
+    """
+
+    kboltzmann = 1.38e-23 #J/K
+
+    gammaD = (2*np.pi)*np.sqrt((1/(wldopp*wldopp)) + (1/(wlrep*wlrep)) - 2*(1/(wldopp*wlrep))*np.cos(alpha))*np.sqrt(kboltzmann*T/(2*mcalcio))
+
+    return gammaD
+

analisis/plots/20221024_BranchingFractionMeasurement/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/20221024_BranchingFractionMeasurement/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/20221024_BranchingFractionMeasurement/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
+
+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, 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 [85.7*1.5]: # esto evalua con distintos radios
+        ax2.errorbar([p/(np.pi*(rad**2)) for p in UVpotVec], Taus, yerr=np.mean(1e6*errsSIM.stdval.values),
+                    fmt='.--', ms=6, lw=.5, \
+                    color="#3949ab",
+                    capsize=2)
+        ax2.errorbar([p/(np.pi*(rad**2)) for p in UVpotVec][4:], Taus_v2[4:], yerr=np.mean(1e6*errsSIM.stdval.values),
+                    fmt='.--', ms=6, lw=.5, \
+                    color="#3949ab",
+                    capsize=2)
+    
+    
+        #ax2.errorbar(potencias/(np.pi*(rad**2)), dataREAL.Tau*1e6, \
+        #             yerr = 1e6*errsSIM.stdval.values, \
+        #             fmt='.--', ms=6, lw=.5, \
+        #             color="#3949ab",
+        #             capsize=3,
+        #             label=f'r={rad}; cambio={cambio}')
+
+fs = 9
+
+ax2.set_xlim([0.47e-4, 0.8e-2])
+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")
+
+

analisis/plots/20221024_BranchingFractionMeasurement/Transitorios_BF.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/20221024_BranchingFractionMeasurement')
+
+Long_files = ['DP_long', 'SP_long']
+Calib_files = ['Fondo_IR_50M', 'Fondo_UV_50M']
+
+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 = 10e-9
+T0 = -0.4e-6
+
+
+BINW_long = 10e-9
+T0_long = -0.4e-6
+
+Long_Heigths = []
+Long_Bins = []
+
+for i, fname in enumerate(Long_files):
+    #print(i)
+    #print(fname)
+    data = h5py.File('Data/Largas/'+str(fname)+'.h5', 'r')
+    counts = np.array(data['datasets']['counts'])
+    bines = np.arange(counts.min(), counts.max()+BINW_long, BINW_long)
+
+    heigs, binsf  = np.histogram(counts, bines[bines>T0_long])
+    
+    Long_Heigths.append(heigs)
+    Long_Bins.append(binsf) 
+
+
+BINW_calib = 10e-9
+T0_calib = -0.4e-6
+
+Calib_Heigths = []
+Calib_Bins = []
+
+for i, fname in enumerate(Calib_files):
+    #print(i)
+    #print(fname)
+    data = h5py.File('Data/Calibrations/'+str(fname)+'.h5', 'r')
+    counts = np.array(data['datasets']['counts'])
+    bines = np.arange(counts.min(), counts.max()+BINW_calib, BINW_calib)
+
+    heigs, binsf  = np.histogram(counts, bines[bines>T0_calib])
+    
+    Calib_Heigths.append(heigs)
+    Calib_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]
+
+
+#%%
+
+RefBins = [t*1e6 for t in Long_Bins[0][:-1]]
+
+plt.figure()
+
+#for i in range(len(Long_Heigths)):
+#    plt.plot(Long_Bins[i][:-1], Long_Heigths[i])
+    
+for i in range(len(Calib_Heigths)):
+    plt.plot(Calib_Bins[i][:-1], Calib_Heigths[i])
+    
+#%%
+"""Calculo branching"""
+
+TotalDP = np.sum(Long_Heigths[0])
+TotalSP = np.sum(Long_Heigths[1])
+
+TotalDPbkg = np.sum(Calib_Heigths[0])
+TotalSPbkg = np.sum(Calib_Heigths[1])
+
+TotalDPfinal = TotalDP-TotalDPbkg
+TotalSPfinal = TotalSP-TotalSPbkg
+
+branch = TotalSPfinal/(TotalSPfinal+TotalDPfinal)
+
+print(branch)
+
+
+
+
+
+