1 | |
---|
2 | ############################################################################ |
---|
3 | #This software was developed by the University of Tennessee as part of the |
---|
4 | #Distributed Data Analysis of Neutron Scattering Experiments (DANSE) |
---|
5 | #project funded by the US National Science Foundation. |
---|
6 | #If you use DANSE applications to do scientific research that leads to |
---|
7 | #publication, we ask that you acknowledge the use of the software with the |
---|
8 | #following sentence: |
---|
9 | #This work benefited from DANSE software developed under NSF award DMR-0520547. |
---|
10 | #copyright 2008, University of Tennessee |
---|
11 | ############################################################################# |
---|
12 | |
---|
13 | """ |
---|
14 | IGOR 2D reduced file reader |
---|
15 | """ |
---|
16 | |
---|
17 | import os |
---|
18 | #import sys |
---|
19 | import numpy |
---|
20 | import math |
---|
21 | #import logging |
---|
22 | from sans.dataloader.data_info import Data2D |
---|
23 | from sans.dataloader.data_info import Detector |
---|
24 | from sans.dataloader.manipulations import reader2D_converter |
---|
25 | |
---|
26 | # Look for unit converter |
---|
27 | has_converter = True |
---|
28 | try: |
---|
29 | from data_util.nxsunit import Converter |
---|
30 | except: |
---|
31 | has_converter = False |
---|
32 | |
---|
33 | class Reader: |
---|
34 | """ Simple data reader for Igor data files """ |
---|
35 | ## File type |
---|
36 | type_name = "IGOR 2D" |
---|
37 | ## Wildcards |
---|
38 | type = ["IGOR 2D files (*.ASC)|*.ASC"] |
---|
39 | ## Extension |
---|
40 | ext=['.ASC', '.asc'] |
---|
41 | |
---|
42 | def read(self,filename=None): |
---|
43 | """ Read file """ |
---|
44 | if not os.path.isfile(filename): |
---|
45 | raise ValueError, \ |
---|
46 | "Specified file %s is not a regular file" % filename |
---|
47 | |
---|
48 | # Read file |
---|
49 | f = open(filename,'r') |
---|
50 | buf = f.read() |
---|
51 | |
---|
52 | # Instantiate data object |
---|
53 | output = Data2D() |
---|
54 | output.filename = os.path.basename(filename) |
---|
55 | detector = Detector() |
---|
56 | if len(output.detector)>0: print str(output.detector[0]) |
---|
57 | output.detector.append(detector) |
---|
58 | |
---|
59 | # Get content |
---|
60 | dataStarted = False |
---|
61 | |
---|
62 | |
---|
63 | lines = buf.split('\n') |
---|
64 | itot = 0 |
---|
65 | x = [] |
---|
66 | y = [] |
---|
67 | |
---|
68 | ncounts = 0 |
---|
69 | |
---|
70 | xmin = None |
---|
71 | xmax = None |
---|
72 | ymin = None |
---|
73 | ymax = None |
---|
74 | |
---|
75 | i_x = 0 |
---|
76 | i_y = -1 |
---|
77 | i_tot_row = 0 |
---|
78 | |
---|
79 | isInfo = False |
---|
80 | isCenter = False |
---|
81 | |
---|
82 | data_conv_q = None |
---|
83 | data_conv_i = None |
---|
84 | |
---|
85 | if has_converter == True and output.Q_unit != '1/A': |
---|
86 | data_conv_q = Converter('1/A') |
---|
87 | # Test it |
---|
88 | data_conv_q(1.0, output.Q_unit) |
---|
89 | |
---|
90 | if has_converter == True and output.I_unit != '1/cm': |
---|
91 | data_conv_i = Converter('1/cm') |
---|
92 | # Test it |
---|
93 | data_conv_i(1.0, output.I_unit) |
---|
94 | |
---|
95 | for line in lines: |
---|
96 | |
---|
97 | # Find setup info line |
---|
98 | if isInfo: |
---|
99 | isInfo = False |
---|
100 | line_toks = line.split() |
---|
101 | # Wavelength in Angstrom |
---|
102 | try: |
---|
103 | wavelength = float(line_toks[1]) |
---|
104 | except: |
---|
105 | msg = "IgorReader: can't read this file, missing wavelength" |
---|
106 | raise ValueError, msg |
---|
107 | |
---|
108 | #Find # of bins in a row assuming the detector is square. |
---|
109 | if dataStarted == True: |
---|
110 | try: |
---|
111 | value = float(line) |
---|
112 | except: |
---|
113 | # Found a non-float entry, skip it |
---|
114 | continue |
---|
115 | |
---|
116 | # Get total bin number |
---|
117 | |
---|
118 | i_tot_row += 1 |
---|
119 | i_tot_row=math.ceil(math.sqrt(i_tot_row))-1 |
---|
120 | #print "i_tot", i_tot_row |
---|
121 | size_x = i_tot_row#192#128 |
---|
122 | size_y = i_tot_row#192#128 |
---|
123 | output.data = numpy.zeros([size_x,size_y]) |
---|
124 | output.err_data = numpy.zeros([size_x,size_y]) |
---|
125 | |
---|
126 | #Read Header and 2D data |
---|
127 | for line in lines: |
---|
128 | # Find setup info line |
---|
129 | if isInfo: |
---|
130 | isInfo = False |
---|
131 | line_toks = line.split() |
---|
132 | # Wavelength in Angstrom |
---|
133 | try: |
---|
134 | wavelength = float(line_toks[1]) |
---|
135 | except: |
---|
136 | msg = "IgorReader: can't read this file, missing wavelength" |
---|
137 | raise ValueError, msg |
---|
138 | # Distance in meters |
---|
139 | try: |
---|
140 | distance = float(line_toks[3]) |
---|
141 | except: |
---|
142 | msg = "IgorReader: can't read this file, missing distance" |
---|
143 | raise ValueError, msg |
---|
144 | |
---|
145 | # Distance in meters |
---|
146 | try: |
---|
147 | transmission = float(line_toks[4]) |
---|
148 | except: |
---|
149 | msg = "IgorReader: can't read this file, " |
---|
150 | msg += "missing transmission" |
---|
151 | raise ValueError, msg |
---|
152 | |
---|
153 | if line.count("LAMBDA")>0: |
---|
154 | isInfo = True |
---|
155 | |
---|
156 | # Find center info line |
---|
157 | if isCenter: |
---|
158 | isCenter = False |
---|
159 | line_toks = line.split() |
---|
160 | |
---|
161 | # Center in bin number: Must substrate 1 because |
---|
162 | #the index starts from 1 |
---|
163 | center_x = float(line_toks[0])-1 |
---|
164 | center_y = float(line_toks[1])-1 |
---|
165 | |
---|
166 | if line.count("BCENT")>0: |
---|
167 | isCenter = True |
---|
168 | |
---|
169 | |
---|
170 | # Find data start |
---|
171 | if line.count("***")>0: |
---|
172 | dataStarted = True |
---|
173 | |
---|
174 | # Check that we have all the info |
---|
175 | if wavelength == None \ |
---|
176 | or distance == None \ |
---|
177 | or center_x == None \ |
---|
178 | or center_y == None: |
---|
179 | msg = "IgorReader:Missing information in data file" |
---|
180 | raise ValueError, msg |
---|
181 | |
---|
182 | if dataStarted == True: |
---|
183 | try: |
---|
184 | value = float(line) |
---|
185 | except: |
---|
186 | # Found a non-float entry, skip it |
---|
187 | continue |
---|
188 | |
---|
189 | # Get bin number |
---|
190 | if math.fmod(itot, i_tot_row)==0: |
---|
191 | i_x = 0 |
---|
192 | i_y += 1 |
---|
193 | else: |
---|
194 | i_x += 1 |
---|
195 | |
---|
196 | output.data[i_y][i_x] = value |
---|
197 | ncounts += 1 |
---|
198 | |
---|
199 | # Det 640 x 640 mm |
---|
200 | # Q = 4pi/lambda sin(theta/2) |
---|
201 | # Bin size is 0.5 cm |
---|
202 | #REmoved +1 from theta = (i_x-center_x+1)*0.5 / distance |
---|
203 | # / 100.0 and |
---|
204 | #REmoved +1 from theta = (i_y-center_y+1)*0.5 / |
---|
205 | # distance / 100.0 |
---|
206 | #ToDo: Need complete check if the following |
---|
207 | # covert process is consistent with fitting.py. |
---|
208 | theta = (i_x-center_x)*0.5 / distance / 100.0 |
---|
209 | qx = 4.0*math.pi/wavelength * math.sin(theta/2.0) |
---|
210 | |
---|
211 | if has_converter == True and output.Q_unit != '1/A': |
---|
212 | qx = data_conv_q(qx, units=output.Q_unit) |
---|
213 | |
---|
214 | if xmin==None or qx<xmin: |
---|
215 | xmin = qx |
---|
216 | if xmax==None or qx>xmax: |
---|
217 | xmax = qx |
---|
218 | |
---|
219 | theta = (i_y-center_y)*0.5 / distance / 100.0 |
---|
220 | qy = 4.0*math.pi/wavelength * math.sin(theta/2.0) |
---|
221 | |
---|
222 | if has_converter == True and output.Q_unit != '1/A': |
---|
223 | qy = data_conv_q(qy, units=output.Q_unit) |
---|
224 | |
---|
225 | if ymin==None or qy<ymin: |
---|
226 | ymin = qy |
---|
227 | if ymax==None or qy>ymax: |
---|
228 | ymax = qy |
---|
229 | |
---|
230 | if not qx in x: |
---|
231 | x.append(qx) |
---|
232 | if not qy in y: |
---|
233 | y.append(qy) |
---|
234 | |
---|
235 | itot += 1 |
---|
236 | |
---|
237 | |
---|
238 | theta = 0.25 / distance / 100.0 |
---|
239 | xstep = 4.0*math.pi/wavelength * math.sin(theta/2.0) |
---|
240 | |
---|
241 | theta = 0.25 / distance / 100.0 |
---|
242 | ystep = 4.0*math.pi/wavelength * math.sin(theta/2.0) |
---|
243 | |
---|
244 | # Store all data ###################################### |
---|
245 | # Store wavelength |
---|
246 | if has_converter==True and output.source.wavelength_unit != 'A': |
---|
247 | conv = Converter('A') |
---|
248 | wavelength = conv(wavelength, units=output.source.wavelength_unit) |
---|
249 | output.source.wavelength = wavelength |
---|
250 | |
---|
251 | # Store distance |
---|
252 | if has_converter==True and detector.distance_unit != 'm': |
---|
253 | conv = Converter('m') |
---|
254 | distance = conv(distance, units=detector.distance_unit) |
---|
255 | detector.distance = distance |
---|
256 | |
---|
257 | # Store transmission |
---|
258 | output.sample.transmission = transmission |
---|
259 | |
---|
260 | # Store pixel size |
---|
261 | pixel = 5.0 |
---|
262 | if has_converter==True and detector.pixel_size_unit != 'mm': |
---|
263 | conv = Converter('mm') |
---|
264 | pixel = conv(pixel, units=detector.pixel_size_unit) |
---|
265 | detector.pixel_size.x = pixel |
---|
266 | detector.pixel_size.y = pixel |
---|
267 | |
---|
268 | # Store beam center in distance units |
---|
269 | detector.beam_center.x = center_x*pixel |
---|
270 | detector.beam_center.y = center_y*pixel |
---|
271 | |
---|
272 | |
---|
273 | # Store limits of the image (2D array) |
---|
274 | xmin =xmin-xstep/2.0 |
---|
275 | xmax =xmax+xstep/2.0 |
---|
276 | ymin =ymin-ystep/2.0 |
---|
277 | ymax =ymax+ystep/2.0 |
---|
278 | if has_converter == True and output.Q_unit != '1/A': |
---|
279 | xmin = data_conv_q(xmin, units=output.Q_unit) |
---|
280 | xmax = data_conv_q(xmax, units=output.Q_unit) |
---|
281 | ymin = data_conv_q(ymin, units=output.Q_unit) |
---|
282 | ymax = data_conv_q(ymax, units=output.Q_unit) |
---|
283 | output.xmin = xmin |
---|
284 | output.xmax = xmax |
---|
285 | output.ymin = ymin |
---|
286 | output.ymax = ymax |
---|
287 | |
---|
288 | # Store x and y axis bin centers |
---|
289 | output.x_bins = x |
---|
290 | output.y_bins = y |
---|
291 | |
---|
292 | # Units |
---|
293 | if data_conv_q is not None: |
---|
294 | output.xaxis("\\rm{Q_{x}}", output.Q_unit) |
---|
295 | output.yaxis("\\rm{Q_{y}}", output.Q_unit) |
---|
296 | else: |
---|
297 | output.xaxis("\\rm{Q_{x}}", 'A^{-1}') |
---|
298 | output.yaxis("\\rm{Q_{y}}", 'A^{-1}') |
---|
299 | |
---|
300 | if data_conv_i is not None: |
---|
301 | output.zaxis("\\rm{Intensity}", output.I_unit) |
---|
302 | else: |
---|
303 | output.zaxis("\\rm{Intensity}","cm^{-1}") |
---|
304 | |
---|
305 | # Store loading process information |
---|
306 | output.meta_data['loader'] = self.type_name |
---|
307 | output = reader2D_converter(output) |
---|
308 | |
---|
309 | return output |
---|
310 | |
---|
311 | if __name__ == "__main__": |
---|
312 | reader = Reader() |
---|
313 | print reader.read("../test/MAR07232_rest.ASC") |
---|
314 | |
---|