#!/usr/bin/env python3
#
'''
.. module:: fitcor
:synopsis: fit correlation of two parameters by rotating the correlation plot and fitting for the minimum width of the data on to the projected y-axis
.. moduleauthor:: Cor de Vries <c.p.de.vries@sron.nl>
'''
import numpy
import matplotlib
import matplotlib.pyplot as plt
def getdata(filename='crossdat.txt'):
'''
routine to read some test data for fitcor routine
'''
sf=open(filename,'r')
n=int(sf.readline().strip(' \n'))
c1=numpy.zeros(n,dtype=float)
c2=numpy.zeros(n,dtype=float)
for i,ll in enumerate(sf.readlines()):
tt=ll.strip(' \n').split()
c1[i]=float(tt[0])
c2[i]=float(tt[1])
return (c1,c2)
def limdist(cc):
'''
compute range of main data distribution in array of values
input:
cc = array of values
output:
xmin = minimum of range
xmax = maximum of range
'''
hh,xx=numpy.histogram(cc,bins=2500) # make histogram
imax=numpy.argmax(hh)
i=hh.size-1
for i in numpy.arange((imax+1),hh.size): # find position of first zero bin below histogram maximum
if hh[i] == 0:
break
xmax=xx[i]
for i in ((imax-1)-numpy.arange(imax)): # find position of first zero bin above histogram maximum
if hh[i] == 0:
break
xmin=xx[i]
ww=xmax-xmin
xmin=xmin-0.5*ww # extend range; add half the distribution width
xmax=xmax+0.5*ww
return (xmin,xmax)
def showcorr(c1,c2,xlabel='X',ylabel='Y',title=''):
'''
routine to plot debug information
'''
f, ax = plt.subplots(1)
ax.plot(c1,c2,'b.')
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_title(title)
plt.show()
plt.close('all')
[docs]def fitcor(c1,c2,nangles=20,dangle=0.02,discard=0.05,fftb=2,debug=False,nolim=False):
'''
fit correlation of two parameters by rotating the correlation plot and fitting for the
minimum width of the data on to the projected y-axis
Args:
* `c1` = first parameter (pulse baseline)
* `c2` = second parameter (pulse optimal filter fit)
Kwargs:
* `nangles` = number of rotation angles to probe
* `dangle` = delta of angle probe (default 0.01 ~ 0.5 degree )
* `discard` = fraction of tail of distribution to be discarded
* `fftb` = number of fft bins from lowest frequency to consider for rotation (1 or 2)
* `nolim` = if **True**, do not limit c2
* `debug` = if **True**, plot detailed fit results [default: False]
Returns:
* `cnew` = c2 corrected for fitted correlation
* `center` = center of distribution
* `ytilt` = slope of correlation
'''
if nolim:
xmax=numpy.amax(c2)
xmin=numpy.amin(c2)
else:
xmin,xmax=limdist(c2)
hh,yy=numpy.histogram(c1,bins=10000)
cc=numpy.cumsum(hh)
efrac=discard # fraction of tail of distribution to be discarded
cndx,=numpy.where( (cc > int(efrac*c1.size)) & (cc < int((1.0-efrac)*c1.size)) )
if cndx.size < 10:
cndx=numpy.arange(cc.size)
xmin=numpy.amin(c2)
xmax=numpy.amax(c2)
ymin=yy[cndx[0]]
ymax=yy[cndx[-1]+1]
indx,=numpy.where( ( c2 > xmin ) & ( c2 < xmax ) &
( c1 > ymin ) & ( c1 < ymax ) ) # index of data to be used
if debug:
showcorr(c1[indx],c2[indx],xlabel='baseline',ylabel='e-fit',title='corr. data')
iter=True
niter=nangles*2+1
angles=numpy.arange(niter,dtype=float)-nangles # rotation angles to probe
ww=numpy.zeros(niter,dtype=float)
xint=float(xmax-xmin)
yint=float(ymax-ymin)
center=numpy.array([ymin+0.5*yint,xmin+0.5*xint],dtype=float) # center of distributions
dy=dangle*xint # delta angle (about 0.5 degree default))
xminb=(xmin+xmax)/2.0-(xmax-xmin)
xmaxb=(xmin+xmax)/2.0+(xmax-xmin)
for i,angle in enumerate(angles): # do for all angles
cs=c2[indx]+(c1[indx]-ymin)/yint*dy*angle # rotate
hh,hx=numpy.histogram(cs,range=(xminb,xmaxb),bins=2048) # project data on y-axis
fh=numpy.absolute(numpy.fft.fft(hh)[1:512]) # power spectrum of projection
if i == 0:
fh0=fh[0] # normalization factor
else:
fh=fh0/fh[0]*fh # normalize
if fftb <= 1:
ww[i]=fh[0]-fh[1] # slope of power spectrum distribution at lowest frequency
else:
ww[i]=fh[0]-(fh[1]+fh[2])/2. # slope of power spectrum distribution at lowest frequency
if ( angle == 0 ) and debug:
f, ax = plt.subplots(2)
ax[0].plot(((hx[0:-1]+hx[1:])/2.0),hh,'b-')
ax[0].set_xlabel('energy coord.')
ax[0].set_ylabel('counts')
ax[0].set_title('angle=%d' % angle)
ax[1].plot(numpy.arange(1,20),fh[0:19],'b-')
ax[1].set_xlabel('frequency')
ax[1].set_ylabel('power')
plt.show()
plt.close('all')
lxww=ww.argsort() # take lowest half of data points
whalf=ww.size//2
angles=angles[lxww[:whalf]]
ww=ww[lxww[:whalf]]
sxangles=angles.argsort()
angles=angles[sxangles]
ww=ww[sxangles]
pp=numpy.polyfit(angles,ww,2) # fit power spectrum slopes as function of angle
wfit=pp[0]*angles**2+pp[1]*angles+pp[2]
fita=-0.5*pp[1]/pp[0] # fitted minimum slope
ytilt=-1.0/yint*dy*fita # compute correlation
cnew=c2+(c1-ymin)/yint*dy*fita # correct c2 for fitted correlation
if debug:
f, ax = plt.subplots(2)
ax[0].plot(angles,ww,'b-',angles,wfit,'r-')
rr=ax[0].axis()
ax[0].plot([fita,fita],[rr[2],rr[3]],'g-')
ax[0].set_xlabel('angle')
ax[0].set_ylabel('narrowness')
ax[0].set_title('result, angle=%3.1f' % fita)
ax[1].plot(c1,c2,'b.')
ax[1].set_xlabel('baseline')
ax[1].set_ylabel('e-fit')
ax[1].set_title('correlation')
rr=ax[1].axis()
y1=center[1]+(rr[0]-center[0])*ytilt
y2=center[1]+(rr[1]-center[0])*ytilt
ax[1].plot([rr[0],rr[1]],[y1,y2],'r-',label=('%5.3f' % (ytilt*1000.0)))
ax[1].legend()
plt.show()
plt.close('all')
return (cnew,center,ytilt)
# ================================================================================
if __name__ == "__main__" : # test program
c1,c2=getdata()
cnew=fitcor(c1,c2,debug=True)