The goal is to digitize the data in a plot of transmission vs. wavelength so that it can be used in the ALICE FIT simulation.
In [1]:
%pylab inline
import matplotlib.pyplot as plt
import numpy as np
In [2]:
import scipy.ndimage as ndi
In [3]:
img = ndi.imread("Cr_quartz2_absorption.png")
print img.shape
In [4]:
#Find indices where the image is red
c = (255, 2, 2)
indices = numpy.where(numpy.all(img == c, axis=-1))
rr = zip(indices[1],297-indices[0]) #zip them together as (x,y)
rr.sort() #sort by x value
ss = np.array(zip(*rr)) #unzip into numpy array for plotting
#Find top line pixel location
c2 = (48, 48, 48)
ind2 = numpy.where(numpy.all(img == c2, axis=-1))
#print zip(ind2[0], ind2[1])
In [5]:
#Plot the line at the locations of the indices, inverting y and starting at 0
plt.plot(ss[0], ss[1],'r-')
plt.plot(ind2[1],297-ind2[0],'ko')
Out[5]:
In [6]:
plt.figure(figsize=(20,10))
plt.imshow(img)
plt.plot(ss[0][82],297-ss[1][82],'bo') #Note, to plot on the image, have to invert y values
plt.plot(ss[0][0:304], 297-ss[1][0:304],'b+')
plt.plot(ss[0][313],297-ss[1][313],'mo')
plt.plot(ss[0][305:327],297-ss[1][305:327],'m+')
plt.plot(ss[0][328],297-ss[1][328],'go')
plt.plot(ss[0][328:405],297-ss[1][328:405],'g+')
plt.plot(ss[0][415],297-ss[1][415],'co')
plt.plot(ss[0][406:423],297-ss[1][406:423],'c+')
plt.plot(ss[0][700],297-ss[1][700],'ko')
plt.plot(ss[0][424:765],297-ss[1][424:765],'k+')
Out[6]:
In [7]:
print "(%d,%d) = (143,1)"%(ss[0][0],ss[1][0])
print "(%d,%d) = (240,88)"%(ss[0][304],ss[1][304])
x1b=float(ss[0][0]); y1b=float(ss[1][0])
x2b=float(ss[0][304]); y2b=float(ss[1][304])
xab=143.; yab=0.
xbb=240.; ybb=88.
sbx=(xbb-xab)/(x2b-x1b)
sby=(ybb-yab)/(y2b-y1b)
#print sb,sby
x = (x1b-x1b)*sbx + xab
y = (y1b-y1b)*sby + yab
print x1b,y1b," = ",x,y
x = (x2b-x1b)*sbx + xab
y = (y2b-y1b)*sby + yab
print x2b,y2b," = ",x,y
In [8]:
print (ss[0][82],ss[1][82])
print (ss[0][82]-x1b)*sbx + xab,(ss[1][82]-y1b)*sby + yab
In [9]:
print "(%d,%d) = (240,88)"%(ss[0][305],ss[1][305])
print "(%d,%d) = (300,89)"%(ss[0][327],ss[1][327])
x1m=float(ss[0][305]); y1m=float(ss[1][305])
x2m=float(ss[0][327]); y2m=float(ss[1][327])
xam=241.; yam=88.
xbm=300.; ybm=89.
smx=(xbm-xam)/(x2m-x1m)
smy=(ybm-yam)/(y2m-y1m)
x = (x1m-x1m)*smx + xam
y = (y1m-y1m)*smy + yam
print x1m,y1m," = ",x,y
x = (x2m-x1m)*smx + xam
y = (y2m-y1m)*smy + yam
print x2m,y2m," = ",x,y
In [10]:
print (ss[0][313],ss[1][313])
print (ss[0][313]-x1m)*smx + xam,(ss[1][313]-y1m)*smy + yam
In [11]:
print "(%d,%d) = (300,89)"%(ss[0][328],ss[1][328])
print "(%d,%d) = (700,90)"%(ss[0][405],ss[1][405])
x1g=float(ss[0][328]); y1g=float(ss[1][328])
x2g=float(ss[0][405]); y2g=float(ss[1][405])
xag=301.; yag=89.
xbg=700.; ybg=90.
sgx=(xbg-xag)/(x2g-x1g)
sgy=(ybg-yag)/(y2g-y1g)
x = (x1g-x1g)*sgx + xag
y = (y1g-y1g)*sgy + yag
print x1g,y1g," = ",x,y
x = (x2g-x1g)*sgx + xag
y = (y2g-y1g)*sgy + yag
print x2g,y2g," = ",x,y
In [12]:
print (ss[0][328],ss[1][328])
print (ss[0][328]-x1g)*sgx + xag,(ss[1][328]-y1g)*sgy + yag
In [13]:
print "(%d,%d) = (700,90)"%(ss[0][406],ss[1][406])
print "(%d,%d) = (1000,91)"%(ss[0][423],ss[1][406])
x1c=float(ss[0][406]); y1c= float(ss[1][406])
x2c=float(ss[0][423]); y2c= float(ss[1][423])
xac=701.; yac=90.
xbc=1000.; ybc=91.
scx=(xbc-xac)/(x2c-x1c)
scy=(ybc-yac)/(y2c-y1c)
x = (x1c-x1c)*scx + xac
y = (y1c-y1c)*scy + yac
print x1c,y1c," = ",x,y
x = (x2c-x1c)*scx + xac
y = (y2c-y1c)*scy + yac
print x2c,y2c," = ",x,y
In [14]:
print (ss[0][415],ss[1][415])
print (ss[0][415]-x1c)*scx + xac,(ss[1][415]-y1c)*scy + yac
In [15]:
print "(%d,%d) = (1000,91)"%(ss[0][424],ss[1][424])
print "(%d,%d) = (4750,0)"%(ss[0][765],ss[1][765])
x1k=float(ss[0][424]); y1k=float(ss[1][424])
x2k=float(ss[0][765]); y2k=float(ss[1][765])
xak=1001.; yak=91.
xbk=4750.; ybk=0.
skx=(xbk-xak)/(x2k-x1k)
sky=(ybk-yak)/(y2k-y1k)
x = (x1k-x1k)*skx + xak
y = (y1k-y1k)*sky + yak
print x1k,y1k," = ",x,y
x = (x2k-x1k)*skx + xak
y = (y2k-y1k)*sky + yak
print x2k,y2k," = ",x,y
In [16]:
print (ss[0][524],ss[1][524])
print (ss[0][524]-x1k)*skx + xak,(ss[1][524]-y1k)*sky + yak
In [110]:
newx = np.zeros(780,dtype=float64)
newy = np.zeros(780,dtype=float64)
for i in range(0,305):
x1b=float(ss[0][0]); y1b=float(ss[1][0])
x2b=float(ss[0][304]); y2b=float(ss[1][304])
xab=143.; yab=0.
xbb=240.; ybb=88.
sbx=(xbb-xab)/(x2b-x1b)
sby=(ybb-yab)/(y2b-y1b)
newx[i] = (ss[0][i]-x1b)*sbx + xab
newy[i] = (ss[1][i]-y1b)*sby + yab
for i in range(305,328):
x1m=float(ss[0][305]); y1m=float(ss[1][305])
x2m=float(ss[0][327]); y2m=float(ss[1][327])
xam=241.; yam=88.
xbm=300.; ybm=89.
smx=(xbm-xam)/(x2m-x1m)
smy=(ybm-yam)/(y2m-y1m)
newx[i] = (ss[0][i]-x1m)*smx + xam
newy[i] = (ss[1][i]-y1m)*smy + yam
for i in range(328,406):
x1g=float(ss[0][328]); y1g=float(ss[1][328])
x2g=float(ss[0][405]); y2g=float(ss[1][405])
xag=301.; yag=89.
xbg=700.; ybg=90.
sgx=(xbg-xag)/(x2g-x1g)
sgy=(ybg-yag)/(y2g-y1g)
newx[i] = (ss[0][i]-x1g)*sgx + xag
newy[i] = (ss[1][i]-y1g)*sgy + yag
for i in range(406,424):
x1c=float(ss[0][406]); y1c= float(ss[1][406])
x2c=float(ss[0][423]); y2c= float(ss[1][423])
xac=701.; yac=90.
xbc=1000.; ybc=91.
scx=(xbc-xac)/(x2c-x1c)
scy=(ybc-yac)/(y2c-y1c)
newx[i] = (ss[0][i]-x1c)*scx + xac
newy[i] = (ss[1][i]-y1c)*scy + yac
for i in range(424,766):
x1k=float(ss[0][424]); y1k=float(ss[1][424])
x2k=float(ss[0][765]); y2k=float(ss[1][765])
xak=1001.; yak=91.
xbk=4750.; ybk=0.
skx=(xbk-xak)/(x2k-x1k)
sky=(ybk-yak)/(y2k-y1k)
newx[i] = (ss[0][i]-x1k)*skx + xak
newy[i] = (ss[1][i]-y1k)*sky + yak
for i in range(766,780):
newx[i] = newx[i-1]+142.
newy[i] = 0.
In [112]:
plt.plot(newx,newy,'r-');
plt.xlim(10,5000)
plt.ylim(-1,100)
plt.xlabel("Wavelength (nm)",fontsize=20)
plt.ylabel("Transmission (%)",fontsize=20)
Out[112]:
In [113]:
plt.plot(newx,newy,'r-');
plt.xlim(2500,5000)
plt.ylim(-1,100)
plt.xlabel("Wavelength (nm)",fontsize=20)
plt.ylabel("Transmission (%)",fontsize=20)
Out[113]:
In [23]:
#Import the stats package
from scipy import stats
from scipy.optimize import curve_fit
In [76]:
print "Range of indices to fit over:"
print (np.where(newx > 4247.))
In [87]:
fitfun = lambda x, a, b, c, d: a + b*x + c*x**2 + d*x**3
funx = np.linspace(0,5000,5000)
In [88]:
p0 = np.array([0.,1.,1.,1.])
popt, pcov = curve_fit(fitfun, newx[561:689], newy[561:689], p0)
In [94]:
plt.plot(newx,newy,'r-');
plt.plot(funx[3390:3780],fitfun(funx[3390:3780],popt[0],popt[1],popt[2],popt[3]),'b-')
plt.xlim(2500,5000)
plt.ylim(-1,100)
plt.xlabel("Wavelength (nm)",fontsize=20)
plt.ylabel("Transmission (%)",fontsize=20)
Out[94]:
In [66]:
p0b = np.array([0.,1.,1.,1.])
poptb, pcovb = curve_fit(fitfun, newx[689:745], newy[689:745], p0b)
In [114]:
plt.plot(newx,newy,'r-');
plt.plot(funx[3390:3780],fitfun(funx[3390:3780],popt[0],popt[1],popt[2],popt[3]),'b-')
plt.plot(funx[3810:4200],fitfun(funx[3810:4200],poptb[0],poptb[1],poptb[2],poptb[3]),'g-')
plt.xlim(3500,5000)
plt.ylim(-1,30)
plt.xlabel("Wavelength (nm)",fontsize=20)
plt.ylabel("Transmission (%)",fontsize=20)
Out[114]:
In [118]:
p0c = np.array([0.,1.,1.,1.])
poptc, pcovc = curve_fit(fitfun, newx[760:], newy[760:], p0c)
In [119]:
plt.plot(newx,newy,'r-');
plt.plot(funx[4200:],fitfun(funx[4200:],poptc[0],poptc[1],poptc[2],poptc[3]),'b-')
plt.xlim(4000,5000)
plt.ylim(-1,20)
plt.xlabel("Wavelength (nm)",fontsize=20)
plt.ylabel("Transmission (%)",fontsize=20)
Out[119]:
In [86]:
funx = np.linspace(0,5000,5000)
plt.plot(newx,newy,'r-');
plt.plot(funx,fitfun(funx,poptc[0],poptc[1],poptc[2],poptc[3]),'b-')
plt.xlim(4000,5000)
plt.ylim(-1,10)
plt.xlabel("Wavelength (nm)",fontsize=20)
plt.ylabel("Transmission (%)",fontsize=20)
Out[86]:
In [85]:
np.where(fitfun(funx,poptc[0],poptc[1],poptc[2],poptc[3])<0)
Out[85]:
In [280]:
x,y = np.loadtxt("DefaultDataset.csv",delimiter=',',dtype=float64,unpack=True)
In [281]:
plt.plot(x,y,'k-')
Out[281]:
In [ ]: