manipulations.py @ e575db9

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 e575db9 was 095ab1b, checked in by Jae Cho <jhjcho@…>, 15 years ago
2d average changed accordingly to the changes in 2d data input system: no proportional method nor sub pixel method used
Property mode set to `100644`
File size: 28.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	def flip_phi(phi):
54	"""
55	Correct phi to within the 0 <= to <= 2pi range
56	@return: phi in >=0 and <=2Pi
57	"""
58	Pi = math.pi
59	if phi < 0:
60	phi_out = phi + 2*Pi
61	elif phi > 2*Pi:
62	phi_out = phi - 2*Pi
63	else:
64	phi_out = phi
65	return phi_out
66
67	def reader2D_converter(data2d=None):
68	"""
69	convert old 2d format opened by IhorReader or danse_reader to new Data2D format
70	@param data2d: 2d array of Data2D object
71	@return: 1d arrays of Data2D object
72	"""
73	if data2d.data==None or data2d.x_bins==None or data2d.y_bins==None:
74	raise ValueError,"Can't convert this data: data=None..."
75
76	from DataLoader.data_info import Data2D
77
78	new_x = numpy.tile(data2d.x_bins, (len(data2d.y_bins),1))
79	new_y = numpy.tile(data2d.y_bins, (len(data2d.x_bins),1))
80	new_y = new_y.swapaxes(0,1)
81
82	new_data = data2d.data.flatten()
83	new_err_data = data2d.err_data.flatten()
84	qx_data = new_x.flatten()
85	qy_data = new_y.flatten()
86	q_data = numpy.sqrt(qx_dataqx_data+qy_dataqy_data)
87	mask = numpy.ones(len(new_data), dtype = bool)
88
89	output = Data2D()
90	output = data2d
91	output.data = new_data
92	output.err_data = new_err_data
93	output.qx_data = qx_data
94	output.qy_data = qy_data
95	output.q_data = q_data
96	output.mask = mask
97
98	return output
99
100	class _Slab(object):
101	"""
102	Compute average I(Q) for a region of interest
103	"""
104	def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0, bin_width=0.001):
105	# Minimum Qx value [A-1]
106	self.x_min = x_min
107	# Maximum Qx value [A-1]
108	self.x_max = x_max
109	# Minimum Qy value [A-1]
110	self.y_min = y_min
111	# Maximum Qy value [A-1]
112	self.y_max = y_max
113	# Bin width (step size) [A-1]
114	self.bin_width = bin_width
115	# If True, I(\|Q\|) will be return, otherwise, negative q-values are allowed
116	self.fold = False
117
118	def __call__(self, data2D): return NotImplemented
119
120	def _avg(self, data2D, maj):
121	"""
122	Compute average I(Q_maj) for a region of interest.
123	The major axis is defined as the axis of Q_maj.
124	The minor axis is the axis that we average over.
125
126	@param data2D: Data2D object
127	@param maj_min: min value on the major axis
128	@return: Data1D object
129	"""
130	if len(data2D.detector) != 1:
131	raise RuntimeError, "_Slab._avg: invalid number of detectors: %g" % len(data2D.detector)
132
133	# Get data
134	data = data2D.data
135	q_data = data2D.q_data
136	err_data = data2D.err_data
137	qx_data = data2D.qx_data
138	qy_data = data2D.qy_data
139
140	# Build array of Q intervals
141	if maj=='x':
142	if self.fold: x_min = 0
143	else: x_min = self.x_min
144	nbins = int(math.ceil((self.x_max-x_min)/self.bin_width))
145	qbins = self.bin_width*numpy.arange(nbins)+ x_min
146	elif maj=='y':
147	if self.fold: y_min = 0
148	else: y_min = self.y_min
149	nbins = int(math.ceil((self.y_max-y_min)/self.bin_width))
150	qbins = self.bin_width*numpy.arange(nbins)+ y_min
151	else:
152	raise RuntimeError, "_Slab._avg: unrecognized axis %s" % str(maj)
153
154	x = numpy.zeros(nbins)
155	y = numpy.zeros(nbins)
156	err_y = numpy.zeros(nbins)
157	y_counts = numpy.zeros(nbins)
158
159	# Average pixelsize in q space
160	for npts in range(len(data)):
161	# default frac
162	frac_x = 0
163	frac_y = 0
164	# get ROI
165	if self.x_min <= qx_data[npts] and self.x_max > qx_data[npts]:
166	frac_x = 1
167	if self.y_min <= qy_data[npts] and self.y_max > qy_data[npts]:
168	frac_y = 1
169
170	frac = frac_x * frac_y
171
172	if frac == 0: continue
173
174	# binning: find axis of q
175	if maj=='x':
176	q_value = qx_data[npts]
177	min = x_min
178	if maj=='y':
179	q_value = qy_data[npts]
180	min = y_min
181	if self.fold and q_value<0: q_value = -q_value
182
183	# bin
184	i_q = int(math.ceil((q_value-min)/self.bin_width)) - 1
185
186	# skip outside of max bins
187	if i_q<0 or i_q>=nbins: continue
188
189	# give it full weight
190	#frac = 1
191
192	#TODO: find better definition of x[i_q] based on q_data
193	x[i_q] = min +(i_q+1)*self.bin_width/2.0
194	y[i_q] += frac * data[npts]
195
196	if err_data == None or err_data[npts]==0.0:
197	err_y[i_q] += frac * frac * math.fabs(data[npts])
198	else:
199	err_y[i_q] += frac * frac * err_data[npts] * err_data[npts]
200	y_counts[i_q] += frac
201
202	# Average the sums
203	for n in range(nbins):
204	err_y[n] = math.sqrt(err_y[n])
205
206	err_y = err_y/y_counts
207	y = y/y_counts
208
209	idx = (numpy.isfinite(y)& numpy.isfinite(x))
210
211	if not idx.any():
212	raise ValueError, "Average Error: No points inside ROI to average..."
213	elif len(y[idx])!= nbins:
214	print "resulted",nbins- len(y[idx]),"empty bin(s) due to tight binning..."
215	return Data1D(x=x[idx], y=y[idx], dy=err_y[idx])
216
217	class SlabY(_Slab):
218	"""
219	Compute average I(Qy) for a region of interest
220	"""
221	def __call__(self, data2D):
222	"""
223	Compute average I(Qy) for a region of interest
224
225	@param data2D: Data2D object
226	@return: Data1D object
227	"""
228	return self._avg(data2D, 'y')
229
230	class SlabX(_Slab):
231	"""
232	Compute average I(Qx) for a region of interest
233	"""
234	def __call__(self, data2D):
235	"""
236	Compute average I(Qx) for a region of interest
237
238	@param data2D: Data2D object
239	@return: Data1D object
240	"""
241	return self._avg(data2D, 'x')
242
243	class Boxsum(object):
244	"""
245	Perform the sum of counts in a 2D region of interest.
246	"""
247	def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0):
248	# Minimum Qx value [A-1]
249	self.x_min = x_min
250	# Maximum Qx value [A-1]
251	self.x_max = x_max
252	# Minimum Qy value [A-1]
253	self.y_min = y_min
254	# Maximum Qy value [A-1]
255	self.y_max = y_max
256
257	def __call__(self, data2D):
258	"""
259	Perform the sum in the region of interest
260
261	@param data2D: Data2D object
262	@return: number of counts, error on number of counts
263	"""
264	y, err_y, y_counts = self._sum(data2D)
265
266	# Average the sums
267	counts = 0 if y_counts==0 else y
268	error = 0 if y_counts==0 else math.sqrt(err_y)
269
270	return counts, error
271
272	def _sum(self, data2D):
273	"""
274	Perform the sum in the region of interest
275	@param data2D: Data2D object
276	@return: number of counts, error on number of counts, number of entries summed
277	"""
278	if len(data2D.detector) != 1:
279	raise RuntimeError, "Circular averaging: invalid number of detectors: %g" % len(data2D.detector)
280
281	# Get data
282	data = data2D.data
283	q_data = data2D.q_data
284	err_data = data2D.err_data
285	qx_data = data2D.qx_data
286	qy_data = data2D.qy_data
287
288	y = 0.0
289	err_y = 0.0
290	y_counts = 0.0
291
292	# Average pixelsize in q space
293	for npts in range(len(data)):
294	# default frac
295	frac_x = 0
296	frac_y = 0
297
298	# get min and max at each points
299	qx = qx_data[npts]
300	qy = qy_data[npts]
301
302	# get the ROI
303	if self.x_min <= qx and self.x_max > qx:
304	frac_x = 1
305	if self.y_min <= qy and self.y_max > qy:
306	frac_y = 1
307
308	#Find the fraction along each directions
309	frac = frac_x * frac_y
310	if frac == 0: continue
311
312	y += frac * data[npts]
313	if err_data == None or err_data[npts]==0.0:
314	err_y += frac * frac * math.fabs(data[npts])
315	else:
316	err_y += frac * frac * err_data[npts] * err_data[npts]
317	y_counts += frac
318
319	return y, err_y, y_counts
320
321
322
323	class Boxavg(Boxsum):
324	"""
325	Perform the average of counts in a 2D region of interest.
326	"""
327	def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0):
328	super(Boxavg, self).__init__(x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max)
329
330	def __call__(self, data2D):
331	"""
332	Perform the sum in the region of interest
333
334	@param data2D: Data2D object
335	@return: average counts, error on average counts
336	"""
337	y, err_y, y_counts = self._sum(data2D)
338
339	# Average the sums
340	counts = 0 if y_counts==0 else y/y_counts
341	error = 0 if y_counts==0 else math.sqrt(err_y)/y_counts
342
343	return counts, error
344
345	def get_pixel_fraction_square(x, xmin, xmax):
346	"""
347	Return the fraction of the length
348	from xmin to x.
349
350	A B
351	+-----------+---------+
352	xmin x xmax
353
354	@param x: x-value
355	@param xmin: minimum x for the length considered
356	@param xmax: minimum x for the length considered
357	@return: (x-xmin)/(xmax-xmin) when xmin < x < xmax
358
359	"""
360	if x<=xmin:
361	return 0.0
362	if x>xmin and x<xmax:
363	return (x-xmin)/(xmax-xmin)
364	else:
365	return 1.0
366
367
368	class CircularAverage(object):
369	"""
370	Perform circular averaging on 2D data
371
372	The data returned is the distribution of counts
373	as a function of Q
374	"""
375	def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.0005):
376	# Minimum radius included in the average [A-1]
377	self.r_min = r_min
378	# Maximum radius included in the average [A-1]
379	self.r_max = r_max
380	# Bin width (step size) [A-1]
381	self.bin_width = bin_width
382
383	def __call__(self, data2D):
384	"""
385	Perform circular averaging on the data
386
387	@param data2D: Data2D object
388	@return: Data1D object
389	"""
390	# Get data
391	data = data2D.data
392	err_data = data2D.err_data
393	q_data = data2D.q_data
394
395	q_data_max = numpy.max(q_data)
396
397	if len(data2D.q_data) == None:
398	raise RuntimeError, "Circular averaging: invalid q_data: %g" % data2D.q_data
399
400	# Build array of Q intervals
401	nbins = int(math.ceil((self.r_max-self.r_min)/self.bin_width))
402	qbins = self.bin_width*numpy.arange(nbins)+self.r_min
403
404	x = numpy.zeros(nbins)
405	y = numpy.zeros(nbins)
406	err_y = numpy.zeros(nbins)
407	y_counts = numpy.zeros(nbins)
408
409	for npt in range(len(data)):
410	frac = 0
411
412	# q-value at the pixel (j,i)
413	q_value = q_data[npt]
414
415	data_n = data[npt]
416
417	## No need to calculate the frac when all data are within range
418	if self.r_min >= self.r_max:
419	raise ValueError, "Limit Error: min > max ???"
420
421	if self.r_min <= q_value and q_value <= self.r_max: frac = 1
422
423	if frac == 0: continue
424
425	i_q = int(math.floor((q_value-self.r_min)/self.bin_width))
426
427	# Take care of the edge case at phi = 2pi.
428	if i_q == nbins:
429	i_q = nbins -1
430
431	y[i_q] += frac * data_n
432
433	if err_data == None or err_data[npt]==0.0:
434	err_y[i_q] += frac * frac * math.fabs(data_n)
435	else:
436	err_y[i_q] += frac * frac * err_data[npt] * err_data[npt]
437	y_counts[i_q] += frac
438
439	## x should be the center value of each bins
440	x = qbins+self.bin_width/2
441
442	# Average the sums
443	for n in range(nbins):
444	err_y[n] = math.sqrt(err_y[n])
445
446	err_y = err_y/y_counts
447	y = y/y_counts
448	idx = (numpy.isfinite(y))&(numpy.isfinite(x))
449
450	if not idx.any():
451	raise ValueError, "Average Error: No points inside ROI to average..."
452	elif len(y[idx])!= nbins:
453	print "resulted",nbins- len(y[idx]),"empty bin(s) due to tight binning..."
454
455	return Data1D(x=x[idx], y=y[idx], dy=err_y[idx])
456
457
458	class Ring(object):
459	"""
460	Defines a ring on a 2D data set.
461	The ring is defined by r_min, r_max, and
462	the position of the center of the ring.
463
464	The data returned is the distribution of counts
465	around the ring as a function of phi.
466
467	Phi_min and phi_max should be defined between 0 and 2*pi
468	in anti-clockwise starting from the x- axis on the left-hand side
469	"""
470	#Todo: remove center.
471	def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0,nbins=20 ):
472	# Minimum radius
473	self.r_min = r_min
474	# Maximum radius
475	self.r_max = r_max
476	# Center of the ring in x
477	self.center_x = center_x
478	# Center of the ring in y
479	self.center_y = center_y
480	# Number of angular bins
481	self.nbins_phi = nbins
482
483	def __call__(self, data2D):
484	"""
485	Apply the ring to the data set.
486	Returns the angular distribution for a given q range
487
488	@param data2D: Data2D object
489	@return: Data1D object
490	"""
491	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
492	raise RuntimeError, "Ring averaging only take plottable_2D objects"
493
494	Pi = math.pi
495
496	# Get data
497	data = data2D.data
498	err_data = data2D.err_data
499	q_data = data2D.q_data
500	qx_data = data2D.qx_data
501	qy_data = data2D.qy_data
502
503	q_data_max = numpy.max(q_data)
504
505	# Set space for 1d outputs
506	phi_bins = numpy.zeros(self.nbins_phi)
507	phi_counts = numpy.zeros(self.nbins_phi)
508	phi_values = numpy.zeros(self.nbins_phi)
509	phi_err = numpy.zeros(self.nbins_phi)
510
511	for npt in range(len(data)):
512	frac = 0
513
514	# q-value at the point (npt)
515	q_value = q_data[npt]
516
517	data_n = data[npt]
518
519	# phi-value at the point (npt)
520	phi_value=math.atan2(qy_data[npt],qx_data[npt])+Pi
521
522	if self.r_min <= q_value and q_value <= self.r_max: frac = 1
523
524	if frac == 0: continue
525
526	# binning
527	i_phi = int(math.floor((self.nbins_phi)phi_value/(2Pi)))
528
529	# Take care of the edge case at phi = 2pi.
530	if i_phi == self.nbins_phi:
531	i_phi = self.nbins_phi -1
532
533	phi_bins[i_phi] += frac * data[npt]
534
535	if err_data == None or err_data[npt] ==0.0:
536	phi_err[i_phi] += frac * frac * math.fabs(data_n)
537	else:
538	phi_err[i_phi] += frac * frac err_data[npt]err_data[npt]
539	phi_counts[i_phi] += frac
540
541	for i in range(self.nbins_phi):
542	phi_bins[i] = phi_bins[i] / phi_counts[i]
543	phi_err[i] = math.sqrt(phi_err[i]) / phi_counts[i]
544	phi_values[i] = 2.0math.pi/self.nbins_phi(1.0*i + 0.5)
545
546	idx = (numpy.isfinite(phi_bins))
547
548	if not idx.any():
549	raise ValueError, "Average Error: No points inside ROI to average..."
550	elif len(phi_bins[idx])!= self.nbins_phi:
551	print "resulted",self.nbins_phi- len(phi_bins[idx]),"empty bin(s) due to tight binning..."
552	return Data1D(x=phi_values[idx], y=phi_bins[idx], dy=phi_err[idx])
553
554	def get_pixel_fraction(qmax, q_00, q_01, q_10, q_11):
555	"""
556	Returns the fraction of the pixel defined by
557	the four corners (q_00, q_01, q_10, q_11) that
558	has q < qmax.
559
560	q_01 q_11
561	y=1 +--------------+
562	\| \|
563	\| \|
564	\| \|
565	y=0 +--------------+
566	q_00 q_10
567
568	x=0 x=1
569
570	"""
571
572	# y side for x = minx
573	x_0 = get_intercept(qmax, q_00, q_01)
574	# y side for x = maxx
575	x_1 = get_intercept(qmax, q_10, q_11)
576
577	# x side for y = miny
578	y_0 = get_intercept(qmax, q_00, q_10)
579	# x side for y = maxy
580	y_1 = get_intercept(qmax, q_01, q_11)
581
582	# surface fraction for a 1x1 pixel
583	frac_max = 0
584
585	if x_0 and x_1:
586	frac_max = (x_0+x_1)/2.0
587
588	elif y_0 and y_1:
589	frac_max = (y_0+y_1)/2.0
590
591	elif x_0 and y_0:
592	if q_00 < q_10:
593	frac_max = x_0*y_0/2.0
594	else:
595	frac_max = 1.0-x_0*y_0/2.0
596
597	elif x_0 and y_1:
598	if q_00 < q_10:
599	frac_max = x_0*y_1/2.0
600	else:
601	frac_max = 1.0-x_0*y_1/2.0
602
603	elif x_1 and y_0:
604	if q_00 > q_10:
605	frac_max = x_1*y_0/2.0
606	else:
607	frac_max = 1.0-x_1*y_0/2.0
608
609	elif x_1 and y_1:
610	if q_00 < q_10:
611	frac_max = 1.0 - (1.0-x_1)*(1.0-y_1)/2.0
612	else:
613	frac_max = (1.0-x_1)*(1.0-y_1)/2.0
614
615	# If we make it here, there is no intercept between
616	# this pixel and the constant-q ring. We only need
617	# to know if we have to include it or exclude it.
618	elif (q_00+q_01+q_10+q_11)/4.0 < qmax:
619	frac_max = 1.0
620
621	return frac_max
622
623	def get_intercept(q, q_0, q_1):
624	"""
625	Returns the fraction of the side at which the
626	q-value intercept the pixel, None otherwise.
627	The values returned is the fraction ON THE SIDE
628	OF THE LOWEST Q.
629
630
631
632	A B
633	+-----------+--------+
634	0 1 <--- pixel size
635
636	Q_0 -------- Q ----- Q_1 <--- equivalent Q range
637
638
639	if Q_1 > Q_0, A is returned
640	if Q_1 < Q_0, B is returned
641
642	if Q is outside the range of [Q_0, Q_1], None is returned
643
644	"""
645	if q_1 > q_0:
646	if (q > q_0 and q <= q_1):
647	return (q-q_0)/(q_1 - q_0)
648	else:
649	if (q > q_1 and q <= q_0):
650	return (q-q_1)/(q_0 - q_1)
651	return None
652
653	class _Sector:
654	"""
655	Defines a sector region on a 2D data set.
656	The sector is defined by r_min, r_max, phi_min, phi_max,
657	and the position of the center of the ring
658	where phi_min and phi_max are defined by the right and left lines wrt central line
659	and phi_max could be less than phi_min.
660
661	Phi is defined between 0 and 2*pi in anti-clockwise starting from the x- axis on the left-hand side
662	"""
663	def __init__(self, r_min, r_max, phi_min=0, phi_max=2*math.pi,nbins=20):
664	self.r_min = r_min
665	self.r_max = r_max
666	self.phi_min = phi_min
667	self.phi_max = phi_max
668	self.nbins = nbins
669
670
671	def _agv(self, data2D, run='phi'):
672	"""
673	Perform sector averaging.
674
675	@param data2D: Data2D object
676	@param run: define the varying parameter ('phi' , 'q' , or 'q2')
677	@return: Data1D object
678	"""
679	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
680	raise RuntimeError, "Ring averaging only take plottable_2D objects"
681	Pi = math.pi
682
683	# Get the all data & info
684	data = data2D.data
685	err_data = data2D.err_data
686	qx_data = data2D.qx_data
687	qy_data = data2D.qy_data
688	q_data = data2D.q_data
689
690	#set space for 1d outputs
691	x = numpy.zeros(self.nbins)
692	y = numpy.zeros(self.nbins)
693	y_err = numpy.zeros(self.nbins)
694	y_counts = numpy.zeros(self.nbins)
695
696	# Get the min and max into the region: 0 <= phi < 2Pi
697	phi_min = flip_phi(self.phi_min)
698	phi_max = flip_phi(self.phi_max)
699
700	q_data_max = numpy.max(q_data)
701
702	for n in range(len(data)):
703	frac = 0
704
705	# q-value at the pixel (j,i)
706	q_value = q_data[n]
707
708
709	data_n = data[n]
710
711	# Is pixel within range?
712	is_in = False
713
714	# phi-value of the pixel (j,i)
715	phi_value=math.atan2(qy_data[n],qx_data[n])+Pi
716
717	## No need to calculate the frac when all data are within range
718	if self.r_min <= q_value and q_value <= self.r_max: frac = 1
719
720	if frac == 0: continue
721
722	#In case of two ROIs (symmetric major and minor regions)(for 'q2')
723	if run.lower()=='q2':
724	## For minor sector wing
725	# Calculate the minor wing phis
726	phi_min_minor = flip_phi(phi_min-Pi)
727	phi_max_minor = flip_phi(phi_max-Pi)
728	# Check if phis of the minor ring is within 0 to 2pi
729	if phi_min_minor > phi_max_minor:
730	is_in = (phi_value > phi_min_minor or phi_value < phi_max_minor)
731	else:
732	is_in = (phi_value > phi_min_minor and phi_value < phi_max_minor)
733
734	#For all cases(i.e.,for 'q', 'q2', and 'phi')
735	#Find pixels within ROI
736	if phi_min > phi_max:
737	is_in = is_in or (phi_value > phi_min or phi_value < phi_max)
738	else:
739	is_in = is_in or (phi_value>= phi_min and phi_value <phi_max)
740
741	if not is_in: frac = 0
742	if frac == 0: continue
743
744	# Check which type of averaging we need
745	if run.lower()=='phi':
746	i_bin = int(math.floor((self.nbins)*(phi_value-self.phi_min)\
747	/(self.phi_max-self.phi_min)))
748	else:
749	i_bin = int(math.floor((self.nbins)*(q_value-self.r_min)/(self.r_max-self.r_min)))
750
751	# Take care of the edge case at phi = 2pi.
752	if i_bin == self.nbins:
753	i_bin = self.nbins -1
754
755	## Get the total y
756	y[i_bin] += frac * data_n
757
758	if err_data == None or err_data[n] ==0.0:
759	y_err[i_bin] += frac * frac * math.fabs(data_n)
760	else:
761	y_err[i_bin] += frac * frac * err_data[n]*err_data[n]
762	y_counts[i_bin] += frac
763
764	# Organize the results
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
769	# The type of averaging: phi,q2, or q
770	# Calculate x[i]should be at the center of the bin
771	if run.lower()=='phi':
772	x[i] = (self.phi_max-self.phi_min)/self.nbins(1.0i + 0.5)+self.phi_min
773	else:
774	x[i] = (self.r_max-self.r_min)/self.nbins(1.0i + 0.5)+self.r_min
775
776	idx = (numpy.isfinite(y)& numpy.isfinite(y_err))
777
778	if not idx.any():
779	raise ValueError, "Average Error: No points inside sector of ROI to average..."
780	elif len(y[idx])!= self.nbins:
781	print "resulted",self.nbins- len(y[idx]),"empty bin(s) due to tight binning..."
782	return Data1D(x=x[idx], y=y[idx], dy=y_err[idx])
783
784	class SectorPhi(_Sector):
785	"""
786	Sector average as a function of phi.
787	I(phi) is return and the data is averaged over Q.
788
789	A sector is defined by r_min, r_max, phi_min, phi_max.
790	The number of bin in phi also has to be defined.
791	"""
792	def __call__(self, data2D):
793	"""
794	Perform sector average and return I(phi).
795
796	@param data2D: Data2D object
797	@return: Data1D object
798	"""
799	return self._agv(data2D, 'phi')
800
801	class SectorQ(_Sector):
802	"""
803	Sector average as a function of Q for both symatric wings.
804	I(Q) is return and the data is averaged over phi.
805
806	A sector is defined by r_min, r_max, phi_min, phi_max.
807	r_min, r_max, phi_min, phi_max >0.
808	The number of bin in Q also has to be defined.
809	"""
810	def __call__(self, data2D):
811	"""
812	Perform sector average and return I(Q).
813
814	@param data2D: Data2D object
815	@return: Data1D object
816	"""
817	return self._agv(data2D, 'q2')
818	if __name__ == "__main__":
819
820	from loader import Loader
821
822
823	d = Loader().load('test/MAR07232_rest.ASC')
824	#d = Loader().load('test/MP_New.sans')
825
826
827	r = SectorQ(r_min=.000001, r_max=.01, phi_min=0.0, phi_max=2*math.pi)
828	o = r(d)
829
830	s = Ring(r_min=.000001, r_max=.01)
831	p = s(d)
832
833	for i in range(len(o.x)):
834	print o.x[i], o.y[i], o.dy[i]
835
836
837

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

SasView

source: sasview/DataLoader/manipulations.py @ e575db9

Download in other formats: