source: trunk/imports/G2pwd_CIF.py @ 2347

Last change on this file since 2347 was 2347, checked in by vondreele, 7 years ago

implement Steven Weigand's fixes for Pilatus 'I' for 'L'
do a 'fix' for powderCif files with value(esd) for powder profile points
add a couple of comments for tutorial stuff
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 15.3 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2016-06-27 14:48:51 +0000 (Mon, 27 Jun 2016) $
4# $Author: vondreele $
5# $Revision: 2347 $
6# $URL: trunk/imports/G2pwd_CIF.py $
7# $Id: G2pwd_CIF.py 2347 2016-06-27 14:48:51Z vondreele $
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: 2347 $")
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            v2 = []
264            for val in cf[blk].get(ycf,'?'):
265                v,e = cif.get_number_with_esd(val)
266                if v is None: # not parsed
267                    vl.append(np.NaN)
268                    v2.append(np.NaN)
269                else:
270                    vl.append(v)
271                    v2.append(max(e,1.0))
272            y = np.array(vl)
273            w = 1./np.array(v2)**2
274            # weights
275            if sui == -1:
276                # no weights
277                vl = np.zeros(len(x)) + 1.
278            else:
279                vl = []
280                sucf = such[sui]
281                if sucf ==  '_pd_proc_ls_weight':
282                    for val in cf[blk].get(sucf,'?'):
283                        v,e = cif.get_number_with_esd(val)
284                        if v is None: # not parsed
285                            vl.append(0.)
286                        else:
287                            vl.append(v)
288                elif sucf ==  '_pd_meas_counts_total':
289                    for val in cf[blk].get(sucf,'?'):
290                        v,e = cif.get_number_with_esd(val)
291                        if v is None: # not parsed
292                            vl.append(0.)
293                        elif v <= 0:
294                            vl.append(1.)
295                        else:
296                            vl.append(1./v)
297                else:
298                    for val in cf[blk].get(sucf,'?'):
299                        v,e = cif.get_number_with_esd(val)
300                        if v is None or e is None: # not parsed or no ESD
301                            vl.append(np.NaN)
302                        elif e <= 0:
303                            vl.append(1.)
304                        else:
305                            vl.append(1./(e*e))
306#            w = np.array(vl)
307            # intensity modification factor
308            if modi >= 1:
309                ycf = modch[modi-1]
310                vl = []
311                for val in cf[blk].get(ycf,'?'):
312                    v,e = cif.get_number_with_esd(val)
313                    if v is None: # not parsed
314                        vl.append(np.NaN)
315                    else:
316                        vl.append(v)
317                y /= np.array(vl)
318                w /= np.array(vl)
319            N = len(x)
320        except Exception as detail:
321            self.errors = "Error reading from selected block"
322            self.errors += "\n  Read exception: "+str(detail)
323            print self.formatName+' read error:'+str(detail) # for testing
324            import traceback
325            traceback.print_exc(file=sys.stdout)
326            #self.errors += "\n  Traceback info:\n"+str(traceback.format_exc())
327            return False
328        finally:
329            self.DoneBusy()
330            print "CIF file, read from selected block"
331
332        self.errors = "Error while storing read values"
333        self.powderdata = [
334                np.array(x), # x-axis values
335                np.array(y), # powder pattern intensities
336                np.array(w), # 1/sig(intensity)^2 values (weights)
337                np.zeros(N), # calc. intensities (zero)
338                np.zeros(N), # calc. background (zero)
339                np.zeros(N), # obs-calc profiles
340            ]
341        self.powderentry[0] = filename
342        #self.powderentry[1] = pos # bank offset (N/A here)
343        self.powderentry[2] = 1 # xye file only has one bank
344        self.idstring = os.path.basename(filename) + ': ' + blk
345        if cf[blk].get('_diffrn_radiation_probe'):
346            if cf[blk]['_diffrn_radiation_probe'] == 'neutron':
347                self.instdict['type'] = 'PNC'
348                #if cf[blk].get('_pd_meas_time_of_flight'): self.instdict['type'] = 'PNT' # not supported yet
349            else:
350                self.instdict['type'] = 'PXC'
351        if cf[blk].get('_diffrn_radiation_wavelength'):
352            val = cf[blk]['_diffrn_radiation_wavelength']
353            wl = []
354            if type(val) is list:
355                for v in val:
356                    w,e = cif.get_number_with_esd(v)
357                    if w: wl.append(w)
358            else:
359                w,e = cif.get_number_with_esd(val)
360                if w: wl.append(w)
361            if wl: self.instdict['wave'] = wl
362        if cf[blk].get('_diffrn_ambient_temperature'):
363            val = cf[blk]['_diffrn_ambient_temperature']
364            w,e = cif.get_number_with_esd(val)
365            if w: 
366                self.Sample['Temperature'] = w
367        return True
Note: See TracBrowser for help on using the repository browser.