[3cd95c8] | 1 | """ |
---|
| 2 | TXT/IGOR 2D Q Map file reader |
---|
| 3 | """ |
---|
| 4 | |
---|
| 5 | """ |
---|
| 6 | This software was developed by the University of Tennessee as part of the |
---|
| 7 | Distributed Data Analysis of Neutron Scattering Experiments (DANSE) |
---|
| 8 | project funded by the US National Science Foundation. |
---|
| 9 | |
---|
| 10 | If you use DANSE applications to do scientific research that leads to |
---|
| 11 | publication, we ask that you acknowledge the use of the software with the |
---|
| 12 | following sentence: |
---|
| 13 | |
---|
| 14 | "This work benefited from DANSE software developed under NSF award DMR-0520547." |
---|
| 15 | |
---|
| 16 | copyright 2008, University of Tennessee |
---|
| 17 | """ |
---|
| 18 | |
---|
| 19 | import os, sys |
---|
| 20 | import numpy |
---|
| 21 | import math, logging |
---|
| 22 | from DataLoader.data_info import Data2D, Detector |
---|
| 23 | |
---|
| 24 | # Look for unit converter |
---|
| 25 | has_converter = True |
---|
| 26 | try: |
---|
| 27 | from data_util.nxsunit import Converter |
---|
| 28 | except: |
---|
| 29 | has_converter = False |
---|
| 30 | |
---|
| 31 | class Reader: |
---|
| 32 | """ Simple data reader for Igor data files """ |
---|
| 33 | ## File type |
---|
| 34 | type_name = "IGOR/DAT 2D Q_map" |
---|
| 35 | ## Wildcards |
---|
| 36 | type = ["IGOR/DAT 2D file in Q_map (*.dat)|*.DAT"] |
---|
| 37 | ## Extension |
---|
| 38 | ext=['.DAT', '.dat'] |
---|
| 39 | |
---|
| 40 | def read(self,filename=None): |
---|
| 41 | """ Read file """ |
---|
| 42 | if not os.path.isfile(filename): |
---|
| 43 | raise ValueError, \ |
---|
| 44 | "Specified file %s is not a regular file" % filename |
---|
| 45 | |
---|
| 46 | # Read file |
---|
| 47 | f = open(filename,'r') |
---|
| 48 | buf = f.read() |
---|
| 49 | |
---|
| 50 | # Instantiate data object |
---|
| 51 | output = Data2D() |
---|
| 52 | output.filename = os.path.basename(filename) |
---|
| 53 | detector = Detector() |
---|
| 54 | if len(output.detector)>0: print str(output.detector[0]) |
---|
| 55 | output.detector.append(detector) |
---|
| 56 | |
---|
| 57 | # Get content |
---|
| 58 | dataStarted = False |
---|
| 59 | |
---|
| 60 | ## Defaults |
---|
| 61 | lines = buf.split('\n') |
---|
| 62 | itot = 0 |
---|
| 63 | x = [] |
---|
| 64 | y = [] |
---|
| 65 | |
---|
| 66 | ncounts = 0 |
---|
| 67 | |
---|
| 68 | wavelength = None |
---|
| 69 | distance = None |
---|
| 70 | transmission = None |
---|
| 71 | |
---|
| 72 | pixel_x = None |
---|
| 73 | pixel_y = None |
---|
| 74 | |
---|
| 75 | i_x = 0 |
---|
| 76 | i_y = -1 |
---|
| 77 | pixels = 0 |
---|
| 78 | |
---|
| 79 | isInfo = False |
---|
| 80 | isCenter = False |
---|
| 81 | |
---|
| 82 | data_conv_q = None |
---|
| 83 | data_conv_i = None |
---|
| 84 | |
---|
| 85 | # Set units: This is the unit assumed for Q and I in the data file. |
---|
| 86 | if has_converter == True and output.Q_unit != '1/A': |
---|
| 87 | data_conv_q = Converter('1/A') |
---|
| 88 | # Test it |
---|
| 89 | data_conv_q(1.0, output.Q_unit) |
---|
| 90 | |
---|
| 91 | if has_converter == True and output.I_unit != '1/cm': |
---|
| 92 | data_conv_i = Converter('1/cm') |
---|
| 93 | # Test it |
---|
| 94 | data_conv_i(1.0, output.I_unit) |
---|
| 95 | |
---|
| 96 | #Set the space for data |
---|
| 97 | data = numpy.zeros(0) |
---|
| 98 | qx_data = numpy.zeros(0) |
---|
| 99 | qy_data = numpy.zeros(0) |
---|
| 100 | q_data = numpy.zeros(0) |
---|
| 101 | dqx_data = numpy.zeros(0) |
---|
| 102 | dqy_data = numpy.zeros(0) |
---|
| 103 | mask = numpy.zeros(0,dtype=bool) |
---|
| 104 | |
---|
| 105 | #Read Header and 2D data |
---|
| 106 | for line in lines: |
---|
| 107 | ## Reading the header applies only to IGOR/NIST 2D q_map data files |
---|
| 108 | # Find setup info line |
---|
| 109 | if isInfo: |
---|
| 110 | isInfo = False |
---|
| 111 | line_toks = line.split() |
---|
| 112 | # Wavelength in Angstrom |
---|
| 113 | try: |
---|
| 114 | wavelength = float(line_toks[1]) |
---|
| 115 | # Units |
---|
| 116 | if has_converter==True and output.source.wavelength_unit != 'A': |
---|
| 117 | conv = Converter('A') |
---|
| 118 | wavelength = conv(wavelength, units=output.source.wavelength_unit) |
---|
| 119 | except: |
---|
| 120 | #Not required |
---|
| 121 | pass |
---|
| 122 | # Distance in mm |
---|
| 123 | try: |
---|
| 124 | distance = float(line_toks[3]) |
---|
| 125 | # Units |
---|
| 126 | if has_converter==True and detector.distance_unit != 'm': |
---|
| 127 | conv = Converter('m') |
---|
| 128 | distance = conv(distance, units=detector.distance_unit) |
---|
| 129 | except: |
---|
| 130 | #Not required |
---|
| 131 | pass |
---|
| 132 | |
---|
| 133 | # Distance in meters |
---|
| 134 | try: |
---|
| 135 | transmission = float(line_toks[4]) |
---|
| 136 | except: |
---|
| 137 | #Not required |
---|
| 138 | pass |
---|
| 139 | |
---|
| 140 | if line.count("LAMBDA")>0: |
---|
| 141 | isInfo = True |
---|
| 142 | |
---|
| 143 | # Find center info line |
---|
| 144 | if isCenter: |
---|
| 145 | isCenter = False |
---|
| 146 | line_toks = line.split() |
---|
| 147 | # Center in bin number |
---|
| 148 | center_x = float(line_toks[0]) |
---|
| 149 | center_y = float(line_toks[1]) |
---|
| 150 | |
---|
| 151 | if line.count("BCENT")>0: |
---|
| 152 | isCenter = True |
---|
| 153 | |
---|
| 154 | # Find data start |
---|
| 155 | if line.count("Data columns") or line.count("ASCII data")>0: |
---|
| 156 | dataStarted = True |
---|
| 157 | continue |
---|
| 158 | |
---|
| 159 | ## Read and get data. |
---|
| 160 | if dataStarted == True: |
---|
| 161 | line_toks = line.split() |
---|
| 162 | if len(line_toks) == 0: |
---|
| 163 | #empty line |
---|
| 164 | continue |
---|
[32e8c78] | 165 | elif len(line_toks) < 3: |
---|
[3cd95c8] | 166 | #Q-map 2d require 3 or more columns of data |
---|
| 167 | raise ValueError,"Igor_Qmap_Reader: Can't read this file: Not a proper file format" |
---|
| 168 | |
---|
| 169 | # defaults |
---|
| 170 | dqx_value = None |
---|
| 171 | dqy_value = None |
---|
| 172 | qz_value = None |
---|
| 173 | mask_value = 0 |
---|
| 174 | |
---|
| 175 | ##Read and get data values |
---|
| 176 | qx_value = float(line_toks[0]) |
---|
| 177 | qy_value = float(line_toks[1]) |
---|
| 178 | value = float(line_toks[2]) |
---|
| 179 | |
---|
| 180 | try: |
---|
| 181 | #Read qz_value if exist on 4th column |
---|
| 182 | qz_value = float(line_toks[3]) |
---|
| 183 | except: |
---|
| 184 | # Found a non-float entry, skip it: Not required. |
---|
| 185 | pass |
---|
| 186 | try: |
---|
| 187 | #Read qz_value if exist on 5th column |
---|
| 188 | dqx_value = float(line_toks[4]) |
---|
| 189 | except: |
---|
| 190 | # Found a non-float entry, skip it: Not required. |
---|
| 191 | pass |
---|
| 192 | try: |
---|
| 193 | #Read qz_value if exist on 6th column |
---|
| 194 | dqy_value = float(line_toks[5]) |
---|
| 195 | except: |
---|
| 196 | # Found a non-float entry, skip it: Not required. |
---|
| 197 | pass |
---|
| 198 | try: |
---|
| 199 | #Read beam block mask if exist on 7th column |
---|
| 200 | mask_value = float(line_toks[6]) |
---|
| 201 | except: |
---|
| 202 | # Found a non-float entry, skip it |
---|
| 203 | pass |
---|
| 204 | |
---|
| 205 | # get data |
---|
| 206 | data = numpy.append(data, value) |
---|
| 207 | qx_data = numpy.append(qx_data, qx_value) |
---|
| 208 | qy_data = numpy.append(qy_data, qy_value) |
---|
| 209 | |
---|
| 210 | # optional data |
---|
| 211 | if dqx_value != None and numpy.isfinite(dqx_value): |
---|
| 212 | dqx_data = numpy.append(dqx_data, dqx_value) |
---|
| 213 | if dqy_value != None and numpy.isfinite(dqy_value): |
---|
| 214 | dqy_data = numpy.append(dqy_data, dqy_value) |
---|
| 215 | |
---|
| 216 | # default data |
---|
| 217 | if qz_value == None or not numpy.isfinite(qz_value): |
---|
| 218 | qz_value = 0 |
---|
[32e8c78] | 219 | if not numpy.isfinite(mask_value): |
---|
[3cd95c8] | 220 | mask_value = 1 |
---|
| 221 | q_data = numpy.append(q_data,numpy.sqrt(qx_value**2+qy_value**2+qz_value**2)) |
---|
[32e8c78] | 222 | # Note: For convenience, mask = False stands for masked, while mask = True for unmasked |
---|
[3cd95c8] | 223 | mask = numpy.append(mask,(mask_value>=1)) |
---|
| 224 | |
---|
[c7c5ef8] | 225 | f.close() |
---|
[32e8c78] | 226 | # If all mask elements are False, put all True |
---|
| 227 | if not mask.any(): mask[mask==False] = True |
---|
| 228 | |
---|
[3cd95c8] | 229 | # Store limits of the image in q space |
---|
| 230 | xmin = numpy.min(qx_data) |
---|
| 231 | xmax = numpy.max(qx_data) |
---|
| 232 | ymin = numpy.min(qy_data) |
---|
| 233 | ymax = numpy.max(qy_data) |
---|
| 234 | |
---|
| 235 | # units |
---|
| 236 | if has_converter == True and output.Q_unit != '1/A': |
---|
| 237 | xmin = data_conv_q(xmin, units=output.Q_unit) |
---|
| 238 | xmax = data_conv_q(xmax, units=output.Q_unit) |
---|
| 239 | ymin = data_conv_q(ymin, units=output.Q_unit) |
---|
| 240 | ymax = data_conv_q(ymax, units=output.Q_unit) |
---|
| 241 | |
---|
| 242 | ## calculate the range of the qx and qy_data |
---|
| 243 | x_size = math.fabs(xmax - xmin) |
---|
| 244 | y_size = math.fabs(ymax - ymin) |
---|
| 245 | |
---|
| 246 | # calculate the number of pixels in the each axes |
---|
| 247 | npix_y = math.floor(math.sqrt(len(data))) |
---|
| 248 | npix_x = math.floor(len(data)/npix_y) |
---|
| 249 | |
---|
| 250 | # calculate the size of bins |
---|
| 251 | xstep = x_size/(npix_x-1) |
---|
| 252 | ystep = y_size/(npix_y-1) |
---|
| 253 | |
---|
| 254 | # store x and y axis bin centers in q space |
---|
| 255 | x_bins = numpy.arange(xmin,xmax+xstep,xstep) |
---|
| 256 | y_bins = numpy.arange(ymin,ymax+ystep,ystep) |
---|
| 257 | |
---|
| 258 | # get the limits of q values |
---|
| 259 | xmin = xmin - xstep/2 |
---|
| 260 | xmax = xmax + xstep/2 |
---|
| 261 | ymin = ymin - ystep/2 |
---|
| 262 | ymax = ymax + ystep/2 |
---|
| 263 | |
---|
| 264 | #Store data in outputs |
---|
| 265 | #TODO: Check the lengths |
---|
| 266 | output.data = data |
---|
| 267 | output.err_data = numpy.sqrt(numpy.abs(data)) |
---|
| 268 | output.qx_data = qx_data |
---|
| 269 | output.qy_data = qy_data |
---|
| 270 | output.q_data = q_data |
---|
| 271 | output.mask = mask |
---|
| 272 | |
---|
| 273 | output.x_bins = x_bins |
---|
| 274 | output.y_bins = y_bins |
---|
| 275 | |
---|
| 276 | output.xmin = xmin |
---|
| 277 | output.xmax = xmax |
---|
| 278 | output.ymin = ymin |
---|
| 279 | output.ymax = ymax |
---|
| 280 | |
---|
| 281 | output.source.wavelength = wavelength |
---|
| 282 | |
---|
| 283 | # Store pixel size in mm |
---|
| 284 | detector.pixel_size.x = pixel_x |
---|
| 285 | detector.pixel_size.y = pixel_y |
---|
| 286 | |
---|
| 287 | # Store the sample to detector distance |
---|
| 288 | detector.distance = distance |
---|
| 289 | |
---|
| 290 | # optional data: if any of dq data == 0, do not pass to output |
---|
| 291 | if len(dqx_data) == len(qx_data): |
---|
| 292 | output.dqx_data = dqx_data |
---|
| 293 | if len(dqy_data) == len(qy_data): |
---|
| 294 | output.dqy_data = dqy_data |
---|
| 295 | |
---|
| 296 | # Units of axes |
---|
| 297 | if data_conv_q is not None: |
---|
| 298 | output.xaxis("\\rm{Q_{x}}", output.Q_unit) |
---|
| 299 | output.yaxis("\\rm{Q_{y}}", output.Q_unit) |
---|
| 300 | else: |
---|
| 301 | output.xaxis("\\rm{Q_{x}}", 'A^{-1}') |
---|
| 302 | output.yaxis("\\rm{Q_{y}}", 'A^{-1}') |
---|
| 303 | if data_conv_i is not None: |
---|
| 304 | output.zaxis("\\rm{Intensity}", output.I_unit) |
---|
| 305 | else: |
---|
| 306 | output.zaxis("\\rm{Intensity}","cm^{-1}") |
---|
| 307 | |
---|
| 308 | # Store loading process information |
---|
| 309 | output.meta_data['loader'] = self.type_name |
---|
| 310 | |
---|
| 311 | return output |
---|
| 312 | |
---|
| 313 | if __name__ == "__main__": |
---|
| 314 | reader = Reader() |
---|
| 315 | print reader.read("../test/exp18_14_igor_2dqxqy.dat") |
---|
| 316 | |
---|