manipulations.py @ 39a365d

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 39a365d was 7983c731, checked in by Jae Cho <jhjcho@…>, 16 years ago
reodered Annulus average data so that the line w/o error bar is not broken.
Property mode set to `100644`
File size: 42.9 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)
	107
	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
[76e2369]	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
[2e83ff3]	380	center_x = data2D.detector[0].beam_center.x/pixel_width_x
	381	center_y = data2D.detector[0].beam_center.y/pixel_width_y
[76e2369]	382
	383	# Find out the maximum Q range
[2e83ff3]	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
[2e83ff3]	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
	443
	444	x[i_q] = q_value
	445	y[i_q] += frac * data2D.data[j][i]
	446	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
	447	err_y[i_q] += frac * frac * math.fabs(data2D.data[j][i])
	448	else:
	449	err_y[i_q] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
	450	y_counts[i_q] += frac
	451
	452	# Average the sums
	453	for i in range(nbins):
	454	if y_counts[i]>0:
	455	err_y[i] = math.sqrt(err_y[i])/y_counts[i]
	456	y[i] = y[i]/y_counts[i]
	457
	458	return Data1D(x=x, y=y, dy=err_y)
	459
	460
	461	class Ring(object):
	462	"""
	463	Defines a ring on a 2D data set.
	464	The ring is defined by r_min, r_max, and
	465	the position of the center of the ring.
	466
	467	The data returned is the distribution of counts
	468	around the ring as a function of phi.
	469
	470	"""
	471
	472	def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0):
	473	# Minimum radius
	474	self.r_min = r_min
	475	# Maximum radius
	476	self.r_max = r_max
	477	# Center of the ring in x
	478	self.center_x = center_x
	479	# Center of the ring in y
	480	self.center_y = center_y
	481	# Number of angular bins
	482	self.nbins_phi = 20
	483
	484	def __call__(self, data2D):
	485	"""
	486	Apply the ring to the data set.
	487	Returns the angular distribution for a given q range
	488
	489	@param data2D: Data2D object
	490	@return: Data1D object
	491	"""
	492	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
	493	raise RuntimeError, "Ring averaging only take plottable_2D objects"
	494
	495	data = data2D.data
	496	qmin = self.r_min
	497	qmax = self.r_max
	498
	499	if len(data2D.detector) != 1:
	500	raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector)
[2e83ff3]	501	pixel_width_x = data2D.detector[0].pixel_size.x
	502	pixel_width_y = data2D.detector[0].pixel_size.y
[76e2369]	503	det_dist = data2D.detector[0].distance
	504	wavelength = data2D.source.wavelength
[2e83ff3]	505	#center_x = self.center_x/pixel_width_x
	506	#center_y = self.center_y/pixel_width_y
	507	center_x = data2D.detector[0].beam_center.x/pixel_width_x
	508	center_y = data2D.detector[0].beam_center.y/pixel_width_y
	509
[76e2369]	510
	511	phi_bins = numpy.zeros(self.nbins_phi)
	512	phi_counts = numpy.zeros(self.nbins_phi)
	513	phi_values = numpy.zeros(self.nbins_phi)
	514	phi_err = numpy.zeros(self.nbins_phi)
	515
[2e83ff3]	516	for i in range(numpy.size(data,1)):
	517	dx = pixel_width_x*(i+0.5 - center_x)
[76e2369]	518
	519	# Min and max x-value for the pixel
[2e83ff3]	520	minx = pixel_width_x*(i - center_x)
	521	maxx = pixel_width_x*(i+1.0 - center_x)
[76e2369]	522
[2e83ff3]	523	for j in range(numpy.size(data,0)):
	524	dy = pixel_width_y*(j+0.5 - center_y)
[76e2369]	525
	526	q_value = get_q(dx, dy, det_dist, wavelength)
	527
	528	# Min and max y-value for the pixel
[2e83ff3]	529	miny = pixel_width_y*(j - center_y)
	530	maxy = pixel_width_y*(j+1.0 - center_y)
[76e2369]	531
	532	# Calculate the q-value for each corner
	533	# q_[x min or max][y min or max]
	534	q_00 = get_q(minx, miny, det_dist, wavelength)
	535	q_01 = get_q(minx, maxy, det_dist, wavelength)
	536	q_10 = get_q(maxx, miny, det_dist, wavelength)
	537	q_11 = get_q(maxx, maxy, det_dist, wavelength)
	538
	539	# Look for intercept between each side of the pixel
	540	# and the constant-q ring for qmax
	541	frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11)
	542
	543	# Look for intercept between each side of the pixel
	544	# and the constant-q ring for qmin
	545	frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11)
	546
	547	# We are interested in the region between qmin and qmax
	548	# therefore the fraction of the surface of the pixel
	549	# that we will use to calculate the number of counts to
	550	# include is given by:
	551
	552	frac = frac_max - frac_min
	553
[2e83ff3]	554	i_phi = int(math.ceil(self.nbins_phi(math.atan2(dy, dx)+math.pi)/(2.0math.pi))) - 1
[76e2369]	555
	556	phi_bins[i_phi] += frac * data[j][i]
	557
	558	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
	559	phi_err[i_phi] += frac * frac * math.fabs(data2D.data[j][i])
	560	else:
	561	phi_err[i_phi] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
	562	phi_counts[i_phi] += frac
	563
	564	for i in range(self.nbins_phi):
	565	phi_bins[i] = phi_bins[i] / phi_counts[i]
	566	phi_err[i] = math.sqrt(phi_err[i]) / phi_counts[i]
[7983c731]	567	phi_values[i] = 2.0math.pi/self.nbins_phi(1.0*i + 0.5)-math.pi # move the pi back to -pi <-->+pi
[76e2369]	568
	569	return Data1D(x=phi_values, y=phi_bins, dy=phi_err)
	570
	571	def get_pixel_fraction(qmax, q_00, q_01, q_10, q_11):
	572	"""
	573	Returns the fraction of the pixel defined by
	574	the four corners (q_00, q_01, q_10, q_11) that
	575	has q < qmax.
	576
	577	q_01 q_11
	578	y=1 +--------------+
	579	\| \|
	580	\| \|
	581	\| \|
	582	y=0 +--------------+
	583	q_00 q_01
	584
	585	x=0 x=1
	586
	587	"""
	588
	589	# y side for x = minx
	590	x_0 = get_intercept(qmax, q_00, q_01)
	591	# y side for x = maxx
	592	x_1 = get_intercept(qmax, q_10, q_11)
	593
	594	# x side for y = miny
	595	y_0 = get_intercept(qmax, q_00, q_10)
	596	# x side for y = maxy
	597	y_1 = get_intercept(qmax, q_01, q_11)
	598
	599	# surface fraction for a 1x1 pixel
	600	frac_max = 0
	601
	602	if x_0 and x_1:
	603	frac_max = (x_0+x_1)/2.0
	604
	605	elif y_0 and y_1:
	606	frac_max = (y_0+y_1)/2.0
	607
	608	elif x_0 and y_0:
	609	if q_00 < q_10:
	610	frac_max = x_0*y_0/2.0
	611	else:
	612	frac_max = 1.0-x_0*y_0/2.0
	613
	614	elif x_0 and y_1:
	615	if q_00 < q_10:
	616	frac_max = x_0*y_1/2.0
	617	else:
	618	frac_max = 1.0-x_0*y_1/2.0
	619
	620	elif x_1 and y_0:
	621	if q_00 > q_10:
	622	frac_max = x_1*y_0/2.0
	623	else:
	624	frac_max = 1.0-x_1*y_0/2.0
	625
	626	elif x_1 and y_1:
	627	if q_00 < q_10:
	628	frac_max = 1.0 - (1.0-x_1)*(1.0-y_1)/2.0
	629	else:
	630	frac_max = (1.0-x_1)*(1.0-y_1)/2.0
	631
	632	# If we make it here, there is no intercept between
	633	# this pixel and the constant-q ring. We only need
	634	# to know if we have to include it or exclude it.
	635	elif (q_00+q_01+q_10+q_11)/4.0 < qmax:
	636	frac_max = 1.0
	637
	638	return frac_max
	639
	640	def get_intercept(q, q_0, q_1):
	641	"""
	642	Returns the fraction of the side at which the
	643	q-value intercept the pixel, None otherwise.
	644	The values returned is the fraction ON THE SIDE
	645	OF THE LOWEST Q.
	646
	647
	648
	649	A B
	650	+-----------+--------+
	651	0 1 <--- pixel size
	652
	653	Q_0 -------- Q ----- Q_1 <--- equivalent Q range
	654
	655
	656	if Q_1 > Q_0, A is returned
	657	if Q_1 < Q_0, B is returned
	658
	659	if Q is outside the range of [Q_0, Q_1], None is returned
	660
	661	"""
	662	if q_1 > q_0:
	663	if (q > q_0 and q <= q_1):
	664	return (q-q_0)/(q_1 - q_0)
	665	else:
	666	if (q > q_1 and q <= q_0):
	667	return (q-q_1)/(q_0 - q_1)
	668	return None
	669
[c2a8523]	670	#This class can be removed.
[fb198a9]	671	class _Sectorold:
[76e2369]	672	"""
	673	Defines a sector region on a 2D data set.
	674	The sector is defined by r_min, r_max, phi_min, phi_max,
[fb198a9]	675	and the position of the center of the ring.
[2e83ff3]	676	Phi is defined between 0 and 2pi
	677	"""
[f2eee4a]	678	def __init__(self, r_min, r_max, phi_min, phi_max,nbins=20):
[2e83ff3]	679	self.r_min = r_min
	680	self.r_max = r_max
	681	self.phi_min = phi_min
	682	self.phi_max = phi_max
[f2eee4a]	683	self.nbins = nbins
[2e83ff3]	684
	685	def _agv(self, data2D, run='phi'):
	686	"""
	687	Perform sector averaging.
	688
	689	@param data2D: Data2D object
	690	@param run: define the varying parameter ('phi' or 'q')
	691	@return: Data1D object
	692	"""
	693	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
	694	raise RuntimeError, "Ring averaging only take plottable_2D objects"
[fb198a9]	695
	696	data = data2D.data
[2e83ff3]	697	qmax = self.r_max
[fb198a9]	698	qmin = self.r_min
[2e83ff3]	699
	700	if len(data2D.detector) != 1:
	701	raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector)
	702	pixel_width_x = data2D.detector[0].pixel_size.x
	703	pixel_width_y = data2D.detector[0].pixel_size.y
	704	det_dist = data2D.detector[0].distance
	705	wavelength = data2D.source.wavelength
	706	center_x = data2D.detector[0].beam_center.x/pixel_width_x
	707	center_y = data2D.detector[0].beam_center.y/pixel_width_y
	708
	709	y = numpy.zeros(self.nbins)
	710	y_counts = numpy.zeros(self.nbins)
	711	x = numpy.zeros(self.nbins)
	712	y_err = numpy.zeros(self.nbins)
	713
	714	for i in range(numpy.size(data,1)):
	715	dx = pixel_width_x*(i+0.5 - center_x)
	716
	717	# Min and max x-value for the pixel
	718	minx = pixel_width_x*(i - center_x)
	719	maxx = pixel_width_x*(i+1.0 - center_x)
	720
	721	for j in range(numpy.size(data,0)):
	722	dy = pixel_width_y*(j+0.5 - center_y)
	723
	724	q_value = get_q(dx, dy, det_dist, wavelength)
	725
	726	# Min and max y-value for the pixel
	727	miny = pixel_width_y*(j - center_y)
	728	maxy = pixel_width_y*(j+1.0 - center_y)
	729
	730	# Calculate the q-value for each corner
	731	# q_[x min or max][y min or max]
	732	q_00 = get_q(minx, miny, det_dist, wavelength)
	733	q_01 = get_q(minx, maxy, det_dist, wavelength)
	734	q_10 = get_q(maxx, miny, det_dist, wavelength)
	735	q_11 = get_q(maxx, maxy, det_dist, wavelength)
	736
	737	# Look for intercept between each side of the pixel
	738	# and the constant-q ring for qmax
	739	frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11)
	740
	741	# Look for intercept between each side of the pixel
	742	# and the constant-q ring for qmin
	743	frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11)
	744
	745	# We are interested in the region between qmin and qmax
	746	# therefore the fraction of the surface of the pixel
	747	# that we will use to calculate the number of counts to
	748	# include is given by:
	749
	750	frac = frac_max - frac_min
	751
	752	# Compute phi and check whether it's within the limits
[fb198a9]	753	phi_value=math.atan2(dy,dx)+math.pi
	754	# if phi_value<self.phi_min or phi_value>self.phi_max:
[2e83ff3]	755	if phi_value<self.phi_min or phi_value>self.phi_max:
	756	continue
	757
	758	# Check which type of averaging we need
	759	if run.lower()=='phi':
	760	i_bin = int(math.ceil(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) - 1
	761	else:
	762	# If we don't need this pixel, skip the rest of the work
	763	#TODO: an improvement here would be to compute the average
	764	# Q for the pixel from the part that is covered by
	765	# the ring defined by q_min/q_max rather than the complete
	766	# pixel
	767	if q_value<self.r_min or q_value>self.r_max:
	768	continue
	769	i_bin = int(math.ceil(self.nbins*(q_value-self.r_min)/(self.r_max-self.r_min))) - 1
	770
	771	try:
	772	y[i_bin] += frac * data[j][i]
	773	except:
	774	import sys
	775	print sys.exc_value
	776	print i_bin, frac
	777
	778	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
	779	y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i])
	780	else:
	781	y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
	782	y_counts[i_bin] += frac
	783
	784	for i in range(self.nbins):
	785	y[i] = y[i] / y_counts[i]
	786	y_err[i] = math.sqrt(y_err[i]) / y_counts[i]
	787	# Check which type of averaging we need
	788	if run.lower()=='phi':
	789	x[i] = (self.phi_max-self.phi_min)/self.nbins(1.0i + 0.5)+self.phi_min
	790	else:
	791	x[i] = (self.r_max-self.r_min)/self.nbins(1.0i + 0.5)+self.r_min
	792
	793	return Data1D(x=x, y=y, dy=y_err)
	794
[fb198a9]	795	class _Sector:
	796	"""
	797	Defines a sector region on a 2D data set.
	798	The sector is defined by r_min, r_max, phi_min, phi_max,
	799	and the position of the center of the ring
	800	where phi_min and phi_max are defined by the right and left lines wrt central line
	801	and phi_max could be less than phi_min.
	802
	803	Phi is defined between 0 and 2pi
	804	"""
	805	def __init__(self, r_min, r_max, phi_min, phi_max,nbins=20):
	806	self.r_min = r_min
	807	self.r_max = r_max
	808	self.phi_min = phi_min
	809	self.phi_max = phi_max
	810	self.nbins = nbins
	811
	812	def _agv(self, data2D, run='phi'):
	813	"""
	814	Perform sector averaging.
	815
	816	@param data2D: Data2D object
	817	@param run: define the varying parameter ('phi' or 'q')
	818	@return: Data1D object
	819	"""
	820	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
	821	raise RuntimeError, "Ring averaging only take plottable_2D objects"
	822
	823	data = data2D.data
	824	qmax = self.r_max
	825	qmin = self.r_min
	826
	827	if len(data2D.detector) != 1:
	828	raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector)
	829	pixel_width_x = data2D.detector[0].pixel_size.x
	830	pixel_width_y = data2D.detector[0].pixel_size.y
	831	det_dist = data2D.detector[0].distance
	832	wavelength = data2D.source.wavelength
	833	center_x = data2D.detector[0].beam_center.x/pixel_width_x
	834	center_y = data2D.detector[0].beam_center.y/pixel_width_y
	835
	836	y = numpy.zeros(self.nbins)
	837	y_counts = numpy.zeros(self.nbins)
	838	x = numpy.zeros(self.nbins)
	839	y_err = numpy.zeros(self.nbins)
	840
[c2a8523]	841	# This If finds qmax within ROI defined by sector lines
[4853c04]	842	temp=0 #to find qmax within ROI or phimax and phimin
	843	temp0=1000000
	844	temp1=0
	845	for i in range(numpy.size(data,1)):
	846	dx = pixel_width_x*(i+0.5 - center_x)
	847	for j in range(numpy.size(data,0)):
[c2a8523]	848
[4853c04]	849	dy = pixel_width_y*(j+0.5 - center_y)
	850	q_value = get_q(dx, dy, det_dist, wavelength)
	851	# Compute phi and check whether it's within the limits
	852	phi_value=math.atan2(dy,dx)+math.pi
	853	if self.phi_max>2*math.pi:
	854	self.phi_max=self.phi_max-2*math.pi
	855	if self.phi_min<0:
	856	self.phi_max=self.phi_max+2*math.pi
[c2a8523]	857
[4853c04]	858	#In case of two ROI (symmetric major and minor regions)(for 'q2')
	859	if run.lower()=='q2':
	860	if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)):
	861	temp_max=self.phi_max+math.pi
	862	temp_min=self.phi_min+math.pi
	863	else:
	864	temp_max=self.phi_max
	865	temp_min=self.phi_min
[c2a8523]	866
[4853c04]	867	if ((temp_max>=math.pi and temp_max<2math.pi)and (temp_min>=math.pi and temp_min<2math.pi)):
	868	if (phi_value<temp_min or phi_value>temp_max):
	869	if (phi_value<temp_min-math.pi or phi_value>temp_max-math.pi):
	870	continue
	871	if (self.phi_max<self.phi_min):
	872	tmp_max=self.phi_max+math.pi
	873	tmp_min=self.phi_min-math.pi
	874	else:
	875	tmp_max=self.phi_max
	876	tmp_min=self.phi_min
	877	if (tmp_min<math.pi and tmp_max>math.pi):
	878	if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)):
	879	continue
	880	#In case of one ROI (major only)(i.e.,for 'q' and 'phi')
	881	else:
	882	if (self.phi_max>=self.phi_min):
	883	if (phi_value<self.phi_min or phi_value>self.phi_max):
	884	continue
	885	else:
	886	if (phi_value<self.phi_min and phi_value>self.phi_max):
	887	continue
[d9629c53]	888	if q_value<qmin or q_value>qmax:
	889	continue
	890
[4853c04]	891	if run.lower()=='phi':
	892	if temp1<phi_value:
	893	temp1=phi_value
	894	if temp0>phi_value:
[d9629c53]	895	temp0=phi_value
	896
[4853c04]	897	elif temp<q_value:
	898	temp=q_value
[d9629c53]	899
[4853c04]	900	if run.lower()=='phi':
	901	self.phi_max=temp1
	902	self.phi_min=temp0
	903	else:
	904	qmax=temp
	905
[fb198a9]	906	for i in range(numpy.size(data,1)):
	907	dx = pixel_width_x*(i+0.5 - center_x)
	908
	909	# Min and max x-value for the pixel
	910	minx = pixel_width_x*(i - center_x)
	911	maxx = pixel_width_x*(i+1.0 - center_x)
	912
	913	for j in range(numpy.size(data,0)):
	914	dy = pixel_width_y*(j+0.5 - center_y)
	915
	916	q_value = get_q(dx, dy, det_dist, wavelength)
	917
	918	# Min and max y-value for the pixel
	919	miny = pixel_width_y*(j - center_y)
	920	maxy = pixel_width_y*(j+1.0 - center_y)
	921
	922	# Calculate the q-value for each corner
	923	# q_[x min or max][y min or max]
	924	q_00 = get_q(minx, miny, det_dist, wavelength)
	925	q_01 = get_q(minx, maxy, det_dist, wavelength)
	926	q_10 = get_q(maxx, miny, det_dist, wavelength)
	927	q_11 = get_q(maxx, maxy, det_dist, wavelength)
	928
[c2a8523]	929	# Compute phi and check whether it's within the limits
	930	phi_value=math.atan2(dy,dx)+math.pi
	931	if self.phi_max>2*math.pi:
	932	self.phi_max=self.phi_max-2*math.pi
	933	if self.phi_min<0:
	934	self.phi_max=self.phi_max+2*math.pi
	935
[fb198a9]	936	# Look for intercept between each side of the pixel
	937	# and the constant-q ring for qmax
	938	frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11)
	939
	940	# Look for intercept between each side of the pixel
	941	# and the constant-q ring for qmin
	942	frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11)
	943
	944	# We are interested in the region between qmin and qmax
	945	# therefore the fraction of the surface of the pixel
	946	# that we will use to calculate the number of counts to
	947	# include is given by:
	948
	949	frac = frac_max - frac_min
	950
[c2a8523]	951	#In case of two ROI (symmetric major and minor regions)(for 'q2')
[fb198a9]	952	if run.lower()=='q2':
	953	if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)):
	954	temp_max=self.phi_max+math.pi
	955	temp_min=self.phi_min+math.pi
	956	else:
	957	temp_max=self.phi_max
	958	temp_min=self.phi_min
	959
	960	if ((temp_max>=math.pi and temp_max<2math.pi)and (temp_min>=math.pi and temp_min<2math.pi)):
	961	if (phi_value<temp_min or phi_value>temp_max):
	962	if (phi_value<temp_min-math.pi or phi_value>temp_max-math.pi):
	963	continue
	964	if (self.phi_max<self.phi_min):
	965	tmp_max=self.phi_max+math.pi
	966	tmp_min=self.phi_min-math.pi
	967	else:
	968	tmp_max=self.phi_max
	969	tmp_min=self.phi_min
	970	if (tmp_min<math.pi and tmp_max>math.pi):
	971	if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)):
	972	continue
[c2a8523]	973	#In case of one ROI (major only)(i.e.,for 'q' and 'phi')
[fb198a9]	974	else:
	975	if (self.phi_max>=self.phi_min):
	976	if (phi_value<self.phi_min or phi_value>self.phi_max):
	977	continue
[d9629c53]	978
[fb198a9]	979	else:
	980	if (phi_value<self.phi_min and phi_value>self.phi_max):
	981	continue
[4853c04]	982	#print "qmax=",qmax,qmin
	983
	984	if q_value<qmin or q_value>qmax:
	985	continue
[fb198a9]	986
	987	# Check which type of averaging we need
	988	if run.lower()=='phi':
	989	i_bin = int(math.ceil(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) - 1
	990	else:
	991	# If we don't need this pixel, skip the rest of the work
	992	#TODO: an improvement here would be to compute the average
	993	# Q for the pixel from the part that is covered by
	994	# the ring defined by q_min/q_max rather than the complete
	995	# pixel
[c2a8523]	996	i_bin = int(math.ceil(self.nbins*(q_value-qmin)/(qmax-qmin))) - 1
	997
[fb198a9]	998	try:
	999	y[i_bin] += frac * data[j][i]
	1000	except:
	1001	import sys
	1002	print sys.exc_value
	1003	print i_bin, frac
	1004
	1005	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
	1006	y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i])
	1007	else:
	1008	y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
	1009	y_counts[i_bin] += frac
[2e83ff3]	1010
[fb198a9]	1011	for i in range(self.nbins):
	1012	y[i] = y[i] / y_counts[i]
	1013	y_err[i] = math.sqrt(y_err[i]) / y_counts[i]
	1014	# Check which type of averaging we need
	1015	if run.lower()=='phi':
[923d926]	1016	#Calculate x[i] and put back the origin of angle back to the right hand side (from the left).
[7983c731]	1017	x[i] = ((self.phi_max-self.phi_min)/self.nbins(1.0i + 0.5)+self.phi_min-2math.pi)180/math.pi
[923d926]	1018	if x[i]<0:
	1019	x[i]=360+x[i]
[fb198a9]	1020	else:
[c2a8523]	1021	x[i] = (qmax-qmin)/self.nbins(1.0i + 0.5)+qmin
[fb198a9]	1022
	1023	return Data1D(x=x, y=y, dy=y_err)
	1024
[2e83ff3]	1025	class SectorPhi(_Sector):
	1026	"""
	1027	Sector average as a function of phi.
	1028	I(phi) is return and the data is averaged over Q.
	1029
	1030	A sector is defined by r_min, r_max, phi_min, phi_max.
	1031	The number of bin in phi also has to be defined.
	1032	"""
	1033	def __call__(self, data2D):
	1034	"""
	1035	Perform sector average and return I(phi).
	1036
	1037	@param data2D: Data2D object
	1038	@return: Data1D object
	1039	"""
	1040	return self._agv(data2D, 'phi')
	1041
[fb198a9]	1042	class SectorQold(_Sector):
[2e83ff3]	1043	"""
	1044	Sector average as a function of Q.
	1045	I(Q) is return and the data is averaged over phi.
	1046
	1047	A sector is defined by r_min, r_max, phi_min, phi_max.
	1048	The number of bin in Q also has to be defined.
[76e2369]	1049	"""
[2e83ff3]	1050	def __call__(self, data2D):
	1051	"""
	1052	Perform sector average and return I(Q).
	1053
	1054	@param data2D: Data2D object
	1055	@return: Data1D object
	1056	"""
	1057	return self._agv(data2D, 'q')
[fb198a9]	1058
	1059	class SectorQ(_Sector):
	1060	"""
	1061	Sector average as a function of Q for both symatric wings.
	1062	I(Q) is return and the data is averaged over phi.
	1063
	1064	A sector is defined by r_min, r_max, phi_min, phi_max.
	1065	r_min, r_max, phi_min, phi_max >0.
	1066	The number of bin in Q also has to be defined.
	1067	"""
	1068	def __call__(self, data2D):
	1069	"""
	1070	Perform sector average and return I(Q).
	1071
	1072	@param data2D: Data2D object
	1073	@return: Data1D object
	1074	"""
	1075	return self._agv(data2D, 'q2')
[76e2369]	1076	if __name__ == "__main__":
	1077
	1078	from loader import Loader
	1079
	1080
[f8d0ee7]	1081	d = Loader().load('test/MAR07232_rest.ASC')
	1082	#d = Loader().load('test/MP_New.sans')
[76e2369]	1083
	1084
[d9629c53]	1085	r = SectorQ(r_min=.000001, r_max=.01, phi_min=0.0, phi_max=2*math.pi)
[f8d0ee7]	1086	o = r(d)
	1087
[d9629c53]	1088	s = Ring(r_min=.000001, r_max=.01)
[2e83ff3]	1089	p = s(d)
[70975f3]	1090
	1091	for i in range(len(o.x)):
	1092	print o.x[i], o.y[i], o.dy[i]
[76e2369]	1093
	1094
	1095

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

SasView

source: sasview/DataLoader/manipulations.py @ 39a365d

Download in other formats: