manipulations.py @ d1dd9d4

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

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

SasView

source: sasview/DataLoader/manipulations.py @ d1dd9d4

Download in other formats: