Changeset 4088


Ignore:
Timestamp:
Aug 11, 2019 10:34:43 PM (2 years ago)
Author:
toby
Message:

Add scripting peak retrieval & document

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIscriptable.py

    r4053 r4088  
    184184:meth:`G2PwdrData.fit_fixed_points`                   Fits background to the specified fixed points.
    185185:meth:`G2PwdrData.getdata`                            Provides access to the diffraction data associated with the histogram.
    186 :meth:`G2PwdrData.Export`                             Writes the diffraction data into a file
     186:meth:`G2PwdrData.Export`                             Writes the diffraction data into a file
     187:meth:`G2PwdrData.add_peak`                           Adds a peak to the peak list. Also see :ref:`PeakRefine`.
     188:meth:`G2PwdrData.set_peakFlags`                      Sets refinement flags for peaks
     189:meth:`G2PwdrData.refine_peaks`                       Starts a peak/background fitting cycle
     190:attr:`G2PwdrData.Peaks`                              Provides access to the peak list data structure
     191:attr:`G2PwdrData.PeakList`                           Provides the peak list parameter values
     192:meth:`G2PwdrData.Export_peaks`                       Writes the peak parameters to a text file
    187193==================================================    ===============================================================================================================
    188194
     
    678684
    679685Note that the parameters must match the object type and method (phase vs. histogram vs. HAP).
     686
     687.. _PeakRefine:
     688
     689=================================
     690Peak Refinement
     691=================================
     692Peak refinement is performed with routines
     693:meth:`G2PwdrData.add_peak`, :meth:`G2PwdrData.set_peakFlags` and
     694:meth:`G2PwdrData.refine_peaks`. Method :meth:`G2PwdrData.Export_peaks` and
     695properties :attr:`G2PwdrData.Peaks` and :attr:`G2PwdrData.PeakList`
     696provide ways to access the results. Note that when peak parameters are
     697refined with :meth:`~G2PwdrData.refine_peaks`, the background may also
     698be refined. Use :meth:`G2PwdrData.set_refinements` to change background
     699settings and the range of data used in the fit. See below for an example
     700peak refinement script, where the data files are taken from the
     701"Rietveld refinement with CuKa lab Bragg-Brentano powder data" tutorial
     702(in https://subversion.xray.aps.anl.gov/pyGSAS/Tutorials/LabData/data/).
     703
     704.. code-block::  python
     705
     706    from __future__ import division, print_function
     707    import os,sys
     708    sys.path.insert(0,'/Users/toby/software/G2/GSASII') # needed to "find" GSAS-II modules
     709    import GSASIIscriptable as G2sc
     710    datadir = os.path.expanduser("~/Scratch/peakfit")
     711    PathWrap = lambda fil: os.path.join(datadir,fil)
     712    gpx = G2sc.G2Project(newgpx=PathWrap('pkfit.gpx'))
     713    hist = gpx.add_powder_histogram(PathWrap('FAP.XRA'), PathWrap('INST_XRY.PRM'),
     714                                    fmthint='GSAS powder')
     715    hist.set_refinements({'Limits': [16.,24.],
     716          'Background': {"no. coeffs": 2,'type': 'chebyschev', 'refine': True}
     717                         })
     718    peak1 = hist.add_peak(1, ttheta=16.8)
     719    peak2 = hist.add_peak(1, ttheta=18.9)
     720    peak3 = hist.add_peak(1, ttheta=21.8)
     721    peak4 = hist.add_peak(1, ttheta=22.9)
     722    hist.set_peakFlags(area=True)
     723    hist.refine_peaks()
     724    hist.set_peakFlags(area=True,pos=True)
     725    hist.refine_peaks()
     726    hist.set_peakFlags(area=True, pos=True, sig=True, gam=True)
     727    hist.refine_peaks()
     728    print('peak positions: ',[i[0] for i in hist.PeakList])
     729    for i in range(len(hist.Peaks['peaks'])):
     730        print('peak',i,'pos=',hist.Peaks['peaks'][i][0],'sig=',hist.Peaks['sigDict']['pos'+str(i)])
     731    hist.Export_peaks('pkfit.txt')
     732    #gpx.save()  # gpx file is not written without this
    680733
    681734.. _CommandlineInterface:
     
    30883141        '''
    30893142        if extension not in exportersByExtension.get('powder',[]):
    3090             print('Known exporters are',exportersByExtension.get('powder',[]))
     3143            print('Defined exporters are:')
     3144            print('  ',list(exportersByExtension.get('powder',[]).keys()))
    30913145            raise G2ScriptException('No Writer for file type = "'+extension+'"')
    30923146        fil = os.path.abspath(os.path.splitext(fileroot)[0]+extension)
     
    32433297        :param float ttheta: peak position as 2Theta (deg)
    32443298
    3245         Note: only one of the parameters dspace, Q or ttheta may be specified
     3299        Note: only one of the parameters: dspace, Q or ttheta may be specified.
     3300        See :ref:`PeakRefine` for an example.
    32463301        '''
    32473302        import GSASIIlattice as G2lat
     
    33083363                                           Parms,Parms2,self.data['data'][1],bxye,[],
    33093364                                           False,controls,None)[0]
     3365
     3366    @property
     3367    def Peaks(self):
     3368        '''Provides a dict with the Peak List parameters
     3369        for this histogram.
     3370
     3371        :returns: dict with two elements where item
     3372          'peaks' is a list of peaks where each element is
     3373          [pos,pos-ref,area,area-ref,sig,sig-ref,gam,gam-ref],
     3374          where the -ref items are refinement flags and item
     3375          'sigDict' is a dict with possible items 'Back;#',
     3376          'pos#', 'int#', 'sig#', 'gam#'
     3377        '''
     3378        return self.data['Peak List']
     3379
     3380    @property
     3381    def PeakList(self):
     3382        '''Provides a list of peaks parameters
     3383        for this histogram.
     3384
     3385        :returns: a list of peaks, where each peak is a list containing
     3386          [pos,area,sig,gam]
     3387          (position, peak area, Gaussian width, Lorentzian width)
     3388           
     3389        '''
     3390        return [i[::2] for i in self.data['Peak List']['peaks']]
     3391   
     3392    def Export_peaks(self,filename):
     3393        '''Write the peaks file. The path is specified by filename
     3394        extension.
     3395       
     3396        :param str filename: name of the file, optionally with a path,
     3397            includes an extension
     3398        :returns: name of file that was written
     3399        '''
     3400        import GSASIIlattice as G2lat
     3401        import math
     3402        nptand = lambda x: np.tan(x*math.pi/180.)
     3403        fil = os.path.abspath(filename)
     3404        fp = open(filename,'w')
     3405        Inst,Inst2 = self.data['Instrument Parameters']
     3406        Type = Inst['Type'][0]
     3407        if 'T' not in Type:
     3408            import GSASIImath as G2mth
     3409            wave = G2mth.getWave(Inst)
     3410        else:
     3411            wave = None
     3412        pkdata = self.data['Peak List']
     3413        peaks = pkdata['peaks']
     3414        sigDict = pkdata['sigDict']
     3415        # code taken from GSASIIdataGUI OnExportPeakList
     3416        fp.write("#%s \n" % (self.name+' Peak List'))
     3417        if wave:
     3418            fp.write('#wavelength = %10.6f\n'%(wave))
     3419        if 'T' in Type:
     3420            fp.write('#%9s %10s %10s %12s %10s %10s %10s %10s %10s\n'%('pos','dsp','esd','int','alp','bet','sig','gam','FWHM'))                                   
     3421        else:
     3422            fp.write('#%9s %10s %10s %12s %10s %10s %10s\n'%('pos','dsp','esd','int','sig','gam','FWHM'))
     3423        for ip,peak in enumerate(peaks):
     3424            dsp = G2lat.Pos2dsp(Inst,peak[0])
     3425            if 'T' in Type:  #TOF - more cols
     3426                esds = {'pos':0.,'int':0.,'alp':0.,'bet':0.,'sig':0.,'gam':0.}
     3427                for name in list(esds.keys()):
     3428                    esds[name] = sigDict.get('%s%d'%(name,ip),0.)
     3429                sig = np.sqrt(peak[8])
     3430                gam = peak[10]
     3431                esddsp = G2lat.Pos2dsp(Inst,esds['pos'])
     3432                FWHM = G2pwd.getgamFW(gam,sig) +(peak[4]+peak[6])*np.log(2.)/(peak[4]*peak[6])     #to get delta-TOF from Gam(peak)
     3433                fp.write("%10.2f %10.5f %10.5f %12.2f %10.3f %10.3f %10.3f %10.3f %10.3f\n" % \
     3434                    (peak[0],dsp,esddsp,peak[2],peak[4],peak[6],peak[8],peak[10],FWHM))
     3435            else:               #CW
     3436                #get esds from sigDict for each peak & put in output - esds for sig & gam from UVWXY?
     3437                esds = {'pos':0.,'int':0.,'sig':0.,'gam':0.}
     3438                for name in list(esds.keys()):
     3439                    esds[name] = sigDict.get('%s%d'%(name,ip),0.)
     3440                sig = np.sqrt(peak[4]) #var -> sig
     3441                gam = peak[6]
     3442                esddsp = 0.5*esds['pos']*dsp/nptand(peak[0]/2.)
     3443                FWHM = G2pwd.getgamFW(gam,sig)      #to get delta-2-theta in deg. from Gam(peak)
     3444                fp.write("%10.4f %10.5f %10.5f %12.2f %10.5f %10.5f %10.5f \n" % \
     3445                    (peak[0],dsp,esddsp,peak[2],np.sqrt(max(0.0001,peak[4]))/100.,peak[6]/100.,FWHM/100.)) #convert to deg
     3446        fp.close()
     3447        return fil
    33103448
    33113449    def SaveProfile(self,filename):
Note: See TracChangeset for help on using the changeset viewer.