manipulations.py @ 9ce41c6

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 9ce41c6 was 8ba103f, checked in by Jae Cho <jhjcho@…>, 16 years ago
fixed circular average bug
Property mode set to `100644`
File size: 43.7 KB

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

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

SasView

source: sasview/DataLoader/manipulations.py @ 9ce41c6

Download in other formats: