source: trunk/imports/G2pwd_CIF.py @ 2342

Last change on this file since 2342 was 2342, checked in by toby, 7 years ago

fix reading of CIF intensities with esds
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 15.1 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2016-06-24 20:36:47 +0000 (Fri, 24 Jun 2016) $
4# $Author: toby $
5# $Revision: 2342 $
6# $URL: trunk/imports/G2pwd_CIF.py $
7# $Id: G2pwd_CIF.py 2342 2016-06-24 20:36:47Z toby $
8########### SVN repository information ###################
9'''
10*Module G2pwd_CIF: CIF powder data*
11------------------------------------
12
13Routine to read in powder data from a CIF.
14
15'''
16import sys
17import numpy as np
18import os.path
19import GSASIIIO as G2IO
20import CifFile as cif # PyCifRW from James Hester
21import GSASIIpath
22GSASIIpath.SetVersionNumber("$Revision: 2342 $")
23
24class CIFpwdReader(G2IO.ImportPowderData):
25    'Routines to import powder data from a CIF file'
26    def __init__(self):
27        super(self.__class__,self).__init__( # fancy way to self-reference
28            extensionlist=('.CIF','.cif'),
29            strictExtension=False,
30            formatName = 'CIF',
31            longFormatName = 'Powder data from CIF'
32            )
33    # Validate the contents
34    def ContentsValidator(self, filepointer):
35        'Use standard CIF validator'
36        return self.CIFValidator(filepointer)
37
38    def Reader(self,filename,filepointer, ParentFrame=None, **kwarg):
39        '''Read powder data from a CIF.
40        If multiple datasets are requested, use self.repeat and buffer caching.
41        '''
42
43        # Define lists of data names used for holding powder diffraction data
44        # entries of a type that are not implemented are commented out in case
45        # we will want them later.
46        xDataItems = (  # "x-axis" data names"
47            ('_pd_meas_2theta_range_min', '_pd_meas_2theta_range_max', '_pd_meas_2theta_range_inc'),
48            ('_pd_proc_2theta_range_min', '_pd_proc_2theta_range_max', '_pd_proc_2theta_range_inc'),
49            '_pd_meas_2theta_scan',
50            '_pd_meas_time_of_flight',
51            '_pd_proc_2theta_corrected',
52            #'_pd_proc_d_spacing',
53            #'_pd_proc_energy_incident',
54            #'_pd_proc_energy_detection',
55            #'_pd_proc_recip_len_q',
56            #'_pd_proc_wavelength',
57        )
58        intDataItems = ( # intensity axis data names
59            '_pd_meas_counts_total',
60            #'_pd_meas_counts_background',
61            #'_pd_meas_counts_container',
62            '_pd_meas_intensity_total',
63            #'_pd_meas_intensity_background',
64            #'_pd_meas_intensity_container',
65            '_pd_proc_intensity_net',
66            '_pd_proc_intensity_total',
67            #'_pd_proc_intensity_bkg_calc',
68            #'_pd_proc_intensity_bkg_fix',
69            '_pd_calc_intensity_net',   # allow computed patterns as input data?
70            '_pd_calc_intensity_total',
71            )
72
73        ESDDataItems = ( # items that can be used to compute uncertainties for obs data
74            '_pd_proc_ls_weight',
75            '_pd_meas_counts_total'
76            )
77       
78        ModDataItems = ( # items that modify the use of the data
79            '_pd_meas_step_count_time',
80            '_pd_meas_counts_monitor',
81            '_pd_meas_intensity_monitor',
82            '_pd_proc_intensity_norm',
83            '_pd_proc_intensity_incident',
84            )
85        rdbuffer = kwarg.get('buffer')
86        cf = None
87        choicelist = None
88        selections = None
89        # reload previously saved values
90        if self.repeat and rdbuffer is not None:
91            cf = rdbuffer.get('lastcif')
92            choicelist = rdbuffer.get('choicelist')
93            print 'debug: Reuse previously parsed CIF'
94            selections = rdbuffer.get('selections')
95        if cf is None:
96            self.ShowBusy() # this can take a while
97            print "Starting parse of CIF file"
98            try:
99                cf = G2IO.ReadCIF(filename)
100            except Exception as detail:
101                self.errors = "Parse or reading of file failed in pyCifRW; check syntax of file in enCIFer or CheckCIF"
102                return False
103            finally:
104                self.DoneBusy()
105            print "CIF file parsed"
106        # scan all blocks for sets of data
107        if choicelist is None:
108            try:
109                choicelist = []
110                for blk in cf.keys():
111                    blkkeys = [k.lower() for k in cf[blk].keys()] # save a list of the data items, since we will use it often
112                    # scan through block for x items
113                    xldict = {}
114                    for x in xDataItems:
115                        if type(x) is tuple: # check for the presence of all three items that define a range of data
116                            if not all([i in blkkeys for i in x]): continue
117                            try:
118                                items = [float(cf[blk][xi]) for xi in x]
119                                l = 1 + int(0.5 + (items[1]-items[0])/items[2])
120                            except:
121                                continue
122                        else:
123                            if x not in blkkeys: continue
124                            l = len(cf[blk][x])
125                        if xldict.get(l) is None:
126                            xldict[l] = [x]
127                        else:
128                            xldict[l].append(x)
129                    # now look for matching intensity items
130                    yldict = {}
131                    suldict = {}
132                    for y in intDataItems:
133                        if y in blkkeys:
134                            l = len(cf[blk][y])
135                            if yldict.get(l) is None:
136                                yldict[l] = [y]
137                            else:
138                                yldict[l].append(y)
139                            # now check if the first item has an uncertainty
140                            if cif.get_number_with_esd(cf[blk][y][0])[1] is None: continue
141                            if suldict.get(l) is None:
142                                suldict[l] = [y]
143                            else:
144                                suldict[l].append(y)
145                    for y in ESDDataItems:
146                        if y in blkkeys:
147                            l = len(cf[blk][y])
148                            if suldict.get(l) is None:
149                                suldict[l] = [y]
150                            else:
151                                suldict[l].append(y)
152                    modldict = {}
153                    for y in ModDataItems:
154                        if y in blkkeys:
155                            l = len(cf[blk][y])
156                            if modldict.get(l) is None:
157                                modldict[l] = [y]
158                            else:
159                                modldict[l].append(y)
160                    for l in xldict:
161                        if yldict.get(l) is None: continue
162                        choicelist.append([blk,l,xldict[l],yldict[l],suldict.get(l,[]),modldict.get(l,[])])
163                        #print blk,l,xldict[l],yldict[l],suldict.get(l,[]),modldict.get(l,[])
164            except Exception as detail:
165                self.errors = "Error scanning blocks"
166                self.errors += "\n  Read exception: "+str(detail)
167                print self.formatName+' read error:'+str(detail) # for testing
168                import traceback
169                traceback.print_exc(file=sys.stdout)
170                #self.errors += "\n  Traceback info:\n"+str(traceback.format_exc())
171                return False
172            print "CIF file scanned for blocks with data"
173        if not choicelist:
174            selblk = None # no block to choose
175            self.errors = "No powder diffraction blocks found"
176            return False
177        elif len(choicelist) == 1: # only one choice
178            selblk = 0
179        elif self.repeat and selections is not None:
180            # we were called to repeat the read
181            print 'debug: repeat #',self.repeatcount,'selection',selections[self.repeatcount]
182            selblk = selections[self.repeatcount]
183            self.repeatcount += 1
184            if self.repeatcount >= len(selections): self.repeat = False
185        else:                       # choose from options
186            # compile a list of choices for the user
187            choices = []
188            for blk,l,x,y,su,mod in choicelist:
189                sx = x[0]
190                if len(x) > 1: sx += '...'
191                sy = y[0]
192                if len(y) > 1: sy += '...'
193                choices.append(
194                    'Block '+str(blk)+', '+str(l)+' points. X='+sx+' & Y='+sy
195                    )
196            selections = self.MultipleBlockSelector(
197                choices,
198                ParentFrame=ParentFrame,
199                title='Select dataset(s) to read from the list below',
200                size=(600,100),
201                header='Dataset Selector')
202            if len(selections) == 0:
203                self.errors = "Abort: block not selected"
204                return False
205            selblk = selections[0] # select first in list
206            if len(selections) > 1: # prepare to loop through again
207                self.repeat = True
208                self.repeatcount = 1
209                if rdbuffer is not None:
210                    rdbuffer['selections'] = selections
211                    rdbuffer['lastcif'] = cf # save the parsed cif for the next loop
212                    rdbuffer['choicelist'] = choicelist # save the parsed choices for the future
213
214        # got a selection, now read it
215        # do we need to ask which fields to read?
216        blk,l,xch,ych,such,modch = choicelist[selblk]
217        xi,yi,sui,modi = 0,0,0,0
218        if len(xch) > 1 or len(ych) > 1 or len(such) > 1 or len(modch) > 0:
219            choices = []
220            chlbls = []
221            chlbls.append('Select the scanned data item')
222            xchlst = []
223            for x in xch:
224                if type(x) is tuple:
225                    xchlst.append(x[0])
226                else:
227                    xchlst.append(x)
228            choices.append(xchlst)
229            chlbls.append('Select the intensity data item')
230            choices.append(ych)
231            chlbls.append('Select the data item to be used for weighting')
232            choices.append(such)
233            chlbls.append('Divide intensities by data item')
234            choices.append(['none']+modch)
235            res = self.MultipleChoicesDialog(choices,chlbls)
236            if not res:
237                self.errors = "Abort: data items not selected"
238                return False
239            xi,yi,sui,modi = res
240
241        # now read in the values
242        try:
243            self.ShowBusy() # this can also take a while
244            # x-values
245            xcf = xch[xi]
246            if type(xcf) is tuple:
247                vals = [float(cf[blk][xi]) for xi in xcf]
248                x = np.array(
249                    [(i*vals[2] + vals[0]) for i in range(1 + int(0.5 + (vals[1]-vals[0])/vals[2]))]
250                    )
251            else:
252                vl = []
253                for val in cf[blk].get(xcf,'?'):
254                    v,e = cif.get_number_with_esd(val)
255                    if v is None: # not parsed
256                        vl.append(np.NaN)
257                    else:
258                        vl.append(v)
259                x = np.array(vl)
260            # y-values
261            ycf = ych[yi]
262            vl = []
263            for val in cf[blk].get(ycf,'?'):
264                v,e = cif.get_number_with_esd(val)
265                if v is None: # not parsed
266                    vl.append(np.NaN)
267                else:
268                    vl.append(v)
269            y = np.array(vl)
270            # weights
271            if sui == -1:
272                # no weights
273                vl = np.zeros(len(x)) + 1.
274            else:
275                vl = []
276                sucf = such[sui]
277                if sucf ==  '_pd_proc_ls_weight':
278                    for val in cf[blk].get(sucf,'?'):
279                        v,e = cif.get_number_with_esd(val)
280                        if v is None: # not parsed
281                            vl.append(0.)
282                        else:
283                            vl.append(v)
284                elif sucf ==  '_pd_meas_counts_total':
285                    for val in cf[blk].get(sucf,'?'):
286                        v,e = cif.get_number_with_esd(val)
287                        if v is None: # not parsed
288                            vl.append(0.)
289                        elif v <= 0:
290                            vl.append(1.)
291                        else:
292                            vl.append(1./v)
293                else:
294                    for val in cf[blk].get(sucf,'?'):
295                        v,e = cif.get_number_with_esd(val)
296                        if v is None or e is None: # not parsed or no ESD
297                            vl.append(np.NaN)
298                        elif e <= 0:
299                            vl.append(1.)
300                        else:
301                            vl.append(1./(e*e))
302            w = np.array(vl)
303            # intensity modification factor
304            if modi >= 1:
305                ycf = modch[modi-1]
306                vl = []
307                for val in cf[blk].get(ycf,'?'):
308                    v,e = cif.get_number_with_esd(val)
309                    if v is None: # not parsed
310                        vl.append(np.NaN)
311                    else:
312                        vl.append(v)
313                y /= np.array(vl)
314                w /= np.array(vl)
315            N = len(x)
316        except Exception as detail:
317            self.errors = "Error reading from selected block"
318            self.errors += "\n  Read exception: "+str(detail)
319            print self.formatName+' read error:'+str(detail) # for testing
320            import traceback
321            traceback.print_exc(file=sys.stdout)
322            #self.errors += "\n  Traceback info:\n"+str(traceback.format_exc())
323            return False
324        finally:
325            self.DoneBusy()
326            print "CIF file, read from selected block"
327
328        self.errors = "Error while storing read values"
329        self.powderdata = [
330                np.array(x), # x-axis values
331                np.array(y), # powder pattern intensities
332                np.array(w), # 1/sig(intensity)^2 values (weights)
333                np.zeros(N), # calc. intensities (zero)
334                np.zeros(N), # calc. background (zero)
335                np.zeros(N), # obs-calc profiles
336            ]
337        self.powderentry[0] = filename
338        #self.powderentry[1] = pos # bank offset (N/A here)
339        self.powderentry[2] = 1 # xye file only has one bank
340        self.idstring = os.path.basename(filename) + ': ' + blk
341        if cf[blk].get('_diffrn_radiation_probe'):
342            if cf[blk]['_diffrn_radiation_probe'] == 'neutron':
343                self.instdict['type'] = 'PNC'
344                #if cf[blk].get('_pd_meas_time_of_flight'): self.instdict['type'] = 'PNT' # not supported yet
345            else:
346                self.instdict['type'] = 'PXC'
347        if cf[blk].get('_diffrn_radiation_wavelength'):
348            val = cf[blk]['_diffrn_radiation_wavelength']
349            wl = []
350            if type(val) is list:
351                for v in val:
352                    w,e = cif.get_number_with_esd(v)
353                    if w: wl.append(w)
354            else:
355                w,e = cif.get_number_with_esd(val)
356                if w: wl.append(w)
357            if wl: self.instdict['wave'] = wl
358        if cf[blk].get('_diffrn_ambient_temperature'):
359            val = cf[blk]['_diffrn_ambient_temperature']
360            w,e = cif.get_number_with_esd(val)
361            if w: 
362                self.Sample['Temperature'] = w
363        return True
Note: See TracBrowser for help on using the repository browser.