cambios sobre el codigo para fitear las mediciones variando la temperatura:...

cambios sobre el codigo para fitear las mediciones variando la temperatura: agregue meterla con la velocidad y con sidebands. La estructura de fiteo es distinta a la de nico, está todo en cpt plotter

analisis/plots/20230920_CPT_TemperatureSens_v2/CPT_plotter_20230920.py

@@ -110,7 +110,7 @@ plt.xlabel('Frecuencia (MHz)')
 plt.ylabel('counts')
 plt.grid()
 plt.legend()
-
+plt.savefig('Ejemplo_CPT_temperaturas_distintas.png')
 #%%
 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

analisis/plots/20230920_CPT_TemperatureSens_v2/Data/EITfit/MM_eightLevel_2repumps_AnalysisFunctions.py

@@ -17,20 +17,63 @@ from MM_eightLevel_2repumps_python_scripts import CPTspectrum8levels,CPTspectrum
 from scipy.optimize import curve_fit
 import random
 from scipy.signal import savgol_filter as sf
+from scipy.stats import norm
 
-
-def PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, beta, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None):
+def PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, betag, betap, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None, dephasing  = False, boltzmann = False):
     """
     solvemode=1: resuelve con np.linalg.solve
     solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
     """
-    tinicial = time.time()
-    ProbeDetuningVectorL, Fluovector = CPTspectrum8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, circularityprobe, beta, drivefreq, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=False, solvemode=1)
-    tfinal = time.time()
+    kboltz = 1.380649e-23
+    m = 6.6551079e-26
+    kg = np.pi * 2/397 *1e9
+    kp = np.pi * 2/866 *1e9
+    if detpvec is None:
+        detpvec = np.arange(freqMin,freqMax,freqStep)    
+
+    if boltzmann == True:
+        FluorescencesVel = []
+        sigmaT = np.sqrt(kboltz*T/m)
+        
+        velvec = np.linspace(-4*sigmaT,4*sigmaT,100)
+        MBprobVec = norm.pdf(velvec,loc = 0, scale = sigmaT)
+        for i in range(len(velvec)):
+            detVel = (detpvec - kp*velvec[i])/(2*np.pi*1e6)
+            _, Fluorescence = CPTspectrum8levels( sg, sp, gPS, gPD, DetDoppler-kg*velvec[i]/(2*np.pi*1e6),u, DopplerLaserLinewidth, ProbeLaserLinewidth, 0,alpha, phidoppler, titadoppler, phiprobe, titaprobe,circularityprobe, betag, betap, drivefreq, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=False,detpvec=detVel)
+            
+            FluorescencesVel.append(Fluorescence)    
+        FluorescencesVel = np.array(FluorescencesVel)
+        MBprobMat = np.tile(MBprobVec,(FluorescencesVel.shape[1],1)).T
+        Fluorescence = np.trapz(FluorescencesVel*MBprobMat,velvec,axis = 0)
+        ProbeDetuningVectorL = np.arange(freqMin,freqMax,freqStep)
+
+    elif dephasing == True:  
+        tinicial = time.time()
+        ProbeDetuningVectorL, Fluorescence = CPTspectrum8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, circularityprobe, betag, betap, drivefreq, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, detpvec = detpvec, plot=False, solvemode=1)
+        tfinal = time.time()
+        
+
+    else:
+        drivefreq =2*np.pi* 1e3
+        FluorescencesVel = []
+        sigmaT = np.sqrt(T*kboltz/m)
+        velvec = np.linspace(-4*sigmaT,4*sigmaT,100)
+        MBprobVec = norm.pdf(velvec,loc = 0, scale = sigmaT)
+        betagVec = velvec*kg/drivefreq
+        betapVec = velvec*kp/drivefreq
+        for i in range(len(betagVec)): 
+            betag,betap = betagVec[i],betapVec[i]
+                                                         
+            Frequencies, Fluorescence = CPTspectrum8levels( sg, sp, gPS, gPD, DetDoppler,u, DopplerLaserLinewidth, ProbeLaserLinewidth, 0,alpha, phidoppler, titadoppler, phiprobe, titaprobe,circularityprobe, betag, betap, drivefreq, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=False,detpvec=None)
+            FluorescencesVel.append(Fluorescence)
+        FluorescencesVel = np.array(FluorescencesVel)
+        MBprobMat = np.tile(MBprobVec,(FluorescencesVel.shape[1],1)).T
+        Fluorescence = np.trapz(FluorescencesVel*MBprobMat,velvec,axis = 0)
+        ProbeDetuningVectorL = Frequencies
     
     #print('Done, Total time: ', round((tfinal-tinicial), 2), "s")   
           
