Source Code for Module reflectometry.reduction.icpformat

  1  # This program is public domain 
  2   
  3  # Author: Paul Kienzle 
  4  # Initial version: William Ratcliff 
  5   
  6  """ 
  7  ICP data reader. 
  8   
  9  summary(filename)  - reads the header information 
 10  read(filename) - reads header information and data 
 11  """ 
 12   
 13  import numpy as N 
 14  import datetime,sys 
 15  from reflectometry.reduction import _reduction 
 16   
 17 -def readdata(fh): 
 18      """ 
 19      Read ICP data, including PSD data if lines contain commas. 
 20      """ 
 21      rows = [] 
 22      blocks = [] 
 23   
 24      line = fh.readline().rstrip() 
 25      linenum = 1 
 26      while line != '': 
 27          # While it might be easy to check for a comment mark on the beginning 
 28          # of the line, supporting this is ill-adviced.  First, users should 
 29          # be strongly discouraged from modifying the original data. 
 30          # Second, sequencing against the automatically generated motor 
 31          # columns will become more complicated.  Let's make life easier 
 32          # and put the masking in the application rather than the data reader. 
 33   
 34          # Process the instrument configuration line and move to the next line 
 35          rows.append([float(val) for val in line.split()]) 
 36          line = fh.readline().rstrip() 
 37          linenum += 1 
 38   
 39          # Build up a multiline detector block by joining all lines that 
 40          # contain a comma 
 41          b = [] 
 42          while ',' in line: 
 43              b.append(line) 
 44              line = fh.readline() 
 45              linenum += 1 
 46   
 47          # If last line ended with a comma then the last number for the 
 48          # the current block is on the current line. 
 49          if b != [] and b[-1].rstrip()[-1] == ",": 
 50              b.append(line) 
 51              line = fh.readline() 
 52              linenum += 1 
 53   
 54          if b != []: 
 55              # Have a detector block so add it 
 56              s = "".join(b) 
 57              # You can use numpy's string to matrix converter if for some 
 58              # reason _reduction doesn't compile, but it is horribly slow. 
 59              #z = N.matrix(s,'i').A 
 60              if blocks != []: 
 61                  # Have an existing block, so we know what size to allocate 
 62                  z = N.empty(blocks[0].shape,'i') 
 63                  i,j = _reduction.str2imat(s,z) 
 64                  if i*j != z.size: 
 65                      raise IOError,"Inconsistent dims at line %d"%linenum 
 66              else: 
 67                  # No existing block.  Worst case is 2 bytes per int. 
 68                  n = int(len(s)/2+1) 
 69                  z = N.empty(n,'i') 
 70                  i,j = _reduction.str2imat(s,z) 
 71                  # Keep the actual size 
 72                  if i==1 or j==1: 
 73                      z = z[:i*j].reshape(i*j) 
 74                  else: 
 75                      z = z[:i*j].reshape(i,j) 
 76              blocks.append(z) 
 77   
 78          elif blocks != []: 
 79              # Oops...missing a detector block.  Set it to zero counts 
 80              # of the same size as the last block 
 81              blocks.append(N.zeros(blocks[-1].shape,'i')) 
 82          # Otherwise no detector block and don't need one 
 83          # Note: this strategy fails to identify that the first 
 84          # detector block is missing; those will be filled in later. 
 85   
 86      # recover from missing leading detector blocks 
 87      if blocks != [] and len(blocks) < len(rows): 
 88          blank = N.zeros(blocks[0].shape,'i') 
 89          blocks = [blank]*(len(blocks)-len(rows)) + blocks 
 90   
 91      # Convert data to arrays 
 92      X = N.array(rows, 'd') 
 93      Z = N.array(blocks) 
 94      return X,Z 
 95   
 96   
 97 -def get_tokenized_line(file): 
 98      """ 
 99      Read the next line of text into a set of words. 
100      """ 
101      line=file.readline() 
102      return line.split() 
103   
104 -def get_quoted_tokens(file): 
105      """ 
106      Build a token list from a line which can be a mix of quoted strings 
107      and unquoted values separated by spaces.  Uses single quotes only. 
108      Does not test for escaped single quotes. 
109      """ 
110      line = file.readline() 
111      tokens = [] 
112      curtoken=None 
113      inquote = False 
114      for c in line: 
115          if c == "'": 
116              if inquote: 
117                  tokens.append("".join(curtoken)) 
118                  curtoken = None 
119                  inquote = False 
120              else: 
121                  curtoken = [] 
122                  inquote = True 
123          elif inquote: 
124              curtoken.append(c) 
125          elif c.isspace(): 
126              if curtoken != None: 
127                  tokens.append("".join(curtoken)) 
128                  curtoken = None 
129          else: 
130              if curtoken == None: 
131                  curtoken = [c] 
132              else: 
133                  curtoken.append(c) 
134   
135      return tokens 
136   
137 -class Lattice(object): pass 
138 -class Motor(object): pass 
139 -class MotorSet(object): pass 
140 -class ColumnSet(object): 
141 -    def __getitem__(self, k): 
142          return getattr(self,k) 
143   
144 -class ICP(object): 
145 -    def __init__(self, path): 
146          self.path = path 
147   
148 -    def readheader1(self): 
149          """ 
150          Read the tow line summary at the start of the ICP data files. 
151          """ 
152   
153          file = self.file 
154          tokens = get_quoted_tokens(file) 
155          self.filename=tokens[0] 
156          self.timestamp = datetime.datetime(2000,1,1) 
157          self.date=self.timestamp.strptime(tokens[1],'%b %d %Y %H:%M') 
158          self.scantype = tokens[2] 
159          self.prefactor = float(tokens[3]) 
160          self.monitor=float(tokens[4]) 
161          self.count_type=tokens[5] 
162          self.points=int(tokens[6]) 
163          self.data_type=tokens[7] 
164   
165          #skip over names of fields 
166          file.readline() 
167   
168          #comment and polarization 
169          line = file.readline() 
170          polarized_index = line.find("F1: O", 52) 
171          if polarized_index > 0: 
172              self.comment = line[:polarized_index].rstrip() 
173              F1 = '+' if line.find("F1: ON", 52)>0 else '-' 
174              F2 = '+' if line.find("F2: ON", 52)>0 else '-' 
175              self.polarization = F1+F2 
176          else: 
177              self.comment = line.rstrip() 
178              self.polarization = "" 
179   
180   
181 -    def readiheader(self): 
182          """ 
183          Read I-buffer structure, excluding motors. 
184          """ 
185          file = self.file 
186   
187          # Read in fields and field names 
188          tokenized=get_tokenized_line(file) 
189          fieldnames = file.readline() 
190          #print tokenized 
191          #print fieldnames 
192   
193          #Collimation    Mosaic    Wavelength   T-Start   Incr.   H-field #Det 
194          self.collimations = [float(s) for s in tokenized[0:4]] 
195          self.mosaic = [float(s) for s in tokenized[4:7]] 
196          self.wavelength=float(tokenized[7]) 
197          self.Tstart=float(tokenized[8]) 
198          self.Tstep=float(tokenized[9]) 
199          self.Hfield=float(tokenized[10]) 
200   
201 -    def readrheader(self): 
202          """ 
203          Read R-buffer structure, excluding motors. 
204          """ 
205          file = self.file 
206          # Read in fields and field names 
207          tokenized=get_tokenized_line(file) 
208          fieldnames = file.readline() 
209          #print tokenized 
210          #print fieldnames 
211   
212          #Mon1    Exp   Dm      Wavel  T-Start  Incr. Hf(Tesla) #Det SclFac 
213          self.Mon1=float(tokenized[0]) 
214          self.Exp=float(tokenized[1]) 
215          self.Dm=float(tokenized[2]) 
216          self.wavelength=float(tokenized[3]) 
217          self.Tstart=float(tokenized[4]) 
218          self.Tstep=float(tokenized[5]) 
219          self.Hfield=float(tokenized[6]) 
220          self.numDet=float(tokenized[7]) 
221          self.SclFac=float(tokenized[8]) 
222   
223 -    def readqheader(self): 
224          """ 
225          Read Q-buffer structure. 
226          """ 
227          file = self.file 
228          #experiment info 
229          tokenized=get_tokenized_line(file) 
230          self.collimations=[float(s) for s in tokenized[0:4]] 
231          self.mosaic=[float(s) for s in tokenized[4:7]] 
232          orient1=[float(s) for s in tokenized[7:10]] 
233          #ignore the "angle" field 
234          orient2=[float(s) for s in tokenized[11:14]] 
235          #skip line with field names 
236          file.readline() 
237          tokenized=get_tokenized_line(file) 
238          lattice=Lattice() 
239          lattice.a=float(tokenized[0]) 
240          lattice.b=float(tokenized[1]) 
241          lattice.c=float(tokenized[2]) 
242          lattice.alpha=float(tokenized[3]) 
243          lattice.beta=float(tokenized[4]) 
244          lattice.gamma=float(tokenized[5]) 
245          self.lattice=lattice 
246          #skip line with field names 
247          file.readline() 
248          tokenized=get_tokenized_line(file) 
249          self.ecenter=float(tokenized[0]) 
250          self.deltae=float(tokenized[1]) 
251          self.ef=float(tokenized[2]) 
252          self.monochromator_dspacing=float(tokenized[3]) 
253          self.analyzer_dspacing=float(tokenized[4]) 
254          self.tstart=float(tokenized[5]) 
255          self.tstep=float(tokenized[6]) 
256          tokenized=get_tokenized_line(file) 
257          self.Efixed=tokenized[4] 
258          tokenized=get_tokenized_line(file) 
259          self.qcenter=[float(s) for s in tokenized[0:3]] 
260          self.qstep=[float(s) for s in tokenized[3:6]] 
261          self.hfield=float(tokenized[6]) 
262          #skip line describing fields 
263          file.readline() 
264   
265 -    def check_wavelength(self, default, overrides): 
266          """ 
267          ICP sometimes records the incorrect wavelength in the file.  Make 
268          sure the right value is being used.  Be annoying about it so that 
269          if the wavelength was changed for a legitimate reason the user can 
270          override.  L is the value in the file.  dectector.wavelength should 
271          already be set to the default for the instrument. 
272          """ 
273          dataset = self.filename[:5] 
274          wavelength = self.wavelength 
275          if dataset in overrides: 
276              # yuck! If already overridden for a particular file in 
277              # a dataset, override for all files in the dataset. 
278              wavelength = overrides[dataset] 
279              message("Using wavelength %s for %s"%(wavelength,dataset)) 
280          elif wavelength == 0: 
281              # yuck! If stored value is 0, use the default 
282              wavelength = default 
283              message("Using default wavelength %s for %s"\ 
284                      %(wavelength,self.path)) 
285          elif abs(default-wavelength)/default > 0.01: 
286              # yuck! Value differs significantly from the default 
287              if question("ICP recorded a wavelength of %s in %s. \ 
288      Do you want to use the default wavelength %s instead?"\ 
289                %(wavelength,self.path,default)): 
290                  wavelength = default 
291          # Regardless of how the value was obtained, use that value for 
292          # the entire dataset 
293          return wavelength 
294   
295   
296 -    def readmotors(self): 
297          """ 
298          Read the 6 motor lines, returning a dictionary of 
299          motor names and start-step-stop values. 
300          E.g., 
301   
302          M = _readmotors(file) 
303          print M['a1'].start 
304          """ 
305          self.motor = MotorSet() 
306          while True:  # read until 'Mot:' line 
307              words=get_tokenized_line(self.file) 
308              if words[0] == 'Mot:': break 
309              motor = Motor() 
310              motor.start=float(words[1]) 
311              motor.step=float(words[2]) 
312              motor.stop=float(words[3]) 
313              name = words[0] if not words[0].isdigit() else 'a'+words[0] 
314              setattr(self.motor,name,motor) 
315   
316 -    def readcolumnheaders(self): 
317          """ 
318          Get a list of column names. Transform the names of certain 
319          columns to make our lives easier elsewhere: 
320                #1 COUNTS -> counts 
321                #2 COUNTS -> counts2 
322                MON -> monitor 
323                MIN -> time 
324                Q(x) -> qx, Q(y) -> qy, Q(z) -> qz 
325          All column names are uppercase. 
326          """ 
327          line = self.file.readline() 
328          line = line.lower() 
329          for (old,new) in (('#1 counts','counts'), 
330                            ('#2 counts','counts2'), 
331                            (' mon ',' monitor '), 
332                            (' min ',' time '), 
333                            ('(',''), 
334                            (')',''), 
335                            ): 
336              line = line.replace(old,new) 
337          self.columnnames = line.split() 
338   
339   
340 -    def readcolumns(self): 
341          ''' 
342          Read and parse ICP data columns listed in columns.  Return a dict of 
343          column name: vector.  If using a position sensitive detector, return 
344          an array of detector values x scan points. 
345          ''' 
346          values,detector = readdata(self.file) 
347          self.column = ColumnSet() 
348          for (c,v) in zip(self.columnnames,values.T): 
349              setattr(self.column,c,v) 
350          self.detector = detector 
351          self.counts = detector if detector.size > 0 else self.column.counts 
352          self.points = len(self.column.counts) 
353   
354 -    def genmotorcolumns(self): 
355          """ 
356          Generate vectors for each of the motors if a vector is not 
357          already stored in the file. 
358          """ 
359          for (M,R) in self.motor.__dict__.iteritems(): 
360              if not hasattr(self.column,M): 
361                  if R.step != 0.: 
362                      vector = N.arange(R.start,R.step,R.stop) 
363                      # truncate to number of points measured 
364                      vector = vector[:self.points]+0 
365                  else: 
366                      vector = R.start * N.ones(self.points) 
367                  setattr(self.column,M,vector) 
368          pass 
369   
370 -    def parseheader(self): 
371          """ 
372          Read and parse ICP header information 
373          """ 
374          # Determine FileType 
375          file = self.file 
376          self.readheader1() 
377          if self.scantype=='I': 
378              self.readiheader() 
379              self.readmotors() 
380          elif self.scantype=='Q': 
381              self.readqheader() 
382          elif self.scantype=='B': 
383              self.readqheader() 
384              self.readmotors() 
385          elif self.scantype=='R': 
386              self.readrheader() 
387              self.readmotors() 
388          else: 
389              raise ValueError, "Unknown scantype %s in ICP file"%scantype 
390          self.readcolumnheaders() 
391   
392 -    def summary(self): 
393          """ 
394          Read header from file, returning a dict of fields. 
395          """ 
396          self.file = gzopen(self.path) 
397          self.parseheader() 
398          data1 = self.file.readline() 
399          data2 = self.file.readline() 
400          self.PSD = (',' in data2) 
401          self.file.close() 
402   
403 -    def read(self): 
404          """ 
405          Read header and data from file, returning a dict of fields. 
406          """ 
407          self.file = gzopen(self.path) 
408          self.parseheader() 
409   
410          #read columns and detector images if available 
411          self.readcolumns() 
412          self.PSD = (self.detector.size>0) 
413   
414          # fill in missing motor columns 
415          self.genmotorcolumns() 
416   
417          self.file.close() 
418   
419 -    def __contains__(self, column): 
420          return hasattr(self.column,column) 
421   
422 -    def counts(self): 
423          if self.detector.size > 1: 
424              return self.detector 
425          else: 
426              return self.column.counts 
427   
428 -def read(filename): 
429      """Read an ICP file and return the corresponding ICP file object""" 
430      icp = ICP(filename) 
431      icp.read() 
432      return icp 
433   
434 -def summary(filename): 
435      """Read an ICP file header and return the corresponding ICP file object""" 
436      icp = ICP(filename) 
437      icp.summary() 
438      return icp 
439   
440 -def gzopen(filename,mode='r'): 
441      """ 
442      Open file or gzip file 
443      """ 
444      if filename.endswith('.gz'): 
445          import gzip 
446          file = gzip.open(filename, mode) 
447      else: 
448          file = open(filename, mode) 
449      return file 
450   
451 -def asdata(icp): 
452      import numpy 
453      import data 
454      d = data.Data() 
455      d.vlabel = 'Counts' 
456      d.v = icp.counts 
457      d.xlabel = icp.columnnames[0].capitalize() 
458      d.x = icp.column[icp.columnnames[0]] 
459      if len(d.v.shape) > 1: 
460          d.ylabel = 'Pixel' 
461          d.y = numpy.arange(d.v.shape[0]) 
462      return d 
463   
464 -def data(filename): 
465      icp = ICP(filename) 
466      icp.read() 
467      return asdata(icp) 
468   
469  # TODO: need message/question functions 
470 -def message(text): pass 
471 -def question(text): return True 
472   
473 -def demo(): 
474      """ 
475      Read and print all command line arguments 
476      """ 
477      import sys 
478      if len(sys.argv) < 2: 
479          print "usage: python icpformat.py file*" 
480      for file in sys.argv[1:]: 
481          fields = read(file) 
482          keys = fields.__dict__.keys() 
483          keys.sort() 
484          for k in keys: print k,getattr(fields,k) 
485   
486 -def plot(filename): 
487      """ 
488      Read and print all command line arguments 
489      """ 
490      import pylab 
491   
492      canvas = pylab.gcf().canvas 
493      d = data(filename) 
494      if len(d.v.shape) > 1: 
495          pylab.gca().pcolorfast(d.xedges,d.yedges,d.v) 
496          pylab.xlabel(d.xlabel) 
497          pylab.ylabel(d.ylabel) 
498      else: 
499          pylab.plot(d.x,d.v) 
500          pylab.xlabel(d.xlabel) 
501          pylab.ylabel(d.vlabel) 
502      pylab.show() 
503   
504 -def plot_demo(): 
505      import sys 
506      if len(sys.argv) != 2: 
507          print "usage: python icpformat.py file" 
508      else: 
509          plot(sys.argv[1]) 
510   
511  if __name__=='__main__': 
512      plot_demo() 
513      #demo() 
514