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

analisis/plots/20231123_CPTconmicromocion3/Data/EITfit/lolo_modelo_full_3niveles.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+4 de ene 2024
+
+@author: lolo
+
+reingenieria del código que anda
+
+Surge de fusionar
+Data/EITfit/threeLevel_2repumps_AnalysisFunctions.py
+Data/EITfit/threeLevel_2repumps_linealpol_python_scripts.py
+
+MAPA de FUNCIONES
+
+
+GenerateNoisyCPT_fit
+ |--> PerformExperiment_8levels_fixedRabi
+   |--> CPTspectrum8levels_fixedRabi --> ndarray,ndarray
+     |--> FullL_efficient --> ndarray
+       |--> dopplerBroadening --> float
+       |--> EffectiveL --> ndarray
+       |--> H0matrix --> ndarray
+       |--> HImatrix --> ndarray
+       |--> CalculateSingleMmatrix --> ndarray
+       |--> Lplusminus --> ndarray,ndarray,ndarray
+       |--> GetL1 --> ndarray
+
+"""
+
+# pylint: disable=C0301,R0913,R0914,W0621
+
+
+
+import os
+import numpy as np
+import time
+import matplotlib.pyplot as plt
+from scipy.signal import argrelextrema
+# from EITfit.threeLevel_2repumps_linealpol_python_scripts import CPTspectrum8levels, CPTspectrum8levels_fixedRabi
+import random
+from scipy.signal import savgol_filter as sf
+
+from  numba import jit,njit
+
+
+
+#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+@njit
+def make_diag(vec):
+    "Construye matris diagonal desde una lista o vector"
+    return np.eye(len(vec))*np.array(vec).reshape(-1,1)
+
+
+
+#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+"""
+Scripts para el calculo de la curva CPT
+"""
+
+@njit
+def H0matrix(Detg, Detp, u):
+    """
+    Calcula la matriz H0 en donde dr es el detuning del doppler, dp es el retuning
+    del repump y u es el campo magnético en Hz/Gauss.
+    Para esto se toma la energía del nivel P como 0
+    """
+
+    eigenEnergies = (Detg-u, Detg+u, -u/3, u/3, Detp-6*u/5,
+                    Detp-2*u/5, Detp+2*u/5, Detp+6*u/5)
+    #pagina 26 de Oberst. los lande del calcio son iguales a Bario.
+    # H0 = np.diag(eigenEnergies)
+    # H0 = np.eye(len(eigenEnergies))*np.array(eigenEnergies).reshape(-1,1)
+    H0 = make_diag(eigenEnergies)
+    return H0
+
+@njit
+def HImatrix(rabG, rabP, phidoppler, titadoppler, phiprobe, titaprobe):
+    """
+    Calcula la matriz de interacción Hsp + Hpd, en donde rabR es la frecuencia de rabi de la transición Doppler SP,
+    rabP es la frecuencia de rabi de la transición repump DP, y las componentes ei_r y ei_p son las componentes de la polarización
+    del campo eléctrico incidente de doppler y repump respectivamente. Deben estar normalizadas a 1
+    """
+    HI = np.zeros((8, 8), dtype=np.complex_)
+
+    i, j = 1, 3
+    HI[i-1, j-1] = (rabG/np.sqrt(3)) * np.cos(titadoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 1, 4
+    HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.sin(titadoppler)*np.exp(1j*phidoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 2, 3
+    HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.sin(titadoppler)*np.exp(-1j*phidoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 2, 4
+    HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.cos(titadoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+
+    i, j = 3, 5
+    HI[i-1, j-1] = -(rabP/2) * np.sin(titaprobe)*np.exp(-1j*phiprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 3, 6
+    HI[i-1, j-1] = -(rabP/np.sqrt(3)) * np.cos(titaprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 3, 7
+    HI[i-1, j-1] = rabP/np.sqrt(12) * np.sin(titaprobe)*np.exp(1j*phiprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+
+    i, j = 4, 6
+    HI[i-1, j-1] = -(rabP/np.sqrt(12)) * np.sin(titaprobe)*np.exp(-1j*phiprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 4, 7
+    HI[i-1, j-1] = -(rabP/np.sqrt(3)) * np.cos(titaprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 4, 8
+    HI[i-1, j-1] = (rabP/2) * np.sin(titaprobe)*np.exp(1j*phiprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    return HI
+
+
+@njit
+def Lplusminus(detr, detp, phirepump, titarepump, forma=1):
+    Hintplus = np.zeros((8, 8), dtype=np.complex_)
+    Hintminus = np.zeros((8, 8), dtype=np.complex_)
+
+
+    Hintplus[4, 2] = (-1/2)*np.sin(titarepump)*np.exp(1j*phirepump)
+    Hintplus[5, 2] = (-1/np.sqrt(3))*np.cos(titarepump)
+    Hintplus[6, 2] = (1/(2*np.sqrt(3)))*np.sin(titarepump)*np.exp(-1j*phirepump)
+
+    Hintplus[5, 3] = (-1/(2*np.sqrt(3)))*np.sin(titarepump)*np.exp(1j*phirepump)
+    Hintplus[6, 3] = (-1/np.sqrt(3))*np.cos(titarepump)
+    Hintplus[7, 3] = (1/2)*np.sin(titarepump)*np.exp(-1j*phirepump)
+
+
+    Hintminus[2, 4] = (-1/2)*np.sin(titarepump)*np.exp(-1j*phirepump)
+    Hintminus[2, 5] = (-1/np.sqrt(3))*np.cos(titarepump)
+    Hintminus[2, 6] = (1/(2*np.sqrt(3)))*np.sin(titarepump)*np.exp(1j*phirepump)
+
+    Hintminus[3, 5] = (-1/(2*np.sqrt(3)))*np.sin(titarepump)*np.exp(-1j*phirepump)
+    Hintminus[3, 6] = (-1/np.sqrt(3))*np.cos(titarepump)
+    Hintminus[3, 7] = (1/2)*np.sin(titarepump)*np.exp(1j*phirepump)
+
+    # if forma==1:
+    if True:
+        Lplus = np.zeros((64, 64), dtype=np.complex_)
+        Lminus = np.zeros((64, 64), dtype=np.complex_)
+        DeltaBar = np.zeros((64, 64), dtype=np.complex_)
+
+        for r  in range(8):
+            for q in range(8):
+                for k in range(8):
+                    for j in range(8):
+                        if j==q:
+                            if (k==2 or k==3) and r > 3:
+                                Lplus[r*8+q][k*8+j] = (-1j)*(Hintplus[r,k])
+                            if (r==2 or r==3) and k > 3:
+                                Lminus[r*8+q][k*8+j] = (-1j)*(Hintminus[r,k])
+                        elif r==k:
+                            if (q==2 or q==3) and j > 3:
+                                Lplus[r*8+q][k*8+j] = (-1j)*(- Hintplus[j,q])
+                            if (j==2 or j==3) and q > 3:
+                                Lminus[r*8+q][k*8+j] = (-1j)*(- Hintminus[j,q])
+    # if forma==2:
+    #     deltaKro = np.diag([1, 1, 1, 1, 1, 1, 1, 1])
+    #     Lplus = (-1j)*(np.kron(Hintplus, deltaKro) - np.kron(deltaKro, Hintplus))
+    #     Lminus = (-1j)*(np.kron(Hintminus, deltaKro) - np.kron(deltaKro, Hintminus))
+
+    DeltaBar = np.zeros((64, 64), dtype=np.complex_)
+    for i in range(64):
+
+        DeltaBar[i, i] = (1j)*(detr - detp)
+
+    # return np.matrix(Lminus), np.matrix(Lplus), np.matrix(DeltaBar)
+    return Lminus, Lplus, DeltaBar
+
+@njit
+def GetL1(Lplus, Lminus, DeltaBar, L0, rabR, nmax):
+    """
+    Devuelve Splus0 y Sminus0
+    """
+
+    # Sp = (-1)*(0.5*rabR)*(np.matrix(np.linalg.inv(L0 - (nmax+1)*DeltaBar))*np.matrix(Lplus))
+    # Sm = (-1)*(0.5*rabR)*(np.matrix(np.linalg.inv(L0 + (nmax+1)*DeltaBar))*np.matrix(Lminus))
+
+    Sp = (-1)*(0.5*rabR)*(np.linalg.inv(L0 - (nmax+1)*DeltaBar)).dot(Lplus)
+    Sm = (-1)*(0.5*rabR)*(np.linalg.inv(L0 + (nmax+1)*DeltaBar)).dot(Lminus)
+
+    for n in list(range(nmax+1))[(nmax+1)::-1][0:len(list(range(nmax+1))[(nmax+1)::-1])-1]: #jaja esto solo es para que vaya de nmax a 1 bajando. debe haber algo mas facil pero kcio
+        # Sp = (-1)*(rabR)*(np.matrix(np.linalg.inv(L0 - n*DeltaBar + rabR*(Lminus*np.matrix(Sp))))*np.matrix(Lplus))
+        # Sm = (-1)*(rabR)*(np.matrix(np.linalg.inv(L0 + n*DeltaBar + rabR*(Lplus*np.matrix(Sm))))*np.matrix(Lminus))
+
+        Sp = (-1)*(rabR)*((np.linalg.inv(L0 - n*DeltaBar + rabR*(Lminus.dot(Sp)))).dot(Lplus))
+        Sm = (-1)*(rabR)*((np.linalg.inv(L0 + n*DeltaBar + rabR*(Lplus.dot(Sm)))).dot(Lminus))
+
+    L1 = 0.5*rabR*(Lminus.dot(Sp) + Lplus.dot(Sm))
+
+    return L1
+
+@njit
+def EffectiveL(gPS, gPD, lwg, lwr, lwp):
+
+    """
+    Siendo Heff = H + EffectiveL, calcula dicho EffectiveL que es (-0.5j)*sumatoria(CmDaga*Cm) que luego sirve para calcular el Liouvilliano
+    """
+    Leff = np.zeros((8, 8), dtype=np.complex_)
+
+    Leff[0, 0] = 2*lwg
+    Leff[1, 1] = 2*lwg
+    Leff[2, 2] = ((2/3)+(1/3))*gPS + ((1/2) + (1/6) + (1/3))*gPD
+    Leff[3, 3] = ((2/3)+(1/3))*gPS + ((1/2) + (1/6) + (1/3))*gPD
+    Leff[4, 4] = 2*(lwr + lwp)
+    Leff[5, 5] = 2*(lwr + lwp)
+    Leff[6, 6] = 2*(lwr + lwp)
+    Leff[7, 7] = 2*(lwr + lwp)
+
+    return (-0.5j)*Leff
+
+@njit
+def CalculateSingleMmatrix(gPS, gPD, lwg, lwr, lwp):
+    """
+    Si tomamos el Liuvilliano como L = (-j)*(Heff*deltak - Heffdaga*deltak) + sum(Mm),
+    esta funcion calcula dichos Mm, que tienen dimensión 64x64 ya que esa es la dimensión del L. Estas componentes
+    salen de hacer la cuenta a mano conociendo los Cm y considerando que Mm[8*(r-1)+s, 8*(k-1)+j] = Cm[r,l] + Cmdaga[j,s] = Cm[r,l] + Cm[s,j]
+    ya que los componentes de Cm son reales.
+    Esta M es la suma de las 8 matrices M.
+
+    """
+
+    M =np.zeros((64, 64), dtype=np.complex_)
+
+    M[0,27] = (2/3)*gPS
+
+    M[9,18] = (2/3)*gPS
+
+    M[0,18] = (1/3)*gPS
+    M[1,19] = -(1/3)*gPS
+    M[8,26] = -(1/3)*gPS
+    M[9,27] = (1/3)*gPS
+
+    M[36,18] = (1/2)*gPD
+    M[37,19] = (1/np.sqrt(12))*gPD
+    M[44,26] = (1/np.sqrt(12))*gPD
+    M[45,27] = (1/6)*gPD
+
+    M[54,18] = (1/6)*gPD
+    M[55,19] = (1/np.sqrt(12))*gPD
+    M[62,26] = (1/np.sqrt(12))*gPD
+    M[63,27] = (1/2)*gPD
+
+    M[45,18] = (1/3)*gPD
+    M[46,19] = (1/3)*gPD
+    M[53,26] = (1/3)*gPD
+    M[54,27] = (1/3)*gPD
+
+    M[0,0] = 2*lwg
+    M[1,1] = 2*lwg
+    M[8,8] = 2*lwg
+    M[9,9] = 2*lwg
+
+    factor1 = 1
+    factor2 = 1
+    factor3 = 1
+    factor4 = 1
+
+
+    #M[36, 45] = lwp
+
+
+    M[36,36] = 2*(lwr + factor1*lwp)
+    M[37,37] = 2*(lwr + factor1*lwp)
+    M[38,38] = 2*(lwr + factor1*lwp)
+    M[39,39] = 2*(lwr + factor1*lwp)
+    M[44,44] = 2*(lwr + factor2*lwp)
+    M[45,45] = 2*(lwr + factor2*lwp)
+    M[46,46] = 2*(lwr + factor2*lwp)
+    M[47,47] = 2*(lwr + factor2*lwp)
+    M[52,52] = 2*(lwr + factor3*lwp)
+    M[53,53] = 2*(lwr + factor3*lwp)
+    M[54,54] = 2*(lwr + factor3*lwp)
+    M[55,55] = 2*(lwr + factor3*lwp)
+    M[60,60] = 2*(lwr + factor4*lwp)
+    M[61,61] = 2*(lwr + factor4*lwp)
+    M[62,62] = 2*(lwr + factor4*lwp)
+    M[63,63] = 2*(lwr + factor4*lwp)
+
+    return M
+
+
+@njit
+def dopplerBroadening(wlg, wlp, alpha, T, mcalcio = 6.655e-23*1e-3):
+    """
+    Calcula el broadening extra semiclásico por temperatura considerando que el ion atrapado se mueve.
+    wlg 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/(wlg*wlg)) + (1/(wlp*wlp)) - 2*(1/(wlg*wlp))*np.cos(alpha))*np.sqrt(kboltzmann*T/(2*mcalcio))
+
+    return gammaD
+
+@njit
+def FullL_efficient(rabG, rabR, rabP, gPS = 0, gPD = 0, Detg = 0, Detr = 0, Detp = 0, u = 0, lwg = 0, lwr=0, lwp = 0,
+          phidoppler=0, titadoppler=0, phiprobe=0, titaprobe=0, phirepump=0, titarepump=0, T = 0, alpha = 0):
+
+    """
+    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)
+
+    #lwr = np.sqrt(lwr**2 + dopplerBroadening(0.397e-6, 0.866e-6, alpha, T)**2)
+    lwg = np.sqrt(lwg**2 + db**2)
+    lwr = np.sqrt(lwr**2 + db**2)
+
+    CC = EffectiveL(gPS, gPD, lwg, lwr, lwp)
+
+    Heff = H0matrix(Detg, Detp, u) + HImatrix(rabG, rabP, phidoppler, titadoppler, phiprobe, titaprobe) + CC
+
+    # Heffdaga = np.matrix(Heff).getH()
+    Heffdaga = Heff.transpose().conj()
+
+    Lfullpartial = np.zeros((64, 64), dtype=np.complex_)
+
+    for r  in range(8):
+        for q in range(8):
+            for k in range(8):
+                for j in range(8):
+                    if j!=q and r!=k:
+                        pass
+                    elif j==q and r!=k:
+                        if (r < 2 and k > 3) or (k < 2 and r > 3) or (r > 3 and k > 3) or (r==0 and k==1) or (r==1 and k==0) or (r==2 and k==3) or (r==3 and k==2): #todo esto sale de analizar explicitamente la matriz y tratar de no calcular cosas de más que dan cero
+                            pass
+                        else:
+                            Lfullpartial[r*8+q][k*8+j] = (-1j)*(Heff[r,k])
+                    elif j!=q and r==k:
+                        if (j < 2 and q > 3) or (q < 2 and j > 3) or (j > 3 and q > 3) or (j==0 and q==1) or (j==1 and q==0) or (j==2 and q==3) or (j==3 and q==2):
+                            pass
+                        else:
+                            Lfullpartial[r*8+q][k*8+j] = (-1j)*(-Heffdaga[j,q])
+                    else:
+                        if Heff[r,k] == Heffdaga[j,q]:
+                            pass
+                        else:
+                            Lfullpartial[r*8+q][k*8+j] = (-1j)*(Heff[r,k]-Heffdaga[j,q])
+
+    M = CalculateSingleMmatrix(gPS, gPD, lwg, lwr, lwp)
+    L0 = Lfullpartial + M
+
+    nmax = 1
+    Lminus, Lplus, DeltaBar = Lplusminus(Detr, Detp, phirepump, titarepump)
+    factor1 = np.exp(1j*0.2*np.pi)
+    factor2 = np.exp(-1j*0.2*np.pi)
+    #print(factor)
+    L1 = GetL1(factor1*Lplus, factor2*Lminus, DeltaBar, L0, rabR, nmax)
+    Lfull = L0 + L1
+
+
+    #NORMALIZACION DE RHO
+    i = 0
+    while i < 64:
+        if i%9 == 0:
+            Lfull[0, i] = 1
+        else:
+            Lfull[0, i] = 0
+        i = i + 1
+
+    return Lfull
+
+
+"""
+Scripts para correr un experimento y hacer el análisis de los datos
+"""
+
+def CalculoTeoricoDarkResonances(u, titadoppler):
+    if titadoppler==0:
+        NegativeDR = [(-7/5)*u, (-3/5)*u, (-1/5)*u, (1/5)*u, (3/5)*u, (7/5)*u]
+    elif titadoppler==90:
+        NegativeDR = [(-11/5)*u, (-7/5)*u, (-3/5)*u, (3/5)*u, (7/5)*u, (11/5)*u]
+
+    PositiveDR = [(-8/5)*u, (-4/5)*u, 0, (4/5)*u, (8/5)*u]
+
+    return NegativeDR, PositiveDR
+
+
+def CPTspectrum8levels(rabG, rabR, rabP, gPS, gPD, Detg, Detr, u, lwg, lwr, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump,
+                      freqMin=-100, freqMax=100, freqStep=1e-1, plot=False, solvemode=1):
+    """
+    Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
+
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+
+    phidoppler, titadoppler = phidoppler*(np.pi/180), titadoppler*(np.pi/180)
+    phiprobe, titaprobe = phiprobe*(np.pi/180), titaprobe*(np.pi/180)
+    phirepump, titarepump = phirepump*(np.pi/180), titarepump*(np.pi/180)
+
+    DetProbeVector = 2*np.pi*np.arange(freqMin*1e6, freqMax*1e6, freqStep*1e6)
+    Detg, Detr = 2*np.pi*Detg*1e6, 2*np.pi*Detr*1e6
+    lwg, lwr, lwp = 2*np.pi*lwg*1e6, 2*np.pi*lwr*1e6, 2*np.pi*lwp*1e6
+
+    #u = 2*np.pi*u*1e6
+
+    Fluovector = []
+
+    tinicial = time.time()
+
+    for Detp in DetProbeVector:
+        L = FullL_efficient(rabG, rabR, rabP, gPS, gPD, Detg, Detr, Detp, u, lwg, lwr, lwp, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, Temp, alpha)
+
+        # if solvemode == 1:
+        if True:
+            rhovectorized = np.linalg.solve(L, np.array([int(i==0) for i in range(64)],dtype=np.complex_))
+            Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
+            Fluovector.append(Fluo)
+
+        # if solvemode == 2:
+        #     Linv = np.linalg.inv(L)
+        #     rhovectorized = [Linv[j][0] for j in range(len(Linv))]
+        #     Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
+        #     Fluovector.append(Fluo)
+
+    tfinal = time.time()
+
+    print('Done, Total time: ', round((tfinal-tinicial), 2), "s")
+
+    DetProbeVectorMHz = np.arange(freqMin, freqMax, freqStep)
+
+    if plot:
+        plt.xlabel('Probe detuning (MHz)')
+        plt.ylabel('Fluorescence (A.U.)')
+        plt.plot(DetProbeVectorMHz, [100*f for f in Fluovector], label=str(titaprobe) + 'º, T: ' + str(Temp*1e3) + ' mK')
+        plt.legend()
+
+    return DetProbeVectorMHz, Fluovector
+
+@njit
+def CPTspectrum8levels_fixedRabi(sg, sr, sp, gPS, gPD, Detg, Detr, u, lwg, lwr, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump,
+                      freqMin=-100, freqMax=100, freqStep=1e-1, plot=False, solvemode=1):
+    """
+    Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
+
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+
+
+
+
+    phidoppler, titadoppler = phidoppler*(np.pi/180), titadoppler*(np.pi/180)
+    phiprobe, titaprobe = phiprobe*(np.pi/180), titaprobe*(np.pi/180)
+    phirepump, titarepump = phirepump*(np.pi/180), titarepump*(np.pi/180)
+
+    DetProbeVector = 2*np.pi*np.arange(freqMin*1e6, freqMax*1e6, freqStep*1e6)
+    Detg, Detr = 2*np.pi*Detg*1e6, 2*np.pi*Detr*1e6
+    #lwg, lwr, lwp = 2*np.pi*lwg*1e6, 2*np.pi*lwr*1e6, 2*np.pi*lwp*1e6
+    lwg, lwr, lwp = lwg*1e6, lwr*1e6, lwp*1e6
+
+    rabG = sg*gPS
+    rabR = sr*gPD
+    rabP = sp*gPD
+    #u = 2*np.pi*u*1e6
+
+    Fluovector = []
+
+    # tinicial = time.time()
+
+    for Detp in DetProbeVector:
+
+        L = FullL_efficient(rabG, rabR, rabP, gPS, gPD, Detg, Detr, Detp, u, lwg, lwr, lwp, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, Temp, alpha)
+
+        # if solvemode == 1:
+        if True:
+            coh = 5
+            rhovectorized = np.linalg.solve(L, np.array([int(i==0) for i in range(64)],dtype=np.complex_))
+            #Fluo = np.abs(rhovectorized[coh])
+            Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
+            Fluovector.append(Fluo)
+
+        # if solvemode == 2:
+        #     Linv = np.linalg.inv(L)
+        #     rhovectorized = [Linv[j][0] for j in range(len(Linv))]
+        #     Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
+        #     Fluovector.append(Fluo)
+
+    # tfinal = time.time()
+
+    # print('Done, Total time: ', round((tfinal-tinicial), 2), "s")
+
+    DetProbeVectorMHz = np.arange(freqMin, freqMax, freqStep)
+
+    # if plot:
+    #     plt.xlabel('Probe detuning (MHz)')
+    #     plt.ylabel('Fluorescence (A.U.)')
+    #     plt.plot(DetProbeVectorMHz, [100*f for f in Fluovector], label=str(titaprobe) + 'º, T: ' + str(Temp*1e3) + ' mK')
+    #     plt.legend()
+
+    return DetProbeVectorMHz, Fluovector
+#%%
+
+if __name__ == "__main__":
+    ub = 9.27e-24
+    h = 6.63e-34
+    c = (ub/h)*1e-4 #en unidades de MHz/G
+
+    B = 25 #campo magnetico en gauss
+
+    u = c*B
+
+
+    sg, sr, sp = 0.5, 1.5, 4 #parámetros de saturación del doppler y repump
+    gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 #anchos de linea de las transiciones
+    rabG, rabR, rabP = sg*gPS, sr*gPD, sp*gPD #frecuencias de rabi
+
+    lwg, lwr, lwp = 0.3, 0.3, 0.3 #ancho de linea de los laseres
+
+    Detg = -25
+    Detr = 20 #detuning del doppler y repump
+
+    Temp = 0.0e-3 #temperatura en K
+    alpha = 0*(np.pi/180) #angulo entre los láseres
+
+    phidoppler, titadoppler = 0, 90
+    phirepump, titarepump = 0, 90
+
+    phiprobe, titaprobe = 0, 90
+
+
+    plotCPT = False
+
+    freqMin = -50
+    freqMax = 50
+    freqStep = 5e-2
+
+    Frequencyvector, Fluovector = CPTspectrum8levels(rabG, rabR, rabP, gPS, gPD, Detg, Detr, u, lwg, lwr, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=plotCPT, solvemode=1)
+
+    NegativeDR, PositiveDR = CalculoTeoricoDarkResonances(u/(2*np.pi*1e6), titadoppler)
+
+    plt.plot(Frequencyvector, [100*f for f in Fluovector], label=str(titaprobe) + 'º, T: ' + str(Temp*1e3) + ' mK')
+    plt.xlabel('Probe detuning (MHz)')
+    plt.ylabel('Fluorescence (A.U.)')
+
+    for PDR in PositiveDR:
+        plt.axvline(Detr+PDR, linestyle='--', linewidth=0.5, color='red')
+    for NDR in NegativeDR:
+        plt.axvline(Detg+NDR, linestyle='--', linewidth=0.5, color='blue')
+
+
+
+
+    #parametros que andan piola:
+
+
+    """
+    ub = 9.27e-24
+    h = 6.63e-34
+    c = (ub/h)*1e-4 #en unidades de MHz/G
+
+    B = 17 #campo magnetico en gauss
+
+    u = c*B
+    #u = 80e6
+
+    sr, sp = 0.53, 4.2
+    gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6
+    rabR, rabP = sr*gPS, sp*gPD
+    lw = 2*np.pi * 0.33e6
+    lwr, lwp = lw, lw #ancho de linea de los laseres
+    dr_spec = - 2*np.pi* 26e6
+
+
+    freqSteps = 500
+    freqMin = -100e6
+    freqMax = 100e6
+
+    dps = 2*np.pi*np.linspace(freqMin, freqMax, freqSteps)
+
+    #dps = [-30e6]
+
+    alfar = 90*(np.pi/180)
+    ex_r, ey_r, ez_r = np.sin(alfar)*np.cos(0), np.sin(alfar)*np.sin(0), np.cos(alfar)
+
+
+    alfap = 90*(np.pi/180)
+    ex_p, ey_p, ez_p = np.sin(alfap)*np.cos(0), np.sin(alfap)*np.sin(0), np.cos(alfap)
+
+    """
+
+
+
+
+
+
+
+
+
+#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+
+
+
+
+def CalculoTeoricoDarkResonances_8levels(u, titadoppler, detuningdoppler, detuningrepump):
+
+    if titadoppler==0:
+        NegativeDR = [(-7/5)*u, (-3/5)*u, (-1/5)*u, (1/5)*u, (3/5)*u, (7/5)*u]
+    elif titadoppler==90:
+        NegativeDR = [(-11/5)*u, (-7/5)*u, (-3/5)*u, (3/5)*u, (7/5)*u, (11/5)*u]
+    else:
+        NegativeDR = [(-11/5)*u, (-7/5)*u, (-3/5)*u, (-1/5)*u, (1/5)*u, (3/5)*u, (7/5)*u, (11/5)*u]
+
+    PositiveDR = [(-8/5)*u, (-4/5)*u, 0, (4/5)*u, (8/5)*u]
+
+    return [detuningdoppler + dr for dr in NegativeDR], [detuningrepump + dr for dr in PositiveDR]
+
+
+def GetClosestIndex(Vector, value, tolerance=1e-3):
+
+    i = 0
+    while i<len(Vector):
+        if abs(Vector[i] - value) < tolerance:
+            return i
+        else:
+            i = i + 1
+    return GetClosestIndex(Vector, value, tolerance=2*tolerance)
+
+
+def FindDRFrequencies(Freq, Fluo, TeoDR, entorno=3):
+    """
+    Busca los indices y la frecuencia de los minimos en un entorno cercano al de la DR.
+    Si no encuentra, devuelve el valor teórico.
+    """
+
+    IndiceDRteo1, IndiceEntornoinicialDRteo1, IndiceEntornofinalDRteo1 = GetClosestIndex(Freq, TeoDR[0]), GetClosestIndex(Freq, TeoDR[0]-entorno), GetClosestIndex(Freq, TeoDR[0]+entorno)
+    IndiceDRteo2, IndiceEntornoinicialDRteo2, IndiceEntornofinalDRteo2 = GetClosestIndex(Freq, TeoDR[1]), GetClosestIndex(Freq, TeoDR[1]-entorno), GetClosestIndex(Freq, TeoDR[1]+entorno)
+    IndiceDRteo3, IndiceEntornoinicialDRteo3, IndiceEntornofinalDRteo3 = GetClosestIndex(Freq, TeoDR[2]), GetClosestIndex(Freq, TeoDR[2]-entorno), GetClosestIndex(Freq, TeoDR[2]+entorno)
+    IndiceDRteo4, IndiceEntornoinicialDRteo4, IndiceEntornofinalDRteo4 = GetClosestIndex(Freq, TeoDR[3]), GetClosestIndex(Freq, TeoDR[3]-entorno), GetClosestIndex(Freq, TeoDR[3]+entorno)
+    IndiceDRteo5, IndiceEntornoinicialDRteo5, IndiceEntornofinalDRteo5 = GetClosestIndex(Freq, TeoDR[4]), GetClosestIndex(Freq, TeoDR[4]-entorno), GetClosestIndex(Freq, TeoDR[4]+entorno)
+    IndiceDRteo6, IndiceEntornoinicialDRteo6, IndiceEntornofinalDRteo6 = GetClosestIndex(Freq, TeoDR[5]), GetClosestIndex(Freq, TeoDR[5]-entorno), GetClosestIndex(Freq, TeoDR[5]+entorno)
+
+    EntornoFreqDR1, EntornoFreqDR2 = Freq[IndiceEntornoinicialDRteo1:IndiceEntornofinalDRteo1], Freq[IndiceEntornoinicialDRteo2:IndiceEntornofinalDRteo2]
+    EntornoFreqDR3, EntornoFreqDR4 = Freq[IndiceEntornoinicialDRteo3:IndiceEntornofinalDRteo3], Freq[IndiceEntornoinicialDRteo4:IndiceEntornofinalDRteo4]
+    EntornoFreqDR5, EntornoFreqDR6 = Freq[IndiceEntornoinicialDRteo5:IndiceEntornofinalDRteo5], Freq[IndiceEntornoinicialDRteo6:IndiceEntornofinalDRteo6]
+
+    EntornoFluoDR1, EntornoFluoDR2 = Fluo[IndiceEntornoinicialDRteo1:IndiceEntornofinalDRteo1], Fluo[IndiceEntornoinicialDRteo2:IndiceEntornofinalDRteo2]
+    EntornoFluoDR3, EntornoFluoDR4 = Fluo[IndiceEntornoinicialDRteo3:IndiceEntornofinalDRteo3], Fluo[IndiceEntornoinicialDRteo4:IndiceEntornofinalDRteo4]
+    EntornoFluoDR5, EntornoFluoDR6 = Fluo[IndiceEntornoinicialDRteo5:IndiceEntornofinalDRteo5], Fluo[IndiceEntornoinicialDRteo6:IndiceEntornofinalDRteo6]
+
+    IndiceFluoMinimaEntorno1, IndiceFluoMinimaEntorno2 = argrelextrema(np.array(EntornoFluoDR1), np.less)[0], argrelextrema(np.array(EntornoFluoDR2), np.less)[0]
+    IndiceFluoMinimaEntorno3, IndiceFluoMinimaEntorno4 = argrelextrema(np.array(EntornoFluoDR3), np.less)[0], argrelextrema(np.array(EntornoFluoDR4), np.less)[0]
+    IndiceFluoMinimaEntorno5, IndiceFluoMinimaEntorno6 = argrelextrema(np.array(EntornoFluoDR5), np.less)[0], argrelextrema(np.array(EntornoFluoDR6), np.less)[0]
+
+    try:
+        FreqDR1 = EntornoFreqDR1[int(IndiceFluoMinimaEntorno1)]
+        IndiceDR1 = GetClosestIndex(Freq, FreqDR1)
+    except:
+        FreqDR1 = TeoDR[0]
+        IndiceDR1 = IndiceDRteo1
+    try:
+        FreqDR2 = EntornoFreqDR2[int(IndiceFluoMinimaEntorno2)]
+        IndiceDR2 = GetClosestIndex(Freq, FreqDR2)
+    except:
+        FreqDR2 = TeoDR[1]
+        IndiceDR2 = IndiceDRteo2
+    try:
+        FreqDR3 = EntornoFreqDR3[int(IndiceFluoMinimaEntorno3)]
+        IndiceDR3 = GetClosestIndex(Freq, FreqDR3)
+    except:
+        FreqDR3 = TeoDR[2]
+        IndiceDR3 = IndiceDRteo3
+    try:
+        FreqDR4 = EntornoFreqDR4[int(IndiceFluoMinimaEntorno4)]
+        IndiceDR4 = GetClosestIndex(Freq, FreqDR4)
+    except:
+        FreqDR4 = TeoDR[3]
+        IndiceDR4 = IndiceDRteo4
+    try:
+        FreqDR5 = EntornoFreqDR5[int(IndiceFluoMinimaEntorno5)]
+        IndiceDR5 = GetClosestIndex(Freq, FreqDR5)
+    except:
+        FreqDR5 = TeoDR[4]
+        IndiceDR5 = IndiceDRteo5
+    try:
+        FreqDR6 = EntornoFreqDR6[int(IndiceFluoMinimaEntorno6)]
+        IndiceDR6 = GetClosestIndex(Freq, FreqDR6)
+    except:
+        FreqDR6 = TeoDR[5]
+        IndiceDR6 = IndiceDRteo6
+
+    return [IndiceDR1, IndiceDR2, IndiceDR3, IndiceDR4, IndiceDR5, IndiceDR6], [FreqDR1, FreqDR2, FreqDR3, FreqDR4, FreqDR5, FreqDR6]
+
+
+def FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=1, frecuenciareferenciacriterioasintotico=-100, getindices=False):
+    """
+    Toma los indices donde estan las DR y evalua su fluorescencia. Esos indices son minimos locales en un entorno
+    cercano a las DR teoricas y, si no hay ningun minimo, toma la teorica.
+    Luego, hace el cociente de esa fluorescencia y un factor de normalización segun NormalizationCriterium:
+        1: Devuelve la fluorescencia absoluta de los minimos
+        2: Devuelve el cociente entre la fluorescencia del minimo y un valor medio entre dos puntos lejanos, como si no
+        hubiera una resonancia oscura y hubiera una recta. Ese valor esta a  DistanciaFrecuenciaCociente del detuning del azul (el punto medio entre las dos DR en este caso)
+        3: Devuelve el cociente entre la fluorescencia del minimo y el valor a -100 MHz (si se hizo de -100 a 100),
+        o el valor limite por izquierda de la curva
+        4: Deuelve el cociente entre la fluorescencia del minimo y el valor de fluorescencia a detuning 0 MHz
+
+    """
+    IndiceDR1, IndiceDR2, IndiceDR3, IndiceDR4, IndiceDR5, IndiceDR6 = IndicesDR[0], IndicesDR[1], IndicesDR[2], IndicesDR[3], IndicesDR[4], IndicesDR[5]
+    FluorescenceOfMinimums = [Fluo[IndiceDR1], Fluo[IndiceDR2], Fluo[IndiceDR3], Fluo[IndiceDR4], Fluo[IndiceDR5], Fluo[IndiceDR6]]
+    FrequencyOfMinimums = [Freq[IndiceDR1], Freq[IndiceDR2], Freq[IndiceDR3], Freq[IndiceDR4], Freq[IndiceDR5], Freq[IndiceDR6]]
+    DistanciaFrecuenciaCociente = 25
+
+    if NormalizationCriterium==0:
+        print('che')
+        return FrequencyOfMinimums, FluorescenceOfMinimums
+
+    if NormalizationCriterium==1:
+        Fluorescenciacerodetuning = Fluo[GetClosestIndex(Freq, 0)]
+        Fluorescenciaasintotica = Fluo[GetClosestIndex(Freq, frecuenciareferenciacriterioasintotico)]
+        return FrequencyOfMinimums, np.array([Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica])
+
+    if NormalizationCriterium==2:
+        k = 0
+        while k < len(Freq):
+            if Freq[k] < detuningdoppler-DistanciaFrecuenciaCociente + 2 and Freq[k] > detuningdoppler-DistanciaFrecuenciaCociente - 2:
+                FluoIzquierda = Fluo[k]
+                indiceizquierda = k
+                print('Izq:', Freq[k])
+                break
+            else:
+                k = k + 1
+
+        l = 0
+        while l < len(Freq):
+            if Freq[l] < detuningdoppler+DistanciaFrecuenciaCociente + 2 and Freq[l] > detuningdoppler+DistanciaFrecuenciaCociente - 2:
+                FluoDerecha = Fluo[l]
+                indicederecha = l
+                print('Der: ', Freq[l])
+                break
+            else:
+                l = l + 1
+        FluoNormDivisor = 0.5*(FluoDerecha+FluoIzquierda)
+        print(FluoNormDivisor)
+
+    if NormalizationCriterium==3:
+        #asintotico
+        FluoNormDivisor = Fluo[GetClosestIndex(Freq, frecuenciareferenciacriterioasintotico)]
+
+    if NormalizationCriterium==4:
+        #este te tira la fluorescencia de detuning 0
+        FluoNormDivisor = Fluo[GetClosestIndex(Freq, 0)]
+
+
+    RelativeFluorescenceOfMinimums = np.array([Fluore/FluoNormDivisor for Fluore in FluorescenceOfMinimums])
+    print('Esto: ', RelativeFluorescenceOfMinimums)
+
+    if NormalizationCriterium==2 and getindices==True:
+        return FrequencyOfMinimums, RelativeFluorescenceOfMinimums, indiceizquierda, indicederecha
+
+    return FrequencyOfMinimums, RelativeFluorescenceOfMinimums
+
+
+def GetFinalMaps(MapasDR1, MapasDR2, MapasDR3, MapasDR4, MapasDR5, MapasDR6):
+    """
+    Nota: esto vale para polarizacion del 397 sigma+ + sigma-. Sino hay que cambiar los coeficientes.
+    La estructura es:
+        MapasDRi = [MapaMedido_criterio1_DRi, MapaMedido_criterio2_DRi, MapaMedido_criterio3_DRi, MapaMedido_criterio4_DRi]
+    """
+    Mapa1 = MapasDR1[0]
+
+    Mapa2pi = np.sqrt(3)*(MapasDR2[1] + MapasDR5[1])
+    Mapa2smas = np.sqrt(12/2)*MapasDR3[1] + (2/np.sqrt(2))*MapasDR6[1]
+    Mapa2smenos = (2/np.sqrt(2))*MapasDR1[1] + np.sqrt(12/2)*MapasDR4[1]
+
+    Mapa3pi = np.sqrt(3)*(MapasDR2[2] + MapasDR5[2])
+    Mapa3smas = np.sqrt(12/2)*MapasDR3[2] + (2/np.sqrt(2))*MapasDR6[2]
+    Mapa3smenos = (2/np.sqrt(2))*MapasDR1[2] + np.sqrt(12/2)*MapasDR4[2]
+
+    return Mapa1, [Mapa2pi, Mapa2smas, Mapa2smenos], [Mapa3pi, Mapa3smas, Mapa3smenos]
+
+
+def CombinateDRwithCG(RelMinMedido1, RelMinMedido2, RelMinMedido3, RelMinMedido4):
+    Fluo1 = RelMinMedido1[0]
+
+    Fluo2pi = np.sqrt(3)*(RelMinMedido2[1] + RelMinMedido2[4])
+    Fluo2smas = np.sqrt(12/2)*RelMinMedido2[2] + (2/np.sqrt(2))*RelMinMedido2[5]
+    Fluo2smenos = (2/np.sqrt(2))*RelMinMedido2[0] + np.sqrt(12/2)*RelMinMedido2[3]
+
+    Fluo3pi = np.sqrt(3)*(RelMinMedido3[1] + RelMinMedido3[4])
+    Fluo3smas = np.sqrt(12/2)*RelMinMedido3[2] + (2/np.sqrt(2))*RelMinMedido3[5]
+    Fluo3smenos = (2/np.sqrt(2))*RelMinMedido3[0] + np.sqrt(12/2)*RelMinMedido3[3]
+
+    return Fluo1, [Fluo2pi, Fluo2smas, Fluo2smenos], [Fluo3pi, Fluo3smas, Fluo3smenos]
+
+
+
+def IdentifyPolarizationCoincidences(theoricalmap, target, tolerance=1e-1):
+    """
+    Busca en un mapa 2D la presencia de un valor target (medido) con tolerancia tolerance.
+    Si lo encuentra, pone un 1. Sino, un 0. Al plotear con pcolor se verá
+    en blanco la zona donde el valor medido se puede hallar.
+    """
+    CoincidenceMatrix = np.zeros((len(theoricalmap), len(theoricalmap[0])))
+
+    i = 0
+    while i<len(theoricalmap):
+        j = 0
+        while j<len(theoricalmap[0]):
+            if abs(theoricalmap[i][j]-target) < tolerance:
+                CoincidenceMatrix[i][j] = 1
+            j=j+1
+        i=i+1
+    return CoincidenceMatrix
+
+def RetrieveAbsoluteCoincidencesBetweenMaps(MapsVectors):
+
+    MatrixSum = np.zeros((len(MapsVectors[0]), len(MapsVectors[0][0])))
+    AbsoluteCoincidencesMatrix = np.zeros((len(MapsVectors[0]), len(MapsVectors[0][0])))
+
+
+    MatrixMapsVectors = []
+
+    for i in range(len(MapsVectors)):
+        MatrixMapsVectors.append(np.matrix(MapsVectors[i]))
+
+    for i in range(len(MatrixMapsVectors)):
+        MatrixSum = MatrixSum + MatrixMapsVectors[i]
+
+    MaxNumberOfCoincidences = np.max(MatrixSum)
+
+    ListMatrixSum = [list(i) for i in list(np.array(MatrixSum))]
+
+    for i in range(len(ListMatrixSum)):
+        for j in range(len(ListMatrixSum[0])):
+            if ListMatrixSum[i][j] == MaxNumberOfCoincidences:
+                AbsoluteCoincidencesMatrix[i][j] = 1
+
+    return AbsoluteCoincidencesMatrix, MaxNumberOfCoincidences
+
+
+def MeasureMeanValueOfEstimatedArea(AbsoluteCoincidencesMap, X, Y):
+    NonZeroIndices = np.nonzero(AbsoluteCoincidencesMap)
+
+    Xsum = 0
+    Xvec = []
+    Ysum = 0
+    Yvec = []
+    N = len(NonZeroIndices[0])
+    for i in range(N):
+        Xsum = Xsum + X[NonZeroIndices[1][i]]
+        Xvec.append(X[NonZeroIndices[1][i]])
+        Ysum = Ysum + Y[NonZeroIndices[0][i]]
+        Yvec.append(Y[NonZeroIndices[0][i]])
+    Xaverage = Xsum/N
+    Yaverage = Ysum/N
+    Xspread = np.std(Xvec)
+    Yspread = np.std(Yvec)
+    return Xaverage, Yaverage, N, Xspread, Yspread
+
+
+def MeasureRelativeFluorescenceFromCPT(Freq, Fluo, u, titadoppler, detuningrepump, detuningdoppler, frefasint=-100, entorno=3):
+    ResonanciasTeoricas, ResonanciasPositivas = CalculoTeoricoDarkResonances_8levels(u, titadoppler, detuningdoppler, detuningrepump)
+    IndicesDR, FreqsDR = FindDRFrequencies(Freq, Fluo, ResonanciasTeoricas, entorno=entorno)
+    FrequencyOfMinimums, RelativeFluorescenceOfMinimums0 = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=0, frecuenciareferenciacriterioasintotico=frefasint)
+    FrequencyOfMinimums, RelativeFluorescenceOfMinimums1 = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=1, frecuenciareferenciacriterioasintotico=frefasint)
+    FrequencyOfMinimums, RelativeFluorescenceOfMinimums2, indiceizquierda, indicederecha = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=2, frecuenciareferenciacriterioasintotico=frefasint, getindices=True)
+    FrequencyOfMinimums, RelativeFluorescenceOfMinimums3 = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=3, frecuenciareferenciacriterioasintotico=frefasint)
+    FrequencyOfMinimums, RelativeFluorescenceOfMinimums4 = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=4, frecuenciareferenciacriterioasintotico=frefasint)
+    print('hola')
+    print(RelativeFluorescenceOfMinimums0)
+    return RelativeFluorescenceOfMinimums0, RelativeFluorescenceOfMinimums1, RelativeFluorescenceOfMinimums2, RelativeFluorescenceOfMinimums3, RelativeFluorescenceOfMinimums4, IndicesDR, [indiceizquierda, indicederecha]
+
+
+def GenerateNoisyCPT(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
+    Frequencyvector, Fluovector = PerformExperiment_8levels(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None)
+    NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
+    return Frequencyvector, NoisyFluovector
+
+
+def GenerateNoisyCPT_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
+    Frequencyvector, Fluovector = PerformExperiment_8levels_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None)
+    NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
+    return Frequencyvector, NoisyFluovector
+
+def GenerateNoisyCPT_fit(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqs, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
+    Frequencyvector, Fluovector = PerformExperiment_8levels_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, min(freqs), max(freqs) + freqs[1]-freqs[0], freqs[1]-freqs[0], plot=False, solvemode=1, detpvec=None)
+    NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
+    return Frequencyvector, NoisyFluovector
+
+
+def AddNoiseToCPT(Fluo, noisefactor):
+    return [f+noisefactor*(2*random.random()-1) for f in Fluo]
+
+
+def SmoothNoisyCPT(Fluo, window=11, poly=3):
+    SmoothenFluo = sf(Fluo, window, poly)
+    return SmoothenFluo
+
+
+def GetMinimaInfo(Freq, Fluo, u, titadoppler, detuningdoppler, detuningrepump, MinimumCriterium=2, NormalizationCriterium=1):
+    """
+    FUNCION VIEJA
+
+    Esta funcion devuelve valores de frecuencias y fluorescencia relativa de los minimos.
+    Minimumcriterion:
+        1: Saca los minimos con funcion argelextrema
+        2: Directamente con las frecuencias teoricas busca las fluorescencias
+
+    Normalizationcriterium:
+        1: Devuelve la fluorescencia absoluta de los minimos
+        2: Devuelve el cociente entre la fluorescencia del minimo y un valor medio entre dos puntos lejanos, como si no
+        hubiera una resonancia oscura y hubiera una recta. Ese valor esta a  DistanciaFrecuenciaCociente del detuning del azul (el punto medio entre las dos DR en este caso)
+        3: Devuelve el cociente entre la fluorescencia del minimo y el valor a -100 MHz (si se hizo de -100 a 100),
+        o el valor limite por izquierda de la curva
+
+    """
+    FluorescenceOfMaximum = max(Fluo)
+    FrequencyOfMaximum = Freq[Fluo.index(FluorescenceOfMaximum)]
+
+
+    #criterio para encontrar los minimos
+    #criterio usando minimos de la fluorescencia calculados con la curva
+    if MinimumCriterium == 1:
+        LocationOfMinimums = argrelextrema(np.array(Fluo), np.less)[0]
+        FluorescenceOfMinimums = np.array([Fluo[i] for i in LocationOfMinimums])
+        FrequencyOfMinimums = np.array([Freq[j] for j in LocationOfMinimums])
+
+
+    #criterio con las DR teoricas
+    if MinimumCriterium == 2:
+        FrecuenciasDRTeoricas, FrecuenciasDRTeoricasPositivas = [darkresonance for darkresonance in CalculoTeoricoDarkResonances_8levels(u, titadoppler, detuningdoppler, detuningrepump)[0]]
+        FrequencyOfMinimums = []
+        FluorescenceOfMinimums =[]
+        print(FrecuenciasDRTeoricas)
+        k=0
+
+        ventanita = 0.001
+        while k < len(Freq):
+            if Freq[k] < FrecuenciasDRTeoricas[0] + ventanita and Freq[k] > FrecuenciasDRTeoricas[0] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            elif Freq[k] < FrecuenciasDRTeoricas[1] + ventanita and Freq[k] > FrecuenciasDRTeoricas[1] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            elif Freq[k] < FrecuenciasDRTeoricas[2] + ventanita and Freq[k] > FrecuenciasDRTeoricas[2] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            elif Freq[k] < FrecuenciasDRTeoricas[3] + ventanita and Freq[k] > FrecuenciasDRTeoricas[3] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            elif Freq[k] < FrecuenciasDRTeoricas[4] + ventanita and Freq[k] > FrecuenciasDRTeoricas[4] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            elif Freq[k] < FrecuenciasDRTeoricas[5] + ventanita and Freq[k] > FrecuenciasDRTeoricas[5] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            k = k + 1
+        print(FrequencyOfMinimums)
+        if len(FrequencyOfMinimums) !=  len(FrecuenciasDRTeoricas):
+            print('NO ANDA BIEN ESTO PAPI, revisalo')
+
+
+    #esto es para establecer un criterio para la fluorescencia relativa
+    DistanciaFrecuenciaCociente = 15
+
+    if NormalizationCriterium==1:
+        FluoNormDivisor = 1
+
+    if NormalizationCriterium==2:
+        k = 0
+        while k < len(Freq):
+            if Freq[k] < detuningdoppler-DistanciaFrecuenciaCociente + 2 and Freq[k] > detuningdoppler-DistanciaFrecuenciaCociente - 2:
+                FluoIzquierda = Fluo[k]
+                print('Izq:', Freq[k])
+                break
+            else:
+                k = k + 1
+
+        l = 0
+        while l < len(Freq):
+            if Freq[l] < detuningdoppler+DistanciaFrecuenciaCociente + 2 and Freq[l] > detuningdoppler+DistanciaFrecuenciaCociente - 2:
+                FluoDerecha = Fluo[l]
+                print('Der: ', Freq[l])
+                break
+            else:
+                l = l + 1
+        FluoNormDivisor = 0.5*(FluoDerecha+FluoIzquierda)
+        print(FluoNormDivisor)
+
+    if NormalizationCriterium==3:
+        FluoNormDivisor = Fluo[0]
+
+    RelativeFluorescenceOfMinimums = np.array([Fluore/FluoNormDivisor for Fluore in FluorescenceOfMinimums])
+
+    return FrequencyOfMinimums, RelativeFluorescenceOfMinimums
+
+
+def GetPlotsofFluovsAngle_8levels(FrequencyOfMinimumsVector, RelativeFluorescenceOfMinimumsVector, u, titadoppler, detuningdoppler, detuningrepump, ventana=0.25, taketheoricalDR=False):
+
+    #primero buscamos las frecuencias referencia que se parezcan a las 6:
+    i = 0
+    FrecuenciasReferenciaBase = FrequencyOfMinimumsVector[0]
+
+    FrecuenciasDRTeoricas = [darkresonance for darkresonance in CalculoTeoricoDarkResonances_8levels(u, titadoppler, detuningdoppler, detuningrepump)[0]]
+
+    while i < len(FrequencyOfMinimumsVector):
+        if len(FrequencyOfMinimumsVector[i])==len(FrecuenciasDRTeoricas):
+            FrecuenciasReferenciaBase = FrequencyOfMinimumsVector[i]
+            print('Cool! Taking the DR identified with any curve')
+            break
+        else:
+            i = i + 1
+            if i==len(FrequencyOfMinimumsVector):
+                print('No hay ningun plot con 5 resonancias oscuras. Tomo las teóricas')
+                FrecuenciasReferenciaBase = FrecuenciasDRTeoricas
+
+    if taketheoricalDR:
+        FrecuenciasReferenciaBase = FrecuenciasDRTeoricas
+
+
+
+    Ventana = abs(ventana*(FrecuenciasReferenciaBase[1] - FrecuenciasReferenciaBase[0])) #ventana separadora de resonancias
+
+    print('Ventana = ', Ventana)
+    DarkResonance1Frequency = []
+    DarkResonance1Fluorescence = []
+
+    DarkResonance2Frequency = []
+    DarkResonance2Fluorescence = []
+
+    DarkResonance3Frequency = []
+    DarkResonance3Fluorescence = []
+
+    DarkResonance4Frequency = []
+    DarkResonance4Fluorescence = []
+
+    DarkResonance5Frequency = []
+    DarkResonance5Fluorescence = []
+
+    DarkResonance6Frequency = []
+    DarkResonance6Fluorescence = []
+
+    i = 0
+
+
+    while i < len(FrequencyOfMinimumsVector):
+        j = 0
+        FrecuenciasReferencia = [i for i in FrecuenciasReferenciaBase]
+
+
+        while j < len(FrequencyOfMinimumsVector[i]):
+
+
+            if abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[0])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[0])-Ventana):
+                DarkResonance1Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance1Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[0] = 0
+
+            elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[1])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[1])-Ventana):
+                DarkResonance2Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance2Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[1] = 0
+
+            elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[2])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[2])-Ventana):
+                DarkResonance3Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance3Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[2] = 0
+
+            elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[3])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[3])-Ventana):
+                DarkResonance4Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance4Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[3] = 0
+
+            elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[4])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[4])-Ventana):
+                DarkResonance5Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance5Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[4] = 0
+
+            elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[5])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[5])-Ventana):
+                DarkResonance6Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance6Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[5] = 0
+
+            else:
+                #print('Algo anduvo mal, por ahi tenes que cambiar la ventana che')
+                pass
+            j = j + 1
+
+        if np.count_nonzero(FrecuenciasReferencia) > 0:
+            if FrecuenciasReferencia[0] != 0:
+                DarkResonance1Frequency.append(FrecuenciasReferencia[0])
+                DarkResonance1Fluorescence.append()
+
+            if FrecuenciasReferencia[1] != 0:
+                DarkResonance2Frequency.append(FrecuenciasReferencia[1])
+                DarkResonance2Fluorescence.append(0)
+
+            if FrecuenciasReferencia[2] != 0:
+                DarkResonance3Frequency.append(FrecuenciasReferencia[2])
+                DarkResonance3Fluorescence.append(0)
+
+            if FrecuenciasReferencia[3] != 0:
+                DarkResonance4Frequency.append(FrecuenciasReferencia[3])
+                DarkResonance4Fluorescence.append(0)
+
+            if FrecuenciasReferencia[4] != 0:
+                DarkResonance5Frequency.append(FrecuenciasReferencia[4])
+                DarkResonance5Fluorescence.append(0)
+
+            if FrecuenciasReferencia[5] != 0:
+                DarkResonance6Frequency.append(FrecuenciasReferencia[5])
+                DarkResonance6Fluorescence.append(0)
+
+        i = i + 1
+
+    return DarkResonance1Frequency, DarkResonance1Fluorescence, DarkResonance2Frequency, DarkResonance2Fluorescence, DarkResonance3Frequency, DarkResonance3Fluorescence, DarkResonance4Frequency, DarkResonance4Fluorescence, DarkResonance5Frequency, DarkResonance5Fluorescence, DarkResonance6Frequency, DarkResonance6Fluorescence, FrecuenciasReferenciaBase
+
+
+
+def PerformExperiment_8levels(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None):
+
+    """
+    Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
+
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+
+    Fluovectors = []
+
+    for titaprobe in titaprobeVec:
+
+        tinicial = time.time()
+
+        ProbeDetuningVectorL, Fluovector = CPTspectrum8levels(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=False, solvemode=1)
+
+        tfinal = time.time()
+
+        print('Done angle ', titarepump, ' Total time: ', round((tfinal-tinicial), 2), "s")
+
+        if plot:
+            plt.figure()
+            plt.xlabel('Repump detuning (MHz')
+            plt.ylabel('Fluorescence (A.U.)')
+            plt.plot(ProbeDetuningVectorL, Fluovector, label=str(titarepump)+'º tita repump, T: ' + str(T*1e3) + ' mK')
+            plt.legend()
+
+        Fluovectors.append(Fluovector)
+
+        if len(titaprobeVec) == 1: #esto es para que no devuelva un vector de vectores si solo fijamos un angulo
+            Fluovectors = Fluovector
+
+    return ProbeDetuningVectorL, Fluovectors
+
+
+
+def PerformExperiment_8levels_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None):
+
+    """
+    Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
+
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+
+    Fluovectors = []
+
+    for titaprobe in titaprobeVec:
+
+        tinicial = time.time()
+
+        ProbeDetuningVectorL, Fluovector = CPTspectrum8levels_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=False, solvemode=1)
+
+        tfinal = time.time()
+
+        print('Done angle ', titarepump, ' Total time: ', round((tfinal-tinicial), 2), "s")
+
+        if plot:
+            plt.figure()
+            plt.xlabel('Repump detuning (MHz')
+            plt.ylabel('Fluorescence (A.U.)')
+            plt.plot(ProbeDetuningVectorL, Fluovector, label=str(titarepump)+'º tita repump, T: ' + str(T*1e3) + ' mK')
+            plt.legend()
+
+        Fluovectors.append(Fluovector)
+
+        if len(titaprobeVec) == 1: #esto es para que no devuelva un vector de vectores si solo fijamos un angulo
+            Fluovectors = Fluovector
+
+    return ProbeDetuningVectorL, Fluovectors

analisis/plots/20231123_CPTconmicromocion3/Data/EITfit/lolo_modelo_full_8niveles.py

@@ -7,6 +7,11 @@
 
 reingenieria del código que anda
 
+Surge de fusionar
+Data/EITfit/MM_eightLevel_2repumps_AnalysisFunctions.py
+Data/EITfit/MM_eightLevel_2repumps_python_scripts.py
+
+
 MAPA de FUNCIONES
 
 CPTspectrum8levels_MM

analisis/plots/20231123_CPTconmicromocion3/lolo_analisis_superajuste.py

 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """
-Ploteo de datos y ajustes
+Ploteo de datos y ajustes del barrido en voltaje del Endcap
+equivalente a pasear el ion en algun entorno del punto de compensación ideal
+y ver los efectos de micromoción (y tal vez temperatura)
+
+
+
 @author: lolo
 """
 
@@ -24,6 +29,7 @@ from time import time
 
 # /home/lolo/Dropbox/marce/LIAF/Trampa_anular/artiq_experiments/analisis/plots/20231123_CPTconmicromocion3/Data/EITfit/MM_eightLevel_2repumps_AnalysisFunctions.py
 
+
 from Data.EITfit.lolo_modelo_full_8niveles import PerformExperiment_8levels_MM
 
 

analisis/plots/20231123_CPTconmicromocion3/lolo_graficos_pub.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Ploteo de datos y ajustes del barrido en voltaje del Endcap
+equivalente a pasear el ion en algun entorno del punto de compensación ideal
+y ver los efectos de micromoción (y tal vez temperatura)
+
+
+
+@author: lolo
+"""
+
+
+import h5py
+import matplotlib.pyplot as plt
+import numpy as np
+
+from scipy.optimize import curve_fit
+from time import time
+from numba import jit,njit
+
+
+#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+#%% Carga de datos y librerías auxiliares
+
+
+## Librerias ##################################################################
+from Data.EITfit.lolo_modelo_full_8niveles import PerformExperiment_8levels_MM
+
+
+# PARAMETROS = np.load('PARAMETROS.npz',allow_pickle=True)
+# for var_name in PARAMETROS.keys():
+#     globals()[var_name] = PARAMETROS[var_name]
+#     print(f'loaded: {var_name}')
+
+
+# Funciones auxiliares
+
+from scipy.stats.distributions import  t,chi2
+
+def estadistica(datos_x,datos_y,modelo,pcov,parametros,nombres=None,alpha=0.05):
+    
+    if nombres is None:
+        nombres = [ f'{j}' for j in range(len(parametros)) ]
+    
+    # Cantidad de parámetros
+    P = len(parametros)
+    
+    # Número de datos
+    N = len(datos_x)
+    
+    # Grados de libertas (Degrees Of Freedom)
+    dof = N-P-1
+
+    # Cauculamos coordenadas del modelo
+    # modelo_x    = datos_x if modelo_x_arr is None else modelo_x_arr
+    # modelo_y    = modelo( modelo_x, *parametros )
+    
+    # Predicción del modelo para los datos_x medidos
+    prediccion_modelo = modelo( datos_x, *parametros )
+    
+    # Calculos de cantidades estadísticas relevantes
+    COV       = pcov                                      # Matriz de Covarianza
+    SE        = np.sqrt(np.diag( COV  ))                        # Standar Error / Error estandar de los parámetros
+    residuos  = datos_y - prediccion_modelo               # diferencia enrte el modelo y los datos
+    
+    SSE       = sum(( residuos )**2 )                     # Resitual Sum of Squares
+    SST       = sum(( datos_y - np.mean(datos_y))**2)        # Total Sum of Squares
+    
+    # http://en.wikipedia.org/wiki/Coefficient_of_determination
+    # Expresa el porcentaje de la varianza que logra explicar el modelos propuesto
+    Rsq       =  1 - SSE/SST                               # Coeficiente de determinación
+    Rsq_adj   = 1-(1-Rsq) * (N-1)/(N-P-1)                  # Coeficiente de determinación Ajustado   
+    
+    # https://en.wikipedia.org/wiki/Pearson_correlation_coefficient#In_least_squares_regression_analysis
+    # Expresa la correlación que hay entre los datos y la predicción del modelo
+    r_pearson = np.corrcoef( datos_y ,  prediccion_modelo )[0,1]
+    
+    # Reduced chi squared
+    # https://en.wikipedia.org/wiki/Reduced_chi-squared_statistic
+    chi2_     = sum( residuos**2 )/N
+    chi2_red  = sum( residuos**2 )/(N-P)
+    
+    # Chi squared test
+    chi2_test = sum( residuos**2 / abs(prediccion_modelo) )
+    # p-value del ajuste
+    p_val  = chi2(dof).cdf( chi2_test )
+    
+    
+    sT = t.ppf(1.0 - alpha/2.0, N - P ) # student T multiplier
+    CI = sT * SE                        # Confidence Interval
+    
+    print('R-squared    ',Rsq)
+    print('R-sq_adjusted',Rsq_adj)
+    print('chi2         ',chi2_)
+    print('chi2_reduced ',chi2_red)
+    print('chi2_test    ',chi2_test)
+    print('r-pearson    ',r_pearson)
+    print('p-value      ',p_val)
+    print('')
+    print('Error Estandard (SE):')
+    for i in range(P):
+        print(f'parametro[{nombres[i]:>5s}]: ' , parametros[i], ' ± ' , SE[i])
+    print('')
+    print('Intervalo de confianza al '+str((1-alpha)*100)+'%:')
+    for i in range(P):
+        print(f'parametro[{nombres[i]:>5s}]: ' , parametros[i], ' ± ' , CI[i])
+    
+    return dict(R2=Rsq,R2_adj=Rsq_adj,chi2=chi2_,chi2_red=chi2_red,
+                chi2_test=chi2_test,r=r_pearson,pvalue=p_val,
+                SE=SE,CI=CI)
+
+
+
+
+
+
+## Datos ######################################################################
+
+"""
+Primero tengo mediciones de espectros cpt de un ion variando la tension dc_A
+"""
+
+#C:\Users\Usuario\Documents\artiq\artiq_experiments\analisis\plots\20220106_CPT_DosLaseres_v08_TISA_DR\Data
+# os.chdir('../20231123_CPTconmicromocion3/Data/')
+
+folder = '../20231123_CPTconmicromocion3/Data/'
+CPT_FILES = f"""
+{folder}/000016262-IR_Scan_withcal_optimized
+{folder}/000016239-IR_Scan_withcal_optimized
+{folder}/000016240-IR_Scan_withcal_optimized
+{folder}/000016241-IR_Scan_withcal_optimized
+{folder}/000016244-IR_Scan_withcal_optimized
+{folder}/000016255-IR_Scan_withcal_optimized
+{folder}/000016256-IR_Scan_withcal_optimized
+{folder}/000016257-IR_Scan_withcal_optimized
+"""
+
+
+def SeeKeys(files):
+    for i, fname in enumerate(files.split()):
+        data = h5py.File(fname+'.h5', 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
+        print(fname)
+        print(list(data['datasets'].keys()))
+
+print(SeeKeys(CPT_FILES))
+
+#carpeta pc nico labo escritorio:
+#C:\Users\Usuario\Documents\artiq\artiq_experiments\analisis\plots\20211101_CPT_DosLaseres_v03\Data
+
+Counts = []
+Freqs = []
+
+AmpTisa = []
+UVCPTAmp = []
+No_measures = []
+Voltages = []
+
+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']['data_array']))
+    #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']))
+    Voltages.append(np.array(data['datasets']['scanning_voltages']))
+
+def Split(array,n):
+    length=len(array)/n
+    splitlist = []
+    jj = 0
+    while jj<length:
+        partial = []
+        ii = 0
+        while ii < n:
+            partial.append(array[jj*n+ii])
+            ii = ii + 1
+        splitlist.append(partial)
+        jj = jj + 1
+    return splitlist
+
+
+CountsSplit = []
+CountsSplit.append(Split(Counts[0],len(Freqs[0])))
+
+
+CountsSplit_2ions = []
+CountsSplit_2ions.append(Split(Counts[4],len(Freqs[4])))
+
+
+## Cargo parámetros fiteados de antes #########################################
+
+
+PARAMETROS = np.load('analisis_superajuste_PARAMETROS.npz',allow_pickle=True)
+for var_name in PARAMETROS.keys():
+    globals()[var_name] = PARAMETROS[var_name]
+    print(f'loaded: {var_name}')
+
+
+if False:
+    # Esto es para correr en caso de necesidad de limpiar todos los vectores de parametros
+    print('Limpio los vectores de parámetros')
+    for var in [ kk for kk in globals().keys() if kk.endswith('_vec') ]:
+        print(f'del {var}')
+        del(globals()[var])
+
+
+## Definiciones de Numba ######################################################
+
+@jit
+def FitEIT_MM_single_plot(Freqs, offset, DetDoppler, SG, SP, SCALE1, OFFSET, BETA1, TEMP):
+    #BETA = 1.8
+    # SG = 0.6
+    # SP = 8.1
+    # TEMP = 0.2e-3
+
+    freqs = [2*f*1e-6-offset for f in Freqs]
+
+    Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, 
+                                                            DopplerLaserLinewidth, ProbeLaserLinewidth, 
+                                                            TEMP, alpha, phidoppler, titadoppler, 
+                                                            phiprobe, titaprobe,  BETA1, drivefreq, 
+                                                            min(freqs), max(freqs)+(freqs[1]-freqs[0]), 
+                                                            freqs[1]-freqs[0], circularityprobe=CircPr, 
+                                                            plot=False, solvemode=1, detpvec=None)
+
+    ScaledFluo1 = np.array([f*SCALE1 + OFFSET for f in Fluorescence1])
+    return ScaledFluo1, Detunings
+
+
+
+@jit
+def FitEIT_MM_single(Freqs, offset, DetDoppler, SG, SP, SCALE1, OFFSET, BETA1, TEMP):
+    "Esta verison de la función devuelve sólo el eje y, para usar de modelo en un ajuste"
+    return FitEIT_MM_single_plot(Freqs, offset, DetDoppler, SG, SP, SCALE1, OFFSET, BETA1, TEMP)[0]
+
+
+
+param_names = 'offset DetDoppler SG SP SCALE1 OFFSET BETA1 TEMP'.split()
+
+#%%
+"""
+AHORA INTENTO SUPER AJUSTES O SEA CON OFFSETXPI Y DETDOPPLER INCLUIDOS
+
+La 0 no ajusta bien incluso con todos los parametros libres
+De la 1 a la 11 ajustan bien
+"""
+
+
+#from EITfit.MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels
+from scipy.optimize import curve_fit
+
+"""
+SUPER AJUSTE (SA)
+"""
+
+phidoppler, titadoppler = 0, 90
+phirepump, titarepump = 0,  0
+phiprobe = 0
+titaprobe = 90
+
+Temp = 0.5e-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 = 32.5e6
+
+#B = (u/(2*np.pi))/c
+
+correccion = 13
+
+#DetDoppler = -11.5-correccion
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6
+alpha = 0
+
+
+drivefreq = 2*np.pi*22.135*1e6
+
+
+SelectedCurveVec = [1,2,3,4,5,6,7,8,9,10,11]
+#SelectedCurveVec = [10]
+
+CircPr = 1
+alpha = 0
+
+
+
+
+
+def hiperbola(x,a,y0,b,x0):
+    """
+    Hiperbola de ecuación:
+        1 =(y-y0)²/a² - (x-x0)²/b² 
+    """
+    return a*np.sqrt(((x-x0)**2+b**2))+y0
+
+#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+#%% Graficamos todos los fiteos
+
+# tmp_datos=(Detuningsshort_vec,Counts_vec,Detuningslong_vec,FittedCounts_vec,SelectedCurveVec)
+
+# for Detunings_3_SA_short,CountsDR,Detunings_3_SA_long,FittedEITpi_3_SA_long,selectedcurve in zip(*tmp_datos):
+#     plt.figure()
+#     plt.errorbar(Detunings_3_SA_short, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
+#     plt.plot(Detunings_3_SA_long, FittedEITpi_3_SA_long, color='darkolivegreen', linewidth=3, label=f'med {selectedcurve}')
+#     #plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
+#     plt.xlabel('Detuning (MHz)')
+#     plt.ylabel('Counts')
+#     plt.legend(loc='upper left', fontsize=20)
+#     plt.grid()
+
+#     print(f'listo med {selectedcurve}')
+#     print(popt_3_SA)
+
+
+fig, axx = plt.subplots( 3,4, figsize=(13,8) ,  constrained_layout=True, sharex=True , sharey=True )
+fig.set_constrained_layout_pads(w_pad=2/72, h_pad=2/72, hspace=0, wspace=0)
+
+tmp_datos=(Detuningsshort_vec,Counts_vec,Detuningslong_vec,FittedCounts_vec,SelectedCurveVec,axx.flatten())
+
+for Detunings_3_SA_short,CountsDR,Detunings_3_SA_long,FittedEITpi_3_SA_long,selectedcurve,ax in zip(*tmp_datos):
+    ax.errorbar(Detunings_3_SA_short, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.3, capsize=2, markersize=2)
+    ax.plot(Detunings_3_SA_long, FittedEITpi_3_SA_long, color='black', linewidth=2, label=f'med {selectedcurve}', alpha=0.7)
+    #plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
+    # ax.set_xlabel('Detuning (MHz)')
+    # ax.set_ylabel('Counts')
+    ax.legend(loc='upper left', fontsize=12)
+    ax.grid(True, ls=":")
+
+    print(f'listo med {selectedcurve}')
+    # print(popt_3_SA)
+
+
+for ax in axx[:,0]:
+    ax.set_ylabel('Counts')
+
+for ax in axx[-1,:]:
+    ax.set_xlabel('Detuning (MHz)')
+
+
+
+
+
+
+
+
+#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+#%% Gráfico central 
+
+import matplotlib.pyplot as plt
+plt.rcParams.update({'font.size': 8})
+
+# Para esto hace falta:
+# sudo apt install dvipng
+
+plt.rcParams['text.usetex']=True
+# plt.rcParams['text.latex.unicode']=True
+# params= {'text.latex.preamble' : [r'\usepackage{amsmath}']}
+# plt.rcParams.update(params)
+
+
+# fig, axx = plt.subplots( 3,4, figsize=(8.6/2.54*2,4) ,  constrained_layout=True, sharex=True , sharey=True )
+# fig.set_constrained_layout_pads(w_pad=2/72, h_pad=2/72, hspace=0, wspace=0)
+
+
+## Layout ###########################
+from matplotlib import gridspec
+figsize=(8.6/2.54*2,4)
+width_ratios=[1,2,1]
+
+fig, axx = plt.subplots(ncols=3, nrows=3, sharex=True, sharey=False,   constrained_layout=True,
+                        figsize=figsize, gridspec_kw=dict(width_ratios=width_ratios))
+fig.set_constrained_layout_pads(w_pad=1/72, h_pad=1/72, hspace=0, wspace=0)
+fig.set_constrained_layout_pads(w_pad=0, h_pad=0, hspace=0, wspace=0)
+
+gs = axx[0, 0].get_gridspec()
+# remove the underlying axes
+for ax in axx[:,1]:
+    ax.remove()
+    
+
+gs2 = gridspec.GridSpec(4,3, width_ratios=width_ratios,left=0.16,right=0.9,bottom=0.12)
+
+ax_central = fig.add_subplot(gs2[:-1, 1])
+
+ax_res = fig.add_subplot(gs2[-1, 1], sharex=ax_central,zorder=-5)
+plt.setp(ax_central.get_xticklabels(), visible=False)
+
+for ax in axx[:,2]:
+    ax.yaxis.tick_right()
+    ax.yaxis.set_label_position("right")
+
+
+
+## CPTs ########################################
+selection= (0,1,3,4,7,8)
+axes_vec = (axx[0,0],axx[1,0],axx[2,0],axx[2,2],axx[1,2],axx[0,2])
+
+tmp_datos=(Detuningsshort_vec[selection,:],
+           Counts_vec[selection,:],
+           Detuningslong_vec[selection,:],
+           FittedCounts_vec[selection,:],
+           np.array(selection)+1,
+           axes_vec,'abcfed')
+
+for Detunings_3_SA_short,CountsDR,Detunings_3_SA_long,FittedEITpi_3_SA_long,selectedcurve,ax,le in zip(*tmp_datos):
+    ax.errorbar(Detunings_3_SA_short, CountsDR, yerr=2*np.sqrt(CountsDR), 
+                fmt='o', color='darkgreen', alpha=0.2, capsize=2, 
+                markersize=2, label='measured data')
+    ax.plot(Detunings_3_SA_long, FittedEITpi_3_SA_long, 
+            color='black', linewidth=2, label=f'micromotion model', alpha=0.7)
+    
+    ax.grid(True, ls=":", color='lightgray')
+    
+    ax.set_yticklabels([ f"{int(y/1000)}k" for y in ax.get_yticks() ])
+    
+    ax.set_title(f'({le})', x=0.1, y=0.78, color='gray')
+    # ax.legend()
+    # ax.legend(ax.get_legend_handles_labels()[0], f'{selectedcurve}')
+    print(f'listo med {selectedcurve}')
+
+
+axx[2,0].set_xlabel('Detuning [MHz]')
+axx[2,2].set_xlabel('Detuning [MHz]')
+
+axx[1,0].set_ylabel('Counts')
+axx[1,2].set_ylabel('Counts')
+
+
+
+
+
+## Grafico central ###############################
+I = slice(None,9)
+par_inicial = (12,0.1,1,-0.13)
+param,pcov = curve_fit(hiperbola,voltages_dcA[I],Betas_vec[I],p0=par_inicial)
+x_hip = np.linspace(-0.23,0.005,200)
+
+
+ax = ax_central
+ax.errorbar(voltages_dcA[I],Betas_vec[I],yerr=ErrorBetas_vec[I],fmt='o',
+            capsize=4,markersize=3,color='C3', label=r'fitted $\beta$')
+ax.plot(x_hip,hiperbola(xhip,*popthip),color='C0', label=r'hyperbola model')
+ax.set_ylabel(r'Modulation factor $\beta$', labelpad=-5)
+ax.set_ylim(-0.05,3)
+ax.set_xlim(-0.22,0)
+ax.set_title(f'(g)', x=0.95, y=0.006, color='gray')
+
+ax = ax_res
+ax.errorbar(voltages_dcA[I],Betas_vec[I]-hiperbola(voltages_dcA[I],*popthip),
+            yerr=ErrorBetas_vec[I],fmt='o',capsize=4,markersize=3,color='C3')
+ax.axhline( 0 , color='C0')
+ax.set_ylabel('Res.', labelpad=-5)
+ax.set_xlabel('Endcap voltage [V]')
+ax.set_title(f'(h)', x=0.95, y=0.72, color='gray')
+
+
+ax_res.get_xticklabels()[-1].set_visible(False)
+
+for ax in [ax_central,ax_res]:
+    ax.grid(True, ls=":", color='lightgray')
+
+
+# print([t*1e3 for t in Temp_vec])
+
+
+# Anotaciones ##########################3
+for jj,ax in zip(selection,axes_vec):
+    
+    Axes_x = 1 if axx.flatten().tolist().index(ax)%3==0 else 0
+    ax_central.annotate("",
+                         xy=(ax_central.get_lines()[0].get_xdata()[jj], ax_central.get_lines()[0].get_ydata()[jj]), 
+                         xycoords=ax_central.transData,
+                         xytext=(Axes_x, 0.5), textcoords=ax.transAxes,
+                         arrowprops=dict(arrowstyle="<-",connectionstyle="arc3,rad=-0.2", color='gray', alpha=0.5),
+                         zorder=-2)
+    
+
+
+# Leyenda ##############################
+h1, l1 = ax_central.get_legend_handles_labels()
+h2, l2 = axx[0,0].get_legend_handles_labels()
+ax_central.legend(h1+h2, l1+l2, loc='upper left')
+
+
+
+fig.align_ylabels([ax_central,ax_res])
+fig.tight_layout()
+
+
+fig.savefig('grafico_central_opcion_A.png', dpi=300)
+fig.savefig('grafico_central_opcion_A.pdf')
+
+
+
+
+
+
+
+# import matplotlib.pyplot as plt
+# from matplotlib import gridspec
+
+# fig = plt.figure()
+# gs = gridspec.GridSpec(1,2)
+# ax1 = fig.add_subplot(gs[0])
+# ax2 = fig.add_subplot(gs[1], sharey=ax1)
+# plt.setp(ax2.get_yticklabels(), visible=False)
+
+# plt.setp([ax1, ax2], title='Test')
+# fig.suptitle('An overall title', size=20)
+# gs.tight_layout(fig, rect=[0, 0, 1, 0.97])
+
+# plt.show()
+
+#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+#%% Endcap hiperbola (con residuos)
+"""
+Veamos cómo varía el Beta con el voltaje del endcap
+"""
+
+import seaborn as sns
+paleta = sns.color_palette("rocket")
+
+
+I = slice(None,9)
+
+voltages_dcA  = Voltages[0][SelectedCurveVec]
+
+
+
+def lineal(x,a,b):
+    return a*x+b
+
+def hiperbola(x,a,y0,b,x0):
+    """
+    Hiperbola de ecuación:
+        1 =(y-y0)²/a² - (x-x0)²/b² 
+    """
+    return a*np.sqrt(((x-x0)**2+b**2))+y0
+
+es_hiperbola = False
+
+    
+# par_inicial = (100,0.1,1,-0.15)
+#               a   y0   b  x0
+par_inicial = (12,0.1,1,-0.13)
+
+
+popthip,pcovhip = curve_fit(hiperbola,voltages_dcA[I],Betas_vec[I],p0=par_inicial)
+
+xhip = np.linspace(-0.23,0.005,200)
+
+# plt.figure()
+# plt.errorbar(voltages_dcA[I],Betas_vec[I],yerr=ErrorBetas_vec[I],fmt='o',capsize=5,markersize=5,color=paleta[1])
+# plt.plot(xhip,hiperbola(xhip,*popthip))
+# # plt.plot(xhip,hiperbola(xhip,*par_inicial),'--',color='red')
+# plt.xlabel('Endcap voltage (V)')
+# plt.ylabel('Modulation factor')
+# plt.grid()
+
+
+
+
+fig, axx = plt.subplots( 2, figsize=(10,7) ,  
+                        constrained_layout=True, sharex=True,
+                        gridspec_kw=dict(height_ratios=[10,2]))
+fig.set_constrained_layout_pads(w_pad=2/72, h_pad=2/72, hspace=0, wspace=0)
+
+ax = axx[0]
+
+ax.errorbar(voltages_dcA[I],Betas_vec[I],yerr=ErrorBetas_vec[I],fmt='o',capsize=5,markersize=5,color=paleta[1])
+ax.plot(xhip,hiperbola(xhip,*popthip))
+# plt.plot(xhip,hiperbola(xhip,*par_inicial),'--',color='red')
+
+ax.set_ylabel('Modulation factor')
+
+ax = axx[1]
+
+ax.errorbar(voltages_dcA[I],Betas_vec[I]-hiperbola(voltages_dcA[I],*popthip),
+            yerr=ErrorBetas_vec[I],fmt='o',capsize=5,markersize=5,color=paleta[1])
+
+ax.set_ylabel('Res.')
+ax.set_xlabel('Endcap voltage (V)')
+
+for ax in axx:
+    ax.grid(True, ls=":", color='lightgray')
+
+
+print([t*1e3 for t in Temp_vec])
+
+
+#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+#%% Este hay que armarlo aún
+
+plt.figure()
+plt.errorbar(voltages_dcA,[t*1e3 for t in Temp_vec],yerr=[t*1e3 for t in ErrorTemp_vec],fmt='o',capsize=5,markersize=5,color=paleta[3])
+# plt.axvline(minimum_voltage,linestyle='dashed',color='grey')
+print(f'\n\nTE FALTA DEFINIR LA VARIABLE minimum_voltage\n\n')
+plt.axhline(0.538)
+plt.xlabel('Endcap voltage (V)')
+plt.ylabel('Temperature (mK)')
+plt.grid()
+#plt.ylim(0,2)
+
+
+
+#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+#%% Ajuste de los betas y la temperatura
+
+from scipy.special import jv
+
+
+def expo(x,tau,A,B):
+    return A*np.exp(x/tau)+B
+
+def cuadratica(x,a,c):
+    return a*(x**2)+c
+
+def InverseMicromotionSpectra(beta, A, det, x0, gamma, B):
+    ftrap=22.1
+    #gamma=30
+    P = ((jv(0, beta)**2)/((((det-x0)**2)+(0.5*gamma)**2)**2))*(-2*(det-x0))
+    i = 1
+    #print(P)
+    while i <= 5:
+        P = P + (-2*(det-x0))*((jv(i, beta))**2)/(((((det-x0)+i*ftrap)**2)+(0.5*gamma)**2)**2) + (-2*(det-x0))*(((jv(-i, beta))**2)/((((det-x0)-i*ftrap)**2)+(0.5*gamma)**2)**2)
+        i = i + 1
+        #print(P)
+    #return 1/(A*P+B)
+    return 1/(A*P+B)
+
+
+def InverseMicromotionSpectra_raw(beta, A, det, B):
+    ftrap=22.1
+    gamma=21
+    P = ((jv(0, beta)**2)/((((det)**2)+(0.5*gamma)**2)**2))*(-2*(det))
+    i = 1
+    #print(P)
+    while i <= 3:
+        P = P + (-2*(det))*((jv(i, beta))**2)/(((((det)+i*ftrap)**2)+(0.5*gamma)**2)**2) + (-2*(det))*(((jv(-i, beta))**2)/((((det)-i*ftrap)**2)+(0.5*gamma)**2)**2)
+        i = i + 1
+        #print(P)
+    return A/P+B
+
+
+"""
+Temperatura vs beta con un ajuste exponencial 
+"""
+
+popt_exp, pcov_exp = curve_fit(expo,Betas_vec[:10],[t*1e3 for t in Temp_vec[:10]])
+popt_quad, pcov_quad = curve_fit(cuadratica,Betas_vec[:10],[t*1e3 for t in Temp_vec[:10]],p0=(1,10))
+#popt_rho22, pcov_rho22 = curve_fit(InverseMicromotionSpectra,Betas_vec,[t*1e3 for t in Temp_vec],p0=(10,10,-10,1,20)) #esto ajusta muy bien
+#popt_rho22, pcov_rho22 = curve_fit(InverseMicromotionSpectra,Betas_vec, [t*1e3 for t in Temp_vec],p0=(-10,-10,10,1,20)) #esto ajusta muy bien
+popt_rho22_raw, pcov_rho22_raw = curve_fit(InverseMicromotionSpectra_raw,Betas_vec[:10], [t*1e3 for t in Temp_vec[:10]],p0=(-10, -10, 1)) #esto ajusta muy bien
+
+
+print(popt_rho22_raw)
+
+betaslong = np.arange(0,2*2.7,0.01)
+
+print(f'Min temp predicted: {InverseMicromotionSpectra_raw(betaslong,*popt_rho22_raw)[100]}')
+
+plt.figure()
+plt.errorbar(Betas_vec[:10],[t*1e3 for t in Temp_vec[:10]],xerr=ErrorBetas_vec[:10], yerr=[t*1e3 for t in ErrorTemp_vec[:10]],fmt='o',capsize=5,markersize=5,color=paleta[3])
+#plt.plot(betaslong,expo(betaslong,*popt_exp),label='Ajuste exponencial')
+#plt.plot(betaslong,cuadratica(betaslong,*popt_quad),label='Ajuste cuadratico')
+#plt.plot(betaslong,InverseMicromotionSpectra(betaslong,*popt_rho22),label='Ajuste cuadratico')
+plt.plot(betaslong,InverseMicromotionSpectra_raw(betaslong,*popt_rho22_raw),label='Ajuste cuadratico')
+
+#plt.axvline(minimum_voltage,linestyle='dashed',color='grey')
+#plt.axhline(0.538)
+plt.xlabel('Modulation factor')
+plt.ylabel('Temperature (mK)')
+plt.grid()
+
+
+#%%
+"""
+Esto no es del super ajuste sino de los ajustes anteriores en donde DetDoppler y offset son puestos a mano
+Aca grafico los betas con su error en funcion de la tension variada.
+Ademas, hago ajuste lineal para primeros y ultimos puntos, ya que espero que
+si la tension hace que la posicion del ion varie linealmente, el beta varia proporcional a dicha posicion.
+"""
+
+import seaborn as sns
+
+def lineal(x,a,b):
+    return a*x+b
+
+paleta = sns.color_palette("rocket")
+
+betavector = [beta1,beta2,beta3,beta4,beta5,beta6,beta7,beta8,beta9]
+errorbetavector = [errorbeta1,errorbeta2,errorbeta3,errorbeta4,errorbeta5,errorbeta6,errorbeta7,errorbeta8,errorbeta9]
+voltages_dcA = Voltages[0][1:10]
+
+poptini,pcovini = curve_fit(lineal,voltages_dcA[0:3],betavector[0:3])
+poptfin,pcovfin = curve_fit(lineal,voltages_dcA[4:],betavector[4:])
+
+minimum_voltage = -(poptini[1]-poptfin[1])/(poptini[0]-poptfin[0]) #voltaje donde se intersectan las rectas, es decir, donde deberia estar el minimo de micromocion
+minimum_modulationfactor = lineal(minimum_voltage,*poptini) #es lo mismo si pongo *poptfin
+
+xini = np.linspace(-0.23,-0.13,100)
+xfin = np.linspace(-0.15,0.005,100)
+
+plt.figure()
+plt.errorbar(voltages_dcA,betavector,yerr=errorbetavector,fmt='o',capsize=5,markersize=5,color=paleta[1])
+plt.plot(xini,lineal(xini,*poptini))
+plt.plot(xfin,lineal(xfin,*poptfin))
+plt.axvline(minimum_voltage,linestyle='dashed',color='grey')
+plt.xlabel('Endcap voltage (V)')
+plt.ylabel('Modulation factor')
+plt.grid()
+
+
+#%%
+"""
+Aca veo la temperatura del ion en funcion del voltaje del endcap, ya que
+al cambiar la cantidad de micromocion, cambia la calidad del enfriado
+"""
+tempvector = np.array([temp1,temp2,temp3,temp4,temp5,temp6,temp7,temp8,temp9])*1e3
+errortempvector = np.array([errortemp1,errortemp2,errortemp3,errortemp4,errortemp5,errortemp6,errortemp7,errortemp8,errortemp9])*1e3
+
+voltages_dcA = Voltages[0][1:10]
+
+plt.figure()
+plt.errorbar(voltages_dcA,tempvector,yerr=errortempvector,fmt='o',capsize=5,markersize=5,color=paleta[3])
+plt.axvline(minimum_voltage,linestyle='dashed',color='grey')
+plt.xlabel('Endcap voltage (V)')
+plt.ylabel('Temperature (mK)')
+plt.grid()
+plt.ylim(0,2)
+
+
+#%%
+"""
+Por las dudas, temperatura en funcion de beta
+"""
+plt.figure()
+plt.errorbar(betavector,tempvector,yerr=errortempvector,xerr=errorbetavector,fmt='o',capsize=5,markersize=5)
+plt.xlabel('Modulation factor')
+plt.ylabel('Temperature (mK)')
+plt.grid()
+
+
+
+#%%
+"""
+Si quiero ver algun parametro del ajuste puntual. el orden es: 0:SG, 1:SP, 2:SCALE1, 3:OFFSET
+"""
+ki=2
+plt.errorbar(np.arange(0,9,1),[popt_1[ki],popt_2[ki],popt_3[ki],popt_4[ki],popt_5[ki],popt_6[ki],popt_7[ki],popt_8[ki],popt_9[ki]],yerr=[np.sqrt(pcov_1[ki,ki]),np.sqrt(pcov_2[ki,ki]),np.sqrt(pcov_3[ki,ki]),np.sqrt(pcov_4[ki,ki]),np.sqrt(pcov_5[ki,ki]),np.sqrt(pcov_6[ki,ki]),np.sqrt(pcov_7[ki,ki]),np.sqrt(pcov_8[ki,ki]),np.sqrt(pcov_9[ki,ki])], fmt='o',capsize=3,markersize=3)
+
+
+#%%
+
+if False:
+    GUARDAR = {}
+    # for var in [ kk for kk in globals().keys() if kk.startswith('pop') ]:
+    #     print(var)
+    #     GUARDAR[var] = globals()[var]
+    # print('')
+    # for var in [ kk for kk in globals().keys() if kk.startswith('pcov') ]:
+    #     print(var)
+    #     GUARDAR[var] = globals()[var]
+
+    # print('')
+    # for var in [ kk for kk in globals().keys() if kk.startswith('Fitted') ]:
+    #     print(var)
+    #     GUARDAR[var] = globals()[var]
+    # print('')
+    for var in [ kk for kk in globals().keys() if kk.endswith('_vec') ]:
+        print(var)
+        GUARDAR[var] = globals()[var]
+
+
+    np.savez('analisis_superajuste_PARAMETROS.npz', **GUARDAR )