-    return ProbeDetuningVectorL, Fluovector
+    return ProbeDetuningVectorL, Fluorescence
 
 def PerformExperiment_8levels_vel(velvect,titavect,phivect,velprob,sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, beta, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None):
     """
@@ -63,10 +106,10 @@ def SmoothNoisyCPT(Fluo, window=11, poly=3):
 
 
 
-def fitCPT_8levels_db(Freq,Fluo,sg,sp,gPS,gPD,DetDoppler,u,DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe,beta, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False):
+def fitCPT_8levels(Freq,Fluo,sg,sp,gPS,gPD,DetDoppler,u,DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe,beta, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False,dephasing = False, boltzmann = False ):
 
     def SpectrumForFit(Freq,sg,sp,T,DetDoppler):
-        freq,spectra = PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, beta, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None)
+        freq,spectra = PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, beta, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=Freq)
         return spectra
 
     

analisis/plots/20230920_CPT_TemperatureSens_v2/Data/EITfit/MM_eightLevel_2repumps_CPTPlotter.py

@@ -10,7 +10,7 @@ Created on Tue Sep  1 17:58:39 2020
 
 import os
 import numpy as np
-from MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels, GenerateNoisyCPT_vel, SmoothNoisyCPT, PerformExperiment_8levels_vel,fitCPT_8levels_db
+from MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels, GenerateNoisyCPT_vel, SmoothNoisyCPT, PerformExperiment_8levels_vel
 import matplotlib.pyplot as plt
 import time
 import h5py
@@ -107,61 +107,66 @@ Tmat = np.zeros((len(s12vec),len(s23vec)))
 Tvec_ajuste = np.zeros(len(Tvec))
 phivect = 0
 
+T = 10e-3
 
 titavect = np.linspace(0,2*np.pi,100)    
   
-beta = 0
+betag = 0
+betap = 0 
 alpha = 0
 
-for l in range(len(Tvec)):
-    for i in range(len(s12vec)):
-        for j in range(len(s23vec)):
-            #print(i,j) 
-            sg = s12vec[i]
-            sp = s23vec[j]
-            
-            T = Tvec[l]  
-
-            velvect = np.linspace(-4 * np.sqrt(kboltz*T/m),4 * np.sqrt(kboltz*T/m),100)
-            
-        
-            MB_prob_vec = np.array(MB_prob_1d(velvect,T))
-            
-
-            freq, Fluo = GenerateNoisyCPT_vel(velvect,titavect,phivect,MB_prob_vec,sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, 0, alpha, phidoppler, titadoppler, phiprobe, titaprobe, beta,drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001)
-             
-            plt.show()
-
-            plt.figure()
-
-            
-            a = plt.plot(freq,Fluo,'.',label = 'T = {:.1f} mK'.format(T*1000))
-            
-    
-
-            freq, Fluo_doppler,popt,pcov = fitCPT_8levels_db(freq,Fluo,sg,sp,gPS,gPD,DetDoppler,u,DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, beta,drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False)
-            
-            Tmat[i,j] = popt[2]
-     
-
-            plt.plot(freq,Fluo_doppler,'-')
-            plt.grid()
-            plt.legend()
-            plt.xlabel(' Detuning repump (MHz)')
-            plt.ylabel('Fluorescence')
-            plt.title('s12 = {:.2f}  s23 = {:.2f}  T = {:.1f} mK Det = {:.1f}\n s12 = {:.2f}  s23 = {:.2f}  T = {:.1f} mK Det = {:.1f}'.format(sg,sp,T*1000,DetDoppler,popt[0],popt[1],popt[2]*1000,popt[3]/1e6))
-
-            plt.savefig('ajuste_espectros_s12_{:.2f}_s23_{:.2f}_T_{:.0f}_mK_gamma_un_cuarto.png'.format(sg,sp,T*1000))
-            
-'''
-fig2,ax2 = plt.subplots()
-im = ax2.imshow(diferencias.T*100,extent = (s12vec[0],s12vec[-1],s23vec[0],s23vec[-1]),origin = 'lower',aspect = 'auto')
-#plt.xlabel('Detuning repump (MHz)')
-#plt.ylabel('Fluorescence (A.U.)')
-plt.colorbar(im,ax = ax2,label = 'Porcentaje diferencia' )
-ax2.set_title('T = {:} mK'.format(1000*T))
-ax2.set_xlabel('s12')
-ax2.set_ylabel('s23')
-fig2.savefig('mapa_comparacion_T_{:.2f}_gamma_un_cuarto.png'.format(T*1000))
-'''
-np.save('mapa_estimacion_T_{:.2f}_gamma_un_cuarto'.format(T*1000))
+Freq,Fluo_sb = PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, betag, betap, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None, dephasing  = False, boltzmann = False)
+Freq,Fluo_boltz = PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, betag, betap, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None, dephasing  = False, boltzmann = True)
+Freq,Fluo_deph = PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, betag, betap, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None, dephasing  = True, boltzmann = False)
+
+#%%
+plt.plot(Freq,Fluo_boltz,label = 'Velocities')
+plt.plot(Freq,Fluo_sb,label = 'Sideband')
+plt.plot(Freq,Fluo_deph, label = 'Dephasing')
+
+plt.legend()
+
+#%%
+
+os.chdir('/home/muri/nubeDF/Documents/codigos/artiq_experiments/analisis/plots/20230920_CPT_TemperatureSens_v2/Data/')
+
+CPT_FILES = """000011345-IR_Scan_withcal_optimized
+000011331-IR_Scan_withcal_optimized
+""" 
+
+for i in range(0,9):
+    CPT_FILES = CPT_FILES + f'0000153{19 + i}-IR_Scan_withcal_optimized\n'
+CPT_FILES = CPT_FILES
+
+
+
+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))
+
+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']))
+
+

analisis/plots/20230920_CPT_TemperatureSens_v2/Data/EITfit/MM_eightLevel_2repumps_python_scripts.py

@@ -105,12 +105,12 @@ def HImatrix(rabG, rabP, phidoppler, titadoppler, phiprobe, titaprobe, circulari
     return HI
 
 
-def LtempCalculus(beta, drivefreq, forma=1):
+def LtempCalculus(betag,betap, drivefreq, forma=1):
     Hint = np.zeros((8, 8), dtype=np.complex_)
     
-    ampg=beta*drivefreq
-    ampr=beta*drivefreq        
-    
+    ampg=betag*drivefreq
+    ampr=betap*drivefreq        
+
     Hint[0,0] = ampg
     Hint[1,1] = ampg
     Hint[4,4] = ampr
@@ -256,7 +256,7 @@ def dopplerBroadening(wlg, wlp, alpha, T, mcalcio = 6.655e-23*1e-3):
     
 
 def FullL(rabG, rabP, gPS = 0, gPD = 0, Detg = 0, Detp = 0, u = 0, lwg = 0, lwp = 0, 
-          phidoppler=0, titadoppler=0, phiprobe=0, titaprobe=0, beta=0, drivefreq=2*np.pi*22.135*1e6, T = 0, alpha = 0, circularityprobe=1):
+          phidoppler=0, titadoppler=0, phiprobe=0, titaprobe=0, betag=0,betap = 0, drivefreq=2*np.pi*22.135*1e6, T = 0, alpha = 0, circularityprobe=1):
 
     """
     Calcula el Liouvilliano total de manera explícita índice a índice. Suma aparte las componentes de las matrices M.
@@ -268,7 +268,7 @@ def FullL(rabG, rabP, gPS = 0, gPD = 0, Detg = 0, Detp = 0, u = 0, lwg = 0, lwp
     
     lwg = np.sqrt(lwg**2 + (db/2)**2)/2
     lwp = np.sqrt(lwp**2 + (db/2)**2)/2    
-    
+
 
     
 
@@ -309,7 +309,7 @@ def FullL(rabG, rabP, gPS = 0, gPD = 0, Detg = 0, Detp = 0, u = 0, lwg = 0, lwp
     #ESTA PARTE ES CUANDO AGREGAS MICROMOCION
     nmax = 7
     #print(nmax)
-    Ltemp, Omega = LtempCalculus(beta, drivefreq)
+    Ltemp, Omega = LtempCalculus(betag,betap, drivefreq)
     #print(factor)
     L1 = GetL1(Ltemp, L0, Omega, nmax)
     Lfull = L0 + L1 #ESA CORRECCION ESTA EN L1
@@ -334,8 +334,6 @@ def FullLvel(vel,tita,phi,rabG, rabP, gPS = 0, gPD = 0, Detg = 0, Detp = 0, u =
     Calcula el Liouvilliano total de manera explícita índice a índice. Suma aparte las componentes de las matrices M.
     Es la más eficiente hasta ahora.
     """
-    
-    
 
 
     db = dopplerBroadening(0.397e-6, 0.866e-6, alpha, T)
@@ -406,7 +404,7 @@ Scripts para correr un experimento y hacer el análisis de los datos
 """
 
 
-def CPTspectrum8levels(sg, sp, gPS, gPD, Detg, u, lwg, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, Circularityprobe, beta, drivefreq, freqMin=-100, freqMax=100, freqStep=1e-1, detpvec = None,plot=False, solvemode=1):
+def CPTspectrum8levels(sg, sp, gPS, gPD, Detg, u, lwg, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, Circularityprobe, betag,betap, drivefreq, freqMin=-100, freqMax=100, freqStep=1e-1, detpvec = None,plot=False, solvemode=1):
 
     """
     ESTA ES LA FUNCION QUE ESTAMOS USANDO
@@ -418,12 +416,13 @@ def CPTspectrum8levels(sg, sp, gPS, gPD, Detg, u, lwg, lwp, Temp, alpha, phidopp
     phidoppler, titadoppler = phidoppler*(np.pi/180), titadoppler*(np.pi/180) 
     phiprobe, titaprobe = phiprobe*(np.pi/180), titaprobe*(np.pi/180)
     
-    if detpvec == None:
+    if detpvec is None:
     
         DetProbeVector = 2*np.pi*np.arange(freqMin*1e6, freqMax*1e6, freqStep*1e6)
     else: 
         DetProbeVector = detpvec*1e6 * 2*np.pi
     
+   
     Detg = 2*np.pi*Detg*1e6
     #lwg, lwr, lwp = 2*np.pi*lwg*1e6, 2*np.pi*lwr*1e6, 2*np.pi*lwp*1e6
     lwg, lwp = lwg*1e6, lwp*1e6
@@ -440,7 +439,7 @@ def CPTspectrum8levels(sg, sp, gPS, gPD, Detg, u, lwg, lwp, Temp, alpha, phidopp
     for Detp in DetProbeVector:
         
         
-        L = FullL(rabG, rabP, gPS, gPD, Detg, Detp, u, lwg, lwp, phidoppler, titadoppler, phiprobe, titaprobe, beta, drivefreq, Temp, alpha, Circularityprobe)
+        L = FullL(rabG, rabP, gPS, gPD, Detg, Detp, u, lwg, lwp, phidoppler, titadoppler, phiprobe, titaprobe, betag,betap, drivefreq, Temp, alpha, Circularityprobe)
         
         if solvemode == 1:
             rhovectorized = np.linalg.solve(L, np.array([int(i==0) for i in range(64)]))