import numpy as np import matplotlib.pyplot as plt import pyfits from scipy.ndimage import median_filter def waverage(x,y): return (sum(x))/sum(y) def standev(x,xbar,N): return np.sqrt(np.sum((x-xbar)**2)/(N)) f= pyfits.open("Calib_cd_6.50.fit") d = f[0].data # 2d array data h = f[0].header # fits header data px = d.sum(0) #converts to 1D array as well as sum the vlaues py = d.sum(1) # same a px but does the sum at differet axis Gain = h.get("EGAIN") # gets te gain value fgx=px*Gain # gives the number of photoelectrons #fgx= median_filter(gx, size=4)# uses a median filter to effectively remove hot pixels secdiffy = np.diff(fgx,2) xint = np.arange(0, 763, 1) bgsubmaxnorm = (fgx[0:-2] - fgx[0:-2].mean())/ fgx[0:-2].max() secdiffynorm = secdiffy / secdiffy.max() xmer=xint[secdiffynorm < -0.1] list_of_lists=[] #creates an empty array sub_list=[] sub_list.append(xmer[0]) #include the first elememt of the array for i in range(0,len(xmer)-1): # loops over from 0 to the lenght of xmer -1, main function of this loop is to form a list of lists which groups of peak values if xmer[i]+1==xmer[i+1]: # if the current element +1 is equal to the next element goes into this statement sub_list.append(xmer[i+1]) # adds the value of the next element to a list else: # however if this element's value +1 does not equal the next element then goes through this statement if len(sub_list)>0: list_of_lists.append(sub_list)# adds all the values of the sublist into list_of_lists sub_list = [] # starts a new sub_list sub_list.append(xmer[i+1]) list_of_lists.append(sub_list)#adds all of the sub_list into list_of_lists print list_of_lists xy=[] y=[] for i in list_of_lists:# goes through groups of peak values and finds the weighted average for each group. #print 'i------>', i #print 'px------>', px[i] multiple= bgsubmaxnorm[i]*i xy.append(multiple) y.append(bgsubmaxnorm[i]) num=[] denom=[] for i in list_of_lists:#This is a loop for trying to find the rms of each of the group to give the width/std of the curve. multiple= bgsubmaxnorm[i]*(np.square(i)) num.append(multiple) denom.append(bgsubmaxnorm[i]) for i in range(0,len(num)): print waverage(num[i],denom[i])# gives a value too large to be the width so different equation used. x=[] meana=[] totN=[] for i in list_of_lists:# this obtains the variables to use for the standard deviation equation for the width. x.append(i) meana.append(np.mean(i)) #totN.append(np.sum(bgsubmaxnorm[i])) totN.append(len(i)-1) print 'x--->', list_of_lists for i in range(0,len(x)):# goes through the list of values obtained in the pervious for loop and calculated the standard deviation of each. if totN[i]!=0: if standev(x[i],meana[i],totN[i])!=0: print 'std---->', 2*2.355*standev(x[i],meana[i],totN[i]) #print 'x---->', x #print 'm---->', meana #print 't---->', totN plt.figure(1) plt.title('Camera angle for 6.5mm(without median filter, normalised values)') plt.plot(bgsubmaxnorm) plt.plot(xint, secdiffynorm) plt.xlabel('Pixel values') plt.ylabel('Normalised values') plt.figure(2) #plt.subplot(211) plt.title('Camera angle for 6.5mm with peaks located (without median filter)') plt.xlabel('Pixel values') plt.ylabel('Photoelectrons') plt.plot(fgx, linewidth='2.0') for i in range(0,len(xy)): #if len(xy[i])>1: # this statement is here to get rid of any peaks that could be due to hot pixels as their "width" is too small. #print 'waverage', waverage(xy[i],y[i]) plt.axvline(waverage(xy[i],y[i]), linewidth='2.0', color='red') plt.show()