analisis/plots/20231212_Bstability/lolo_CPT_plotter_20231212.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Acá hay un análisis de la estabilidad del campo magnético (presuntamente).
+Se midieron N veces el mismo espectro y se analiza la variación del spliting
+de picos CPT que se debe al efecto Zeeman.
+
+El gráfico final muestra el corrimiento en procentual
+
+
+@author: lolo
+"""
+
+
+
+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 de una resonancia oscura DD multiples veces a lo largo de una noche para ver estabilidad de B
+"""
+
+#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/20231212_Bstability/Data/')
+
+# CPT_FILES = """000016432-IR_Scan_withcal_optimized
+# 000016433-IR_Scan_withcal_optimized
+# 000016434-IR_Scan_withcal_optimized
+# 000016435-IR_Scan_withcal_optimized
+# 000016436-IR_Scan_withcal_optimized
+# 000016437-IR_Scan_withcal_optimized
+# 000016438-IR_Scan_withcal_optimized
+# 000016439-IR_Scan_withcal_optimized
+# 000016440-IR_Scan_withcal_optimized
+# 000016441-IR_Scan_withcal_optimized
+# 000016442-IR_Scan_withcal_optimized
+# 000016443-IR_Scan_withcal_optimized
+# """ 
+
+folder = '../20231212_Bstability/Data/'
+CPT_FILES = f"""
+{folder}/000016434-IR_Scan_withcal_optimized
+{folder}/000016435-IR_Scan_withcal_optimized
+{folder}/000016436-IR_Scan_withcal_optimized
+{folder}/000016437-IR_Scan_withcal_optimized
+{folder}/000016438-IR_Scan_withcal_optimized
+{folder}/000016439-IR_Scan_withcal_optimized
+{folder}/000016440-IR_Scan_withcal_optimized
+{folder}/000016441-IR_Scan_withcal_optimized
+{folder}/000016442-IR_Scan_withcal_optimized
+{folder}/000016443-IR_Scan_withcal_optimized
+""".strip()
+
+
+CALIB_FILES = f"""{folder}/000016430-IR_Scan_withcal_optimized"""
+
+
+
+def SeeKeys(files):
+    for i, fname in enumerate(files.split()):
+        data = h5py.File(fname+'.h5', 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
+        print(fname)
+        print(list(data['datasets'].keys()))
+
+print(SeeKeys(CPT_FILES))
+
+#carpeta pc nico labo escritorio:
+#C:\Users\Usuario\Documents\artiq\artiq_experiments\analisis\plots\20211101_CPT_DosLaseres_v03\Data
+
+Counts = []
+Freqs = []
+
+CalibCounts = []
+CalibFreqs = []
+
+AmpTisa = []
+UVCPTAmp = []
+No_measures = []
+Voltages = []
+
+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']['data_array']))
+    #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']))
+    Voltages.append(np.array(data['datasets']['scanning_voltages']))
+
+
+for i, fname in enumerate(CALIB_FILES.split()):
+    print(str(i) + ' - ' + fname)
+    data = h5py.File(fname+'.h5', 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
+
+    CalibFreqs.append(np.array(data['datasets']['IR1_Frequencies']))
+    CalibCounts.append(np.array(data['datasets']['counts_spectrum']))
+
+
+def Split(array,n):
+    length=len(array)/n
+    splitlist = []
+    jj = 0
+    while jj<length:
+        partial = []
+        ii = 0
+        while ii < n:
+            partial.append(array[jj*n+ii])
+            ii = ii + 1
+        splitlist.append(partial)
+        jj = jj + 1
+    return splitlist
+
+
+CountsSplit = []
+k=0
+for k in range(len(Counts)):
+    CountsSplit.append(Split(Counts[k],len(Freqs[k])))
+
+#%%
+from scipy.optimize import curve_fit
+
+def lorentzian(x,A,B,x0,g,C):
+    return 2*(A/np.pi)*g/(g**2 + 4*(x-x0)**2)+B+C*(x-x0)
+
+Freqscal = [2*f*1e-6 for f in CalibFreqs[0]]
+Countscal = CalibCounts[0]
+
+popt_dr1, pcov_dr1 = curve_fit(lorentzian,Freqscal[37:47],Countscal[37:47],p0=(-1000,1000,436,1,1))
+popt_dr2, pcov_dr2 = curve_fit(lorentzian,Freqscal[90:120],Countscal[90:120],p0=(-1000,1000,443,1,1))
+
+DeltaFreqs = popt_dr2[2]-popt_dr1[2]
+ZeroFrequency = 0.5*(popt_dr2[2]+popt_dr1[2])
+
+plt.figure()
+plt.plot(Freqscal,Countscal,'o')
+plt.plot(Freqscal,lorentzian(Freqscal,*popt_dr1))
+plt.plot(Freqscal,lorentzian(Freqscal,*popt_dr2))
+plt.axvline(ZeroFrequency)
+
+
+print(DeltaFreqs)
+
+"""
+Estas cuentas estan en el cuaderno SMILE MORE WORRY LESS pag 25.
+La resonancia de la izquierda esta a (-4/5)*u. La de la derecha esta a (4/5)*u. 
+Por ende la diferencia es (8/5)*u.
+Definimos u como 1.4 MHz/G * B. Entonces Despejamos B facilmente.
+"""
+
+ub = 9.27e-24
+h = 6.63e-34
+u = 1e-6*(ub/h)*1e-4 #en unidades de MHz/G
+
+MagneticField = DeltaFreqs/((8/5)*u)
+
+print(f'Magnetic field: {MagneticField}')
+
+
+
+#%%
+
+"""
+Ploteo la cpt de referencia / plotting the reference CPT
+"""
+
+freqs = [2*f*1e-6 for f in Freqs[0]]
+
+def lorentzian(x,A,B,x0,g,C):
+    return 2*(A/np.pi)*g/(g**2 + 4*(x-x0)**2)+B+C*(x-x0)
+
+# ii_plot = 11
+# jj_plot = 0
+
+ii_plot = 9
+jj_plot = 0
+
+ii_problematic = []
+jj_problematic = []
+
+Centers = []
+Widths = []
+test = []
+for ii in range(len(CountsSplit)):
+    for jj in range(len(CountsSplit[0])):    
+#        print(ii)
+#        print(jj)            
+        try:
+            if ii==2 and jj==11:
+                popt_lorentz, pcov_lorentz = curve_fit(lorentzian, freqs[:-10], CountsSplit[ii][jj][:-10],p0=(-1000,1000,436,1,1))
+            
+            elif ii==2 and jj==12:
+                popt_lorentz, pcov_lorentz = curve_fit(lorentzian, freqs[40:], CountsSplit[ii][jj][40:],p0=(-1000,1000,436,1,1))
+            
+            elif ii==4 and jj==1:
+                popt_lorentz, pcov_lorentz = curve_fit(lorentzian, freqs[:-86], CountsSplit[ii][jj][:-86],p0=(-1000,1000,436,1,1))
+            
+            elif ii==4 and jj==2:
+                popt_lorentz = [0,0,0,0,0]
+
+            elif ii==4 and jj==7:
+                popt_lorentz = [0,0,0,0,0]
+
+            elif ii==4 and jj==12:
+                popt_lorentz = [0,0,0,0,0]
+    
+            elif ii==4 and jj==13:
+                popt_lorentz = [0,0,0,0,0]
+
+            elif ii==4 and jj==14:
+                popt_lorentz = [0,0,0,0,0]
+
+            elif ii==11 and jj==2:
+                popt_lorentz = [0,0,0,0,0]                
+
+            elif ii==11 and jj==3:
+                popt_lorentz = [0,0,0,0,0]                
+            
+            else:
+                popt_lorentz, pcov_lorentz = curve_fit(lorentzian, freqs, CountsSplit[ii][jj],p0=(-1000,1000,436,1,1))
+                
+            if popt_lorentz[2]>435.95 or popt_lorentz[2]<435.8:
+                if popt_lorentz[2]==0:
+                    pass
+                else:
+                    ii_problematic.append(ii)
+                    jj_problematic.append(jj)
+                    
+        except:
+            popt_lorentz=[0,0,0,0]
+            print("except")
+        print(ii,jj)
+        if ii == ii_plot and jj == jj_plot:
+            print('MATCH')
+            test.append(popt_lorentz)
+        Centers.append(popt_lorentz[2])
+        Widths.append(popt_lorentz[3])
+
+prob = 4
+
+print(ii_problematic[prob])
+print(jj_problematic[prob])
+
+kk=-83
+
+plt.figure()
+plt.plot(freqs, CountsSplit[ii_problematic[prob]][jj_problematic[prob]])
+plt.plot(freqs[kk], CountsSplit[ii_problematic[prob]][jj_problematic[prob]][kk],'o',markersize=10)
+plt.plot(freqs,lorentzian(freqs,*test[0]))
+#%%
+
+"""
+Usando que la DR de la izquierda esta a (-4/5)u, donde u = 1.4 MHz/G * B, 
+despejo y convierto la posicion de esa resonancia a campo magnetico
+"""
+
+def ConvertFreqsToMagneticField(f,zerofreq,u):
+    return np.abs(f-zerofreq)*(5/4)/(1.4)
+
+
+lentotal = len(CountsSplit)*len(CountsSplit[0])
+medtime=4/60
+timevec = np.linspace(0,medtime*lentotal, lentotal)
+
+plt.figure()
+plt.plot(timevec[4:],ConvertFreqsToMagneticField(Centers,ZeroFrequency,u)[4:],'o')
+plt.ylim(3.670,3.730)
+plt.xlabel('Time (h)')
+plt.ylabel('Magnetic field (G)')
+
+plt.figure()
+plt.plot(timevec[4:],[100*c/3.718 for c in ConvertFreqsToMagneticField(Centers,ZeroFrequency,u)][4:],'o')
+plt.ylim(98.5,100.1)
+plt.xlabel('Time (h)')
+plt.ylabel('Magnetic field variation (percent)')
+
+

analisis/plots/20231214_CPTconmicromocioncristals/lolo_CPT_plotter_20231214.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+NO entiendo bien que es esto
+
+
+@author: lolo
+"""
+
+
+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
+
+"""
+Primero tengo mediciones de espectros cpt de un ion variando la tension dc_A
+"""
+
+#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/20231214_CPTconmicromocioncristals/Data/')
+
+
+
+folder = '../20231214_CPTconmicromocioncristals/Data/Single/'
+
+SINGLECPT_FILES = f"""
+{folder}/000016453-IR_Scan_withcal_optimized
+{folder}/000016454-IR_Scan_withcal_optimized
+{folder}/000016455-IR_Scan_withcal_optimized
+{folder}/000016456-IR_Scan_withcal_optimized
+{folder}/000016457-IR_Scan_withcal_optimized
+{folder}/000016458-IR_Scan_withcal_optimized
+{folder}/000016459-IR_Scan_withcal_optimized
+{folder}/000016461-IR_Scan_withcal_optimized
+""".strip()
+
+
+
+folder = '../20231214_CPTconmicromocioncristals/Data/Multi'
+
+MULTICPT_FILES = f"""
+{folder}/000016460-IR_Scan_withcal_optimized
+{folder}/000016462-IR_Scan_withcal_optimized
+{folder}/000016463-IR_Scan_withcal_optimized
+{folder}/000016464-IR_Scan_withcal_optimized
+""".strip()
+
+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(MULTICPT_FILES))
+
+#carpeta pc nico labo escritorio:
+#C:\Users\Usuario\Documents\artiq\artiq_experiments\analisis\plots\20211101_CPT_DosLaseres_v03\Data
+
+SingleCounts = []
+SingleFreqs = []
+MultiCounts = []
+MultiFreqs = []
+
+AmpTisa = []
+UVCPTAmp = []
+No_measures = []
+Voltages = []
+
+for i, fname in enumerate(MULTICPT_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.
+    MultiFreqs.append(np.array(data['datasets']['IR1_Frequencies']))
+    MultiCounts.append(np.array(data['datasets']['data_array']))
+    #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']))
+    Voltages.append(np.array(data['datasets']['scanning_voltages']))
+    
+for i, fname in enumerate(SINGLECPT_FILES.split()):
+    print(str(i) + ' - ' + fname)
+    #print(fname)
+    data = h5py.File(fname+'.h5', 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
+    SingleFreqs.append(np.array(data['datasets']['IR1_Frequencies']))
+    SingleCounts.append(np.array(data['datasets']['data_array']))
+
+
+def Split(array,n):
+    length=len(array)/n
+    splitlist = []
+    jj = 0
+    while jj<length:
+        partial = []
+        ii = 0
+        while ii < n:
+            partial.append(array[jj*n+ii])
+            ii = ii + 1
+        splitlist.append(partial)
+        jj = jj + 1
+    return splitlist
+
+
+CountsSplit = []
+
+for kk in range(len(MultiCounts)):
+    CountsSplit.append(Split(MultiCounts[kk],len(MultiFreqs[kk])))
+
+
+#%%
+
+"""
+Ploteo la cpt de referencia / plotting the reference CPT
+"""
+
+jvec = [0] # de la 1 a la 9 vale la pena, despues no
+
+
+plt.figure()
+i = 0
+for j in jvec:
+    plt.errorbar([2*f*1e-6 for f in SingleFreqs[j]], SingleCounts[j], yerr=np.sqrt(SingleCounts[j]), fmt='o', capsize=2, markersize=2)
+    i = i + 1
+plt.xlabel('Frecuencia (MHz)')
+plt.ylabel('counts')
+plt.grid()
+#for dr in drs:
+#    plt.axvline(dr)
+    #plt.axvline(dr+drive)
+plt.legend()
+
+
+
+#%%
+
+"""
+Ploteo curvas de la multi1
+
+meds:
+    0: dcA, 11 voltajes
+    1: dcA, 21 voltajes
+    2: compOven, 21 voltajes
+    3: dcA, 31 voltajes
+
+"""
+
+med=0
+# jvec = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21]
+
+jvec = [ ii for ii in range(len(CountsSplit[med])) ]
+
+kk=9
+
+plt.figure()
+i = 0
+for j in jvec:
+    plt.errorbar([2*f*1e-6 for f in MultiFreqs[med]], CountsSplit[med][j], yerr=np.sqrt(CountsSplit[med][j]), fmt='o', capsize=2, markersize=2)
+    #plt.plot([2*f*1e-6 for f in MultiFreqs[med]][kk], CountsSplit[med][j][kk],'o',markersize=10)
+    i = i + 1
+plt.xlabel('Frecuencia (MHz)')
+plt.ylabel('counts')
+plt.grid()
+#for dr in drs:
+#    plt.axvline(dr)
+    #plt.axvline(dr+drive)
+plt.legend()
+
+print(CountsSplit[med][j][9])
+print(CountsSplit[med][j][10])
+print(CountsSplit[med][j][11])
+print(CountsSplit[med][j][12])
+
+
+
+#%%
+#from EITfit.MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels
+from scipy.optimize import curve_fit
+import time
+
+cwd = os.getcwd()
+os.chdir('../20231123_CPTconmicromocion3')
+from Data.EITfit.lolo_modelo_full_8niveles import PerformExperiment_8levels_MM
+os.chdir(cwd)
+
+# from Data.EITfit.lolo_modelo_full_8niveles import PerformExperiment_8levels_MM
+
+
+
+"""
+MEDICION 1: ajusto una curva con dos iones con un modelo que considera solo uno
+"""
+
+phidoppler, titadoppler = 0, 90
+phirepump, titarepump = 0,  0
+phiprobe = 0
+titaprobe = 90
+
+Temp = 0.5e-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 = 32.5e6
+
+#B = (u/(2*np.pi))/c
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 
+alpha = 0
+
+
+drivefreq = 2*np.pi*22.135*1e6
+
+selectedcurve=0
+
+FreqsDR = SingleFreqs[selectedcurve]
+CountsDR = SingleCounts[selectedcurve]
+   
+freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+CircPr = 1
+alpha = 0
+
+def FitEIT_MM_1ion(Freqs, offset, DetDoppler, SG, SP, SCALE1, OFFSET, BETA1, TEMP, plot=False):
+#def FitEIT_MM(freqs, SG, SP, SCALE1, OFFSET, BETA1):
+    #BETA = 1.8
+    # SG = 0.6
+    # SP = 8.1
+    # TEMP = 0.2e-3
+
+    freqs = [2*f*1e-6-offset for f in Freqs]
+  
+    Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+
+    ScaledFluo1 = np.array([f*SCALE1 + OFFSET for f in Fluorescence1])
+    if plot:
+        return ScaledFluo1, Detunings
+    else:
+        return ScaledFluo1
+    #return ScaledFluo1
+
+do_fit = True
+if do_fit:
+    popt_1, pcov_1 = curve_fit(FitEIT_MM_1ion, FreqsDR, CountsDR, p0=[430, -25, 0.9, 6.2, 3e4, 1.34e3, 2, (np.pi**2)*1e-3], bounds=((0, -50, 0, 0, 0, 0, 0, 0), (1000, 0, 2, 20, 5e6, 5e4, 10, (np.pi**2)*10e-3)))
+
+    FittedEITpi_1_short, Detunings_1_short = FitEIT_MM_1ion(FreqsDR, *popt_1, plot=True)
+    freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+    FittedEITpi_1_long, Detunings_1_long = FitEIT_MM_1ion(freqslong, *popt_1, plot=True)
+
+
+plt.figure()
+plt.errorbar(Detunings_1_short, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
+plt.plot(Detunings_1_long, FittedEITpi_1_long, color='darkolivegreen', linewidth=3, label='med 1')
+#plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
+plt.xlabel('Detuning (MHz)')
+plt.ylabel('Counts')
+plt.legend(loc='upper left', fontsize=20)
+plt.grid()
+
+
+
+#%%
+#from EITfit.MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels
+from scipy.optimize import curve_fit
+import time
+
+"""
+MEDICION 1: ahora la ajusto pero considerando contribucion de dos iones
+"""
+
+phidoppler, titadoppler = 0, 90
+phirepump, titarepump = 0,  0
+phiprobe = 0
+titaprobe = 90
+
+Temp = 0.5e-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 = 32.5e6
+
+#B = (u/(2*np.pi))/c
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 
+alpha = 0
+
+
+drivefreq = 2*np.pi*22.135*1e6
+
+selectedcurve=0
+
+FreqsDR = SingleFreqs[selectedcurve]
+CountsDR = SingleCounts[selectedcurve]
+   
+#freqslong = np.arange(min(FreqsDR)*0.2, max(FreqsDR)*2+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+CircPr = 1
+alpha = 0
+
+def FitEIT_MM_2ion(Freqs, offset, DetDoppler, SG, SP, SCALE1, SCALE2, OFFSET1, OFFSET2, BETA1, BETA2, TEMP, plot=False):
+#def FitEIT_MM(freqs, SG, SP, SCALE1, OFFSET, BETA1):
+    #BETA = 1.8
+    # SG = 0.6
+    # SP = 8.1
+    # TEMP = 0.2e-3
+
+    freqs = [2*f*1e-6-offset for f in Freqs]
+  
+    Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+    Detunings, Fluorescence2 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA2, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+
+
+    ScaledFluo1 = np.array([f*SCALE1 + OFFSET1 for f in Fluorescence1])
+    ScaledFluo2 = np.array([f*SCALE2 + OFFSET2 for f in Fluorescence2])
+    
+    if plot:
+        return ScaledFluo1+ScaledFluo2, Detunings
+    else:
+        return ScaledFluo1+ScaledFluo2
+    #return ScaledFluo1
+
+do_fit = True
+if do_fit:
+    popt_1_2ion, pcov_1_2ion = curve_fit(FitEIT_MM_2ion, FreqsDR, CountsDR, p0=[445, -32, 0.5, 7, 2e4, 1e4, 2e3, 1.5e3, 2, 1, 0.5e-3], bounds=((0, -50, 0, 0, 0, 0, 0,0, 0,0, 0), (1000, 0, 2, 20, 5e6, 5e6, 5e4,5e4, 10, 10,20e-3)))
+
+    FittedEITpi_1_short_2ion, Detunings_1_short_2ion = FitEIT_MM_2ion(FreqsDR, *popt_1_2ion, plot=True)
+freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+FittedEITpi_1_long_2ion, Detunings_1_long_2ion = FitEIT_MM_2ion(freqslong, *popt_1_2ion, plot=True)
+
+plt.figure()
+plt.errorbar(Detunings_1_short_2ion, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
+plt.plot(Detunings_1_long_2ion, FittedEITpi_1_long_2ion, color='darkolivegreen', linewidth=3, label='med 1')
+#plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
+plt.xlabel('Detuning (MHz)')
+plt.ylabel('Counts')
+#plt.xlim(-80,50)
+plt.legend(loc='upper left', fontsize=20)
+plt.grid()
+
+
+
+#%%
+#from EITfit.MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels
+from scipy.optimize import curve_fit
+import time
+
+"""
+MEDICION MULTI INDIVIDUAL
+
+VEO EL AJUSTE DE UNA DEL AS CURVAS MULTI PARA VER COMO AJUSTA
+
+"""
+
+phidoppler, titadoppler = 0, 90
+phirepump, titarepump = 0,  0
+phiprobe = 0
+titaprobe = 90
+
+Temp = 0.5e-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 = 32.5e6
+
+#B = (u/(2*np.pi))/c
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 
+alpha = 0
+
+
+drivefreq = 2*np.pi*22.135*1e6
+
+measurement = 2
+#selectedcurve=10
+selectedcurvevec=[10]
+
+#popt_vecs = []
+#pcov_vecs = []
+
+for selectedcurve in selectedcurvevec:
+    FreqsDR = MultiFreqs[measurement]
+    CountsDR = CountsSplit[measurement][selectedcurve]
+       
+    if selectedcurve==9 and measurement==1:
+        CountsDR[10]=4132+89
+        CountsDR[11]=4132+2*89
+
+    freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+    CircPr = 1
+    alpha = 0
+    
+    
+    
+    def FitEIT_MM_2ion(Freqs, offset, DetDoppler, SG, SP, SCALE1, SCALE2, OFFSET, BETA1, BETA2, TEMP, U, plot=False):
+    #def FitEIT_MM(freqs, SG, SP, SCALE1, OFFSET, BETA1):
+        #BETA = 1.8
+        # SG = 0.6
+        # SP = 8.1
+        # TEMP = 0.2e-3
+    
+        freqs = [2*f*1e-6-offset for f in Freqs]
+      
+        Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, U, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+        Detunings, Fluorescence2 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, U, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA2, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+    
+    
+        ScaledFluo1 = np.array([f*SCALE1 + 0.5*OFFSET for f in Fluorescence1])
+        ScaledFluo2 = np.array([f*SCALE2 + 0.5*OFFSET for f in Fluorescence2])
+        
+        if plot:
+            return ScaledFluo1+ScaledFluo2, Detunings
+        else:
+            return ScaledFluo1+ScaledFluo2
+        #return ScaledFluo1
+
+
+
+    def FitEIT_MM_1ion(Freqs, offset1, offset2, DetDoppler, SG, SP, SCALE1, OFFSET, BETA1, TEMP, plot=False):
+    #def FitEIT_MM(freqs, SG, SP, SCALE1, OFFSET, BETA1):
+        #BETA = 1.8
+        # SG = 0.6
+        # SP = 8.1
+        # TEMP = 0.2e-3
+    
+        freqs = [2*f*1e-6-offset1 for f in Freqs]
+      
+        Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+        #Detunings, Fluorescence2 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA2, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+    
+    
+        ScaledFluo1 = np.array([f*SCALE1 + 1*OFFSET for f in Fluorescence1])
+        #ScaledFluo2 = np.array([f*SCALE2 + 0.5*OFFSET for f in Fluorescence2])
+        
+        if plot:
+            return ScaledFluo1, Detunings
+        else:
+            return ScaledFluo1
+
+    
+    do_fit = True
+    if do_fit:
+        try:
+            #popt_multi1_1ion_test, pcov_multi1_1ion_test = curve_fit(FitEIT_MM_1ion, FreqsDR, CountsDR, p0=[448.2, -44.8, 0.5, 6.6, 3.8e4, 1.5e-1, 4.2, 1.4e-3], bounds=((0, -50, 0, 0, 0, 0, 0, 0), (1000, 0, 2, 20, 5e6, 5e4, 10, 20e-3)))
+            popt_multi1_2ion_test, pcov_multi1_2ion_test = curve_fit(FitEIT_MM_2ion, FreqsDR, CountsDR, p0=[448.2, -44.8, 0.5, 6.6, 3.8e4, 1.26e5, 1000, 4.2, 1.3, 1.4e-3, 32e6], bounds=((0, -50, 0, 0, 0, 0, 800, 0,0, 0, 28e6), (1000, 0, 2, 20, 5e6, 5e6, 2000, 10, 10,20e-3,40e6)))
+        except:
+            popt_multi1_2ion_test = [0,0,0,0,0,0,0,0,0,0]
+            pcov_multi1_2ion_test = [0]
+    
+
+        # FittedEITpi_multi1_short_1ion, Detunings_multi1_short_1ion = FitEIT_MM_1ion(FreqsDR, *popt_multi1_1ion_test, plot=True)
+        # freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+        # FittedEITpi_multi1_long_1ion, Detunings_multi1_long_1ion = FitEIT_MM_1ion(freqslong, *popt_multi1_1ion_test, plot=True)
+
+
+        FittedEITpi_multi1_short_2ion, Detunings_multi1_short_2ion = FitEIT_MM_2ion(FreqsDR, *popt_multi1_2ion_test, plot=True)
+        freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+        FittedEITpi_multi1_long_2ion, Detunings_multi1_long_2ion = FitEIT_MM_2ion(freqslong, *popt_multi1_2ion_test, plot=True)
+
+
+        #popt_vecs.append(popt_multi1_2ion)
+        #pcov_vecs.append(pcov_multi1_2ion)
+
+    print(f'Listo {selectedcurve}')
+
+    # plt.figure()
+    # plt.errorbar(Detunings_multi1_short_1ion, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
+    # plt.plot(Detunings_multi1_long_1ion, FittedEITpi_multi1_long_1ion, color='red', linewidth=3, label=f'selcurve: {selectedcurve}')
+    # plt.title('1 ion model')
+    # plt.xlabel('Detuning (MHz)')
+    # plt.ylabel('Counts')
+    # plt.legend(loc='upper left', fontsize=20)
+    # plt.grid()
+
+
+    plt.figure()
+    plt.errorbar(Detunings_multi1_short_2ion, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
+    plt.plot(Detunings_multi1_long_2ion, FittedEITpi_multi1_long_2ion, color='darkolivegreen', linewidth=3, label=f'selcurve: {selectedcurve}')
+    plt.title('2 ion model')
+    plt.xlabel('Detuning (MHz)')
+    plt.ylabel('Counts')
+    plt.legend(loc='upper left', fontsize=20)
+    plt.grid()
+    
+
+# plt.plot(detunings,'o')
+
+
+#%%
+#from EITfit.MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels
+from scipy.optimize import curve_fit
+import time
+
+"""
+MEDICION MULTI 1
+
+Cada bloque ajusta un grupo de mediciones porque sino es un lio
+
+"""
+
+phidoppler, titadoppler = 0, 90
+phirepump, titarepump = 0,  0
+phiprobe = 0
+titaprobe = 90
+
+Temp = 0.5e-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 = 32.5e6
+
+#B = (u/(2*np.pi))/c
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 
+alpha = 0
+
+
+drivefreq = 2*np.pi*22.135*1e6
+
+measurement = 1
+#selectedcurve=10
+selectedcurvevec=[7,8,9,10,11,12,13,14]
+
+popt_vecs = []
+pcov_vecs = []
+
+for selectedcurve in selectedcurvevec:
+    FreqsDR = MultiFreqs[measurement]
+    CountsDR = CountsSplit[measurement][selectedcurve]
+       
+    freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+    CircPr = 1
+    alpha = 0
+    
+    
+    
+    def FitEIT_MM_2ion(Freqs, offset, DetDoppler, SG, SP, SCALE1, SCALE2, OFFSET, BETA1, BETA2, TEMP, U, plot=False):
+    #def FitEIT_MM(freqs, SG, SP, SCALE1, OFFSET, BETA1):
+        #BETA = 1.8
+        # SG = 0.6
+        # SP = 8.1
+        # TEMP = 0.2e-3
+    
+        freqs = [2*f*1e-6-offset for f in Freqs]
+      
+        Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, U, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+        Detunings, Fluorescence2 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, U, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA2, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+    
+    
+        ScaledFluo1 = np.array([f*SCALE1 + 0.5*OFFSET for f in Fluorescence1])
+        ScaledFluo2 = np.array([f*SCALE2 + 0.5*OFFSET for f in Fluorescence2])
+        
+        if plot:
+            return ScaledFluo1+ScaledFluo2, Detunings
+        else:
+            return ScaledFluo1+ScaledFluo2
+        #return ScaledFluo1
+    
+    do_fit = True
+    if do_fit:
+        try:
+            popt_multi1_2ion, pcov_multi1_2ion = curve_fit(FitEIT_MM_2ion, FreqsDR, CountsDR, p0=[448.2, -44.8, 0.5, 6.6, 3.8e4, 1.26e5, 1.5e-1, 4.2, 1.3, 1.4e-3, 32e6], bounds=((0, -50, 0, 0, 0, 0, 0, 0,0, 0,25e6), (1000, 0, 2, 20, 5e6, 5e6, 5e4, 10, 10,20e-3,40e6)))
+        except:
+            popt_multi1_2ion = [0,0,0,0,0,0,0,0,0,0,0]
+            pcov_multi1_2ion = [0]
+    
+        FittedEITpi_multi1_short_2ion, Detunings_multi1_short_2ion = FitEIT_MM_2ion(FreqsDR, *popt_multi1_2ion, plot=True)
+        freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+        FittedEITpi_multi1_long_2ion, Detunings_multi1_long_2ion = FitEIT_MM_2ion(freqslong, *popt_multi1_2ion, plot=True)
+
+        popt_vecs.append(popt_multi1_2ion)
+        pcov_vecs.append(pcov_multi1_2ion)
+
+    print(f'Listo {selectedcurve}')
+
+    plt.figure()
+    plt.errorbar(Detunings_multi1_short_2ion, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
+    plt.plot(Detunings_multi1_long_2ion, FittedEITpi_multi1_long_2ion, color='darkolivegreen', linewidth=3, label=f'selcurve: {selectedcurve}')
+    #plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
+    plt.xlabel('Detuning (MHz)')
+    plt.ylabel('Counts')
+    plt.legend(loc='upper left', fontsize=20)
+    plt.grid()
+
+# plt.plot(detunings,'o')
+
+
+#from EITfit.MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels
+from scipy.optimize import curve_fit
+import time
+
+"""
+MEDICION MULTI 1
+
+otras
+"""
+
+phidoppler, titadoppler = 0, 90
+phirepump, titarepump = 0,  0
+phiprobe = 0
+titaprobe = 90
+
+Temp = 0.5e-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 = 32.5e6
+
+#B = (u/(2*np.pi))/c
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 
+alpha = 0
+
+
+drivefreq = 2*np.pi*22.135*1e6
+
+measurement = 1
+#selectedcurve=10
+selectedcurvevec=[15,16,17,18]
+
+popt_vecs2 = []
+pcov_vecs2 = []
+
+for selectedcurve in selectedcurvevec:
+    FreqsDR = MultiFreqs[measurement]
+    CountsDR = CountsSplit[measurement][selectedcurve]
+       
+    freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+    CircPr = 1
+    alpha = 0
+    
+    
+    
+    def FitEIT_MM_2ion(Freqs, offset, DetDoppler, SG, SP, SCALE1, SCALE2, OFFSET, BETA1, BETA2, TEMP, U, plot=False):
+    #def FitEIT_MM(freqs, SG, SP, SCALE1, OFFSET, BETA1):
+        #BETA = 1.8
+        # SG = 0.6
+        # SP = 8.1
+        # TEMP = 0.2e-3
+    
+        freqs = [2*f*1e-6-offset for f in Freqs]
+      
+        Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, U, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+        Detunings, Fluorescence2 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, U, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA2, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+    
+    
+        ScaledFluo1 = np.array([f*SCALE1 + 0.5*OFFSET for f in Fluorescence1])
+        ScaledFluo2 = np.array([f*SCALE2 + 0.5*OFFSET for f in Fluorescence2])
+        
+        if plot:
+            return ScaledFluo1+ScaledFluo2, Detunings
+        else:
+            return ScaledFluo1+ScaledFluo2
+        #return ScaledFluo1
+    
+    do_fit = True
+    if do_fit:
+        try:
+            popt_multi1_2ion, pcov_multi1_2ion = curve_fit(FitEIT_MM_2ion, FreqsDR, CountsDR, p0=[447.5, -44.0, 0.7, 10, 5.0e4, 6e4, 1.5e-10, 3.7, 1.3, 1.1e-3, 32e6], bounds=((0, -50, 0, 0, 0, 0, 0, 0,0, 0, 25e6), (1000, 0, 2, 20, 5e6, 5e6, 5e4, 10, 10,20e-3, 40e6)))
+        except:
+            popt_multi1_2ion = [0,0,0,0,0,0,0,0,0,0,0]
+            pcov_multi1_2ion = [0]
+    
+        FittedEITpi_multi1_short_2ion, Detunings_multi1_short_2ion = FitEIT_MM_2ion(FreqsDR, *popt_multi1_2ion, plot=True)
+        freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+        FittedEITpi_multi1_long_2ion, Detunings_multi1_long_2ion = FitEIT_MM_2ion(freqslong, *popt_multi1_2ion, plot=True)
+
+        popt_vecs2.append(popt_multi1_2ion)
+        pcov_vecs2.append(pcov_multi1_2ion)
+
+    print(f'Listo {selectedcurve}')
+
+    plt.figure()
+    plt.errorbar(Detunings_multi1_short_2ion, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
+    plt.plot(Detunings_multi1_long_2ion, FittedEITpi_multi1_long_2ion, color='darkolivegreen', linewidth=3, label=f'selcurve: {selectedcurve}')
+    #plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
+    plt.xlabel('Detuning (MHz)')
+    plt.ylabel('Counts')
+    plt.legend(loc='upper left', fontsize=20)
+    plt.grid()
+
+#%%
+
+#from EITfit.MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels
+from scipy.optimize import curve_fit
+import time
+
+"""
+MEDICION MULTI 1
+
+otras
+"""
+
+phidoppler, titadoppler = 0, 90
+phirepump, titarepump = 0,  0
+phiprobe = 0
+titaprobe = 90
+
+Temp = 0.5e-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 = 32.5e6
+
+#B = (u/(2*np.pi))/c
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 
+alpha = 0
+
+
+drivefreq = 2*np.pi*22.135*1e6
+
+measurement = 1
+#selectedcurve=10
+selectedcurvevec=[6]
+
+popt_vecs3 = []
+pcov_vecs3 = []
+
+for selectedcurve in selectedcurvevec:
+    FreqsDR = MultiFreqs[measurement]
+    CountsDR = CountsSplit[measurement][selectedcurve]
+       
+    freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+    CircPr = 1
+    alpha = 0
+    
+    
+    
+    def FitEIT_MM_2ion(Freqs, offset, DetDoppler, SG, SP, SCALE1, SCALE2, OFFSET, BETA1, BETA2, TEMP, U, plot=False):
+    #def FitEIT_MM(freqs, SG, SP, SCALE1, OFFSET, BETA1):
+        #BETA = 1.8
+        # SG = 0.6
+        # SP = 8.1
+        # TEMP = 0.2e-3
+    
+        freqs = [2*f*1e-6-offset for f in Freqs]
+      
+        Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, U, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+        Detunings, Fluorescence2 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, U, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA2, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+    
+    
+        ScaledFluo1 = np.array([f*SCALE1 + 0.5*OFFSET for f in Fluorescence1])
+        ScaledFluo2 = np.array([f*SCALE2 + 0.5*OFFSET for f in Fluorescence2])
+        
+        if plot:
+            return ScaledFluo1+ScaledFluo2, Detunings
+        else:
+            return ScaledFluo1+ScaledFluo2
+        #return ScaledFluo1
+    
+    do_fit = True
+    if do_fit:
+        try:
+            popt_multi1_2ion, pcov_multi1_2ion = curve_fit(FitEIT_MM_2ion, FreqsDR, CountsDR, p0=[448.2, -45.8, 0.6, 10, 7.3e4, 2.9e4, 1.3e3, 3.7, 1.1, 3.4e-3,32e6], bounds=((0, -50, 0, 0, 0, 0, 0, 0,0, 0,25e6), (1000, 0, 2, 20, 5e6, 5e6, 5e4, 10, 10,20e-3,40e6)))
+        except:
+            popt_multi1_2ion = [0,0,0,0,0,0,0,0,0,0,0]
+            pcov_multi1_2ion = [0]
+    
+        FittedEITpi_multi1_short_2ion, Detunings_multi1_short_2ion = FitEIT_MM_2ion(FreqsDR, *popt_multi1_2ion, plot=True)
+        freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+        FittedEITpi_multi1_long_2ion, Detunings_multi1_long_2ion = FitEIT_MM_2ion(freqslong, *popt_multi1_2ion, plot=True)
+
+        popt_vecs3.append(popt_multi1_2ion)
+        pcov_vecs3.append(pcov_multi1_2ion)
+
+    print(f'Listo {selectedcurve}')
+
+    plt.figure()
+    plt.errorbar(Detunings_multi1_short_2ion, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
+    plt.plot(Detunings_multi1_long_2ion, FittedEITpi_multi1_long_2ion, color='darkolivegreen', linewidth=3, label=f'selcurve: {selectedcurve}')
+    #plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
+    plt.xlabel('Detuning (MHz)')
+    plt.ylabel('Counts')
+    plt.legend(loc='upper left', fontsize=20)
+    plt.grid()
+
+#%%
+betas1s = []
+betas2s = []
+temps = []
+detunings = []
+
+# for i in range(len(popt_vecs3)):
+#     betas1s.append(popt_vecs3[i][7])
+#     betas2s.append(popt_vecs3[i][8])
+#     temps.append(popt_vecs3[i][9])
+#     detunings.append(popt_vecs3[i][6])
+
+for i in range(len(popt_vecs)):
+    betas1s.append(popt_vecs[i][7])
+    betas2s.append(popt_vecs[i][8])
+    temps.append(popt_vecs[i][9])
+    detunings.append(popt_vecs[i][])
+
+for i in range(len(popt_vecs2)):
+    betas1s.append(popt_vecs2[i][7])
+    betas2s.append(popt_vecs2[i][8])
+    temps.append(popt_vecs2[i][9])
+    detunings.append(popt_vecs2[i][5])
+
+# plt.figure()
+# plt.plot(betas1s,'o')
+# plt.plot(betas2s,'o')
+
+plt.figure()
+plt.plot(detunings,'o')
+
+
+

analisis/plots/20231218_CPT_DosLaseres_Reflotoajustes/lolo_CPT_plotter_20231218.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+??
+
+
+@author: lolo
+"""
+
+
+
+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
+
+"""
+CPT con tres laseres pero lso dos IR son el mismo entonces las DD son mas finas
+"""
+
+#C:\Users\Usuario\Documents\artiq\artiq_experiments\analisis\plots\20211223_CPT_DosLaseres_v07_ChristmasSpecial\Data
+
+folder = '../20231218_CPT_DosLaseres_Reflotoajustes/Data'
+
+ALL_FILES = f"""
+{folder}/000016420-IR_Scan_withcal_optimized
+""" 
+
+
+def SeeKeys(files):
+    for i, fname in enumerate(files.split()):
+        data = h5py.File(fname+'.h5', 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
+        print(fname)
+        print(list(data['datasets'].keys()))
+
+print(SeeKeys(ALL_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(ALL_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']))
+
+
+
+#%%
+
+#Barriendo angulo del IR con tisa apagado
+
+jvec = [0]
+
+jselected = jvec
+
+plt.figure()
+i = 0
+for j in jvec:
+    if j in jselected:
+        plt.errorbar([2*f*1e-6 for f in Freqs[j]], Counts[j], yerr=np.sqrt(Counts[j]), fmt='o', capsize=2, markersize=2)
+    #plt.plot([2*f*1e-6 for f in Freqs[j]], Counts[j], 'o-', label=f'Amp Tisa: {AmpTisa[i]}', mb  arkersize=3)
+    i = i + 1
+plt.xlabel('Frecuencia (MHz)')
+plt.ylabel('counts')
+plt.grid()
+plt.legend()
+
+#%%
+from scipy.optimize import curve_fit
+import time
+
+cwd = os.getcwd()
+os.chdir('../20231123_CPTconmicromocion3')
+from Data.EITfit.lolo_modelo_full_3niveles import GenerateNoisyCPT_fit
+os.chdir(cwd)
+
+# from Data.EITfit.lolo_modelo_full_8niveles import PerformExperiment_8levels_MM
+
+
+phidoppler, titadoppler = 0, 90
+phirepump, titarepump = 0,  90
+phiprobe = 0
+titaprobe = 0.1
+
+Temp = 0.5e-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 = 32.5e6
+
+#B = (u/(2*np.pi))/c
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 
+alpha = 0
+
+
+drivefreq = 2*np.pi*22.135*1e6
+noiseamplitude = 0
+
+selectedcurve=0
+
+FreqsDR = Freqs[selectedcurve]
+CountsDR = Counts[selectedcurve]
+   
+freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+CircPr = 1
+alpha = 0
+
+def FitEIT_MM_1ion(Freqs, offset, DetDoppler, DetRepump, SG, SP, SR, SCALE1, OFFSET, TEMP, U, plot=False):
+#def FitEIT_MM(freqs, SG, SP, SCALE1, OFFSET, BETA1):
+    # BETA1 = 0
+    # SG = 0.6
+    # SP = 8.1
+    # TEMP = 0.2e-3
+    # U = 32.5e6
+    freqs = [2*f*1e-6-offset for f in Freqs]
+  
+    #Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+    Detunings, Fluorescence1 = 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)
+
+    ScaledFluo1 = np.array([f*SCALE1 + OFFSET for f in Fluorescence1])
+    if plot:
+        return ScaledFluo1, Detunings
+    else:
+        return ScaledFluo1
+    #return ScaledFluo1
+
+do_fit = True
+if do_fit:
+    popt_1, pcov_1 = curve_fit(FitEIT_MM_1ion, FreqsDR, CountsDR, p0=[430, -25, 12, 0.9, 6.2, 3, 3e4, 2e3, 0.5e-3, 32e6], bounds=((0, -100, -20, 0, 0, 0, 0, 0, 0,20e6), (1000, 0, 50, 2, 20, 20, 5e6, 5e4, 15e-3,40e6)))
+
+    FittedEITpi_1_short, Detunings_1_short = FitEIT_MM_1ion(FreqsDR, *popt_1, plot=True)
+    freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+    FittedEITpi_1_long, Detunings_1_long = FitEIT_MM_1ion(freqslong, *popt_1, plot=True)
+
+#%%
+
+plt.figure()
+plt.errorbar(Detunings_1_short, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='red', alpha=0.5, capsize=2, markersize=2)
+plt.plot(Detunings_1_long, FittedEITpi_1_long, color='darkolivegreen', linewidth=3, label='med 1')
+#plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
+plt.xlabel('Detuning (MHz)')
+plt.ylabel('Counts')
+#plt.xlim(-20,0)
+plt.legend(loc='upper left', fontsize=20)
+plt.grid()
+
+#%%
+# u = 32.5e6
+
+# B = (u/(2*np.pi))/c
+
+# correccion = 8 #con 8 fitea bien
+
+# offsetxpi = 440+1+correccion
+# DetDoppler = -5.0-correccion
+
+# FreqsDRpi_3 = [2*f*1e-6-offsetxpi+14 for f in Freqs_B[5]]
+# CountsDRpi_3 = Counts_B[5]
+
+# freqslongpi_3 = np.arange(min(FreqsDRpi_3), max(FreqsDRpi_3)+FreqsDRpi_3[1]-FreqsDRpi_3[0], 0.1*(FreqsDRpi_3[1]-FreqsDRpi_3[0]))
+
+# #[1.71811842e+04 3.34325038e-17]
+
+# def FitEITpi(freqs, SG, SP):
+#     temp = 2e-3
+#     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)
+#     FinalFluo = [f*6.554e4 + 1.863e3 for f in MeasuredFluo]
+#     return FinalFluo
+
+# popt_tisaoff, pcov_tisaoff = curve_fit(FitEITpi, FreqsDRpi_3, CountsDRpi_3, p0=[0.5, 4.5], bounds=((0, 0), (2, 10)))
+
+# print(popt_tisaoff)
+
+# Sat_3 = popt_tisaoff[0]
+# Det_3 = popt_tisaoff[1]
+
+# FittedEITpi_3 = FitEITpi(freqslongpi_3, *popt_tisaoff)
+
+# plt.figure()
+# plt.errorbar(FreqsDRpi_3, CountsDRpi_3, yerr=2*np.sqrt(CountsDRpi_3), fmt='o', capsize=2, markersize=2)
+# plt.plot(freqslongpi_3, FittedEITpi_3)
+# #plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
+
+# FreqsCalibradas_B = FreqsDRpi_3
+
+
+
+
+
+
+
+
+
+

analisis/plots/20231226_CPTconmicromocion4/lolo_CPT_plotter_20231226.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+??
+
+
+@author: lolo
+"""
+
+
+
+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
+
+"""
+Primero tengo mediciones de espectros cpt de un ion variando la tension dc_A
+"""
+
+#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/20231226_CPTconmicromocion4/Data/')
+
+# CPT_FILES = """000016531-IR_Scan_withcal_optimized
+# 000016532-IR_Scan_withcal_optimized
+# """
+
+
+folder = '../20231226_CPTconmicromocion4/Data'
+
+CPT_FILES = f"""
+{folder}/000016531-IR_Scan_withcal_optimized
+{folder}/000016532-IR_Scan_withcal_optimized
+""" 
+
+
+def SeeKeys(files):
+    for i, fname in enumerate(files.split()):
+        data = h5py.File(fname+'.h5', 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
+        print(fname)
+        print(list(data['datasets'].keys()))
+
+print(SeeKeys(CPT_FILES))
+
+#carpeta pc nico labo escritorio:
+#C:\Users\Usuario\Documents\artiq\artiq_experiments\analisis\plots\20211101_CPT_DosLaseres_v03\Data
+
+Counts = []
+Freqs = []
+
+AmpTisa = []
+UVCPTAmp = []
+No_measures = []
+Voltages = []
+
+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']['data_array']))
+    #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']))
+    Voltages.append(np.array(data['datasets']['scanning_voltages']))
+
+def Split(array,n):
+    length=len(array)/n
+    splitlist = []
+    jj = 0
+    while jj<length:
+        partial = []
+        ii = 0
+        while ii < n:
+            partial.append(array[jj*n+ii])
+            ii = ii + 1
+        splitlist.append(partial)
+        jj = jj + 1
+    return splitlist
+
+
+CountsSplit = []
+CountsSplit.append(Split(Counts[0],len(Freqs[0])))
+CountsSplit.append(Split(Counts[1],len(Freqs[1])))
+
+#%%
+
+"""
+Ploteo la cpt de referencia / plotting the reference CPT
+"""
+
+medic = 1 #puede ser 0 o 1
+
+jvec = [28] # de la 1 a la 9 vale la pena, despues no
+
+drs = [390.5, 399.5, 406, 413.5]
+
+drive=22.1
+
+Frequencies = Freqs[medic]
+
+plt.figure()
+i = 0
+for j in jvec:
+    plt.errorbar([2*f*1e-6 for f in Frequencies], CountsSplit[medic][j], yerr=np.sqrt(CountsSplit[0][j]), fmt='o', capsize=2, markersize=2)
+    i = i + 1
+plt.xlabel('Frecuencia (MHz)')
+plt.ylabel('counts')
+plt.grid()
+#for dr in drs:
+#    plt.axvline(dr)
+    #plt.axvline(dr+drive)
+plt.legend()
+
+#%%
+"""
+SUPERAJUSTE: TESTEO FITEOS
+"""
+
+
+cwd = os.getcwd()
+os.chdir('../20231123_CPTconmicromocion3')
+from Data.EITfit.lolo_modelo_full_8niveles import PerformExperiment_8levels_MM
+os.chdir(cwd)
+
+
+
+from scipy.optimize import curve_fit
+import time
+
+"""
+SUPER AJUSTE (SA)
+"""
+
+phidoppler, titadoppler = 0, 90
+phirepump, titarepump = 0,  0
+phiprobe = 0
+titaprobe = 90
+
+Temp = 0.5e-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 = 32.5e6
+
+#B = (u/(2*np.pi))/c
+
+# correccion = 13
+
+#DetDoppler = -11.5-correccion
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6
+alpha = 0
+
+
+drivefreq = 2*np.pi*22.135*1e6
+
+measurement = 1
+SelectedCurveVec = [30]
+
+
+for selectedcurve in SelectedCurveVec:
+
+    FreqsDR = Freqs[measurement]
+    CountsDR = CountsSplit[measurement][selectedcurve]
+
+    freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+    CircPr = 1
+    alpha = 0
+
+
+    def FitEIT_MM_single(Freqs, offset, DetDoppler, SG, SP, SCALE1, OFFSET, BETA1, TEMP, U, DL, PL, plot=False):
+    #def FitEIT_MM(freqs, SG, SP, SCALE1, OFFSET, BETA1):
+        #BETA = 1.8
+        # SG = 0.6
+        # SP = 8.1
+        # TEMP = 0.2e-3
+
+        freqs = [2*f*1e-6-offset for f in Freqs]
+
+        #Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, U, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+        Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, U, DL, PL, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+
+        ScaledFluo1 = np.array([f*SCALE1 + OFFSET for f in Fluorescence1])
+        if plot:
+            return ScaledFluo1, Detunings
+        else:
+            return ScaledFluo1
+        #return ScaledFluo1
+
+    do_fit = True
+    if do_fit:
+        t1 = time.time()
+        print(f'Beginning {selectedcurve}')
+        popt_3_SA, pcov_3_SA = curve_fit(FitEIT_MM_single, FreqsDR, CountsDR, p0=[446, -36, 0.67, 8.5, 6e4, 250, 1, 0.3e-3, 32e6, 0.1, 0.15], bounds=((0, -50, 0, 0, 0, 0, 0, 0, 31e6, 0, 0), (1000, -20, 2, 20, 5e5, 5e4, 10, 100e-3, 34e6, 10, 10)))
+        print(f'Ended {selectedcurve}, time elapsed: {round(time.time()-t1)} s')
+
+        FittedEITpi_3_SA_short, Detunings_3_SA_short = FitEIT_MM_single(FreqsDR, *popt_3_SA, plot=True)
+        freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+        FittedEITpi_3_SA_long, Detunings_3_SA_long = FitEIT_MM_single(freqslong, *popt_3_SA, plot=True)
+
+
+    plt.figure()
+    plt.errorbar(Detunings_3_SA_short, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
+    plt.plot(Detunings_3_SA_long, FittedEITpi_3_SA_long, color='darkolivegreen', linewidth=3, label=f'med {selectedcurve}')
+    #plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
+    plt.xlabel('Detuning (MHz)')
+    plt.ylabel('Counts')
+    plt.legend(loc='upper left', fontsize=20)
+    plt.grid()
+
+    print(popt_3_SA)
+
+    print(f'Beta: {popt_3_SA[6]} +- {np.sqrt(pcov_3_SA[6,6])}')
+
+#%%
+"""
+AHORA INTENTO SUPER AJUSTES O SEA CON OFFSETXPI Y DETDOPPLER INCLUIDOS
+"""
+
+
+from scipy.optimize import curve_fit
+import time
+
+"""
+SUPER AJUSTE (SA)
+"""
+
+phidoppler, titadoppler = 0, 90
+phirepump, titarepump = 0,  0
+phiprobe = 0
+titaprobe = 90
+
+Temp = 0.5e-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 = 32.5e6
+
+#B = (u/(2*np.pi))/c
+
+correccion = 13
+
+#DetDoppler = -11.5-correccion
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6
+alpha = 0
+
+
+drivefreq = 2*np.pi*22.135*1e6
+
+measurement = 1
+SelectedCurveVec = np.arange(0,31,1)
+#SelectedCurveVec = [17]
+
+popt_SA_vec = []
+pcov_SA_vec = []
+Detuningsshort_vec = []
+Counts_vec = []
+Detuningslong_vec = []
+FittedCounts_vec = []
+
+Betas_vec = []
+ErrorBetas_vec = []
+Temp_vec = []
+ErrorTemp_vec = []
+
+DetuningsUV_vec = []
+ErrorDetuningsUV_vec = []
+
+for selectedcurve in SelectedCurveVec:
+
+    FreqsDR = Freqs[measurement]
+    CountsDR = CountsSplit[measurement][selectedcurve]
+
+    freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+    CircPr = 1
+    alpha = 0
+
+
+    def FitEIT_MM_single(Freqs, offset, DetDoppler, SG, SP, SCALE1, OFFSET, BETA1, TEMP, U, plot=False):
+    #def FitEIT_MM(freqs, SG, SP, SCALE1, OFFSET, BETA1):
+        #BETA = 1.8
+        # SG = 0.6
+        # SP = 8.1
+        # TEMP = 0.2e-3
+
+        freqs = [2*f*1e-6-offset for f in Freqs]
+
+        Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, U, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+
+        ScaledFluo1 = np.array([f*SCALE1 + OFFSET for f in Fluorescence1])
+        if plot:
+            return ScaledFluo1, Detunings
+        else:
+            return ScaledFluo1
+        #return ScaledFluo1
+
+    do_fit = True
+    if do_fit:
+        t1 = time.time()
+        print(f'Beginning {selectedcurve}')
+        popt_3_SA, pcov_3_SA = curve_fit(FitEIT_MM_single, FreqsDR, CountsDR, p0=[443, -32, 0.67, 8, 4e4, 350, 1, 0.5e-3, 32e6], bounds=((0, -50, 0, 0, 0, 0, 0, 0, 31e6), (1000, -20, 2, 20, 5e5, 5e4, 10, 100e-3, 34e6)))
+        print(f'Ended {selectedcurve}, time elapsed: {round(time.time()-t1)} s')
+
+        popt_SA_vec.append(popt_3_SA)
+        pcov_SA_vec.append(pcov_3_SA)
+
+        FittedEITpi_3_SA_short, Detunings_3_SA_short = FitEIT_MM_single(FreqsDR, *popt_3_SA, plot=True)
+        freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+        FittedEITpi_3_SA_long, Detunings_3_SA_long = FitEIT_MM_single(freqslong, *popt_3_SA, plot=True)
+
+    DetuningsUV_vec.append(popt_3_SA[1])
+    ErrorDetuningsUV_vec.append(np.sqrt(pcov_3_SA[1,1]))
+
+    Betas_vec.append(popt_3_SA[6])
+    ErrorBetas_vec.append(np.sqrt(pcov_3_SA[6,6]))
+    Temp_vec.append(popt_3_SA[7])
+    ErrorTemp_vec.append(np.sqrt(pcov_3_SA[7,7]))
+
+    Detuningsshort_vec.append(Detunings_3_SA_short)
+    Counts_vec.append(CountsDR)
+    Detuningslong_vec.append(Detunings_3_SA_long)
+    FittedCounts_vec.append(FittedEITpi_3_SA_long)
+
+
+    plt.figure()
+    plt.errorbar(Detunings_3_SA_short, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
+    plt.plot(Detunings_3_SA_long, FittedEITpi_3_SA_long, color='darkolivegreen', linewidth=3, label=f'med {selectedcurve}')
+    #plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
+    plt.xlabel('Detuning (MHz)')
+    plt.ylabel('Counts')
+    plt.legend(loc='upper left', fontsize=20)
+    plt.grid()
+
+    print(f'listo med {selectedcurve}')
+    print(popt_3_SA)
+
+
+#%%
+"""
+Grafico distintas variables que salieron del SUper ajuste
+"""
+
+import seaborn as sns
+paleta = sns.color_palette("rocket")
+
+medfin = 19
+
+
+voltages_dcA = Voltages[0][0:medfin]
+
+def lineal(x,a,b):
+    return a*x+b
+
+def hiperbola(x,a,b,c,x0):
+    return a*np.sqrt(((x-x0)**2+c**2))+b
+
+hiperbola_or_linear = True
+
+if hiperbola_or_linear:
+    #popthip,pcovhip = curve_fit(hiperbola,voltages_dcA[0:20],Betas_vec[0:20],p0=(100,0.1,1,-0.15))
+    
+    #xhip = np.linspace(-0.23,0.055,200)
+    
+    popthip,pcovhip = curve_fit(hiperbola,voltages_dcA,Betas_vec[:medfin],p0=(100,0.1,1,-0.15))
+
+    xhip = np.linspace(-0.16,0.05,200)
+
+    plt.figure()
+    plt.errorbar(voltages_dcA,Betas_vec[0:medfin],yerr=ErrorBetas_vec[:medfin],fmt='o',capsize=5,markersize=5,color=paleta[1])
+    plt.plot(xhip,hiperbola(xhip,*popthip))
+    plt.xlabel('Endcap voltage (V)')
+    plt.ylabel('Modulation factor')
+    plt.ylim(-1,4)
+    plt.grid()
+
+else:
+    poptini,pcovini = curve_fit(lineal,voltages_dcA[0:3],Betas_vec[0:3])
+    poptfin,pcovfin = curve_fit(lineal,voltages_dcA[4:],Betas_vec[4:])
+
+    minimum_voltage = -(poptini[1]-poptfin[1])/(poptini[0]-poptfin[0]) #voltaje donde se intersectan las rectas, es decir, donde deberia estar el minimo de micromocion
+    minimum_modulationfactor = lineal(minimum_voltage,*poptini) #es lo mismo si pongo *poptfin
+
+    xini = np.linspace(-0.23,-0.13,100)
+    xfin = np.linspace(-0.15,0.005,100)
+
+    plt.figure()
+    plt.errorbar(voltages_dcA,Betas_vec,yerr=ErrorBetas_vec,fmt='o',capsize=5,markersize=5,color=paleta[1])
+    plt.plot(xini,lineal(xini,*poptini))
+    plt.plot(xfin,lineal(xfin,*poptfin))
+    plt.y
+    plt.axvline(minimum_voltage,linestyle='dashed',color='grey')
+    plt.xlabel('Endcap voltage (V)')
+    plt.ylabel('Modulation factor')
+    plt.grid()
+
+
+print([t*1e3 for t in Temp_vec])
+
+plt.figure()
+plt.errorbar(voltages_dcA,[t*1e3 for t in Temp_vec[:medfin]],yerr=[t*1e3 for t in ErrorTemp_vec[:medfin]],fmt='o',capsize=5,markersize=5,color=paleta[3])
+#plt.axvline(minimum_voltage,linestyle='dashed',color='grey')
+plt.axhline(0.538)
+#plt.yscale('log')
+plt.xlabel('Endcap voltage (V)')
+plt.ylabel('Temperature (mK)')
+plt.grid()
+#plt.ylim(0,2)
+
+#%%
+"""
+Ahora hago un ajuste con una hiperbola porque tiene mas sentido, por el hecho
+de que en el punto optimo el ion no esta en el centro de la trampa
+sino que esta a una distancia d
+"""
+def hiperbola(x,a,b,c,x0):
+    return a*np.sqrt(((x-x0)**2+c**2))+b
+
+medfin = 11
+voltages_dcA = Voltages[0][1:medfin]
+
+
+popthip,pcovhip = curve_fit(hiperbola,voltages_dcA[:10],Betas_vec[:10],p0=(100,0.1,1,-0.15))
+
+xhip = np.linspace(-0.23,0.055,200)
+
+plt.figure()
+plt.errorbar(voltages_dcA,Betas_vec[:medfin-1],yerr=ErrorBetas_vec[:medfin-1],fmt='o',capsize=5,markersize=5,color=paleta[1])
+plt.plot(xhip,hiperbola(xhip,*popthip))
+plt.xlabel('Endcap voltage (V)')
+plt.ylabel('Modulation factor')
+#plt.yscale('log')
+plt.grid()
+
+
+
+
+#%%
+from scipy.special import jv
+
+
+def expo(x,tau,A,B):
+    return A*np.exp(x/tau)+B
+
+def cuadratica(x,a,c):
+    return a*(x**2)+c
+
+def InverseMicromotionSpectra(beta, A, det, x0, gamma, B):
+    ftrap=22.1
+    #gamma=30
+    P = ((jv(0, beta)**2)/((((det-x0)**2)+(0.5*gamma)**2)**2))*(-2*(det-x0))
+    i = 1
+    #print(P)
+    while i <= 5:
+        P = P + (-2*(det-x0))*((jv(i, beta))**2)/(((((det-x0)+i*ftrap)**2)+(0.5*gamma)**2)**2) + (-2*(det-x0))*(((jv(-i, beta))**2)/((((det-x0)-i*ftrap)**2)+(0.5*gamma)**2)**2)
+        i = i + 1
+        #print(P)
+    #return 1/(A*P+B)
+    return 1/(A*P+B)
+
+
+def MicromotionSpectra(beta,det, gamma):
+    ftrap=22.1
+    #gamma=23
+    P = (jv(0, beta)**2)/(((det)**2)+(0.5*gamma)**2)
+    i = 1
+    #print(P)
+    while i <= 5:
+        P = P + ((jv(i, beta))**2)/((((det)+i*ftrap)**2)+(0.5*gamma)**2) + ((jv(-i, beta))**2)/((((det)-i*ftrap)**2)+(0.5*gamma)**2) 
+        i = i + 1
+        #print(P)
+    return P
+
+
+def polynomial(x,a,b,c,d,e):
+    b=0
+    d=0
+    return a+b*x+c*x*x+d*x*x*x+e*x*x*x*x
+
+def InverseDerivMicromotionSpectra(beta, det, gamma):
+    ftrap=22.1
+    #gamma=23
+    #det = -gamma/2
+    P = ((jv(0, beta)**2)/((((det)**2)+(0.5*gamma)**2)**2))*(-2*(det))
+    i = 1
+    #print(P)
+    while i <= 5:
+        P = P + (-2*(det))*((jv(i, beta))**2)/(((((det)+i*ftrap)**2)+(0.5*gamma)**2)**2) + (-2*(det))*(((jv(-i, beta))**2)/((((det)-i*ftrap)**2)+(0.5*gamma)**2)**2)
+        i = i + 1
+        #print(P)
+    return 1/P
+
+
+def FinalTemp(beta,det, C,D):
+    gamma = 21
+    #det=-11
+    #D=-0.8
+    #C = 1.68656122e-03  
+    #D = 6.64227010e-02
+    #C=0
+    
+    #print(MicromotionSpectra(beta,det,gamma))
+    return (C*MicromotionSpectra(beta,det,gamma)+D*beta**2)*InverseDerivMicromotionSpectra(beta, det, gamma)
+    #return (C*MicromotionSpectra(beta,det,gamma))*InverseDerivMicromotionSpectra(beta, det, gamma)
+
+"""
+Temperatura vs beta con un ajuste exponencial
+"""
+
+medfin = 10
+
+popt_exp, pcov_exp = curve_fit(expo,Betas_vec[:medfin],[t*1e3 for t in Temp_vec[:medfin]])
+#popt_quad, pcov_quad = curve_fit(cuadratica,Betas_vec[:11],[t*1e3 for t in Temp_vec[:11]],p0=(1,10))
+#popt_rho22, pcov_rho22 = curve_fit(InverseMicromotionSpectra,Betas_vec,[t*1e3 for t in Temp_vec],p0=(10,10,-10,1,20)) #esto ajusta muy bien
+#popt_rho22, pcov_rho22 = curve_fit(InverseMicromotionSpectra,Betas_vec, [t*1e3 for t in Temp_vec],p0=(-10,-10,10,1,20)) #esto ajusta muy bien
+#popt_rho22_raw, pcov_rho22_raw = curve_fit(InverseMicromotionSpectra_raw,Betas_vec[:7], [t*1e3 for t in Temp_vec[:7]],p0=(-0.1, -10, 1)) #esto ajusta muy bien
+popt_rho22_balance, pcov_rho22_balance = curve_fit(FinalTemp,Betas_vec[:medfin], [t*1e3 for t in Temp_vec[:medfin]],p0=(-10, 10,1)) #esto ajusta muy bien
+popt_rho22_poly, pcov_rho22_poly = curve_fit(polynomial,Betas_vec[:medfin], [t*1e3 for t in Temp_vec[:medfin]],p0=(1,2,3,4,10)) #esto ajusta muy bien
+
+
+print(popt_rho22_balance)
+
+betaslong = np.arange(0,2.8,0.01)
+
+print(f'Min temp predicted: {FinalTemp(betaslong,*popt_rho22_balance)[0]}')
+
+print(f'Detuning: {popt_rho22_balance[0]} MHz')
+print(f'rho22 coeff: {popt_rho22_balance[1]}')
+print(f'betasquared coeff: {popt_rho22_balance[2]}')
+
+print(f'cociente de los coeff: {popt_rho22_balance[2]/popt_rho22_balance[1]}')
+
+
+print(f'params: {popt_rho22_balance}')
+print(f'errores: {np.sqrt(np.diag(pcov_rho22_balance))}')
+
+k_plot = medfin
+
+plt.figure()
+plt.errorbar(Betas_vec[:k_plot],[t*1e3 for t in Temp_vec[:k_plot]],xerr=ErrorBetas_vec[:k_plot], yerr=[t*1e3 for t in ErrorTemp_vec[:k_plot]],fmt='o',capsize=5,markersize=5,color=paleta[3])
+#plt.plot(betaslong,polynomial(betaslong,*popt_rho22_poly),label='Ajuste exponencial')
+#plt.plot(betaslong,FinalTemp(betaslong,popt_rho22_balance[0],popt_rho22_balance[1],popt_rho22_balance[2]*1),label='Ajuste con espectro modulado')
+plt.xlim(-0,3)
+plt.ylim(0,1)
+
+#plt.axvline(minimum_voltage,linestyle='dashed',color='grey')
+#plt.axhline(0.538)
+plt.xlabel('Modulation factor')
+plt.ylabel('Temperature (mK)')
+plt.legend()
+plt.grid()
+
+#%%
+
+hbar=1.05e-34
+gammita = 22.1e6
+kx = 2*np.pi/(0.397e-6)
+kb = 1.38e-23
+masita = 6.6e-26
+
+coeff = gammita*hbar*hbar*kx*kx/(6*kb*masita)
+print(coeff)
+
+rfheatrate = coeff*popt_rho22_balance[2]/popt_rho22_balance[1]
+
+print(f'heating rate due to rf heating: {rfheatrate*1e3} mK/s')
+
+
+
+#%%
+"""
+Esto no es del super ajuste sino de los ajustes anteriores en donde DetDoppler y offset son puestos a mano
+Aca grafico los betas con su error en funcion de la tension variada.
+Ademas, hago ajuste lineal para primeros y ultimos puntos, ya que espero que
+si la tension hace que la posicion del ion varie linealmente, el beta varia proporcional a dicha posicion.
+"""
+
+import seaborn as sns
+
+def lineal(x,a,b):
+    return a*x+b
+
+paleta = sns.color_palette("rocket")
+
+betavector = [beta1,beta2,beta3,beta4,beta5,beta6,beta7,beta8,beta9]
+errorbetavector = [errorbeta1,errorbeta2,errorbeta3,errorbeta4,errorbeta5,errorbeta6,errorbeta7,errorbeta8,errorbeta9]
+voltages_dcA = Voltages[0][1:10]
+
+poptini,pcovini = curve_fit(lineal,voltages_dcA[0:3],betavector[0:3])
+poptfin,pcovfin = curve_fit(lineal,voltages_dcA[4:],betavector[4:])
+
+minimum_voltage = -(poptini[1]-poptfin[1])/(poptini[0]-poptfin[0]) #voltaje donde se intersectan las rectas, es decir, donde deberia estar el minimo de micromocion
+minimum_modulationfactor = lineal(minimum_voltage,*poptini) #es lo mismo si pongo *poptfin
+
+xini = np.linspace(-0.23,-0.13,100)
+xfin = np.linspace(-0.15,0.005,100)
+
+plt.figure()
+plt.errorbar(voltages_dcA,betavector,yerr=errorbetavector,fmt='o',capsize=5,markersize=5,color=paleta[1])
+plt.plot(xini,lineal(xini,*poptini))
+plt.plot(xfin,lineal(xfin,*poptfin))
+plt.axvline(minimum_voltage,linestyle='dashed',color='grey')
+plt.xlabel('Endcap voltage (V)')
+plt.ylabel('Modulation factor')
+plt.grid()
+
+
+#%%
+"""
+Aca veo la temperatura del ion en funcion del voltaje del endcap, ya que
+al cambiar la cantidad de micromocion, cambia la calidad del enfriado
+"""
+tempvector = np.array([temp1,temp2,temp3,temp4,temp5,temp6,temp7,temp8,temp9])*1e3
+errortempvector = np.array([errortemp1,errortemp2,errortemp3,errortemp4,errortemp5,errortemp6,errortemp7,errortemp8,errortemp9])*1e3
+
+voltages_dcA = Voltages[0][1:10]
+
+plt.figure()
+plt.errorbar(voltages_dcA,tempvector,yerr=errortempvector,fmt='o',capsize=5,markersize=5,color=paleta[3])
+plt.axvline(minimum_voltage,linestyle='dashed',color='grey')
+plt.xlabel('Endcap voltage (V)')
+plt.ylabel('Temperature (mK)')
+plt.grid()
+plt.ylim(0,2)
+
+
+#%%
+"""
+Por las dudas, temperatura en funcion de beta
+"""
+plt.figure()
+plt.errorbar(betavector,tempvector,yerr=errortempvector,xerr=errorbetavector,fmt='o',capsize=5,markersize=5)
+plt.xlabel('Modulation factor')
+plt.ylabel('Temperature (mK)')
+plt.grid()
+
+
+
+#%%
+"""
+Si quiero ver algun parametro del ajuste puntual. el orden es: 0:SG, 1:SP, 2:SCALE1, 3:OFFSET
+"""
+ki=2
+plt.errorbar(np.arange(0,9,1),[popt_1[ki],popt_2[ki],popt_3[ki],popt_4[ki],popt_5[ki],popt_6[ki],popt_7[ki],popt_8[ki],popt_9[ki]],yerr=[np.sqrt(pcov_1[ki,ki]),np.sqrt(pcov_2[ki,ki]),np.sqrt(pcov_3[ki,ki]),np.sqrt(pcov_4[ki,ki]),np.sqrt(pcov_5[ki,ki]),np.sqrt(pcov_6[ki,ki]),np.sqrt(pcov_7[ki,ki]),np.sqrt(pcov_8[ki,ki]),np.sqrt(pcov_9[ki,ki])], fmt='o',capsize=3,markersize=3)
+
+#%%
+"""
+AHORA VAMOS A MEDICIONES CON MAS DE UN ION!!!
+"""
+
+"""
+Ploteo la cpt de referencia / plotting the reference CPT
+
+1: 2 iones, -100 mV dcA
+2: 2 iones, -150 mV dcA
+3: 2 iones, -50 mV dcA
+4: 2 iones, 5 voltajes (el ion se va en la 4ta medicion y en la 5ta ni esta)
+5, 6 y 7: 3 iones en donde el scaneo esta centrado en distintos puntos
+
+"""
+
+jvec = [3] # desde la 1, pero la 4 no porque es un merge de curvitas
+
+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()
+#for dr in drs:
+#    plt.axvline(dr)
+    #plt.axvline(dr+drive)
+plt.legend()
+
+
+#%%
+"""
+Mergeo la 5, 6 y 7
+"""
+
+Freqs5 = [2*f*1e-6 for f in Freqs[5]]
+Freqs6 = [2*f*1e-6 for f in Freqs[6]]
+Freqs7 = [2*f*1e-6 for f in Freqs[7]]
+
+Counts5 = Counts[5]
+Counts6 = Counts[6]
+Counts7 = Counts[7]
+
+i_1_ini = 0
+i_1 = 36
+i_2_ini = 0
+i_2 = 24
+
+f_1 = 18
+f_2 = 30
+
+
+scale_1 = 0.92
+scale_2 = 0.98
+
+#Merged_freqs_test = [f-f_2 for f in Freqs6[i_2_ini:i_2]]+[f-f_1 for f in Freqs5[i_1_ini:i_1]]+Freqs7
+
+#plt.plot(Merged_freqs_test,'o')
+
+
+
+Merged_freqs = [f-f_2 for f in Freqs6[0:i_2]]+[f-f_1 for f in Freqs5[0:i_1]]+Freqs7
+Merged_counts = [scale_2*c for c in Counts6[0:i_2]]+[scale_1*c for c in Counts5[0:i_1]]+list(Counts7)
+
+Merged_freqs_rescaled = np.linspace(np.min(Merged_freqs),np.max(Merged_freqs),len(Merged_freqs))
+
+#drs = [391.5, 399.5, 405.5, 414]
+drs = [370,379,385,391.5]
+
+plt.figure()
+i = 0
+for j in jvec:
+    plt.plot([f-f_1 for f in Freqs5[0:i_1]], [scale_1*c for c in Counts5[0:i_1]],'o')
+    plt.plot([f-f_2 for f in Freqs6[0:i_2]], [scale_2*c for c in Counts6[0:i_2]],'o')
+    plt.plot(Freqs7, Counts7,'o')
+    plt.errorbar(Merged_freqs, Merged_counts, yerr=np.sqrt(Merged_counts), fmt='o', capsize=2, markersize=2)
+
+
+    i = i + 1
+plt.xlabel('Frecuencia (MHz)')
+plt.ylabel('counts')
+plt.grid()
+for dr in drs:
+    plt.axvline(dr)
+    plt.axvline(dr+drive, color='red', linestyle='dashed', alpha=0.3)
+    plt.axvline(dr-drive, color='red', linestyle='dashed', alpha=0.3)
+
+plt.legend()
+
+#%%
+"""
+ajusto la mergeada de 3 iones
+"""
+
+phidoppler, titadoppler = 0, 90
+phirepump, titarepump = 0,  0
+phiprobe = 0
+titaprobe = 90
+
+Temp = 0.5e-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 = 32.5e6
+
+#B = (u/(2*np.pi))/c
+
+correccion = -20
+
+offsetxpi = 438+correccion
+DetDoppler = -35-correccion-22
+
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6
+alpha = 0
+
+
+drivefreq = 2*np.pi*22.135*1e6
+
+FreqsDR = [f-offsetxpi for f in Merged_freqs]
+CountsDR = Merged_counts
+
+freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+
+CircPr = 1
+alpha = 0
+
+
+def FitEIT_MM(freqs, SG, SP, SCALE1, SCALE2, SCALE3, OFFSET, BETA1, BETA2, BETA3):
+#def FitEIT_MM(freqs, SG, SP, SCALE1, OFFSET, BETA1):
+    #BETA = 1.8
+    # SG = 0.6
+    # SP = 8.1
+    TEMP = 0.1e-3
+
+
+    #BETA1, BETA2, BETA3 = 0, 0, 2
+
+    Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+    Detunings, Fluorescence2 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA2, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+    Detunings, Fluorescence3 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA3, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+
+
+    ScaledFluo1 = np.array([f*SCALE1 for f in Fluorescence1])
+    ScaledFluo2 = np.array([f*SCALE2 for f in Fluorescence2])
+    ScaledFluo3 = np.array([f*SCALE3 for f in Fluorescence3])
+    return ScaledFluo1+ScaledFluo2+ScaledFluo3+OFFSET
+    #return ScaledFluo1
+
+
+do_fit = True
+if do_fit:
+    popt_3ions, pcov_3ions = curve_fit(FitEIT_MM, FreqsDR, CountsDR, p0=[0.6, 6.2, 3.5e5, 3.5e5, 3.5e5, 2e3, 1, 1, 1], bounds=((0, 0, 0, 0, 0, 0, 0, 0, 0), (2, 20, 5e8, 5e8, 5e8, 7e3, 10, 10, 10)))
+#popt, pcov = curve_fit(FitEIT_MM, FreqsDR, CountsDR, p0=[0.8, 8, 4e4, 3.5e3, 0], bounds=((0, 0, 0, 0, 0), (2, 15, 1e5, 1e5, 10)))
+
+#array([7.12876797e-01, 7.92474752e+00, 4.29735308e+04, 1.74240582e+04,
+       #1.53401696e+03, 1.17073206e-06, 2.53804151e+00])
+
+
+FittedEITpi_3ions = FitEIT_MM(freqslong, *popt_3ions)
+#FittedEITpi_3ions = FitEIT_MM(freqslong, popt_3ions[0],popt_3ions[1],popt_3ions[2],popt_3ions[3],popt_3ions[4],popt_3ions[5],4,2,0)
+
+#FittedEITpi_3ions = FitEIT_MM(freqslong, *popt_3ions)
+
+
+print(popt_3ions)
+
+plt.figure()
+plt.errorbar(FreqsDR, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
+plt.plot(freqslong, FittedEITpi_3ions, color='darkgreen', linewidth=3)
+#plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
+plt.xlabel('Detuning (MHz)')
+plt.ylabel('Counts')
+plt.title(f'Corr:{correccion},DetD:{DetDoppler}')
+plt.grid()
+
+
+
+
+
+#%%
+"""
+Veo la medicion de varios voltajes uno atras de otro
+Se va en medio de la medicion 4, y en la 5 ni esta
+"""
+
+jvec = [2] # desde la 1, pero la 4 no porque es un merge de curvitas
+
+Freqs
+
+plt.figure()
+i = 0
+for j in jvec:
+    plt.errorbar([2*f*1e-6 for f in Freqs[4]], CountsSplit_2ions[0][j], yerr=np.sqrt(CountsSplit_2ions[0][j]), fmt='o', capsize=2, markersize=2)
+    i = i + 1
+plt.xlabel('Frecuencia (MHz)')
+plt.ylabel('counts')
+plt.grid()
+#for dr in drs:
+#    plt.axvline(dr)
+    #plt.axvline(dr+drive)
+plt.legend()
+
+
+#%%
+#from EITfit.MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels
+from scipy.optimize import curve_fit
+import time
+
+"""
+AJUSTO LA CPT DE 2 IONES CON UN MODELO EN DONDE SUMO DOS ESPECTROS CON BETAS DISTINTOS
+"""
+
+phidoppler, titadoppler = 0, 90
+phirepump, titarepump = 0,  0
+phiprobe = 0
+titaprobe = 90
+
+Temp = 0.5e-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 = 32.5e6
+
+#B = (u/(2*np.pi))/c
+
+correccion = 27
+
+offsetxpi = 421+correccion
+DetDoppler = -16-correccion+5
+
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6
+alpha = 0
+
+
+drivefreq = 2*np.pi*22.135*1e6
+
+FreqsDR = [2*f*1e-6-offsetxpi for f in Freqs[1]]
+CountsDR = Counts[1]
+
+freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+
+CircPr = 1
+alpha = 0
+
+
+def FitEIT_MM(freqs, SG, SP, SCALE1, SCALE2, OFFSET):
+#def FitEIT_MM(freqs, SG, SP, SCALE1, OFFSET, BETA1):
+    #BETA = 1.8
+    # SG = 0.6
+    # SP = 8.1
+    TEMP = 0.1e-3
+
+    BETA1, BETA2 = 3, 0
+
+    Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+    Detunings, Fluorescence2 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA2, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+
+    ScaledFluo1 = np.array([f*SCALE1 + OFFSET for f in Fluorescence1])
+    ScaledFluo2 = np.array([f*SCALE2 + OFFSET for f in Fluorescence2])
+    return ScaledFluo1+ScaledFluo2
+    #return ScaledFluo1
+
+
+do_fit = True
+if do_fit:
+    popt_2ions_1, pcov_2ions_1 = curve_fit(FitEIT_MM, FreqsDR, CountsDR, p0=[0.9, 6.2, 3.5e3, 2.9e3, 3e3], bounds=((0, 0, 0, 0, 0), (2, 20, 5e8, 5e8, 8e3)))
+#popt, pcov = curve_fit(FitEIT_MM, FreqsDR, CountsDR, p0=[0.8, 8, 4e4, 3.5e3, 0], bounds=((0, 0, 0, 0, 0), (2, 15, 1e5, 1e5, 10)))
+
+#array([7.12876797e-01, 7.92474752e+00, 4.29735308e+04, 1.74240582e+04,
+       #1.53401696e+03, 1.17073206e-06, 2.53804151e+00])
+
+FittedEITpi_2sp = FitEIT_MM(freqslong, *popt_2ions_1)
+#FittedEITpi = FitEIT_MM(freqslong, 0.8, 8, 4e4, 3.5e3, 0)
+
+# beta1_2ions = popt_2ions_1[5]
+# beta2_2ions = popt_2ions_1[6]
+
+# errbeta1_2ions = np.sqrt(pcov_2ions_1[5,5])
+# errbeta2_2ions = np.sqrt(pcov_2ions_1[6,6])
+
+"""
+Estos params dan bien poniendo beta2=0 y correccion=0 y son SG, SP, SCALE1, SCALE2, OFFSET, BETA1
+#array([9.03123248e-01, 6.25865542e+00, 3.47684055e+04, 2.92076804e+04, 1.34556420e+03, 3.55045904e+00])
+"""
+
+"""
+Ahora considerando ambos betas, con los parametros iniciales dados por los que se obtuvieron con beta2=0
+y correccion=0 dan estos parametros que son los de antes pero con BETA2 incluido:
+array([8.52685426e-01, 7.42939084e+00, 3.61998310e+04, 3.40160472e+04, 8.62651715e+02, 3.89756335e+00, 7.64867601e-01])
+"""
+
+#arreglito = np.array([8.52685426e-01, 7.42939084e+00, 3.61998310e+04, 3.40160472e+04, 8.62651715e+02, 3.89756335e+00, 7.64867601e-01])
+
+FittedEITpi_2ions_1 = FitEIT_MM(freqslong, *popt_2ions_1)
+
+
+print(popt_2ions_1)
+
+plt.figure()
+plt.errorbar(FreqsDR, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
+plt.plot(freqslong, FittedEITpi_2ions_1, color='darkgreen', linewidth=3)
+#plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
+plt.xlabel('Detuning (MHz)')
+plt.ylabel('Counts')
+plt.title(f'Corr:{correccion},DetD:{DetDoppler}')
+plt.grid()
+
+
+#%%
+"""
+SUPER AJUSTE PARA MED DE 2 IONES
+"""
+
+
+#from EITfit.MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels
+from scipy.optimize import curve_fit
+import time
+
+
+phidoppler, titadoppler = 0, 90
+phirepump, titarepump = 0,  0
+phiprobe = 0
+titaprobe = 90
+
+Temp = 0.5e-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 = 32.5e6
+
+#B = (u/(2*np.pi))/c
+
+correccion = 13
+
+#DetDoppler = -11.5-correccion
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6
+alpha = 0
+
+
+drivefreq = 2*np.pi*22.135*1e6
+
+
+SelectedCurveVec = [3]
+
+popt_SA_vec_2ions = []
+pcov_SA_vec_2ions = []
+
+for selectedcurve in SelectedCurveVec:
+
+    FreqsDR = Freqs[selectedcurve]
+    CountsDR = Counts[selectedcurve]
+
+    freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+
+    CircPr = 1
+    alpha = 0
+
+
+    def FitEIT_MM_single(Freqs, offset, DetDoppler, SG, SP, SCALE1, SCALE2, OFFSET, BETA1, BETA2, TEMP, plot=False):
+    #def FitEIT_MM(freqs, SG, SP, SCALE1, OFFSET, BETA1):
+        #BETA = 1.8
+        #SG = 0.6
+        #SP = 8.1
+        # TEMP = 0.2e-3
+
+        freqs = [2*f*1e-6-offset for f in Freqs]
+
+        Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+        Detunings, Fluorescence2 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA2, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+
+        ScaledFluo1 = np.array([f*SCALE1 + OFFSET for f in Fluorescence1])
+        ScaledFluo2 = np.array([f*SCALE2 for f in Fluorescence2])
+        if plot:
+            return ScaledFluo1+ScaledFluo2, Detunings
+        else:
+            return ScaledFluo1+ScaledFluo2
+        #return ScaledFluo1
+
+    do_fit = True
+    if do_fit:
+        popt_3_SA_2ions, pcov_3_SA_2ions = curve_fit(FitEIT_MM_single, FreqsDR, CountsDR, p0=[448, -42, 0.6, 8.1, 4e4, 4e4, 6e3, 1, 1.2, 0.5e-3], bounds=((0, -100,0, 0, 0,0,0,0,0, 0), (1000, 0, 2, 20,5e6, 5e6,5e4, 10, 10,10e-3)))
+
+        #popt_3_SA_2ions = [448, -42, 8e4, 6e3, 2, 0.5e-3]
+
+
+        popt_SA_vec_2ions.append(popt_3_SA_2ions)
+        pcov_SA_vec_2ions.append(pcov_3_SA_2ions)
+
+        FittedEITpi_3_SA_short, Detunings_3_SA_short = FitEIT_MM_single(FreqsDR, *popt_3_SA_2ions, plot=True)
+        freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+        FittedEITpi_3_SA_long, Detunings_3_SA_long = FitEIT_MM_single(freqslong, *popt_3_SA_2ions, plot=True)
+
+
+    plt.figure()
+    plt.errorbar(Detunings_3_SA_short, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
+    plt.plot(Detunings_3_SA_long, FittedEITpi_3_SA_long, color='darkolivegreen', linewidth=3, label=f'med {selectedcurve}')
+    #plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
+    plt.xlabel('Detuning (MHz)')
+    plt.ylabel('Counts')
+    plt.legend(loc='upper left', fontsize=20)
+    plt.grid()
+
+    print(f'listo med {selectedcurve}')
+    print(popt_3_SA_2ions)
+    #print(f'Detdop:{popt_3_SA[1]},popt_3_SA:{popt[0]}')
+
+
+#%%
+#from EITfit.MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels
+from scipy.optimize import curve_fit
+import time
+
+"""
+AJUSTO LA CPT DE 2 IONES CON UN MODELO EN DONDE SUMO DOS ESPECTROS CON BETAS DISTINTOS
+"""
+
+phidoppler, titadoppler = 0, 90
+phirepump, titarepump = 0,  0
+phiprobe = 0
+titaprobe = 90
+
+Temp = 0.5e-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 = 32.5e6
+
+#B = (u/(2*np.pi))/c
+
+correccion = 27
+
+offsetxpi = 421+correccion
+DetDoppler = -16-correccion+5
+
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6
+alpha = 0
+
+
+drivefreq = 2*np.pi*22.135*1e6
+
+FreqsDR = [2*f*1e-6-offsetxpi for f in Freqs[1]]
+CountsDR = Counts[1]
+
+freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+
+CircPr = 1
+alpha = 0
+
+
+def FitEIT_MM(freqs, SG, SP, SCALE1, SCALE2, OFFSET):
+#def FitEIT_MM(freqs, SG, SP, SCALE1, OFFSET, BETA1):
+    #BETA = 1.8
+    # SG = 0.6
+    # SP = 8.1
+    TEMP = 0.1e-3
+
+    BETA1, BETA2 = 3, 0
+
+    Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+    Detunings, Fluorescence2 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA2, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+
+    ScaledFluo1 = np.array([f*SCALE1 + OFFSET for f in Fluorescence1])
+    ScaledFluo2 = np.array([f*SCALE2 + OFFSET for f in Fluorescence2])
+    return ScaledFluo1+ScaledFluo2
+    #return ScaledFluo1
+
+
+do_fit = True
+if do_fit:
+    popt_2ions_1, pcov_2ions_1 = curve_fit(FitEIT_MM, FreqsDR, CountsDR, p0=[0.9, 6.2, 3.5e3, 2.9e3, 3e3], bounds=((0, 0, 0, 0, 0), (2, 20, 5e8, 5e8, 8e3)))
+#popt, pcov = curve_fit(FitEIT_MM, FreqsDR, CountsDR, p0=[0.8, 8, 4e4, 3.5e3, 0], bounds=((0, 0, 0, 0, 0), (2, 15, 1e5, 1e5, 10)))
+
+#array([7.12876797e-01, 7.92474752e+00, 4.29735308e+04, 1.74240582e+04,
+       #1.53401696e+03, 1.17073206e-06, 2.53804151e+00])
+
+FittedEITpi_2sp = FitEIT_MM(freqslong, *popt_2ions_1)
+#FittedEITpi = FitEIT_MM(freqslong, 0.8, 8, 4e4, 3.5e3, 0)
+
+# beta1_2ions = popt_2ions_1[5]
+# beta2_2ions = popt_2ions_1[6]
+
+# errbeta1_2ions = np.sqrt(pcov_2ions_1[5,5])
+# errbeta2_2ions = np.sqrt(pcov_2ions_1[6,6])
+
+"""
+Estos params dan bien poniendo beta2=0 y correccion=0 y son SG, SP, SCALE1, SCALE2, OFFSET, BETA1
+#array([9.03123248e-01, 6.25865542e+00, 3.47684055e+04, 2.92076804e+04, 1.34556420e+03, 3.55045904e+00])
+"""
+
+"""
+Ahora considerando ambos betas, con los parametros iniciales dados por los que se obtuvieron con beta2=0
+y correccion=0 dan estos parametros que son los de antes pero con BETA2 incluido:
+array([8.52685426e-01, 7.42939084e+00, 3.61998310e+04, 3.40160472e+04, 8.62651715e+02, 3.89756335e+00, 7.64867601e-01])
+"""
+
+#arreglito = np.array([8.52685426e-01, 7.42939084e+00, 3.61998310e+04, 3.40160472e+04, 8.62651715e+02, 3.89756335e+00, 7.64867601e-01])
+
+FittedEITpi_2ions_1 = FitEIT_MM(freqslong, *popt_2ions_1)
+
+
+print(popt_2ions_1)
+
+plt.figure()
+plt.errorbar(FreqsDR, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
+plt.plot(freqslong, FittedEITpi_2ions_1, color='darkgreen', linewidth=3)
+#plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
+plt.xlabel('Detuning (MHz)')
+plt.ylabel('Counts')
+plt.title(f'Corr:{correccion},DetD:{DetDoppler}')
+plt.grid()
+
+
+#%%
+"""
+AHORA INTENTO SUPER AJUSTES O SEA CON OFFSETXPI Y DETDOPPLER INCLUIDOS
+"""
+
+
+#from EITfit.MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels
+from scipy.optimize import curve_fit
+import time
+
+"""
+SUPER AJUSTE (SA)
+
+
+"""
+
+phidoppler, titadoppler = 0, 90
+phirepump, titarepump = 0,  0
+phiprobe = 0
+titaprobe = 90
+
+Temp = 0.5e-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 = 32.5e6
+
+#B = (u/(2*np.pi))/c
+
+correccion = 13
+
+#DetDoppler = -11.5-correccion
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6
+alpha = 0
+
+
+drivefreq = 2*np.pi*22.135*1e6
+
+
+#SelectedCurveVec = [1,2,3,4,5,6,7,8,9]
+SelectedCurveVec = [0]
+
+# popt_SA_vec = []
+# pcov_SA_vec = []
+# Detuningsshort_vec = []
+# Counts_vec = []
+# Detuningslong_vec = []
+# FittedCounts_vec = []
+
+# Betas_vec = []
+# ErrorBetas_vec = []
+# Temp_vec = []
+# ErrorTemp_vec = []
+
+# DetuningsUV_vec = []
+# ErrorDetuningsUV_vec = []
+
+for selectedcurve in SelectedCurveVec:
+
+#selectedcurve = 2 #IMPORTANTE: SELECCIONA LA MEDICION
+
+    FreqsDR = Freqs[0]
+    CountsDR = CountsSplit[0][selectedcurve]
+
+
+    if selectedcurve==1:
+        CountsDR[100]=0.5*(CountsDR[99]+CountsDR[101])
+        CountsDR[105]=0.5*(CountsDR[104]+CountsDR[106])
+    if selectedcurve==2:
+        CountsDR[67]=0.5*(CountsDR[66]+CountsDR[68])
+        CountsDR[71]=0.5*(CountsDR[70]+CountsDR[72])
+    if selectedcurve==6:
+        CountsDR[1]=0.5*(CountsDR[0]+CountsDR[2])
+        CountsDR[76]=0.5*(CountsDR[75]+CountsDR[77])
+    if selectedcurve==7:
+        CountsDR[117]=0.5*(CountsDR[116]+CountsDR[118])
+
+    freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+    CircPr = 1
+    alpha = 0
+
+
+    def FitEIT_MM_single(Freqs, offset, DetDoppler, SG, SP, SCALE1, OFFSET, BETA1, TEMP, U, plot=False):
+    #def FitEIT_MM(freqs, SG, SP, SCALE1, OFFSET, BETA1):
+        #BETA = 1.8
+        # SG = 0.6
+        # SP = 8.1
+        # TEMP = 0.2e-3
+
+        freqs = [2*f*1e-6-offset for f in Freqs]
+
+        Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, U, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+
+        ScaledFluo1 = np.array([f*SCALE1 + OFFSET for f in Fluorescence1])
+        if plot:
+            return ScaledFluo1, Detunings
+        else:
+            return ScaledFluo1
+        #return ScaledFluo1
+
+    do_fit = True
+    if do_fit:
+        popt_3_SA, pcov_3_SA = curve_fit(FitEIT_MM_single, FreqsDR, CountsDR, p0=[430, -25, 0.9, 6.2, 3e4, 1.34e3, 2, (np.pi**2)*1e-3, 32e6], bounds=((0, -50, 0, 0, 0, 0, 0, 0, 25e6), (1000, 0, 2, 20, 5e4, 5e4, 10, (np.pi**2)*10e-3, 40e6)))
+
+        # popt_SA_vec.append(popt_3_SA)
+        # pcov_SA_vec.append(pcov_3_SA)
+
+        FittedEITpi_3_SA_short, Detunings_3_SA_short = FitEIT_MM_single(FreqsDR, *popt_3_SA, plot=True)
+        freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+        FittedEITpi_3_SA_long, Detunings_3_SA_long = FitEIT_MM_single(freqslong, *popt_3_SA, plot=True)
+
+    # DetuningsUV_vec.append(popt_3_SA[1])
+    # ErrorDetuningsUV_vec.append(np.sqrt(pcov_3_SA[1,1]))
+
+    # Betas_vec.append(popt_3_SA[6])
+    # ErrorBetas_vec.append(np.sqrt(pcov_3_SA[6,6]))
+    # Temp_vec.append(popt_3_SA[7])
+    # ErrorTemp_vec.append(np.sqrt(pcov_3_SA[7,7]))
+
+    # Detuningsshort_vec.append(Detunings_3_SA_short)
+    # Counts_vec.append(CountsDR)
+    # Detuningslong_vec.append(Detunings_3_SA_long)
+    # FittedCounts_vec.append(FittedEITpi_3_SA_long)
+
+
+    plt.figure()
+    plt.errorbar(Detunings_3_SA_short, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
+    plt.plot(Detunings_3_SA_long, FittedEITpi_3_SA_long, color='darkolivegreen', linewidth=3, label=f'med {selectedcurve}')
+    #plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
+    plt.xlabel('Detuning (MHz)')
+    plt.ylabel('Counts')
+    plt.legend(loc='upper left', fontsize=20)
+    plt.grid()
+
+    print(f'listo med {selectedcurve}')
+    print(popt_3_SA)