RDS_piezo.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_v5/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


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()))


Piezo1Counts = []
Piezo1Frequencies = []

PIEZO1_FILES = np.arange(791, 826,1)

for i in PIEZO1_FILES:
    #print(str(i) + ' - ' + fname)
    data = h5py.File(f'VaryingBeamlocation/Piezo/000014{i}-IR_Scan_withcal_optimized'+'.h5', 'r')
    Piezo1Counts.append(np.array(data['datasets']['counts_spectrum']))
    Piezo1Frequencies.append(np.array(data['datasets']['IR1_Frequencies']))


Piezo2Counts = []
Piezo2Frequencies = []

PIEZO2_FILES = list(np.arange(834, 841,1))+list(np.arange(842, 872,1))

for i in PIEZO2_FILES:
    #print(str(i) + ' - ' + fname)
    data = h5py.File(f'VaryingBeamlocation/Piezo2/000014{i}-IR_Scan_withcal_optimized'+'.h5', 'r')
    Piezo2Counts.append(np.array(data['datasets']['counts_spectrum']))
    Piezo2Frequencies.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 1*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

TODO EL OAM
"""

palette = sns.color_palette("tab10")

pmlocmedvec = list(np.arange(0,12,1))+[13,12]+list(np.arange(15,len(PIEZO1_FILES),1))

"""
Hay que invertir la 12 con la 13, y la 14 es la misma que la 12, por las dudas
"""


#pmlocmedvec = [21]

idxvecdr = [85,185, 185,182,182,232, 162,162,175,177,217, 217,184,184,188,186,186,195,195,198,195,198,196,192,192,192,178,175,175,170,140,140,190,185]


plt.figure()

bkg = np.min(Piezo1Counts[1])

pmdepthsdr=[]
errorpmdepthsdr=[]

Intensity = []
errorIntensity = []

idxtest = 185
print(idxtest)

jj=0
for med in pmlocmedvec:
    print(med)
    
    if med == 21 or med == 22 or med == 23:
        pmdepthsdr.append(1-(Piezo1Counts[med][1:][idxvecdr[jj]]-bkg)/(np.mean(Piezo1Counts[med][1:][0:20])-bkg))
        errorpmdepthsdr.append(ErrorDRdepth(Piezo1Counts[med][1:][idxvecdr[jj]],np.mean(Piezo1Counts[med][1:][0:20]), bkg))
            
        Intens = np.mean(Piezo1Counts[med][1:][0:10])-bkg
        
        Intensity.append(Intens)
        errorIntensity.append(2*np.sqrt(np.mean(Piezo1Counts[med][1:][0:10]))+np.sqrt(bkg))

    else:
        pmdepthsdr.append(1-(Piezo1Counts[med][1:][idxvecdr[jj]]-bkg)/(np.mean(Piezo1Counts[med][1:][0:20])-bkg))
        errorpmdepthsdr.append(ErrorDRdepth(Piezo1Counts[med][1:][idxvecdr[jj]],np.mean(Piezo1Counts[med][1:][0:20]), bkg))
            
        Intens = np.mean(Piezo1Counts[med][1:][0:20])-bkg
        
        Intensity.append(Intens)
        errorIntensity.append(2*np.sqrt(np.mean(Piezo1Counts[med][1:][0:20]))+np.sqrt(bkg))
        
    
        
    plt.plot([2*f*1e-6 for f in Piezo1Frequencies[med][1:]], [c for c in Piezo1Counts[med][1:]], '-o', markersize=2, alpha=0.7)
    plt.plot([2*f*1e-6 for f in Piezo1Frequencies[med][1:]][idxtest], [c for c in Piezo1Counts[med][1:]][idxtest], 'o',markersize=14)
    jj=jj+1
# plt.xlabel('Frecuencia (MHz)')
# plt.ylabel('Counts')
# plt.xlim(432, 446.5)
# plt.grid()
# plt.legend()
# #plt.title('Espectros para distintas geometrías')


plt.figure()
plt.errorbar(np.arange(0,len(Intensity),1), [i/np.max(Intensity) for i in Intensity], yerr=[i/np.max(Intensity) for i in errorIntensity], fmt='o',capsize=2, markersize=8)
plt.errorbar(np.arange(0,len(Intensity),1), [p for p in pmdepthsdr], yerr=errorpmdepthsdr, fmt='o',capsize=2, markersize=8)
plt.xlabel('Ion position')
plt.ylabel('Intensity / DR Relative depth')
#plt.xticks([1,2,3,4,5])
#plt.xlim(200,3200)
plt.ylim(-0.1,1.1)
plt.grid()
#plt.axvline(3, color='salmon')
plt.legend()

#%%

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

TODO EL OAM
"""

palette = sns.color_palette("tab10")

mmlocmedvec = list(np.arange(0,10,1))+[11,10]+list(np.arange(12,len(PIEZO2_FILES),1))

"""
s
"""

#mmlocmedvec = [18]
#idxvecdr = [126,129,129,129,128,128,216,215,215,215,242,215,162,162,162,162,162,132,138,129,242,182,182,208,205,205,205,181,140,170,140,140,140,138,128,128,126,]
idxvecdr = [126,129,129,129,128,128,216,215,215,215,242,215,162,162,162,162,162,132,132,129,242,182,182,208,205,205,205,181,140,170,140,140,140,138,128,128,126,]
#idxtest = idxvecdr[18]

plt.figure()

bkg = np.min(Piezo2Counts[1])

mmdepthsdr=[]
errormmdepthsdr=[]

Intensity2 = []
errorIntensity2 = []


print(idxtest)

jj=0
for med in mmlocmedvec:
    print(med)
    
    mmdepthsdr.append(1-(Piezo2Counts[med][1:][idxvecdr[jj]]-bkg)/(np.mean(Piezo2Counts[med][1:][0:20])-bkg))
    errormmdepthsdr.append(ErrorDRdepth(Piezo2Counts[med][1:][idxvecdr[jj]],np.mean(Piezo2Counts[med][1:][0:20]), bkg))
        
    Intens = np.mean(Piezo2Counts[med][1:][0:20])-bkg
    
    Intensity2.append(Intens)
    errorIntensity2.append(2*np.sqrt(np.mean(Piezo2Counts[med][1:][0:20]))+np.sqrt(bkg))
      
    
        
    plt.plot([2*f*1e-6 for f in Piezo2Frequencies[med][1:]], [c for c in Piezo2Counts[med][1:]], '-o', markersize=2, alpha=0.7)
    plt.plot([2*f*1e-6 for f in Piezo2Frequencies[med][1:]][idxtest], [c for c in Piezo2Counts[med][1:]][idxtest], 'o',markersize=14)
    jj=jj+1
# plt.xlabel('Frecuencia (MHz)')
# plt.ylabel('Counts')
# plt.xlim(432, 446.5)
# plt.grid()
# plt.legend()
# #plt.title('Espectros para distintas geometrías')


plt.figure()
#plt.errorbar(np.arange(0,len(Intensity2),1), [i/np.max(Intensity2) for i in Intensity2], yerr=[i/np.max(Intensity2) for i in errorIntensity2], fmt='o',capsize=2, markersize=8)
plt.errorbar(np.arange(0,len(Intensity2),1), [p for p in mmdepthsdr], yerr=errormmdepthsdr, fmt='o',capsize=2, markersize=8)
plt.xlabel('Ion position')
plt.ylabel('Intensity / DR Relative depth')
#plt.xticks([1,2,3,4,5])
#plt.xlim(200,3200)
plt.ylim(-0.1,1.1)
plt.grid()
#plt.axvline(3, color='salmon')
plt.legend()


#%%


plt.figure()
plt.errorbar(np.arange(0,len(Intensity),1), [i/np.max(Intensity) for i in Intensity], yerr=[i/np.max(Intensity) for i in errorIntensity], fmt='o',capsize=2, markersize=8)
plt.errorbar(np.arange(0,len(Intensity),1), [p for p in pmdepthsdr], yerr=errorpmdepthsdr, fmt='o',capsize=2, markersize=8)
plt.errorbar(np.arange(0,len(Intensity),1), [p for p in mmdepthsdr[3:]], yerr=[0.5*m for m in errormmdepthsdr[3:]], fmt='o',capsize=2, markersize=8)
plt.xlabel('Ion position')
plt.ylabel('Intensity / DR Relative depth')
#plt.xticks([1,2,3,4,5])
#plt.xlim(200,3200)
plt.ylim(-0.1,1.1)
plt.grid()
#plt.axvline(3, color='salmon')
plt.legend()


#%%
"""
Intento ajustar un modelo para la profundidad
"""

def FunctionTarget(x,x0,A):
    return (A/(1+(x0/x)**2))**1

lim=10

xvec =  np.arange(0,len(Intensity),1)-15.5
xveclong = np.arange(np.min(xvec)-lim, np.max(xvec)+lim,0.01)
popt,pcov=curve_fit(FunctionTarget,xvec,pmdepthsdr)

plt.figure()
plt.errorbar(xvec, pmdepthsdr, yerr=errorpmdepthsdr, fmt='o',capsize=2, markersize=8)
#plt.plot(xveclong,FunctionTarget(xveclong,*popt))
plt.plot(xveclong,FunctionTarget(xveclong,5.93,0.83))
plt.xlabel('Ion position')
plt.ylabel('Intensity / DR Relative depth')
#plt.ylim(-0.1,1.1)
plt.grid()
#plt.axvline(3, color='salmon')
plt.legend()

print(popt)


#%%
"""
Ahora voy a intentar ajustarlas con una lorentziana que es mejor
"""
import seaborn as sns
"""
Resonancias DD configuracion +2/-2 colineal variando la ubicacion del ion en los haces

TODO EL OAM
"""
def Lorentzian( x, A, B, x0, gam ):
    return A * gam**2 / ( gam**2 + ( x - x0 )**2) + B


palette = sns.color_palette("tab10")

pmlocmedvec = list(np.arange(0,12,1))+[13,12]+list(np.arange(15,len(PIEZO1_FILES),1))

"""
Hay que invertir la 12 con la 13, y la 14 es la misma que la 12, por las dudas
"""


#pmlocmedvec = [35]


plt.figure()

bkg = np.min(Piezo1Counts[1])

pmdepthsdr=[]
errorpmdepthsdr=[]

Intensity = []
errorIntensity = []

idxtest = 185
print(idxtest)

jj=0
for med in pmlocmedvec:

    Freqs = [2*f*1e-6 for f in Piezo1Frequencies[med][1:]]
    Counts = [c for c in Piezo1Counts[med][1:]]

    if med==30:
        Freqs = Freqs[100:]
        Counts = Counts[100:]
        popt, pcov = curve_fit(Lorentzian, Freqs, Counts, p0=(-200,2100,435.8,0.05), bounds=((-10000,0,435.7,0),(0,1e4, 436.1, 1)))
        
    else:
        popt, pcov = curve_fit(Lorentzian, Freqs, Counts, p0=(-200,2100,435.8,0.05), bounds=((-10000,0,435.5,0),(0,1e4, 436.1, 1)))

    pmdepthsdr.append(1-(np.min(Lorentzian(Freqs,*popt))-bkg)/(popt[1]-bkg))
    errorpmdepthsdr.append(ErrorDRdepth(np.min(Lorentzian(Freqs,*popt)),popt[1], bkg))
        
    Intens = popt[1]
    
    Intensity.append(Intens)
    # errorIntensity.append(2*np.sqrt(np.mean(Piezo1Counts[med][1:][0:20]))+np.sqrt(bkg))
        
    
    if med in [8,21,25]:
        plt.plot([2*f*1e-6 for f in Piezo1Frequencies[med][1:]], [c for c in Piezo1Counts[med][1:]], '-o', markersize=2, alpha=0.7)
        plt.plot(Freqs,Lorentzian(Freqs,*popt))
    
    jj=jj+1
# plt.xlabel('Frecuencia (MHz)')
# plt.ylabel('Counts')
plt.xlim(435.2, 436.5)
plt.grid()
# plt.legend()
# #plt.title('Espectros para distintas geometrías')


plt.figure()
plt.plot(np.arange(0,len(Intensity),1), [i/np.max(Intensity) for i in Intensity], '-o',markersize=8)
plt.plot(np.arange(0,len(Intensity),1), [p for p in pmdepthsdr], 'o',markersize=8)
plt.xlabel('Ion position')
plt.ylabel('Intensity / DR Relative depth')
#plt.xticks([1,2,3,4,5])
#plt.xlim(200,3200)
plt.ylim(-0.1,1.1)
plt.grid()
#plt.axvline(3, color='salmon')
plt.legend()

#%%
"""
Intento ajustar un modelo para la profundidad
"""

def FunctionTarget(x,x0,A):
    return (A/(1+(x0/x)**2))**1

lim=10

xvec =  np.arange(0,len(Intensity),1)-15.5
xveclong = np.arange(np.min(xvec)-lim, np.max(xvec)+lim,0.01)
popt,pcov=curve_fit(FunctionTarget,xvec,pmdepthsdr)

plt.figure()
plt.errorbar(xvec, pmdepthsdr, yerr=errorpmdepthsdr, fmt='o',capsize=2, markersize=8)
#plt.plot(xveclong,FunctionTarget(xveclong,*popt))
plt.plot(xveclong,FunctionTarget(xveclong,5.93,0.83))
plt.xlabel('Ion position')
plt.ylabel('Intensity / DR Relative depth')
#plt.ylim(-0.1,1.1)
plt.grid()
#plt.axvline(3, color='salmon')
plt.legend()

print(popt)