Commit a253bee2 authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

new way to produce phase models

parent 088a4ad6
......@@ -3,7 +3,8 @@
include $(SOPHYABASE)/include/sophyamake.inc
MINUITDIR = /sps/baoradio/JEC/PAON4/Minuit2-5.34.14/
### CCIN2P3 MINUITDIR = /sps/baoradio/JEC/PAON4/Minuit2-5.34.14/
MINUITDIR = /Users/campagne/Travail/kits/Minuit2-5.34.14/
MINUITLIB = -L$(MINUITDIR)/lib
MINUITINC = -I$(MINUITDIR)/include
......
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
def pcubic(x,a,b,c,d):
return a*x**3 + b*x**2 + c*x + d
def LoadPlotData(rootdir, src, polar):
''' Load freq, ch1, ch2,..., chn data '''
return np.loadtxt(rootdir+'phases-'+src+'-'+polar+'-RE.txt')
def Clean(x,y,thr=2):
'''
Fit a cubic polynomail over x, y.
Set a threshold at <thr>sigmas of residual to clean the x and y
'''
popt, pcov = curve_fit(pcubic,x,y)
res = y-pcubic(x,*popt)
resStd = np.nanstd(res)
cut = np.where( np.fabs(res) < thr*resStd )
# recompute the residual std dev on good points
resStd = np.nanstd(y[cut]-pcubic(x[cut],*popt))
return x[cut], y[cut] , resStd
def FitPolynomial(x,y,degMax=30,iqrMax=0.005):
polDeg = 2
#better conditionnement
xmean = np.mean(x)
xshifted = x-xmean
rc = 1
while polDeg < degMax :
polDeg += 1
coeff = np.polyfit(xshifted, y, polDeg)
residuals = (y - np.poly1d(coeff)(xshifted))
q75, q25 = np.percentile(residuals/(2*np.pi), [75 ,25])
iqrn = (q75 - q25)/1.34898
if iqrn <= iqrMax:
rc = 0
break
return rc, xmean, polDeg, coeff, residuals, iqrn
def FitPolar(data):
fig, axs = plt.subplots(2*(data.shape[1]-1),1,sharey=False,sharex=True,figsize=(7,7))
freq = data[:,0]
infos_list = []
for i in range(1,data.shape[1]):
xclean, yclean, _ = Clean(freq, data[:,i])
rc, fcenter, deg, coeff, res, iqrn = FitPolynomial(xclean,yclean)
info = np.array([i+1])
info = np.append(info,[fcenter])
info = np.append(info,[deg])
info = np.append(info,coeff)
infos_list.append(info)
idx0 = 2*i-2
idx1 = idx0 + 1
axs[idx0].plot(freq, data[:,i],label='Phi'+str(i)+polar)
axs[idx0].plot(xclean, yclean,ls='--',label='Phi'+str(i)+polar+' cleaned')
axs[idx0].plot(freq,np.poly1d(coeff)(freq-fcenter),'g-',label='fit'+str(deg))
axs[idx0].set_ylim(-np.pi,np.pi)
axs[idx0].grid()
axs[idx0].legend(loc=1)
axs[idx1].plot(xclean,res,'g-',label='residuals')
axs[idx1].set_ylim(-0.1,0.1)
axs[idx1].grid()
axs[idx1].legend(loc=1)
ax0 = fig.add_subplot(111, frame_on=False) # creating a single axes
ax0.set_xticks([])
ax0.set_yticks([])
ax0.set_xlabel('Freq (MHz)', labelpad=25)
ax0.set_ylabel('Phase (rad)', labelpad=45)
fig.suptitle(src+' (polar:'+polar+')')
return infos_list
###
rootdir = './CasA/'
src='hv_A00_17juil18'
polars = ['H','V']
for idx, polar in enumerate(polars):
data = LoadPlotData(rootdir, src, polar)
infos = FitPolar(data)
filename = 'PhaseParam-'+src+'-'+polar+'.txt'
with open(filename,"w") as f:
f.write("\n".join(" ".join(map(str,x)) for x in infos))
#for bookkeeping
plt.savefig(src+'-'+polar+'-fit.pdf',bbox_inches = 'tight')
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment