manipulations.py @ ca10d8e

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 ca10d8e was ca10d8e, checked in by Gervaise Alina <gervyh@…>, 16 years ago
reverse code from version 1237 with modification on output axis unit
Property mode set to `100644`
File size: 42.3 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
[ca10d8e]	233	sign=1
[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)
[ca10d8e]	238	if minx>=0:
	239	sign=1
	240	else:
	241	sign=-1
	242	qxmin = sign*get_q(minx, 0.0, det_dist, wavelength)
	243	if maxx>=0:
	244	sign=1
	245	else:
	246	sign=-1
	247	qxmax = sign*get_q(maxx, 0.0, det_dist, wavelength)
[f8d0ee7]	248
	249	# Get the count fraction in x for that pixel
	250	frac_min = get_pixel_fraction_square(self.x_min, qxmin, qxmax)
	251	frac_max = get_pixel_fraction_square(self.x_max, qxmin, qxmax)
	252	frac_x = frac_max - frac_min
	253
[2e83ff3]	254	for j in range(numpy.size(data2D.data,0)):
[f8d0ee7]	255	# Min and max y-value for the pixel
[2e83ff3]	256	miny = pixel_width_y*(j - center_y)
	257	maxy = pixel_width_y*(j+1.0 - center_y)
[ca10d8e]	258	if miny>=0:
	259	sign=1
	260	else:
	261	sign=-1
[f8d0ee7]	262
[ca10d8e]	263	qymin = sign*get_q(0.0, miny, det_dist, wavelength)
	264	if maxy>=0:
	265	sign=1
	266	else:
	267	sign=-1
	268
	269	qymax = sign*get_q(0.0, maxy, det_dist, wavelength)
[f8d0ee7]	270
	271	# Get the count fraction in x for that pixel
	272	frac_min = get_pixel_fraction_square(self.y_min, qymin, qymax)
	273	frac_max = get_pixel_fraction_square(self.y_max, qymin, qymax)
	274	frac_y = frac_max - frac_min
	275
	276	frac = frac_x * frac_y
	277
	278	y += frac * data2D.data[j][i]
	279	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
	280	err_y += frac * frac * math.fabs(data2D.data[j][i])
	281	else:
	282	err_y += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
	283	y_counts += frac
	284
	285	return y, err_y, y_counts
	286
	287	class Boxavg(Boxsum):
	288	"""
	289	Perform the average of counts in a 2D region of interest.
	290	"""
	291	def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0):
	292	super(Boxavg, self).__init__(x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max)
	293
	294	def __call__(self, data2D):
	295	"""
	296	Perform the sum in the region of interest
	297
	298	@param data2D: Data2D object
	299	@return: average counts, error on average counts
	300	"""
	301	y, err_y, y_counts = self._sum(data2D)
	302
	303	# Average the sums
	304	counts = 0 if y_counts==0 else y/y_counts
	305	error = 0 if y_counts==0 else math.sqrt(err_y)/y_counts
	306
	307	return counts, error
	308
	309	def get_pixel_fraction_square(x, xmin, xmax):
	310	"""
	311	Return the fraction of the length
	312	from xmin to x.
	313
	314	A B
	315	+-----------+---------+
	316	xmin x xmax
	317
	318	@param x: x-value
	319	@param xmin: minimum x for the length considered
	320	@param xmax: minimum x for the length considered
	321	@return: (x-xmin)/(xmax-xmin) when xmin < x < xmax
	322
	323	"""
	324	if x<=xmin:
	325	return 0.0
	326	if x>xmin and x<xmax:
	327	return (x-xmin)/(xmax-xmin)
	328	else:
	329	return 1.0
	330
[76e2369]	331
	332	class CircularAverage(object):
	333	"""
	334	Perform circular averaging on 2D data
	335
	336	The data returned is the distribution of counts
	337	as a function of Q
	338	"""
	339	def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.001):
	340	# Minimum radius included in the average [A-1]
	341	self.r_min = r_min
	342	# Maximum radius included in the average [A-1]
	343	self.r_max = r_max
	344	# Bin width (step size) [A-1]
	345	self.bin_width = bin_width
	346
	347	def __call__(self, data2D):
	348	"""
	349	Perform circular averaging on the data
	350
	351	@param data2D: Data2D object
	352	@return: Data1D object
	353	"""
	354	if len(data2D.detector) != 1:
	355	raise RuntimeError, "Circular averaging: invalid number of detectors: %g" % len(data2D.detector)
	356
[2e83ff3]	357	pixel_width_x = data2D.detector[0].pixel_size.x
	358	pixel_width_y = data2D.detector[0].pixel_size.y
[76e2369]	359	det_dist = data2D.detector[0].distance
	360	wavelength = data2D.source.wavelength
[2e83ff3]	361	center_x = data2D.detector[0].beam_center.x/pixel_width_x
	362	center_y = data2D.detector[0].beam_center.y/pixel_width_y
[76e2369]	363
	364	# Find out the maximum Q range
[2e83ff3]	365	xwidth = numpy.size(data2D.data,1)*pixel_width_x
[76e2369]	366	dx_max = xwidth - data2D.detector[0].beam_center.x
	367	if xwidth-dx_max>dx_max:
	368	dx_max = xwidth-dx_max
	369
[2e83ff3]	370	ywidth = numpy.size(data2D.data,0)*pixel_width_y
[76e2369]	371	dy_max = ywidth - data2D.detector[0].beam_center.y
	372	if ywidth-dy_max>dy_max:
	373	dy_max = ywidth-dy_max
	374
	375	qmax = get_q(dx_max, dy_max, det_dist, wavelength)
	376
	377	# Build array of Q intervals
	378	nbins = int(math.ceil((qmax-self.r_min)/self.bin_width))
	379	qbins = self.bin_width*numpy.arange(nbins)+self.r_min
	380
	381	x = numpy.zeros(nbins)
	382	y = numpy.zeros(nbins)
	383	err_y = numpy.zeros(nbins)
	384	y_counts = numpy.zeros(nbins)
	385
[2e83ff3]	386	for i in range(numpy.size(data2D.data,1)):
	387	dx = pixel_width_x*(i+0.5 - center_x)
[76e2369]	388
	389	# Min and max x-value for the pixel
[2e83ff3]	390	minx = pixel_width_x*(i - center_x)
	391	maxx = pixel_width_x*(i+1.0 - center_x)
[76e2369]	392
[2e83ff3]	393	for j in range(numpy.size(data2D.data,0)):
	394	dy = pixel_width_y*(j+0.5 - center_y)
[76e2369]	395
	396	q_value = get_q(dx, dy, det_dist, wavelength)
	397
	398	# Min and max y-value for the pixel
[2e83ff3]	399	miny = pixel_width_y*(j - center_y)
	400	maxy = pixel_width_y*(j+1.0 - center_y)
[76e2369]	401
	402	# Calculate the q-value for each corner
	403	# q_[x min or max][y min or max]
	404	q_00 = get_q(minx, miny, det_dist, wavelength)
	405	q_01 = get_q(minx, maxy, det_dist, wavelength)
	406	q_10 = get_q(maxx, miny, det_dist, wavelength)
	407	q_11 = get_q(maxx, maxy, det_dist, wavelength)
	408
	409	# Look for intercept between each side of the pixel
	410	# and the constant-q ring for qmax
	411	frac_max = get_pixel_fraction(self.r_max, q_00, q_01, q_10, q_11)
	412
	413	# Look for intercept between each side of the pixel
	414	# and the constant-q ring for qmin
	415	frac_min = get_pixel_fraction(self.r_min, q_00, q_01, q_10, q_11)
	416
	417	# We are interested in the region between qmin and qmax
	418	# therefore the fraction of the surface of the pixel
	419	# that we will use to calculate the number of counts to
	420	# include is given by:
	421	frac = frac_max - frac_min
	422
	423	i_q = int(math.ceil((q_value-self.r_min)/self.bin_width)) - 1
	424
	425	x[i_q] = q_value
	426	y[i_q] += frac * data2D.data[j][i]
	427	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
	428	err_y[i_q] += frac * frac * math.fabs(data2D.data[j][i])
	429	else:
	430	err_y[i_q] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
	431	y_counts[i_q] += frac
	432
	433	# Average the sums
	434	for i in range(nbins):
	435	if y_counts[i]>0:
	436	err_y[i] = math.sqrt(err_y[i])/y_counts[i]
	437	y[i] = y[i]/y_counts[i]
	438
	439	return Data1D(x=x, y=y, dy=err_y)
	440
	441
	442	class Ring(object):
	443	"""
	444	Defines a ring on a 2D data set.
	445	The ring is defined by r_min, r_max, and
	446	the position of the center of the ring.
	447
	448	The data returned is the distribution of counts
	449	around the ring as a function of phi.
	450
	451	"""
	452
	453	def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0):
	454	# Minimum radius
	455	self.r_min = r_min
	456	# Maximum radius
	457	self.r_max = r_max
	458	# Center of the ring in x
	459	self.center_x = center_x
	460	# Center of the ring in y
	461	self.center_y = center_y
	462	# Number of angular bins
	463	self.nbins_phi = 20
	464
	465	def __call__(self, data2D):
	466	"""
	467	Apply the ring to the data set.
	468	Returns the angular distribution for a given q range
	469
	470	@param data2D: Data2D object
	471	@return: Data1D object
	472	"""
	473	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
	474	raise RuntimeError, "Ring averaging only take plottable_2D objects"
	475
	476	data = data2D.data
	477	qmin = self.r_min
	478	qmax = self.r_max
	479
	480	if len(data2D.detector) != 1:
	481	raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector)
[2e83ff3]	482	pixel_width_x = data2D.detector[0].pixel_size.x
	483	pixel_width_y = data2D.detector[0].pixel_size.y
[76e2369]	484	det_dist = data2D.detector[0].distance
	485	wavelength = data2D.source.wavelength
[2e83ff3]	486	#center_x = self.center_x/pixel_width_x
	487	#center_y = self.center_y/pixel_width_y
	488	center_x = data2D.detector[0].beam_center.x/pixel_width_x
	489	center_y = data2D.detector[0].beam_center.y/pixel_width_y
	490
[76e2369]	491
	492	phi_bins = numpy.zeros(self.nbins_phi)
	493	phi_counts = numpy.zeros(self.nbins_phi)
	494	phi_values = numpy.zeros(self.nbins_phi)
	495	phi_err = numpy.zeros(self.nbins_phi)
	496
[2e83ff3]	497	for i in range(numpy.size(data,1)):
	498	dx = pixel_width_x*(i+0.5 - center_x)
[76e2369]	499
	500	# Min and max x-value for the pixel
[2e83ff3]	501	minx = pixel_width_x*(i - center_x)
	502	maxx = pixel_width_x*(i+1.0 - center_x)
[76e2369]	503
[2e83ff3]	504	for j in range(numpy.size(data,0)):
	505	dy = pixel_width_y*(j+0.5 - center_y)
[76e2369]	506
	507	q_value = get_q(dx, dy, det_dist, wavelength)
	508
	509	# Min and max y-value for the pixel
[2e83ff3]	510	miny = pixel_width_y*(j - center_y)
	511	maxy = pixel_width_y*(j+1.0 - center_y)
[76e2369]	512
	513	# Calculate the q-value for each corner
	514	# q_[x min or max][y min or max]
	515	q_00 = get_q(minx, miny, det_dist, wavelength)
	516	q_01 = get_q(minx, maxy, det_dist, wavelength)
	517	q_10 = get_q(maxx, miny, det_dist, wavelength)
	518	q_11 = get_q(maxx, maxy, det_dist, wavelength)
	519
	520	# Look for intercept between each side of the pixel
	521	# and the constant-q ring for qmax
	522	frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11)
	523
	524	# Look for intercept between each side of the pixel
	525	# and the constant-q ring for qmin
	526	frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11)
	527
	528	# We are interested in the region between qmin and qmax
	529	# therefore the fraction of the surface of the pixel
	530	# that we will use to calculate the number of counts to
	531	# include is given by:
	532
	533	frac = frac_max - frac_min
	534
[2e83ff3]	535	i_phi = int(math.ceil(self.nbins_phi(math.atan2(dy, dx)+math.pi)/(2.0math.pi))) - 1
[76e2369]	536
	537	phi_bins[i_phi] += frac * data[j][i]
	538
	539	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
	540	phi_err[i_phi] += frac * frac * math.fabs(data2D.data[j][i])
	541	else:
	542	phi_err[i_phi] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
	543	phi_counts[i_phi] += frac
	544
	545	for i in range(self.nbins_phi):
	546	phi_bins[i] = phi_bins[i] / phi_counts[i]
	547	phi_err[i] = math.sqrt(phi_err[i]) / phi_counts[i]
	548	phi_values[i] = 2.0math.pi/self.nbins_phi(1.0*i + 0.5)
	549
	550	return Data1D(x=phi_values, y=phi_bins, dy=phi_err)
	551
	552	def get_pixel_fraction(qmax, q_00, q_01, q_10, q_11):
	553	"""
	554	Returns the fraction of the pixel defined by
	555	the four corners (q_00, q_01, q_10, q_11) that
	556	has q < qmax.
	557
	558	q_01 q_11
	559	y=1 +--------------+
	560	\| \|
	561	\| \|
	562	\| \|
	563	y=0 +--------------+
	564	q_00 q_01
	565
	566	x=0 x=1
	567
	568	"""
	569
	570	# y side for x = minx
	571	x_0 = get_intercept(qmax, q_00, q_01)
	572	# y side for x = maxx
	573	x_1 = get_intercept(qmax, q_10, q_11)
	574
	575	# x side for y = miny
	576	y_0 = get_intercept(qmax, q_00, q_10)
	577	# x side for y = maxy
	578	y_1 = get_intercept(qmax, q_01, q_11)
	579
	580	# surface fraction for a 1x1 pixel
	581	frac_max = 0
	582
	583	if x_0 and x_1:
	584	frac_max = (x_0+x_1)/2.0
	585
	586	elif y_0 and y_1:
	587	frac_max = (y_0+y_1)/2.0
	588
	589	elif x_0 and y_0:
	590	if q_00 < q_10:
	591	frac_max = x_0*y_0/2.0
	592	else:
	593	frac_max = 1.0-x_0*y_0/2.0
	594
	595	elif x_0 and y_1:
	596	if q_00 < q_10:
	597	frac_max = x_0*y_1/2.0
	598	else:
	599	frac_max = 1.0-x_0*y_1/2.0
	600
	601	elif x_1 and y_0:
	602	if q_00 > q_10:
	603	frac_max = x_1*y_0/2.0
	604	else:
	605	frac_max = 1.0-x_1*y_0/2.0
	606
	607	elif x_1 and y_1:
	608	if q_00 < q_10:
	609	frac_max = 1.0 - (1.0-x_1)*(1.0-y_1)/2.0
	610	else:
	611	frac_max = (1.0-x_1)*(1.0-y_1)/2.0
	612
	613	# If we make it here, there is no intercept between
	614	# this pixel and the constant-q ring. We only need
	615	# to know if we have to include it or exclude it.
	616	elif (q_00+q_01+q_10+q_11)/4.0 < qmax:
	617	frac_max = 1.0
	618
	619	return frac_max
	620
	621	def get_intercept(q, q_0, q_1):
	622	"""
	623	Returns the fraction of the side at which the
	624	q-value intercept the pixel, None otherwise.
	625	The values returned is the fraction ON THE SIDE
	626	OF THE LOWEST Q.
	627
	628
	629
	630	A B
	631	+-----------+--------+
	632	0 1 <--- pixel size
	633
	634	Q_0 -------- Q ----- Q_1 <--- equivalent Q range
	635
	636
	637	if Q_1 > Q_0, A is returned
	638	if Q_1 < Q_0, B is returned
	639
	640	if Q is outside the range of [Q_0, Q_1], None is returned
	641
	642	"""
	643	if q_1 > q_0:
	644	if (q > q_0 and q <= q_1):
	645	return (q-q_0)/(q_1 - q_0)
	646	else:
	647	if (q > q_1 and q <= q_0):
	648	return (q-q_1)/(q_0 - q_1)
	649	return None
	650
[c2a8523]	651	#This class can be removed.
[fb198a9]	652	class _Sectorold:
[76e2369]	653	"""
	654	Defines a sector region on a 2D data set.
	655	The sector is defined by r_min, r_max, phi_min, phi_max,
[fb198a9]	656	and the position of the center of the ring.
[2e83ff3]	657	Phi is defined between 0 and 2pi
	658	"""
[f2eee4a]	659	def __init__(self, r_min, r_max, phi_min, phi_max,nbins=20):
[2e83ff3]	660	self.r_min = r_min
	661	self.r_max = r_max
	662	self.phi_min = phi_min
	663	self.phi_max = phi_max
[f2eee4a]	664	self.nbins = nbins
[2e83ff3]	665
	666	def _agv(self, data2D, run='phi'):
	667	"""
	668	Perform sector averaging.
	669
	670	@param data2D: Data2D object
	671	@param run: define the varying parameter ('phi' or 'q')
	672	@return: Data1D object
	673	"""
	674	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
	675	raise RuntimeError, "Ring averaging only take plottable_2D objects"
[fb198a9]	676
	677	data = data2D.data
[2e83ff3]	678	qmax = self.r_max
[fb198a9]	679	qmin = self.r_min
[2e83ff3]	680
	681	if len(data2D.detector) != 1:
	682	raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector)
	683	pixel_width_x = data2D.detector[0].pixel_size.x
	684	pixel_width_y = data2D.detector[0].pixel_size.y
	685	det_dist = data2D.detector[0].distance
	686	wavelength = data2D.source.wavelength
	687	center_x = data2D.detector[0].beam_center.x/pixel_width_x
	688	center_y = data2D.detector[0].beam_center.y/pixel_width_y
	689
	690	y = numpy.zeros(self.nbins)
	691	y_counts = numpy.zeros(self.nbins)
	692	x = numpy.zeros(self.nbins)
	693	y_err = numpy.zeros(self.nbins)
	694
	695	for i in range(numpy.size(data,1)):
	696	dx = pixel_width_x*(i+0.5 - center_x)
	697
	698	# Min and max x-value for the pixel
	699	minx = pixel_width_x*(i - center_x)
	700	maxx = pixel_width_x*(i+1.0 - center_x)
	701
	702	for j in range(numpy.size(data,0)):
	703	dy = pixel_width_y*(j+0.5 - center_y)
	704
	705	q_value = get_q(dx, dy, det_dist, wavelength)
	706
	707	# Min and max y-value for the pixel
	708	miny = pixel_width_y*(j - center_y)
	709	maxy = pixel_width_y*(j+1.0 - center_y)
	710
	711	# Calculate the q-value for each corner
	712	# q_[x min or max][y min or max]
	713	q_00 = get_q(minx, miny, det_dist, wavelength)
	714	q_01 = get_q(minx, maxy, det_dist, wavelength)
	715	q_10 = get_q(maxx, miny, det_dist, wavelength)
	716	q_11 = get_q(maxx, maxy, det_dist, wavelength)
	717
	718	# Look for intercept between each side of the pixel
	719	# and the constant-q ring for qmax
	720	frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11)
	721
	722	# Look for intercept between each side of the pixel
	723	# and the constant-q ring for qmin
	724	frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11)
	725
	726	# We are interested in the region between qmin and qmax
	727	# therefore the fraction of the surface of the pixel
	728	# that we will use to calculate the number of counts to
	729	# include is given by:
	730
	731	frac = frac_max - frac_min
	732
	733	# Compute phi and check whether it's within the limits
[fb198a9]	734	phi_value=math.atan2(dy,dx)+math.pi
	735	# if phi_value<self.phi_min or phi_value>self.phi_max:
[2e83ff3]	736	if phi_value<self.phi_min or phi_value>self.phi_max:
	737	continue
	738
	739	# Check which type of averaging we need
	740	if run.lower()=='phi':
	741	i_bin = int(math.ceil(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) - 1
	742	else:
	743	# If we don't need this pixel, skip the rest of the work
	744	#TODO: an improvement here would be to compute the average
	745	# Q for the pixel from the part that is covered by
	746	# the ring defined by q_min/q_max rather than the complete
	747	# pixel
	748	if q_value<self.r_min or q_value>self.r_max:
	749	continue
	750	i_bin = int(math.ceil(self.nbins*(q_value-self.r_min)/(self.r_max-self.r_min))) - 1
	751
	752	try:
	753	y[i_bin] += frac * data[j][i]
	754	except:
	755	import sys
	756	print sys.exc_value
	757	print i_bin, frac
	758
	759	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
	760	y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i])
	761	else:
	762	y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
	763	y_counts[i_bin] += frac
	764
	765	for i in range(self.nbins):
	766	y[i] = y[i] / y_counts[i]
	767	y_err[i] = math.sqrt(y_err[i]) / y_counts[i]
	768	# Check which type of averaging we need
	769	if run.lower()=='phi':
	770	x[i] = (self.phi_max-self.phi_min)/self.nbins(1.0i + 0.5)+self.phi_min
	771	else:
	772	x[i] = (self.r_max-self.r_min)/self.nbins(1.0i + 0.5)+self.r_min
	773
	774	return Data1D(x=x, y=y, dy=y_err)
	775
[fb198a9]	776	class _Sector:
	777	"""
	778	Defines a sector region on a 2D data set.
	779	The sector is defined by r_min, r_max, phi_min, phi_max,
	780	and the position of the center of the ring
	781	where phi_min and phi_max are defined by the right and left lines wrt central line
	782	and phi_max could be less than phi_min.
	783
	784	Phi is defined between 0 and 2pi
	785	"""
	786	def __init__(self, r_min, r_max, phi_min, phi_max,nbins=20):
	787	self.r_min = r_min
	788	self.r_max = r_max
	789	self.phi_min = phi_min
	790	self.phi_max = phi_max
	791	self.nbins = nbins
	792
	793	def _agv(self, data2D, run='phi'):
	794	"""
	795	Perform sector averaging.
	796
	797	@param data2D: Data2D object
	798	@param run: define the varying parameter ('phi' or 'q')
	799	@return: Data1D object
	800	"""
	801	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
	802	raise RuntimeError, "Ring averaging only take plottable_2D objects"
	803
	804	data = data2D.data
	805	qmax = self.r_max
	806	qmin = self.r_min
	807
	808	if len(data2D.detector) != 1:
	809	raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector)
	810	pixel_width_x = data2D.detector[0].pixel_size.x
	811	pixel_width_y = data2D.detector[0].pixel_size.y
	812	det_dist = data2D.detector[0].distance
	813	wavelength = data2D.source.wavelength
	814	center_x = data2D.detector[0].beam_center.x/pixel_width_x
	815	center_y = data2D.detector[0].beam_center.y/pixel_width_y
	816
	817	y = numpy.zeros(self.nbins)
	818	y_counts = numpy.zeros(self.nbins)
	819	x = numpy.zeros(self.nbins)
	820	y_err = numpy.zeros(self.nbins)
	821
[c2a8523]	822	# This If finds qmax within ROI defined by sector lines
[4853c04]	823	temp=0 #to find qmax within ROI or phimax and phimin
	824	temp0=1000000
	825	temp1=0
	826	for i in range(numpy.size(data,1)):
	827	dx = pixel_width_x*(i+0.5 - center_x)
	828	for j in range(numpy.size(data,0)):
[c2a8523]	829
[4853c04]	830	dy = pixel_width_y*(j+0.5 - center_y)
	831	q_value = get_q(dx, dy, det_dist, wavelength)
	832	# Compute phi and check whether it's within the limits
	833	phi_value=math.atan2(dy,dx)+math.pi
	834	if self.phi_max>2*math.pi:
	835	self.phi_max=self.phi_max-2*math.pi
	836	if self.phi_min<0:
	837	self.phi_max=self.phi_max+2*math.pi
[c2a8523]	838
[4853c04]	839	#In case of two ROI (symmetric major and minor regions)(for 'q2')
	840	if run.lower()=='q2':
	841	if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)):
	842	temp_max=self.phi_max+math.pi
	843	temp_min=self.phi_min+math.pi
	844	else:
	845	temp_max=self.phi_max
	846	temp_min=self.phi_min
[c2a8523]	847
[4853c04]	848	if ((temp_max>=math.pi and temp_max<2math.pi)and (temp_min>=math.pi and temp_min<2math.pi)):
	849	if (phi_value<temp_min or phi_value>temp_max):
	850	if (phi_value<temp_min-math.pi or phi_value>temp_max-math.pi):
	851	continue
	852	if (self.phi_max<self.phi_min):
	853	tmp_max=self.phi_max+math.pi
	854	tmp_min=self.phi_min-math.pi
	855	else:
	856	tmp_max=self.phi_max
	857	tmp_min=self.phi_min
	858	if (tmp_min<math.pi and tmp_max>math.pi):
	859	if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)):
	860	continue
	861	#In case of one ROI (major only)(i.e.,for 'q' and 'phi')
	862	else:
	863	if (self.phi_max>=self.phi_min):
	864	if (phi_value<self.phi_min or phi_value>self.phi_max):
	865	continue
	866	else:
	867	if (phi_value<self.phi_min and phi_value>self.phi_max):
	868	continue
[d9629c53]	869	if q_value<qmin or q_value>qmax:
	870	continue
	871
[4853c04]	872	if run.lower()=='phi':
	873	if temp1<phi_value:
	874	temp1=phi_value
	875	if temp0>phi_value:
[d9629c53]	876	temp0=phi_value
	877
[4853c04]	878	elif temp<q_value:
	879	temp=q_value
[d9629c53]	880
[4853c04]	881	if run.lower()=='phi':
	882	self.phi_max=temp1
	883	self.phi_min=temp0
	884	else:
	885	qmax=temp
	886
[fb198a9]	887	for i in range(numpy.size(data,1)):
	888	dx = pixel_width_x*(i+0.5 - center_x)
	889
	890	# Min and max x-value for the pixel
	891	minx = pixel_width_x*(i - center_x)
	892	maxx = pixel_width_x*(i+1.0 - center_x)
	893
	894	for j in range(numpy.size(data,0)):
	895	dy = pixel_width_y*(j+0.5 - center_y)
	896
	897	q_value = get_q(dx, dy, det_dist, wavelength)
	898
	899	# Min and max y-value for the pixel
	900	miny = pixel_width_y*(j - center_y)
	901	maxy = pixel_width_y*(j+1.0 - center_y)
	902
	903	# Calculate the q-value for each corner
	904	# q_[x min or max][y min or max]
	905	q_00 = get_q(minx, miny, det_dist, wavelength)
	906	q_01 = get_q(minx, maxy, det_dist, wavelength)
	907	q_10 = get_q(maxx, miny, det_dist, wavelength)
	908	q_11 = get_q(maxx, maxy, det_dist, wavelength)
	909
[c2a8523]	910	# Compute phi and check whether it's within the limits
	911	phi_value=math.atan2(dy,dx)+math.pi
	912	if self.phi_max>2*math.pi:
	913	self.phi_max=self.phi_max-2*math.pi
	914	if self.phi_min<0:
	915	self.phi_max=self.phi_max+2*math.pi
	916
[fb198a9]	917	# Look for intercept between each side of the pixel
	918	# and the constant-q ring for qmax
	919	frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11)
	920
	921	# Look for intercept between each side of the pixel
	922	# and the constant-q ring for qmin
	923	frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11)
	924
	925	# We are interested in the region between qmin and qmax
	926	# therefore the fraction of the surface of the pixel
	927	# that we will use to calculate the number of counts to
	928	# include is given by:
	929
	930	frac = frac_max - frac_min
	931
[c2a8523]	932	#In case of two ROI (symmetric major and minor regions)(for 'q2')
[fb198a9]	933	if run.lower()=='q2':
	934	if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)):
	935	temp_max=self.phi_max+math.pi
	936	temp_min=self.phi_min+math.pi
	937	else:
	938	temp_max=self.phi_max
	939	temp_min=self.phi_min
	940
	941	if ((temp_max>=math.pi and temp_max<2math.pi)and (temp_min>=math.pi and temp_min<2math.pi)):
	942	if (phi_value<temp_min or phi_value>temp_max):
	943	if (phi_value<temp_min-math.pi or phi_value>temp_max-math.pi):
	944	continue
	945	if (self.phi_max<self.phi_min):
	946	tmp_max=self.phi_max+math.pi
	947	tmp_min=self.phi_min-math.pi
	948	else:
	949	tmp_max=self.phi_max
	950	tmp_min=self.phi_min
	951	if (tmp_min<math.pi and tmp_max>math.pi):
	952	if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)):
	953	continue
[c2a8523]	954	#In case of one ROI (major only)(i.e.,for 'q' and 'phi')
[fb198a9]	955	else:
	956	if (self.phi_max>=self.phi_min):
	957	if (phi_value<self.phi_min or phi_value>self.phi_max):
	958	continue
[d9629c53]	959
[fb198a9]	960	else:
	961	if (phi_value<self.phi_min and phi_value>self.phi_max):
	962	continue
[4853c04]	963	#print "qmax=",qmax,qmin
	964
	965	if q_value<qmin or q_value>qmax:
	966	continue
[fb198a9]	967
	968	# Check which type of averaging we need
	969	if run.lower()=='phi':
	970	i_bin = int(math.ceil(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) - 1
	971	else:
	972	# If we don't need this pixel, skip the rest of the work
	973	#TODO: an improvement here would be to compute the average
	974	# Q for the pixel from the part that is covered by
	975	# the ring defined by q_min/q_max rather than the complete
	976	# pixel
[c2a8523]	977	i_bin = int(math.ceil(self.nbins*(q_value-qmin)/(qmax-qmin))) - 1
	978
[fb198a9]	979	try:
	980	y[i_bin] += frac * data[j][i]
	981	except:
	982	import sys
	983	print sys.exc_value
	984	print i_bin, frac
	985
	986	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
	987	y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i])
	988	else:
	989	y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
	990	y_counts[i_bin] += frac
[2e83ff3]	991
[fb198a9]	992	for i in range(self.nbins):
	993	y[i] = y[i] / y_counts[i]
	994	y_err[i] = math.sqrt(y_err[i]) / y_counts[i]
	995	# Check which type of averaging we need
	996	if run.lower()=='phi':
[923d926]	997	#Calculate x[i] and put back the origin of angle back to the right hand side (from the left).
	998	x[i] = ((self.phi_max-self.phi_min)/self.nbins(1.0i + 0.5)+self.phi_min-math.pi)*180/math.pi
	999	if x[i]<0:
	1000	x[i]=360+x[i]
[fb198a9]	1001	else:
[c2a8523]	1002	x[i] = (qmax-qmin)/self.nbins(1.0i + 0.5)+qmin
[fb198a9]	1003
	1004	return Data1D(x=x, y=y, dy=y_err)
	1005
[2e83ff3]	1006	class SectorPhi(_Sector):
	1007	"""
	1008	Sector average as a function of phi.
	1009	I(phi) is return and the data is averaged over Q.
	1010
	1011	A sector is defined by r_min, r_max, phi_min, phi_max.
	1012	The number of bin in phi also has to be defined.
	1013	"""
	1014	def __call__(self, data2D):
	1015	"""
	1016	Perform sector average and return I(phi).
	1017
	1018	@param data2D: Data2D object
	1019	@return: Data1D object
	1020	"""
	1021	return self._agv(data2D, 'phi')
	1022
[fb198a9]	1023	class SectorQold(_Sector):
[2e83ff3]	1024	"""
	1025	Sector average as a function of Q.
	1026	I(Q) is return and the data is averaged over phi.
	1027
	1028	A sector is defined by r_min, r_max, phi_min, phi_max.
	1029	The number of bin in Q also has to be defined.
[76e2369]	1030	"""
[2e83ff3]	1031	def __call__(self, data2D):
	1032	"""
	1033	Perform sector average and return I(Q).
	1034
	1035	@param data2D: Data2D object
	1036	@return: Data1D object
	1037	"""
	1038	return self._agv(data2D, 'q')
[fb198a9]	1039
	1040	class SectorQ(_Sector):
	1041	"""
	1042	Sector average as a function of Q for both symatric wings.
	1043	I(Q) is return and the data is averaged over phi.
	1044
	1045	A sector is defined by r_min, r_max, phi_min, phi_max.
	1046	r_min, r_max, phi_min, phi_max >0.
	1047	The number of bin in Q also has to be defined.
	1048	"""
	1049	def __call__(self, data2D):
	1050	"""
	1051	Perform sector average and return I(Q).
	1052
	1053	@param data2D: Data2D object
	1054	@return: Data1D object
	1055	"""
	1056	return self._agv(data2D, 'q2')
[76e2369]	1057	if __name__ == "__main__":
	1058
	1059	from loader import Loader
	1060
	1061
[f8d0ee7]	1062	d = Loader().load('test/MAR07232_rest.ASC')
	1063	#d = Loader().load('test/MP_New.sans')
[76e2369]	1064
	1065
[d9629c53]	1066	r = SectorQ(r_min=.000001, r_max=.01, phi_min=0.0, phi_max=2*math.pi)
[f8d0ee7]	1067	o = r(d)
	1068
[d9629c53]	1069	s = Ring(r_min=.000001, r_max=.01)
[2e83ff3]	1070	p = s(d)
[70975f3]	1071
	1072	for i in range(len(o.x)):
	1073	print o.x[i], o.y[i], o.dy[i]
[76e2369]	1074
	1075
	1076

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

SasView

source: sasview/DataLoader/manipulations.py @ ca10d8e

Download in other formats: