RDS_location4.py

import h5py
import matplotlib.pyplot as plt
import numpy as np
import sys
import re
import ast
from scipy.optimize import curve_fit
import os
from scipy import interpolate

#Mediciones barriendo angulo del TISA y viendo kicking de resonancias oscuras

#C:\Users\Usuario\Documents\artiq\artiq_experiments\analisis\plots\20220106_CPT_DosLaseres_v08_TISA_DR\Data

os.chdir('/home/nico/Documents/artiq_experiments/analisis/plots/20230817_RotationalDopplerShift_v4/Data')



"""
en este codigo ploteo espectros CPT de resonancias D-D para configuracion +2/+2 y +2/-2 (usando pentaprisma)
"""

def find_nearest(array, value):
    array = np.asarray(array)
    idx = (np.abs(array - value)).argmin()
    return idx


PM_FILES = """VaryingBeamlocation/Plusminus/000014461-IR_Scan_withcal_optimized
VaryingBeamlocation/Plusminus/000014462-IR_Scan_withcal_optimized
VaryingBeamlocation/Plusminus/000014463-IR_Scan_withcal_optimized
VaryingBeamlocation/Plusminus/000014464-IR_Scan_withcal_optimized
VaryingBeamlocation/Plusminus/000014465-IR_Scan_withcal_optimized
VaryingBeamlocation/Plusminus/000014466-IR_Scan_withcal_optimized
VaryingBeamlocation/Plusminus/000014467-IR_Scan_withcal_optimized
VaryingBeamlocation/Plusminus/000014468-IR_Scan_withcal_optimized
"""


MM_FILES = """VaryingBeamlocation/Minusminus/000014491-IR_Scan_withcal_optimized
VaryingBeamlocation/Minusminus/000014494-IR_Scan_withcal_optimized
VaryingBeamlocation/Minusminus/000014495-IR_Scan_withcal_optimized
VaryingBeamlocation/Minusminus/000014496-IR_Scan_withcal_optimized
VaryingBeamlocation/Minusminus/000014497-IR_Scan_withcal_optimized
VaryingBeamlocation/Minusminus/000014499-IR_Scan_withcal_optimized
VaryingBeamlocation/Minusminus/000014500-IR_Scan_withcal_optimized
VaryingBeamlocation/Minusminus/000014502-IR_Scan_withcal_optimized
"""


#HS = high statistics
MM_FILES_HS = """VaryingBeamlocation/Minusminus/000014492-IR_Scan_withcal_optimized
VaryingBeamlocation/Minusminus/000014498-IR_Scan_withcal_optimized
VaryingBeamlocation/Minusminus/000014501-IR_Scan_withcal_optimized
"""

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


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(PM_FILES))
print(SeeKeys(MM_FILES))


#carpeta pc nico labo escritorio:
#C:\Users\Usuario\Documents\artiq\artiq_experiments\analisis\plots\20211101_CPT_DosLaseres_v03\Data

PmCounts = []
PmFrequencies = []

for i, fname in enumerate(PM_FILES.split()):
    print(str(i) + ' - ' + fname)
    data = h5py.File(fname+'.h5', 'r')
    #Amplitudes.append(np.array(data['datasets']['amplitudes']))
    PmCounts.append(np.array(data['datasets']['counts_spectrum']))
    PmFrequencies.append(np.array(data['datasets']['IR1_Frequencies']))

MmCounts = []
MmFrequencies = []

for i, fname in enumerate(MM_FILES.split()):
    print(str(i) + ' - ' + fname)
    data = h5py.File(fname+'.h5', 'r')
    #Amplitudes.append(np.array(data['datasets']['amplitudes']))
    MmCounts.append(np.array(data['datasets']['counts_spectrum']))
    MmFrequencies.append(np.array(data['datasets']['IR1_Frequencies']))

HS_MmCounts = []
HS_MmFrequencies = []

for i, fname in enumerate(MM_FILES_HS.split()):
    print(str(i) + ' - ' + fname)
    data = h5py.File(fname+'.h5', 'r')
    #Amplitudes.append(np.array(data['datasets']['amplitudes']))
    HS_MmCounts.append(np.array(data['datasets']['counts_spectrum']))
    HS_MmFrequencies.append(np.array(data['datasets']['IR1_Frequencies']))
    

def ErrorDRdepth(p, f, b):
    ep = np.sqrt(p)
    ef = np.sqrt(f)
    eb = np.sqrt(b)
    derivadap = 1/((f-b)**2)
    derivadaf = ((p-b)/((f-b)**2))**2
    derivadab = ((p-f)/((f-b)**2))**2
    return 2*np.sqrt(derivadap*ep*ep + derivadaf*ef*ef + derivadab*eb*eb)

#%%
import seaborn as sns
"""
Resonancias DD configuracion +2/-2 colineal variando la ubicacion del ion en los haces
"""

palette = sns.color_palette("tab10")

# pmlocmedvec = [0, 4, 3, 5, 1, 7, 6, 2]
# idxvec = [95, 95, 95, 94, 93, 98, 97, 97]

pmlocmedvec = [7,2,0,4,3,5,1]

idxvecdr1 = [98,97, 95,95,95,94,93]
idxvecdr2 = [238, 233, 236, 236, 235, 235,235]

#pmlocmedvec = [2]

LocVecs = ['1', '2', '3', '4', '5', '6', '7', '8']

LocVecsint = np.arange(1,8,1)

plt.figure()

ftrap = 22.1

DR1 = 435.8
DR2 = 444.2

bkg = 120

pmdepthsdr1=[]
errorpmdepthsdr1=[]

pmdepthsdr2=[]
errorpmdepthsdr2=[]


idxtest = 235

jj=0
for med in pmlocmedvec:
    pmdepthsdr1.append(1-(PmCounts[med][1:][idxvecdr1[jj]]-bkg)/(np.mean(PmCounts[med][1:][0:20])-bkg))
    pmdepthsdr2.append(1-(PmCounts[med][1:][idxvecdr2[jj]]-bkg)/(np.mean(PmCounts[med][1:][0:20])-bkg))
    errorpmdepthsdr1.append(ErrorDRdepth(PmCounts[med][1:][idxvecdr1[jj]],np.mean(PmCounts[med][1:][0:20]), bkg))
    errorpmdepthsdr2.append(ErrorDRdepth(PmCounts[med][1:][idxvecdr2[jj]],np.mean(PmCounts[med][1:][0:20]), bkg))
    plt.plot([2*f*1e-6 for f in PmFrequencies[med][1:]], [c for c in PmCounts[med][1:]], '-o', color=palette[med],markersize=2, alpha=0.7, label=f'{LocVecs[jj]}')
    #plt.plot([2*f*1e-6 for f in PmFrequencies[med][1:]][idxtest], [c for c in PmCounts[med][1:]][idxtest], 'o', color=palette[med],markersize=14, label=f'{LocVecs[jj]}')
    jj=jj+1
plt.xlabel('Frecuencia (MHz)')
plt.ylabel('Counts')
#plt.xlim(442,445)
plt.grid()
plt.legend()
#plt.title('Espectros para distintas geometrías')


plt.figure()
plt.errorbar(LocVecsint, pmdepthsdr1, yerr=errorpmdepthsdr1, fmt='o',capsize=2, markersize=5)
plt.errorbar(LocVecsint, pmdepthsdr2, yerr=errorpmdepthsdr2, fmt='o',capsize=2, markersize=5)
plt.ylabel('DR1 Relative depth config -2/+2')
plt.xticks([])
plt.ylim(0,1)


#%%
import seaborn as sns
"""
Resonancias DD configuracion -2/-2 colineal variando la ubicacion del ion en los 
Mido la profundidad de las resonancias oscuras
"""

palette = sns.color_palette("tab10")

# pmlocmedvec = [0, 4, 3, 5, 1, 7, 6, 2]
# idxvec = [95, 95, 95, 94, 93, 98, 97, 97]

mmlocmedvec = [0,1,2,3,4,5,6]

idxvecdr1 = [160, 160, 160, 160, 159,159, 159]
idxvecdr2 = [394, 394, 394, 394, 394, 394, 394]

#mmlocmedvec = [7]

LocVecs = ['1', '2', '3', '4', '5', '6', '7', '8']

LocVecsint = np.arange(1,8,1)

plt.figure()

ftrap = 22.1

DR1 = 435.8
DR2 = 444.2

bkg = 40

mmdepthsdr1=[]
errormmdepthsdr1=[]

mmdepthsdr2=[]
errormmdepthsdr2=[]


idxtest = 394

jj=0
for med in mmlocmedvec:
    mmdepthsdr1.append(1-(MmCounts[med][1:][idxvecdr1[jj]]-bkg)/(np.mean(MmCounts[med][1:][0:20])-bkg))
    mmdepthsdr2.append(1-(MmCounts[med][1:][idxvecdr2[jj]]-bkg)/(np.mean(MmCounts[med][1:][0:20])-bkg))
    errormmdepthsdr1.append(ErrorDRdepth(MmCounts[med][1:][idxvecdr1[jj]],np.mean(MmCounts[med][1:][0:20]), bkg))
    errormmdepthsdr2.append(ErrorDRdepth(MmCounts[med][1:][idxvecdr2[jj]],np.mean(MmCounts[med][1:][0:20]), bkg))
    plt.plot([2*f*1e-6 for f in MmFrequencies[med][1:]], [c for c in MmCounts[med][1:]], '-o', color=palette[med],markersize=2, alpha=0.7, label=f'{LocVecs[jj]}')
    #plt.plot([2*f*1e-6 for f in MmFrequencies[med][1:]][idxtest], [c for c in MmCounts[med][1:]][idxtest], 'o', color=palette[med],markersize=14, label=f'{LocVecs[jj]}')
    jj=jj+1
plt.xlabel('Frecuencia (MHz)')
plt.ylabel('Counts')
#plt.xlim(442,445)
plt.grid()
plt.legend()
plt.title('Espectros para distintas geometrías')


plt.figure()
plt.errorbar(LocVecsint, mmdepthsdr1, yerr=errormmdepthsdr1, fmt='o',capsize=2, markersize=5)
plt.errorbar(LocVecsint, mmdepthsdr2, yerr=errormmdepthsdr2, fmt='o',capsize=2, markersize=5)
plt.ylabel('DR1 Relative depth config -2/-2')
plt.xticks([])
plt.ylim(0,1)

#%%
"""
grafico la profundidad de la resonancia oscura en funcion de la ubicacion del ion
en el perfil del haz. la ubicacion 1 y 7 son horizontales, izquierda y derecha.
la ubicacion 4 es el centro, arriba. el resto, intermedias
"""

plt.figure()
plt.errorbar(LocVecsint, mmdepthsdr1, yerr=errormmdepthsdr1, fmt='o',capsize=2, markersize=6, color='darkred', label='Same OAM')
plt.errorbar(LocVecsint, pmdepthsdr1, yerr=errorpmdepthsdr2, fmt='o',capsize=2, markersize=6, color='coral', label='Opposite OAM')
plt.ylabel('DR1 Relative depth')
plt.xlabel('Ion angular position')
#plt.xticks([])
plt.grid()
plt.ylim(0,1)
plt.legend()

#%%
"""
Grafico ahora dos gráficos que se vean 3 curvas con buena estadística para mostrar
"""
from scipy.optimize import curve_fit
from scipy.signal import savgol_filter as sf

def linealfunc(x,a,b):
    return a*x+b

def remove_outliers(array):
    """
    work in progress...
    """
    array2 = []
    threshold = 300
    for jj in range(len(array)):
        if jj == len(array)-1:
            array2.append(array[jj])
        else:
            if np.abs(array[jj+1]-array[jj])<threshold:
               array2.append(array[jj])
            else:
                array2.append((array[jj+3]+array[jj])*0.5)
    return array2

def normalizeplot(xvec, yvec,skip=0):
    i_1, i_2, i_3, i_4 = 0, 20, -10-skip, -1-skip
    popt, pcov = curve_fit(linealfunc, list(xvec[i_1:i_2])+list(xvec[i_3:i_4]), list(yvec[i_1:i_2])+list(yvec[i_3:i_4]))
    print(popt)
    yvecnorm = []
    for jj in range(len(yvec)):
        print(linealfunc(xvec[jj],*popt))
        yvecnorm.append(yvec[jj]/linealfunc(xvec[jj],*popt))
    return yvecnorm
#pmlocmedvec = [7,2,0,4,3,5,1]



pmlocmedvec = [3,5,7]
mmlocmedvec = [0,1,2]

plt.figure()

plotpm = 1
plotmm = 0


if plotpm:
    jj=0
    for med in pmlocmedvec:
        Freqspm = [2*f*1e-6 for f in PmFrequencies[med][1:]]
        Countspm = [c-120 for c in PmCounts[med][1:]]
        NormCountspm = normalizeplot(Freqspm, Countspm)
        plt.plot(Freqspm, NormCountspm, '-o', color=palette[jj], markersize=2, alpha=1, label=f'{LocVecs[jj]}')
        #plt.plot(Freqs, sf(NormCounts,5,2), 'o', color=palette[jj], markersize=2, label=f'{LocVecs[jj]}')
        #plt.plot([2*f*1e-6 for f in PmFrequencies[med][1:]][idxtest], [c for c in PmCounts[med][1:]][idxtest], 'o', color=palette[med],markersize=14, label=f'{LocVecs[jj]}')
        jj=jj+1

if plotmm:
    jj=0
    for med in mmlocmedvec:
        Freqsmm = [2*f*1e-6 for f in HS_MmFrequencies[med][1:-8]]
        Countsmm = [c-120 for c in HS_MmCounts[med][1:-8]]
        NormCountsmm = normalizeplot(Freqsmm, Countsmm, skip=0)
        plt.plot(Freqsmm, NormCountsmm, '-o', color=palette[jj], markersize=2, alpha=1, label=f'{LocVecs[jj]}')
        #plt.plot(Freqs, sf(NormCounts,5,2), 'o', color=palette[jj], markersize=2, label=f'{LocVecs[jj]}')
        #plt.plot([2*f*1e-6 for f in PmFrequencies[med][1:]][idxtest], [c for c in PmCounts[med][1:]][idxtest], 'o', color=palette[med],markersize=14, label=f'{LocVecs[jj]}')
        jj=jj+1

plt.ylim(0,1.2)
plt.xlabel('Frecuencia (MHz)')
plt.ylabel('Normalized Counts')
#plt.ylim(412,2500)
plt.grid()
#plt.legend()