source: trunk/imports/G2pwd_fxye.py

Last change on this file was 5112, checked in by vondreele, 5 months ago

implement energy dispersive x-ray data as a new type 'PXE'; assume Gaussian peak shape only; no sample broadening considered. Peak resolution is only ~0.5% not good enough to see sample broadening
import, plotting, peak fitting, indexing & cell refinement from peak positions "complete"
data is old (from ~20 yrs ago) so no idea about modern ED data.
Start on use in Pawley/LeBail? refinement - Rietveld excluded; currently not working
Some work on ISODISTORT implementation

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 17.1 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2021-12-15 22:56:41 +0000 (Wed, 15 Dec 2021) $
4# $Author: toby $
5# $Revision: 5112 $
6# $URL: trunk/imports/G2pwd_fxye.py $
7# $Id: G2pwd_fxye.py 5112 2021-12-15 22:56:41Z toby $
8########### SVN repository information ###################
9'''
10*Module G2pwd_fxye: GSAS data files*
11------------------------------------
12Routine to read in powder data in a variety of formats
13that are defined for GSAS.
14
15'''
16from __future__ import division, print_function
17import os.path as ospath
18import platform
19import numpy as np
20import GSASIIobj as G2obj
21import GSASIIpath
22GSASIIpath.SetVersionNumber("$Revision: 5112 $")
23
24class GSAS_ReaderClass(G2obj.ImportPowderData):
25    'Routines to import powder data from a GSAS files'
26    def __init__(self):
27        super(self.__class__,self).__init__( # fancy way to self-reference
28            extensionlist=('.fxye','.raw','.gsas','.gda','.gsa','.gss','.RAW','.GSAS','.GDA','.GSA','.dat'),
29            strictExtension=False,
30            formatName = 'GSAS powder data',
31            longFormatName = 'GSAS powder data files (.fxye, .raw, .gsas...)'
32            )
33        self.clockWd = {}
34        self.TimeMap = {}
35        self.dnames = []
36        self.scriptable = True
37
38    # Validate the contents -- look for a bank line
39    def ContentsValidator(self, filename):
40        'Validrate by checking to see if the file has BANK lines & count them'
41        #print 'ContentsValidator: '+self.formatName
42        nBanks= 0
43        if '2' in platform.python_version_tuple()[0]:
44            fp = open(filename,'Ur')
45        else:
46            fp = open(filename,'r',encoding='latin-1')
47        fname = ospath.basename(fp.name)
48        for i,line in enumerate(fp):
49            self.GSAS = True
50            if i==0: # first line is always a comment
51                continue
52            if i==1 and line[:4].lower() == 'inst':
53                # 2nd line is optional instrument parameter file
54                continue
55            if line[0] == '#': continue
56            if line[:4] == 'BANK':
57                self.dnames.append(fname+' '+' '.join(line.split()[:2]))
58                nBanks += 1
59                continue
60            elif line[:7] == 'Monitor': continue
61            elif line [:8] == 'TIME_MAP':          #LANSCE TOF data
62                continue
63#            else:
64#                if not all(c in ' 0123456789.eE+-' for c in line[:-1]):
65#                    self.errors = 'Unexpected information in line: '+str(i+1)
66#                    self.errors += '  '+str(line)
67#                    return False
68        if nBanks:
69            if not len(self.selections):
70                self.selections = range(nBanks)
71            fp.close()
72            return True
73        self.errors = 'No BANK records found'
74        fp.close()
75        return False # no bank records
76
77    def Reader(self,filename, ParentFrame=None, **kwarg):
78        '''Read a GSAS (old formats) file of type FXY, FXYE, ESD or STD types.
79        If multiple datasets are requested, use self.repeat and buffer caching.
80        '''
81        def GetFXYEdata(File,Pos,Bank):
82            File.seek(Pos)
83            x = []
84            y = []
85            w = []
86            S = File.readline()
87            while S and S[:4] != 'BANK' and S[0] != '#':
88                vals = S.split()
89                x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
90                f = float(vals[1])
91                s = float(vals[2])
92                if f <= 0.0 or s <= 0.0:
93                    y.append(0.0)
94                    w.append(0.0)
95                else:
96                    y.append(float(vals[1]))
97                    w.append(1.0/float(vals[2])**2)
98                S = File.readline()
99            N = len(x)
100            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]   
101           
102        def GetFXYdata(File,Pos,Bank):
103            File.seek(Pos)
104            x = []
105            y = []
106            w = []
107            S = File.readline()
108            while S and S[:4] != 'BANK' and S[0] != '#':
109                vals = S.split()
110                x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
111                f = float(vals[1])
112                if f > 0.0:
113                    y.append(f)
114                    w.append(1.0/f)
115                else:             
116                    y.append(0.0)
117                    w.append(0.0)
118                S = File.readline()
119            N = len(x)
120            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
121           
122        def GetESDdata(File,Pos,Bank):
123            File.seek(Pos)
124            cons = Bank.split()
125            if 'TIME_MAP' == cons[4]:
126                start = 0
127                step = 1
128            else:
129                start = float(cons[5])/100.0               #CW: from centidegrees to degrees
130                step = float(cons[6])/100.0
131            x = []
132            y = []
133            w = []
134            S = File.readline()
135            j = 0
136            while S and S[:4] != 'BANK' and S[0] != '#':
137                if 'TIME_MAP' in S:
138                    break
139                for i in range(0,80,16):
140                    if S[i:i+8] == 8*' ':
141                        break
142                    xi = start+step*j
143                    yi = sfloat(S[i:i+8])
144                    ei = sfloat(S[i+8:i+16])
145                    x.append(xi)
146                    if yi > 0.0:
147                        y.append(yi)
148                        w.append(1.0/ei**2)
149                    else:             
150                        y.append(0.0)
151                        w.append(0.0)
152                    j += 1
153                S = File.readline()
154            N = len(x)
155            if self.clockWd:
156                if len(self.clockWd) == 1:
157                    x = Tmap2TOF(self.TimeMap['1'],self.clockWd['1'])
158                else:
159                    x = Tmap2TOF(self.TimeMap[cons[5]],self.clockWd[cons[5]])
160            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
161       
162        def GetSTDdata(File,Pos,Bank):
163            File.seek(Pos)
164            cons = Bank.split()
165            Nch = int(cons[2])
166            if self.clockWd:
167                start = 0
168                step = 1
169            else:
170                start = float(cons[5])/100.0               #CW: from centidegrees to degrees
171                step = float(cons[6])/100.0                 #NB TOF 0.1*ms!
172            x = []
173            y = []
174            w = []
175            S = File.readline()
176            j = 0
177            Ecoef = np.zeros(4)
178            if 'EDS' in cons[4]:
179               Ecoef = np.fromstring(' '.join(cons[5:9]),sep=' ')
180            while S and S[:4] != 'BANK' and S[0] != '#':
181                if 'TIME_MAP' in S or '\x1a' in S:
182                    break
183                for i in range(0,80,8):
184                    if S[i:i+10] == 10*' ':
185                        break
186                    xi = start+step*j
187                    if 'EDS' in Bank:   #Energy dispersive - convert CN to EkeV
188                        xi = Ecoef[0]+(Ecoef[1]+(Ecoef[2]+Ecoef[3]*j)*j)*j
189                    ni = max(sint(S[i:i+2]),1)
190                    yi = max(sfloat(S[i+2:i+8]),0.0)
191                    if yi:
192                        vi = yi/ni
193                    else:
194                        yi = 0.0
195                        vi = 0.0
196                    j += 1
197                    if j < Nch:
198                        x.append(xi)
199                        if vi <= 0.:
200                            y.append(0.)
201                            w.append(0.)
202                        else:
203                            y.append(yi)
204                            w.append(1.0/vi)
205                S = File.readline()
206            N = len(x)
207            if self.clockWd:
208                if len(self.clockWd) == 1:
209                    x = Tmap2TOF(self.TimeMap['1'],self.clockWd['1'])
210                else:
211                    x = Tmap2TOF(self.TimeMap[cons[5]],self.clockWd[cons[5]])
212            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
213           
214        def GetALTdata(File,Pos,Bank):
215            File.seek(Pos)
216            cons = Bank.split()
217            x = []
218            y = []
219            w = []
220            S = File.readline()
221            j = 0
222            while S and S[:4] != 'BANK' and S[0] != '#':
223                for i in range(0,80,20):
224                    if S[i:i+8] == 8*' ':
225                        break
226                    xi = sfloat(S[i:i+8])/3200.
227                    if not xi:
228                        break
229                    yi = sfloat(S[i+8:i+15])/1000.
230                    ei = sfloat(S[i+15:i+21])/1000.
231                    x.append(xi)
232                    if yi > 0.0:
233                        y.append(yi)
234                        w.append(1.0/ei**2)
235                    else:             
236                        y.append(0.0)
237                        w.append(0.0)
238                    j += 1
239                S = File.readline()
240            N = len(x)
241            if self.clockWd:
242                x = Tmap2TOF(self.TimeMap[cons[5]],self.clockWd[cons[5]])
243            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
244           
245        def GetTimeMap(File,Pos,TimeMap):
246            File.seek(Pos)
247            cons = TimeMap[8:].split()
248            mapNo = cons[0]
249            if 'TIME_MAP10' in TimeMap: mapNo = '1'   #HIPD cluge!
250            Nch = int(cons[1])
251            Nrec = int(cons[2])
252            clockWd = float(cons[4])/1000.          #in mus
253            TMap = np.zeros(Nch+2,dtype=int)
254            ind = 0
255            for i in range(Nrec):
256                S = File.readline().rstrip('\n')
257                vals = S.split()
258                for val in vals:
259                    TMap[ind] = int(val)
260                    ind += 1
261            TMap = np.reshape(TMap,(-1,3))
262            TMax = TMap[-1][0]
263            Nch = TMap[-2][0]+(TMax-TMap[-2][1]+TMap[-2][2]-1)/TMap[-2][2]
264            TMap[-1] = [Nch+1,TMax,0]
265            TMap = TMap.T
266            TMap[0] -= 1
267            return TMap.T,clockWd,mapNo
268           
269        def Tmap2TOF(TMap,clockWd):
270            TOF = []
271            Tch,T,Step = TMap[0]
272            for tmap in TMap[1:]:
273                tch,t,step = tmap
274                TOF += [T+Step*(i-Tch) for i in range(Tch,tch)]
275                Tch,T,Step = tmap
276            TOF = np.array(TOF)*clockWd
277            return TOF
278
279        Banks = []
280        Pos = []
281        rdbuffer = kwarg.get('buffer')
282        title = ''
283        comments = None
284        if '2' in platform.python_version_tuple()[0]:
285            fp = open(filename,'Ur')
286        else:
287            fp = open(filename,'r',encoding='latin-1')
288        # reload previously saved values - used for multibank reads
289        if self.repeat and rdbuffer is not None:
290            Banks = rdbuffer.get('Banks')
291            Pos = rdbuffer.get('Pos')
292            self.selections = rdbuffer.get('selections')
293            comments = rdbuffer.get('comments')
294
295        # read through the file and find the beginning of each bank
296        # Save the offset (Pos), BANK line (Banks), comments for each bank
297        #
298        # This is going to need a fair amount of work to track line numbers
299        # in the input file.
300        if len(Banks) != len(Pos) or len(Banks) == 0:
301            i = -1
302            while True:
303                i += 1
304                S = fp.readline()
305                if len(S) == 0: break
306                   
307                if i==0: # first line is always a comment
308                    self.errors = 'Error reading title'
309                    title = S[:-1]
310                    comments = [[title,]]
311                    continue
312                if i==1 and S[:4].lower() == 'inst' and ':' in S:
313                    # 2nd line is instrument parameter file (optional)
314                    self.errors = 'Error reading instrument parameter filename'
315                    self.instparm = S.split(':')[1].strip('[]').strip()
316                    continue
317                if S[0] == '#': # allow comments anywhere in the file
318                    # comments in fact should only preceed BANK lines
319                    comments[-1].append(S[:-1])
320                    continue
321                if S[:4] == 'BANK':
322                    self.errors = 'Error reading bank:'
323                    self.errors += '  '+str(S)
324                    comments.append([title,])
325                    Banks.append(S)
326                    Pos.append(fp.tell())
327                if S[:8] == 'TIME_MAP':     #assumes one time map; HIPPO has multiple time maps
328                    if len(Banks) == 0:
329                        self.errors = 'Error reading time map before any bank lines'
330                    else:
331                        self.errors = 'Error reading time map after bank:\n  '+str(Banks[-1])
332                    timemap,clockwd,mapNo = GetTimeMap(fp,fp.tell(),S)
333                    self.TimeMap[mapNo] = timemap
334                    self.clockWd[mapNo] = clockwd
335                   
336
337        # Now select the bank to read
338        if not Banks: # use of ContentsValidator should prevent this error
339            print (self.formatName+' scan error: no BANK records')
340            selblk = None # no block to choose
341            self.errors = 'No BANK records found (strange!)'
342            return False
343        elif len(Banks) == 1: # only one Bank, don't ask
344            selblk = 0
345        elif self.repeat and self.selections is not None:
346            # we were called to repeat the read
347            #print 'debug: repeat #',self.repeatcount,'selection',self.selections[self.repeatcount]
348            selblk = self.selections[self.repeatcount]
349            self.repeatcount += 1
350            if self.repeatcount >= len(self.selections): self.repeat = False
351        else:                       # choose from options
352            selblk = self.selections[0] # select first in list
353            if len(self.selections) > 1: # prepare to loop through again
354                self.repeat = True
355                self.repeatcount = 1
356                if rdbuffer is not None:
357                    rdbuffer['Banks'] = Banks
358                    rdbuffer['Pos'] = Pos
359                    rdbuffer['selections'] = self.selections
360                    rdbuffer['comments'] = comments
361
362        # got a selection, now read it
363        Bank = Banks[selblk]
364        try:
365            bnkNo = int(Bank.split()[1])
366        except ValueError:
367            bnkNo = 1
368        if 'FXYE' in Bank:
369            self.errors = 'Error reading FXYE data in Bank\n  '+Banks[selblk]
370            self.powderdata = GetFXYEdata(fp,Pos[selblk],Banks[selblk])
371        elif 'FXY' in Bank:
372            self.errors = 'Error reading FXY data in Bank\n  '+Banks[selblk]
373            self.powderdata = GetFXYdata(fp,Pos[selblk],Banks[selblk])
374        elif 'ESD' in Bank:
375            self.errors = 'Error reading ESD data in Bank\n  '+Banks[selblk]
376            self.powderdata = GetESDdata(fp,Pos[selblk],Banks[selblk])
377        elif 'STD' in Bank:
378            self.errors = 'Error reading STD data in Bank\n  '+Banks[selblk]
379            self.powderdata = GetSTDdata(fp,Pos[selblk],Banks[selblk])
380        elif 'ALT' in Bank:
381            self.errors = 'Error reading ALT data in Bank\n  '+Banks[selblk]
382            self.powderdata = GetALTdata(fp,Pos[selblk],Banks[selblk])
383        else:
384            self.errors = 'Error reading STD data in Bank\n  '+Banks[selblk]
385            self.powderdata = GetSTDdata(fp,Pos[selblk],Banks[selblk])
386
387        self.errors = 'Error processing information after read complete'
388        if comments is not None:
389            self.comments = comments[selblk]
390        self.powderentry[0] = filename
391        self.powderentry[1] = Pos # position offset (never used, I hope)
392        self.powderentry[2] = bnkNo #selblk+1 # bank number
393        self.idstring = ospath.basename(filename) + ' Bank '+str(bnkNo) #selblk+1)
394        self.numbanks=len(Banks)
395        # scan comments for temperature & radius
396        Temperature = 300.
397        for S in self.comments:
398            if 'Temp' in S.split('=')[0]:
399                try:
400                    Temperature = float(S.split('=')[1])
401                except:
402                    pass
403            elif 'Gonio' in S.split('=')[0]:
404                try:
405                    self.Sample['Gonio. radius'] = float(S.split('=')[1])
406                except:
407                    pass
408            elif 'Omega' in S.split('=')[0] or 'Theta' in S.split('=')[0]:  #HIPD weirdness
409                try:
410                    self.Sample['Omega'] = float(S.split('=')[1])
411                except:
412                    pass
413            elif 'Chi' in S.split('=')[0]:
414                try:
415                    self.Sample['Chi'] = float(S.split('=')[1])
416                except:
417                    pass                   
418            elif 'Phi' in S.split('=')[0]:
419                try:
420                    self.Sample['Phi'] = float(S.split('=')[1])
421                except:
422                    pass                   
423        self.Sample['Temperature'] = Temperature
424        fp.close()
425        return True       
426
427def sfloat(S):
428    'convert a string to a float, treating an all-blank string as zero'
429    if S.strip():
430        return float(S)
431    else:
432        return 0.0
433
434def sint(S):
435    'convert a string to an integer, treating an all-blank string as zero'
436    if S.strip():
437        return int(S)
438    else:
439        return 0
440
Note: See TracBrowser for help on using the repository browser.