manipulations.py @ 3ba03a4c

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since 3ba03a4c was 8ba103f, checked in by Jae Cho <jhjcho@…>, 16 years ago
fixed circular average bug
Property mode set to `100644`
File size: 43.7 KB

Rev	Line
[76e2369]	1	"""
[f8d0ee7]	2	Data manipulations for 2D data sets.
	3	Using the meta data information, various types of averaging
	4	are performed in Q-space
[76e2369]	5	"""
	6
	7	"""
	8	This software was developed by the University of Tennessee as part of the
	9	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
	10	project funded by the US National Science Foundation.
	11
	12	See the license text in license.txt
	13
	14	copyright 2008, University of Tennessee
	15	"""
	16	#TODO: copy the meta data from the 2D object to the resulting 1D object
	17
	18	from data_info import plottable_2D, Data1D
	19	import math
	20	import numpy
	21
	22	def get_q(dx, dy, det_dist, wavelength):
	23	"""
	24	@param dx: x-distance from beam center [mm]
	25	@param dy: y-distance from beam center [mm]
	26	@return: q-value at the given position
	27	"""
	28	# Distance from beam center in the plane of detector
	29	plane_dist = math.sqrt(dxdx + dydy)
	30	# Half of the scattering angle
	31	theta = 0.5*math.atan(plane_dist/det_dist)
	32	return (4.0math.pi/wavelength)math.sin(theta)
[acb37d9]	33
	34	def get_q_compo(dx, dy, det_dist, wavelength,compo=None):
	35	#This reduces tiny error at very large q.
	36	#Implementation of this func is not started yet.<--ToDo
	37	if dy==0:
	38	if dx>=0:
	39	angle_xy=0
	40	else:
	41	angle_xy=math.pi
	42	else:
	43	angle_xy=math.atan(dx/dy)
	44
	45	if compo=="x":
	46	out=get_q(dx, dy, det_dist, wavelength)*cos(angle_xy)
	47	elif compo=="y":
	48	out=get_q(dx, dy, det_dist, wavelength)*sin(angle_xy)
	49	else:
	50	out=get_q(dx, dy, det_dist, wavelength)
	51	return out
[76e2369]	52
[70975f3]	53	class _Slab(object):
	54	"""
	55	Compute average I(Q) for a region of interest
	56	"""
	57	def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0, bin_width=0.001):
	58	# Minimum Qx value [A-1]
	59	self.x_min = x_min
	60	# Maximum Qx value [A-1]
	61	self.x_max = x_max
	62	# Minimum Qy value [A-1]
	63	self.y_min = y_min
	64	# Maximum Qy value [A-1]
	65	self.y_max = y_max
	66	# Bin width (step size) [A-1]
	67	self.bin_width = bin_width
	68	# If True, I(\|Q\|) will be return, otherwise, negative q-values are allowed
	69	self.fold = False
	70
	71	def __call__(self, data2D): return NotImplemented
	72
	73	def _avg(self, data2D, maj):
	74	"""
	75	Compute average I(Q_maj) for a region of interest.
	76	The major axis is defined as the axis of Q_maj.
	77	The minor axis is the axis that we average over.
	78
	79	@param data2D: Data2D object
	80	@param maj_min: min value on the major axis
	81	@return: Data1D object
	82	"""
	83	if len(data2D.detector) != 1:
	84	raise RuntimeError, "_Slab._avg: invalid number of detectors: %g" % len(data2D.detector)
	85
	86	pixel_width_x = data2D.detector[0].pixel_size.x
	87	pixel_width_y = data2D.detector[0].pixel_size.y
	88	det_dist = data2D.detector[0].distance
	89	wavelength = data2D.source.wavelength
	90	center_x = data2D.detector[0].beam_center.x/pixel_width_x
	91	center_y = data2D.detector[0].beam_center.y/pixel_width_y
	92
	93	# Build array of Q intervals
	94	if maj=='x':
	95	nbins = int(math.ceil((self.x_max-self.x_min)/self.bin_width))
	96	qbins = self.bin_width*numpy.arange(nbins)+self.x_min
	97	elif maj=='y':
	98	nbins = int(math.ceil((self.y_max-self.y_min)/self.bin_width))
	99	qbins = self.bin_width*numpy.arange(nbins)+self.y_min
	100	else:
	101	raise RuntimeError, "_Slab._avg: unrecognized axis %s" % str(maj)
	102
	103	x = numpy.zeros(nbins)
	104	y = numpy.zeros(nbins)
	105	err_y = numpy.zeros(nbins)
	106	y_counts = numpy.zeros(nbins)
[ee3b2c0c]	107
[70975f3]	108	for i in range(numpy.size(data2D.data,1)):
	109	# Min and max x-value for the pixel
	110	minx = pixel_width_x*(i - center_x)
	111	maxx = pixel_width_x*(i+1.0 - center_x)
	112
	113	qxmin = get_q(minx, 0.0, det_dist, wavelength)
	114	qxmax = get_q(maxx, 0.0, det_dist, wavelength)
	115
	116	# Get the count fraction in x for that pixel
	117	frac_min = get_pixel_fraction_square(self.x_min, qxmin, qxmax)
	118	frac_max = get_pixel_fraction_square(self.x_max, qxmin, qxmax)
	119	frac_x = frac_max - frac_min
	120
	121	if frac_x == 0:
	122	continue
	123
	124	if maj=='x':
	125	dx = pixel_width_x*(i+0.5 - center_x)
	126	q_value = get_q(dx, 0.0, det_dist, wavelength)
	127	if self.fold==False and dx<0:
	128	q_value = -q_value
	129	i_q = int(math.ceil((q_value-self.x_min)/self.bin_width)) - 1
	130
	131	if i_q<0 or i_q>=nbins:
	132	continue
	133
	134	for j in range(numpy.size(data2D.data,0)):
	135	# Min and max y-value for the pixel
	136	miny = pixel_width_y*(j - center_y)
	137	maxy = pixel_width_y*(j+1.0 - center_y)
	138
	139	qymin = get_q(0.0, miny, det_dist, wavelength)
	140	qymax = get_q(0.0, maxy, det_dist, wavelength)
	141
	142	# Get the count fraction in x for that pixel
	143	frac_min = get_pixel_fraction_square(self.y_min, qymin, qymax)
	144	frac_max = get_pixel_fraction_square(self.y_max, qymin, qymax)
	145	frac_y = frac_max - frac_min
	146
	147	frac = frac_x * frac_y
	148
	149	if frac == 0:
	150	continue
[76e2369]	151
[70975f3]	152	if maj=='y':
	153	dy = pixel_width_y*(j+0.5 - center_y)
	154	q_value = get_q(0.0, dy, det_dist, wavelength)
	155	if self.fold==False and dy<0:
	156	q_value = -q_value
	157	i_q = int(math.ceil((q_value-self.y_min)/self.bin_width)) - 1
	158
	159	if i_q<0 or i_q>=nbins:
	160	continue
	161
	162	x[i_q] = q_value
	163	y[i_q] += frac * data2D.data[j][i]
	164	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
	165	err_y[i_q] += frac * frac * math.fabs(data2D.data[j][i])
	166	else:
	167	err_y[i_q] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
	168	y_counts[i_q] += frac
[8ba103f]	169
[70975f3]	170	# Average the sums
	171	for i in range(nbins):
	172	if y_counts[i]>0:
	173	err_y[i] = math.sqrt(err_y[i])/y_counts[i]
	174	y[i] = y[i]/y_counts[i]
	175
	176	return Data1D(x=x, y=y, dy=err_y)
	177
	178	class SlabY(_Slab):
	179	"""
	180	Compute average I(Qy) for a region of interest
	181	"""
	182	def __call__(self, data2D):
	183	"""
	184	Compute average I(Qy) for a region of interest
	185
	186	@param data2D: Data2D object
	187	@return: Data1D object
	188	"""
	189	return self._avg(data2D, 'y')
	190
	191	class SlabX(_Slab):
	192	"""
	193	Compute average I(Qx) for a region of interest
	194	"""
	195	def __call__(self, data2D):
	196	"""
	197	Compute average I(Qx) for a region of interest
	198
	199	@param data2D: Data2D object
	200	@return: Data1D object
	201	"""
	202	return self._avg(data2D, 'x')
[f8d0ee7]	203
	204	class Boxsum(object):
	205	"""
	206	Perform the sum of counts in a 2D region of interest.
	207	"""
	208	def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0):
	209	# Minimum Qx value [A-1]
	210	self.x_min = x_min
	211	# Maximum Qx value [A-1]
	212	self.x_max = x_max
	213	# Minimum Qy value [A-1]
	214	self.y_min = y_min
	215	# Maximum Qy value [A-1]
	216	self.y_max = y_max
	217
	218	def __call__(self, data2D):
	219	"""
	220	Perform the sum in the region of interest
	221
	222	@param data2D: Data2D object
	223	@return: number of counts, error on number of counts
	224	"""
	225	y, err_y, y_counts = self._sum(data2D)
	226
	227	# Average the sums
	228	counts = 0 if y_counts==0 else y
	229	error = 0 if y_counts==0 else math.sqrt(err_y)
	230
	231	return counts, error
	232
	233	def _sum(self, data2D):
	234	"""
	235	Perform the sum in the region of interest
	236	@param data2D: Data2D object
	237	@return: number of counts, error on number of counts, number of entries summed
	238	"""
	239	if len(data2D.detector) != 1:
	240	raise RuntimeError, "Circular averaging: invalid number of detectors: %g" % len(data2D.detector)
	241
[2e83ff3]	242	pixel_width_x = data2D.detector[0].pixel_size.x
	243	pixel_width_y = data2D.detector[0].pixel_size.y
[f8d0ee7]	244	det_dist = data2D.detector[0].distance
	245	wavelength = data2D.source.wavelength
[2e83ff3]	246	center_x = data2D.detector[0].beam_center.x/pixel_width_x
	247	center_y = data2D.detector[0].beam_center.y/pixel_width_y
[f8d0ee7]	248
	249	y = 0.0
	250	err_y = 0.0
	251	y_counts = 0.0
[ca10d8e]	252	sign=1
[2e83ff3]	253	for i in range(numpy.size(data2D.data,1)):
[f8d0ee7]	254	# Min and max x-value for the pixel
[2e83ff3]	255	minx = pixel_width_x*(i - center_x)
	256	maxx = pixel_width_x*(i+1.0 - center_x)
[ca10d8e]	257	if minx>=0:
	258	sign=1
	259	else:
	260	sign=-1
	261	qxmin = sign*get_q(minx, 0.0, det_dist, wavelength)
	262	if maxx>=0:
	263	sign=1
	264	else:
	265	sign=-1
	266	qxmax = sign*get_q(maxx, 0.0, det_dist, wavelength)
[f8d0ee7]	267
	268	# Get the count fraction in x for that pixel
	269	frac_min = get_pixel_fraction_square(self.x_min, qxmin, qxmax)
	270	frac_max = get_pixel_fraction_square(self.x_max, qxmin, qxmax)
	271	frac_x = frac_max - frac_min
	272
[2e83ff3]	273	for j in range(numpy.size(data2D.data,0)):
[f8d0ee7]	274	# Min and max y-value for the pixel
[2e83ff3]	275	miny = pixel_width_y*(j - center_y)
	276	maxy = pixel_width_y*(j+1.0 - center_y)
[ca10d8e]	277	if miny>=0:
	278	sign=1
	279	else:
	280	sign=-1
[f8d0ee7]	281
[ca10d8e]	282	qymin = sign*get_q(0.0, miny, det_dist, wavelength)
	283	if maxy>=0:
	284	sign=1
	285	else:
	286	sign=-1
	287
	288	qymax = sign*get_q(0.0, maxy, det_dist, wavelength)
[f8d0ee7]	289
	290	# Get the count fraction in x for that pixel
	291	frac_min = get_pixel_fraction_square(self.y_min, qymin, qymax)
	292	frac_max = get_pixel_fraction_square(self.y_max, qymin, qymax)
	293	frac_y = frac_max - frac_min
	294
	295	frac = frac_x * frac_y
	296
	297	y += frac * data2D.data[j][i]
	298	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
	299	err_y += frac * frac * math.fabs(data2D.data[j][i])
	300	else:
	301	err_y += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
	302	y_counts += frac
	303
	304	return y, err_y, y_counts
	305
	306	class Boxavg(Boxsum):
	307	"""
	308	Perform the average of counts in a 2D region of interest.
	309	"""
	310	def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0):
	311	super(Boxavg, self).__init__(x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max)
	312
	313	def __call__(self, data2D):
	314	"""
	315	Perform the sum in the region of interest
	316
	317	@param data2D: Data2D object
	318	@return: average counts, error on average counts
	319	"""
	320	y, err_y, y_counts = self._sum(data2D)
	321
	322	# Average the sums
	323	counts = 0 if y_counts==0 else y/y_counts
	324	error = 0 if y_counts==0 else math.sqrt(err_y)/y_counts
	325
	326	return counts, error
	327
	328	def get_pixel_fraction_square(x, xmin, xmax):
	329	"""
	330	Return the fraction of the length
	331	from xmin to x.
	332
	333	A B
	334	+-----------+---------+
	335	xmin x xmax
	336
	337	@param x: x-value
	338	@param xmin: minimum x for the length considered
	339	@param xmax: minimum x for the length considered
	340	@return: (x-xmin)/(xmax-xmin) when xmin < x < xmax
	341
	342	"""
	343	if x<=xmin:
	344	return 0.0
	345	if x>xmin and x<xmax:
	346	return (x-xmin)/(xmax-xmin)
	347	else:
	348	return 1.0
	349
[76e2369]	350
	351	class CircularAverage(object):
	352	"""
	353	Perform circular averaging on 2D data
	354
	355	The data returned is the distribution of counts
	356	as a function of Q
	357	"""
	358	def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.001):
	359	# Minimum radius included in the average [A-1]
	360	self.r_min = r_min
	361	# Maximum radius included in the average [A-1]
	362	self.r_max = r_max
	363	# Bin width (step size) [A-1]
	364	self.bin_width = bin_width
	365
	366	def __call__(self, data2D):
	367	"""
	368	Perform circular averaging on the data
	369
	370	@param data2D: Data2D object
	371	@return: Data1D object
	372	"""
	373	if len(data2D.detector) != 1:
	374	raise RuntimeError, "Circular averaging: invalid number of detectors: %g" % len(data2D.detector)
	375
[2e83ff3]	376	pixel_width_x = data2D.detector[0].pixel_size.x
	377	pixel_width_y = data2D.detector[0].pixel_size.y
[76e2369]	378	det_dist = data2D.detector[0].distance
	379	wavelength = data2D.source.wavelength
[8ba103f]	380	center_x = data2D.detector[0].beam_center.x/pixel_width_x+0.5
	381	center_y = data2D.detector[0].beam_center.y/pixel_width_y+0.5
[76e2369]	382
	383	# Find out the maximum Q range
[8ba103f]	384	xwidth = (numpy.size(data2D.data,1))*pixel_width_x
[76e2369]	385	dx_max = xwidth - data2D.detector[0].beam_center.x
	386	if xwidth-dx_max>dx_max:
	387	dx_max = xwidth-dx_max
	388
[8ba103f]	389	ywidth = (numpy.size(data2D.data,0))*pixel_width_y
[76e2369]	390	dy_max = ywidth - data2D.detector[0].beam_center.y
	391	if ywidth-dy_max>dy_max:
	392	dy_max = ywidth-dy_max
	393
	394	qmax = get_q(dx_max, dy_max, det_dist, wavelength)
	395
	396	# Build array of Q intervals
	397	nbins = int(math.ceil((qmax-self.r_min)/self.bin_width))
	398	qbins = self.bin_width*numpy.arange(nbins)+self.r_min
	399
	400	x = numpy.zeros(nbins)
	401	y = numpy.zeros(nbins)
	402	err_y = numpy.zeros(nbins)
	403	y_counts = numpy.zeros(nbins)
	404
[2e83ff3]	405	for i in range(numpy.size(data2D.data,1)):
	406	dx = pixel_width_x*(i+0.5 - center_x)
[76e2369]	407
	408	# Min and max x-value for the pixel
[2e83ff3]	409	minx = pixel_width_x*(i - center_x)
	410	maxx = pixel_width_x*(i+1.0 - center_x)
[76e2369]	411
[2e83ff3]	412	for j in range(numpy.size(data2D.data,0)):
	413	dy = pixel_width_y*(j+0.5 - center_y)
[76e2369]	414
	415	q_value = get_q(dx, dy, det_dist, wavelength)
	416
	417	# Min and max y-value for the pixel
[2e83ff3]	418	miny = pixel_width_y*(j - center_y)
	419	maxy = pixel_width_y*(j+1.0 - center_y)
[76e2369]	420
	421	# Calculate the q-value for each corner
	422	# q_[x min or max][y min or max]
	423	q_00 = get_q(minx, miny, det_dist, wavelength)
	424	q_01 = get_q(minx, maxy, det_dist, wavelength)
	425	q_10 = get_q(maxx, miny, det_dist, wavelength)
	426	q_11 = get_q(maxx, maxy, det_dist, wavelength)
	427
	428	# Look for intercept between each side of the pixel
	429	# and the constant-q ring for qmax
	430	frac_max = get_pixel_fraction(self.r_max, q_00, q_01, q_10, q_11)
	431
	432	# Look for intercept between each side of the pixel
	433	# and the constant-q ring for qmin
	434	frac_min = get_pixel_fraction(self.r_min, q_00, q_01, q_10, q_11)
	435
	436	# We are interested in the region between qmin and qmax
	437	# therefore the fraction of the surface of the pixel
	438	# that we will use to calculate the number of counts to
	439	# include is given by:
	440	frac = frac_max - frac_min
	441
	442	i_q = int(math.ceil((q_value-self.r_min)/self.bin_width)) - 1
[8ba103f]	443	if q_value > qmax or q_value < self.r_min:
	444	continue
	445
[76e2369]	446	x[i_q] = q_value
	447	y[i_q] += frac * data2D.data[j][i]
	448	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
	449	err_y[i_q] += frac * frac * math.fabs(data2D.data[j][i])
	450	else:
	451	err_y[i_q] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
	452	y_counts[i_q] += frac
[8ba103f]	453
[76e2369]	454	# Average the sums
[8ba103f]	455	nzero = 0
[76e2369]	456	for i in range(nbins):
	457	if y_counts[i]>0:
	458	err_y[i] = math.sqrt(err_y[i])/y_counts[i]
	459	y[i] = y[i]/y_counts[i]
[8ba103f]	460	else:
	461	nzero += 1
	462	## Get rid of NULL points
	463	tx = numpy.zeros(nbins-nzero)
	464	ty = numpy.zeros(nbins-nzero)
	465	terr_y = numpy.zeros(nbins-nzero)
	466	j=0
	467	for i in range(nbins):
	468	if err_y[i] != 0 :
	469	tx[j] = x[i]
	470	ty[j] = y[i]
	471	terr_y[j] = err_y[i]
	472	j+=1
	473
	474	return Data1D(x=tx, y=ty, dy=terr_y)
[76e2369]	475
	476
	477	class Ring(object):
	478	"""
	479	Defines a ring on a 2D data set.
	480	The ring is defined by r_min, r_max, and
	481	the position of the center of the ring.
	482
	483	The data returned is the distribution of counts
	484	around the ring as a function of phi.
	485
	486	"""
	487
[bd89dea]	488	def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0,nbins=20 ):
[76e2369]	489	# Minimum radius
	490	self.r_min = r_min
	491	# Maximum radius
	492	self.r_max = r_max
	493	# Center of the ring in x
	494	self.center_x = center_x
	495	# Center of the ring in y
	496	self.center_y = center_y
	497	# Number of angular bins
[8ba103f]	498	self.nbins_phi = nbins
[76e2369]	499
	500	def __call__(self, data2D):
	501	"""
	502	Apply the ring to the data set.
	503	Returns the angular distribution for a given q range
	504
	505	@param data2D: Data2D object
	506	@return: Data1D object
	507	"""
	508	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
	509	raise RuntimeError, "Ring averaging only take plottable_2D objects"
	510
	511	data = data2D.data
	512	qmin = self.r_min
	513	qmax = self.r_max
	514
	515	if len(data2D.detector) != 1:
	516	raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector)
[2e83ff3]	517	pixel_width_x = data2D.detector[0].pixel_size.x
	518	pixel_width_y = data2D.detector[0].pixel_size.y
[76e2369]	519	det_dist = data2D.detector[0].distance
	520	wavelength = data2D.source.wavelength
[2e83ff3]	521	#center_x = self.center_x/pixel_width_x
	522	#center_y = self.center_y/pixel_width_y
	523	center_x = data2D.detector[0].beam_center.x/pixel_width_x
	524	center_y = data2D.detector[0].beam_center.y/pixel_width_y
	525
[76e2369]	526
	527	phi_bins = numpy.zeros(self.nbins_phi)
	528	phi_counts = numpy.zeros(self.nbins_phi)
	529	phi_values = numpy.zeros(self.nbins_phi)
	530	phi_err = numpy.zeros(self.nbins_phi)
	531
[2e83ff3]	532	for i in range(numpy.size(data,1)):
	533	dx = pixel_width_x*(i+0.5 - center_x)
[76e2369]	534
	535	# Min and max x-value for the pixel
[2e83ff3]	536	minx = pixel_width_x*(i - center_x)
	537	maxx = pixel_width_x*(i+1.0 - center_x)
[76e2369]	538
[2e83ff3]	539	for j in range(numpy.size(data,0)):
	540	dy = pixel_width_y*(j+0.5 - center_y)
[76e2369]	541
	542	q_value = get_q(dx, dy, det_dist, wavelength)
	543
	544	# Min and max y-value for the pixel
[2e83ff3]	545	miny = pixel_width_y*(j - center_y)
	546	maxy = pixel_width_y*(j+1.0 - center_y)
[76e2369]	547
	548	# Calculate the q-value for each corner
	549	# q_[x min or max][y min or max]
	550	q_00 = get_q(minx, miny, det_dist, wavelength)
	551	q_01 = get_q(minx, maxy, det_dist, wavelength)
	552	q_10 = get_q(maxx, miny, det_dist, wavelength)
	553	q_11 = get_q(maxx, maxy, det_dist, wavelength)
	554
	555	# Look for intercept between each side of the pixel
	556	# and the constant-q ring for qmax
	557	frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11)
	558
	559	# Look for intercept between each side of the pixel
	560	# and the constant-q ring for qmin
	561	frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11)
	562
	563	# We are interested in the region between qmin and qmax
	564	# therefore the fraction of the surface of the pixel
	565	# that we will use to calculate the number of counts to
	566	# include is given by:
	567
	568	frac = frac_max - frac_min
	569
[2e83ff3]	570	i_phi = int(math.ceil(self.nbins_phi(math.atan2(dy, dx)+math.pi)/(2.0math.pi))) - 1
[76e2369]	571
	572	phi_bins[i_phi] += frac * data[j][i]
	573
	574	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
	575	phi_err[i_phi] += frac * frac * math.fabs(data2D.data[j][i])
	576	else:
	577	phi_err[i_phi] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
	578	phi_counts[i_phi] += frac
	579
	580	for i in range(self.nbins_phi):
	581	phi_bins[i] = phi_bins[i] / phi_counts[i]
	582	phi_err[i] = math.sqrt(phi_err[i]) / phi_counts[i]
[7983c731]	583	phi_values[i] = 2.0math.pi/self.nbins_phi(1.0*i + 0.5)-math.pi # move the pi back to -pi <-->+pi
[76e2369]	584
	585	return Data1D(x=phi_values, y=phi_bins, dy=phi_err)
	586
	587	def get_pixel_fraction(qmax, q_00, q_01, q_10, q_11):
	588	"""
	589	Returns the fraction of the pixel defined by
	590	the four corners (q_00, q_01, q_10, q_11) that
	591	has q < qmax.
	592
	593	q_01 q_11
	594	y=1 +--------------+
	595	\| \|
	596	\| \|
	597	\| \|
	598	y=0 +--------------+
	599	q_00 q_01
	600
	601	x=0 x=1
	602
	603	"""
	604
	605	# y side for x = minx
	606	x_0 = get_intercept(qmax, q_00, q_01)
	607	# y side for x = maxx
	608	x_1 = get_intercept(qmax, q_10, q_11)
	609
	610	# x side for y = miny
	611	y_0 = get_intercept(qmax, q_00, q_10)
	612	# x side for y = maxy
	613	y_1 = get_intercept(qmax, q_01, q_11)
	614
	615	# surface fraction for a 1x1 pixel
	616	frac_max = 0
	617
	618	if x_0 and x_1:
	619	frac_max = (x_0+x_1)/2.0
	620
	621	elif y_0 and y_1:
	622	frac_max = (y_0+y_1)/2.0
	623
	624	elif x_0 and y_0:
	625	if q_00 < q_10:
	626	frac_max = x_0*y_0/2.0
	627	else:
	628	frac_max = 1.0-x_0*y_0/2.0
	629
	630	elif x_0 and y_1:
	631	if q_00 < q_10:
	632	frac_max = x_0*y_1/2.0
	633	else:
	634	frac_max = 1.0-x_0*y_1/2.0
	635
	636	elif x_1 and y_0:
	637	if q_00 > q_10:
	638	frac_max = x_1*y_0/2.0
	639	else:
	640	frac_max = 1.0-x_1*y_0/2.0
	641
	642	elif x_1 and y_1:
	643	if q_00 < q_10:
	644	frac_max = 1.0 - (1.0-x_1)*(1.0-y_1)/2.0
	645	else:
	646	frac_max = (1.0-x_1)*(1.0-y_1)/2.0
	647
	648	# If we make it here, there is no intercept between
	649	# this pixel and the constant-q ring. We only need
	650	# to know if we have to include it or exclude it.
	651	elif (q_00+q_01+q_10+q_11)/4.0 < qmax:
	652	frac_max = 1.0
	653
	654	return frac_max
	655
	656	def get_intercept(q, q_0, q_1):
	657	"""
	658	Returns the fraction of the side at which the
	659	q-value intercept the pixel, None otherwise.
	660	The values returned is the fraction ON THE SIDE
	661	OF THE LOWEST Q.
	662
	663
	664
	665	A B
	666	+-----------+--------+
	667	0 1 <--- pixel size
	668
	669	Q_0 -------- Q ----- Q_1 <--- equivalent Q range
	670
	671
	672	if Q_1 > Q_0, A is returned
	673	if Q_1 < Q_0, B is returned
	674
	675	if Q is outside the range of [Q_0, Q_1], None is returned
	676
	677	"""
	678	if q_1 > q_0:
	679	if (q > q_0 and q <= q_1):
	680	return (q-q_0)/(q_1 - q_0)
	681	else:
	682	if (q > q_1 and q <= q_0):
	683	return (q-q_1)/(q_0 - q_1)
	684	return None
	685
[c2a8523]	686	#This class can be removed.
[fb198a9]	687	class _Sectorold:
[76e2369]	688	"""
	689	Defines a sector region on a 2D data set.
	690	The sector is defined by r_min, r_max, phi_min, phi_max,
[fb198a9]	691	and the position of the center of the ring.
[2e83ff3]	692	Phi is defined between 0 and 2pi
	693	"""
[f2eee4a]	694	def __init__(self, r_min, r_max, phi_min, phi_max,nbins=20):
[2e83ff3]	695	self.r_min = r_min
	696	self.r_max = r_max
	697	self.phi_min = phi_min
	698	self.phi_max = phi_max
[f2eee4a]	699	self.nbins = nbins
[2e83ff3]	700
	701	def _agv(self, data2D, run='phi'):
	702	"""
	703	Perform sector averaging.
	704
	705	@param data2D: Data2D object
	706	@param run: define the varying parameter ('phi' or 'q')
	707	@return: Data1D object
	708	"""
	709	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
	710	raise RuntimeError, "Ring averaging only take plottable_2D objects"
[fb198a9]	711
	712	data = data2D.data
[2e83ff3]	713	qmax = self.r_max
[fb198a9]	714	qmin = self.r_min
[2e83ff3]	715
	716	if len(data2D.detector) != 1:
	717	raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector)
	718	pixel_width_x = data2D.detector[0].pixel_size.x
	719	pixel_width_y = data2D.detector[0].pixel_size.y
	720	det_dist = data2D.detector[0].distance
	721	wavelength = data2D.source.wavelength
	722	center_x = data2D.detector[0].beam_center.x/pixel_width_x
	723	center_y = data2D.detector[0].beam_center.y/pixel_width_y
	724
	725	y = numpy.zeros(self.nbins)
	726	y_counts = numpy.zeros(self.nbins)
	727	x = numpy.zeros(self.nbins)
	728	y_err = numpy.zeros(self.nbins)
	729
	730	for i in range(numpy.size(data,1)):
	731	dx = pixel_width_x*(i+0.5 - center_x)
	732
	733	# Min and max x-value for the pixel
	734	minx = pixel_width_x*(i - center_x)
	735	maxx = pixel_width_x*(i+1.0 - center_x)
	736
	737	for j in range(numpy.size(data,0)):
	738	dy = pixel_width_y*(j+0.5 - center_y)
	739
	740	q_value = get_q(dx, dy, det_dist, wavelength)
	741
	742	# Min and max y-value for the pixel
	743	miny = pixel_width_y*(j - center_y)
	744	maxy = pixel_width_y*(j+1.0 - center_y)
	745
	746	# Calculate the q-value for each corner
	747	# q_[x min or max][y min or max]
	748	q_00 = get_q(minx, miny, det_dist, wavelength)
	749	q_01 = get_q(minx, maxy, det_dist, wavelength)
	750	q_10 = get_q(maxx, miny, det_dist, wavelength)
	751	q_11 = get_q(maxx, maxy, det_dist, wavelength)
	752
	753	# Look for intercept between each side of the pixel
	754	# and the constant-q ring for qmax
	755	frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11)
	756
	757	# Look for intercept between each side of the pixel
	758	# and the constant-q ring for qmin
	759	frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11)
	760
	761	# We are interested in the region between qmin and qmax
	762	# therefore the fraction of the surface of the pixel
	763	# that we will use to calculate the number of counts to
	764	# include is given by:
	765
	766	frac = frac_max - frac_min
	767
	768	# Compute phi and check whether it's within the limits
[fb198a9]	769	phi_value=math.atan2(dy,dx)+math.pi
	770	# if phi_value<self.phi_min or phi_value>self.phi_max:
[2e83ff3]	771	if phi_value<self.phi_min or phi_value>self.phi_max:
	772	continue
	773
	774	# Check which type of averaging we need
	775	if run.lower()=='phi':
	776	i_bin = int(math.ceil(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) - 1
	777	else:
	778	# If we don't need this pixel, skip the rest of the work
	779	#TODO: an improvement here would be to compute the average
	780	# Q for the pixel from the part that is covered by
	781	# the ring defined by q_min/q_max rather than the complete
	782	# pixel
	783	if q_value<self.r_min or q_value>self.r_max:
	784	continue
	785	i_bin = int(math.ceil(self.nbins*(q_value-self.r_min)/(self.r_max-self.r_min))) - 1
	786
	787	try:
	788	y[i_bin] += frac * data[j][i]
	789	except:
	790	import sys
	791	print sys.exc_value
	792	print i_bin, frac
	793
	794	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
	795	y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i])
	796	else:
	797	y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
	798	y_counts[i_bin] += frac
	799
	800	for i in range(self.nbins):
	801	y[i] = y[i] / y_counts[i]
	802	y_err[i] = math.sqrt(y_err[i]) / y_counts[i]
	803	# Check which type of averaging we need
	804	if run.lower()=='phi':
	805	x[i] = (self.phi_max-self.phi_min)/self.nbins(1.0i + 0.5)+self.phi_min
	806	else:
	807	x[i] = (self.r_max-self.r_min)/self.nbins(1.0i + 0.5)+self.r_min
	808
	809	return Data1D(x=x, y=y, dy=y_err)
	810
[fb198a9]	811	class _Sector:
	812	"""
	813	Defines a sector region on a 2D data set.
	814	The sector is defined by r_min, r_max, phi_min, phi_max,
	815	and the position of the center of the ring
	816	where phi_min and phi_max are defined by the right and left lines wrt central line
	817	and phi_max could be less than phi_min.
	818
	819	Phi is defined between 0 and 2pi
	820	"""
	821	def __init__(self, r_min, r_max, phi_min, phi_max,nbins=20):
	822	self.r_min = r_min
	823	self.r_max = r_max
	824	self.phi_min = phi_min
	825	self.phi_max = phi_max
	826	self.nbins = nbins
	827
	828	def _agv(self, data2D, run='phi'):
	829	"""
	830	Perform sector averaging.
	831
	832	@param data2D: Data2D object
	833	@param run: define the varying parameter ('phi' or 'q')
	834	@return: Data1D object
	835	"""
	836	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
	837	raise RuntimeError, "Ring averaging only take plottable_2D objects"
	838
	839	data = data2D.data
	840	qmax = self.r_max
	841	qmin = self.r_min
	842
	843	if len(data2D.detector) != 1:
	844	raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector)
	845	pixel_width_x = data2D.detector[0].pixel_size.x
	846	pixel_width_y = data2D.detector[0].pixel_size.y
	847	det_dist = data2D.detector[0].distance
	848	wavelength = data2D.source.wavelength
	849	center_x = data2D.detector[0].beam_center.x/pixel_width_x
	850	center_y = data2D.detector[0].beam_center.y/pixel_width_y
	851
	852	y = numpy.zeros(self.nbins)
	853	y_counts = numpy.zeros(self.nbins)
	854	x = numpy.zeros(self.nbins)
	855	y_err = numpy.zeros(self.nbins)
	856
[c2a8523]	857	# This If finds qmax within ROI defined by sector lines
[4853c04]	858	temp=0 #to find qmax within ROI or phimax and phimin
	859	temp0=1000000
	860	temp1=0
	861	for i in range(numpy.size(data,1)):
	862	dx = pixel_width_x*(i+0.5 - center_x)
	863	for j in range(numpy.size(data,0)):
[c2a8523]	864
[4853c04]	865	dy = pixel_width_y*(j+0.5 - center_y)
	866	q_value = get_q(dx, dy, det_dist, wavelength)
	867	# Compute phi and check whether it's within the limits
	868	phi_value=math.atan2(dy,dx)+math.pi
	869	if self.phi_max>2*math.pi:
	870	self.phi_max=self.phi_max-2*math.pi
	871	if self.phi_min<0:
	872	self.phi_max=self.phi_max+2*math.pi
[c2a8523]	873
[8ba103f]	874	#In case of two ROI (symmetric major and minor wings)(for 'q2')
[4853c04]	875	if run.lower()=='q2':
	876	if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)):
	877	temp_max=self.phi_max+math.pi
	878	temp_min=self.phi_min+math.pi
	879	else:
	880	temp_max=self.phi_max
	881	temp_min=self.phi_min
[c2a8523]	882
[4853c04]	883	if ((temp_max>=math.pi and temp_max<2math.pi)and (temp_min>=math.pi and temp_min<2math.pi)):
	884	if (phi_value<temp_min or phi_value>temp_max):
	885	if (phi_value<temp_min-math.pi or phi_value>temp_max-math.pi):
	886	continue
	887	if (self.phi_max<self.phi_min):
	888	tmp_max=self.phi_max+math.pi
	889	tmp_min=self.phi_min-math.pi
	890	else:
	891	tmp_max=self.phi_max
	892	tmp_min=self.phi_min
	893	if (tmp_min<math.pi and tmp_max>math.pi):
	894	if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)):
	895	continue
[8d36420]	896	#In case of one ROI (major only)(i.e.,for 'q')and do nothing for 'phi'.
	897	elif run.lower()=='q':
[4853c04]	898	if (self.phi_max>=self.phi_min):
	899	if (phi_value<self.phi_min or phi_value>self.phi_max):
	900	continue
	901	else:
	902	if (phi_value<self.phi_min and phi_value>self.phi_max):
	903	continue
[d9629c53]	904	if q_value<qmin or q_value>qmax:
	905	continue
	906
[4853c04]	907	if run.lower()=='phi':
	908	if temp1<phi_value:
	909	temp1=phi_value
	910	if temp0>phi_value:
[d9629c53]	911	temp0=phi_value
	912
[4853c04]	913	elif temp<q_value:
	914	temp=q_value
[d9629c53]	915
[4853c04]	916	if run.lower()=='phi':
	917	self.phi_max=temp1
	918	self.phi_min=temp0
	919	else:
	920	qmax=temp
[8d36420]	921	#Beam center is already corrected, but the calculation below assumed it was not.
[8ba103f]	922	# Thus Beam center shifted back to uncorrected value. ToDo: cleanup the mess.
[ee3b2c0c]	923	center_x=center_x+0.5
	924	center_y=center_y+0.5
[fb198a9]	925	for i in range(numpy.size(data,1)):
	926	dx = pixel_width_x*(i+0.5 - center_x)
	927
	928	# Min and max x-value for the pixel
	929	minx = pixel_width_x*(i - center_x)
	930	maxx = pixel_width_x*(i+1.0 - center_x)
	931
	932	for j in range(numpy.size(data,0)):
	933	dy = pixel_width_y*(j+0.5 - center_y)
	934
	935	q_value = get_q(dx, dy, det_dist, wavelength)
	936
	937	# Min and max y-value for the pixel
	938	miny = pixel_width_y*(j - center_y)
	939	maxy = pixel_width_y*(j+1.0 - center_y)
	940
	941	# Calculate the q-value for each corner
	942	# q_[x min or max][y min or max]
	943	q_00 = get_q(minx, miny, det_dist, wavelength)
	944	q_01 = get_q(minx, maxy, det_dist, wavelength)
	945	q_10 = get_q(maxx, miny, det_dist, wavelength)
	946	q_11 = get_q(maxx, maxy, det_dist, wavelength)
	947
[c2a8523]	948	# Compute phi and check whether it's within the limits
	949	phi_value=math.atan2(dy,dx)+math.pi
	950	if self.phi_max>2*math.pi:
	951	self.phi_max=self.phi_max-2*math.pi
	952	if self.phi_min<0:
	953	self.phi_max=self.phi_max+2*math.pi
	954
[fb198a9]	955	# Look for intercept between each side of the pixel
	956	# and the constant-q ring for qmax
	957	frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11)
	958
	959	# Look for intercept between each side of the pixel
	960	# and the constant-q ring for qmin
	961	frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11)
	962
	963	# We are interested in the region between qmin and qmax
	964	# therefore the fraction of the surface of the pixel
	965	# that we will use to calculate the number of counts to
	966	# include is given by:
	967
	968	frac = frac_max - frac_min
	969
[c2a8523]	970	#In case of two ROI (symmetric major and minor regions)(for 'q2')
[fb198a9]	971	if run.lower()=='q2':
	972	if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)):
	973	temp_max=self.phi_max+math.pi
	974	temp_min=self.phi_min+math.pi
	975	else:
	976	temp_max=self.phi_max
	977	temp_min=self.phi_min
	978
	979	if ((temp_max>=math.pi and temp_max<2math.pi)and (temp_min>=math.pi and temp_min<2math.pi)):
	980	if (phi_value<temp_min or phi_value>temp_max):
	981	if (phi_value<temp_min-math.pi or phi_value>temp_max-math.pi):
	982	continue
	983	if (self.phi_max<self.phi_min):
	984	tmp_max=self.phi_max+math.pi
	985	tmp_min=self.phi_min-math.pi
	986	else:
	987	tmp_max=self.phi_max
	988	tmp_min=self.phi_min
	989	if (tmp_min<math.pi and tmp_max>math.pi):
	990	if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)):
	991	continue
[c2a8523]	992	#In case of one ROI (major only)(i.e.,for 'q' and 'phi')
[fb198a9]	993	else:
	994	if (self.phi_max>=self.phi_min):
	995	if (phi_value<self.phi_min or phi_value>self.phi_max):
	996	continue
[d9629c53]	997
[fb198a9]	998	else:
	999	if (phi_value<self.phi_min and phi_value>self.phi_max):
	1000	continue
[4853c04]	1001	#print "qmax=",qmax,qmin
	1002
	1003	if q_value<qmin or q_value>qmax:
	1004	continue
[fb198a9]	1005
	1006	# Check which type of averaging we need
	1007	if run.lower()=='phi':
	1008	i_bin = int(math.ceil(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) - 1
	1009	else:
	1010	# If we don't need this pixel, skip the rest of the work
	1011	#TODO: an improvement here would be to compute the average
	1012	# Q for the pixel from the part that is covered by
	1013	# the ring defined by q_min/q_max rather than the complete
	1014	# pixel
[c2a8523]	1015	i_bin = int(math.ceil(self.nbins*(q_value-qmin)/(qmax-qmin))) - 1
	1016
[fb198a9]	1017	try:
	1018	y[i_bin] += frac * data[j][i]
	1019	except:
	1020	import sys
	1021	print sys.exc_value
	1022	print i_bin, frac
	1023
	1024	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
	1025	y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i])
	1026	else:
	1027	y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
	1028	y_counts[i_bin] += frac
[2e83ff3]	1029
[fb198a9]	1030	for i in range(self.nbins):
	1031	y[i] = y[i] / y_counts[i]
	1032	y_err[i] = math.sqrt(y_err[i]) / y_counts[i]
	1033	# Check which type of averaging we need
	1034	if run.lower()=='phi':
[923d926]	1035	#Calculate x[i] and put back the origin of angle back to the right hand side (from the left).
[7983c731]	1036	x[i] = ((self.phi_max-self.phi_min)/self.nbins(1.0i + 0.5)+self.phi_min-2math.pi)180/math.pi
[923d926]	1037	if x[i]<0:
[77bf7c0]	1038	x[i]=180+x[i]
[fb198a9]	1039	else:
[c2a8523]	1040	x[i] = (qmax-qmin)/self.nbins(1.0i + 0.5)+qmin
[fb198a9]	1041
	1042	return Data1D(x=x, y=y, dy=y_err)
	1043
[2e83ff3]	1044	class SectorPhi(_Sector):
	1045	"""
	1046	Sector average as a function of phi.
	1047	I(phi) is return and the data is averaged over Q.
	1048
	1049	A sector is defined by r_min, r_max, phi_min, phi_max.
	1050	The number of bin in phi also has to be defined.
	1051	"""
	1052	def __call__(self, data2D):
	1053	"""
	1054	Perform sector average and return I(phi).
	1055
	1056	@param data2D: Data2D object
	1057	@return: Data1D object
	1058	"""
	1059	return self._agv(data2D, 'phi')
	1060
[fb198a9]	1061	class SectorQold(_Sector):
[2e83ff3]	1062	"""
	1063	Sector average as a function of Q.
	1064	I(Q) is return and the data is averaged over phi.
	1065
	1066	A sector is defined by r_min, r_max, phi_min, phi_max.
	1067	The number of bin in Q also has to be defined.
[76e2369]	1068	"""
[2e83ff3]	1069	def __call__(self, data2D):
	1070	"""
	1071	Perform sector average and return I(Q).
	1072
	1073	@param data2D: Data2D object
	1074	@return: Data1D object
	1075	"""
	1076	return self._agv(data2D, 'q')
[fb198a9]	1077
	1078	class SectorQ(_Sector):
	1079	"""
	1080	Sector average as a function of Q for both symatric wings.
	1081	I(Q) is return and the data is averaged over phi.
	1082
	1083	A sector is defined by r_min, r_max, phi_min, phi_max.
	1084	r_min, r_max, phi_min, phi_max >0.
	1085	The number of bin in Q also has to be defined.
	1086	"""
	1087	def __call__(self, data2D):
	1088	"""
	1089	Perform sector average and return I(Q).
	1090
	1091	@param data2D: Data2D object
	1092	@return: Data1D object
	1093	"""
	1094	return self._agv(data2D, 'q2')
[76e2369]	1095	if __name__ == "__main__":
	1096
	1097	from loader import Loader
	1098
	1099
[f8d0ee7]	1100	d = Loader().load('test/MAR07232_rest.ASC')
	1101	#d = Loader().load('test/MP_New.sans')
[76e2369]	1102
	1103
[d9629c53]	1104	r = SectorQ(r_min=.000001, r_max=.01, phi_min=0.0, phi_max=2*math.pi)
[f8d0ee7]	1105	o = r(d)
	1106
[d9629c53]	1107	s = Ring(r_min=.000001, r_max=.01)
[2e83ff3]	1108	p = s(d)
[70975f3]	1109
	1110	for i in range(len(o.x)):
	1111	print o.x[i], o.y[i], o.dy[i]
[76e2369]	1112
	1113
	1114

Note: See TracBrowser for help on using the repository browser.

SasView

source: sasview/DataLoader/manipulations.py @ 3ba03a4c

Download in other formats: