manipulations.py @ 1fbb343

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

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

SasView

source: sasview/DataLoader/manipulations.py @ 1fbb343

Download in other formats: