Changeset 5271


Ignore:
Timestamp:
May 5, 2022 9:29:34 PM (19 months ago)
Author:
toby
Message:

add in Laue Fringe peak fitting

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIdataGUI.py

    r5265 r5271  
    46284628        if self.G2plotNB:
    46294629            self.G2plotNB.Destroy()
    4630         if self.undofile:
     4630        if self.undofile and os.path.exists(self.undofile):
    46314631            os.remove(self.undofile)
    46324632        sys.exit()
    … …  
    79447944        #     import importlib as imp
    79457945        #     imp.reload(G2pdG)
     7946        #     imp.reload(G2pwd)
     7947        #     imp.reload(G2plt)
    79467948        #     print('reloading G2pdG')
    79477949        G2pdG.UpdatePeakGrid(G2frame,data)

  • trunk/GSASIImath.py

    r5260 r5271  
    47614761        gam = getTOFgamma(ins,dsp)
    47624762        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
    4763     elif 'C' in Parms['Type'][0]:
     4763    elif 'C' in Parms['Type'][0] or 'LF' in Parms['Type'][0]:
    47644764        for x in ['U','V','W','X','Y','Z']:
    47654765            ins[x] = Parms.get(x,[0.0,0.0])[ind]

  • trunk/GSASIIplot.py

    r5268 r5271  
    37223722                        G2frame.dataWindow.movePeak.Enable(len(selectedPeaks) == 1) # allow peak move from table when one peak is selected
    37233723                        for i,item in enumerate(data['peaks']):
     3724                            if type(item) is dict: continue
    37243725                            if i in selectedPeaks:
    37253726                                Ni = N+1

  • trunk/GSASIIpwd.py

    r5269 r5271  
    5252    print ('pydiffax is not available for this platform')
    5353import GSASIIfiles as G2fil
    54 #import LaueFringe as LF
    5554   
    5655# trig functions in degrees
    … …  
    864863    alpPink = lambda pos,alp0,alp1: alp0+alp1*tand(pos/2.)
    865864    betPink = lambda pos,bet0,bet1: bet0+bet1*tand(pos/2.)
    866     if 'T' in Inst['Type'][0]:
     865    if 'LF' in Inst['Type'][0]:
     866        return 3
     867    elif 'T' in Inst['Type'][0]:
    867868        dsp = pos/Inst['difC'][1]
    868869        alp = alpTOF(dsp,Inst['alpha'][1])
    … …  
    14091410    yb = getBackground('',parmDict,bakType,dataType,xdata,fixback)[0]
    14101411    yc = np.zeros_like(yb)
    1411     # if parmDict.get('LaueFringe',False):
    1412     #     if 'Lam1' in parmDict.keys():
    1413     #         lam = parmDict['Lam1']
    1414     #         lam2 = parmDict['Lam2']
    1415     #         Ka2 = True
    1416     #         lamRatio = 360*(lam2-lam)/(np.pi*lam)
    1417     #         kRatio = parmDict['I(L2)/I(L1)']
    1418     #     else:
    1419     #         lam = parmDict['Lam']
    1420     #         Ka2 = False
    1421     #     shol = 0
    1422     #     # loop over peaks
    1423     #     iPeak = -1
    1424     #     cells =  parmDict['ncell']
    1425     #     while True:
    1426     #         iPeak += 1
    1427     #         try:
    1428     #             pos = parmDict['pos'+str(iPeak)]
    1429     #             #tth = (pos-parmDict['Zero'])
    1430     #             intens = parmDict['int'+str(iPeak)]
    1431     #             damp =  parmDict['damp'+str(iPeak)]
    1432     #             asym =  parmDict['asym'+str(iPeak)]
    1433     #             sig =  parmDict['sig'+str(iPeak)]
    1434     #             gam =  parmDict['gam'+str(iPeak)]
    1435     #             fmin = 4 # for now make peaks 4 degrees wide
    1436     #             fmin = min(0.9*abs(xdata[-1] - xdata[0]),fmin) # unless the data range is smaller
    1437     #             iBeg = np.searchsorted(xdata,pos-fmin/2)
    1438     #             iFin = np.searchsorted(xdata,pos+fmin/2)
    1439     #             if not iBeg+iFin:       # skip peak below low limit
    1440     #                 continue
    1441     #             elif not iBeg-iFin:     # got peak above high limit (peaks sorted, so we can stop)
    1442     #                 break
    1443     #             LF.LaueFringePeakCalc(xdata,yc,lam,pos,intens,damp,asym,sig,gam,shol,cells,fmin)
    1444     #             if Ka2:
    1445     #                 pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
    1446     #                 iBeg = np.searchsorted(xdata,pos2-fmin)
    1447     #                 iFin = np.searchsorted(xdata,pos2+fmin)
    1448     #                 if iBeg-iFin:
    1449     #                     LF.LaueFringePeakCalc(xdata,yc,lam2,pos2,intens*kRatio,damp,asym,sig,gam,shol,cells,fmin)
    1450     #         except KeyError:        #no more peaks to process
    1451     #             return yb+yc
    1452     # elif 'C' in dataType:
    1453     if 'C' in dataType:
     1412    if 'LF' in dataType:
     1413        if 'Lam1' in parmDict.keys():
     1414            lam = parmDict['Lam1']
     1415            lam2 = parmDict['Lam2']
     1416            Ka2 = True
     1417            lamRatio = 360*(lam2-lam)/(np.pi*lam)
     1418            kRatio = parmDict['I(L2)/I(L1)']
     1419        else:
     1420            lam = parmDict['Lam']
     1421            Ka2 = False
     1422        shol = 0
     1423        # loop over peaks
     1424        iPeak = -1
     1425        try:
     1426            ncells =  parmDict['ncell']
     1427            clat =  parmDict['clat']
     1428        except KeyError: # no Laue info must be bkg fit
     1429            print('Laue Fit: no params, assuming bkg fit')
     1430            return yb
     1431        while True:
     1432            iPeak += 1
     1433            try:
     1434                #Qcen = 2 * np.pi * lam * (iPeak+1) / parmDict['clat']
     1435                l = parmDict['l'+str(iPeak)]
     1436                pos = 360 * np.arcsin(0.5 * lam * l / parmDict['clat']) / np.pi
     1437                parmDict['pos'+str(iPeak)] = pos
     1438                #tth = (pos-parmDict['Zero'])
     1439                intens = parmDict['int'+str(iPeak)]
     1440                damp =  parmDict['damp'+str(iPeak)]
     1441                asym =  parmDict['asym'+str(iPeak)]
     1442                sig =  parmDict['sig'+str(iPeak)]
     1443                gam =  parmDict['gam'+str(iPeak)]
     1444                fmin = 8 # for now make peaks 8 degrees wide
     1445                fmin = min(0.9*abs(xdata[-1] - xdata[0]),fmin) # unless the data range is smaller
     1446                iBeg = np.searchsorted(xdata,pos-fmin/2)
     1447                iFin = np.searchsorted(xdata,pos+fmin/2)
     1448                if not iBeg+iFin:       # skip peak below low limit
     1449                    continue
     1450                elif not iBeg-iFin:     # got peak above high limit (peaks sorted, so we can stop)
     1451                    break
     1452                #LF.plotme(fmin,lam,pos,intens,sig,gam,shol,ncells,clat,damp,asym)
     1453                LaueFringePeakCalc(xdata,yc,lam,pos,intens,sig,gam,shol,ncells,clat,damp,asym,fmin)
     1454                if Ka2:
     1455                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
     1456                    iBeg = np.searchsorted(xdata,pos2-fmin)
     1457                    iFin = np.searchsorted(xdata,pos2+fmin)
     1458                    if iBeg-iFin:
     1459                        LaueFringePeakCalc(xdata,yc,lam2,pos2,intens*kRatio,sig,gam,shol,ncells,clat,damp,asym,fmin)
     1460            except KeyError:        #no more peaks to process
     1461                return yb+yc
     1462    elif 'C' in dataType:
    14541463        shl = max(parmDict['SH/L'],0.002)
    14551464        Ka2 = False
    … …  
    16471656            ip = names.index(parm)
    16481657            dMdv[varyList.index(name)] = dMdpk[4*int(Id)+ip]
    1649     # if parmDict.get('LaueFringe',False):
    1650     #     for i,name in enumerate(varyList):
    1651     #         if not np.all(dMdv[i] == 0): continue
    1652     #         deltaParmDict = parmDict.copy()
    1653     #         delta = max(parmDict[name]/1e5,0.001)
    1654     #         deltaParmDict[name] += delta
    1655     #         intArrP = getPeakProfile(dataType,deltaParmDict,xdata,fixback,varyList,bakType)
    1656     #         deltaParmDict[name] -= 2*delta
    1657     #         intArrM = getPeakProfile(dataType,deltaParmDict,xdata,fixback,varyList,bakType)
    1658     #         dMdv[i] = 0.5 * (intArrP - intArrM) / delta
    1659     #     return dMdv
     1658    if 'LF' in dataType:
     1659        for i,name in enumerate(varyList):
     1660            if not np.all(dMdv[i] == 0): continue
     1661            deltaParmDict = parmDict.copy()
     1662            delta = max(parmDict[name]/1e5,0.001)
     1663            deltaParmDict[name] += delta
     1664            #print('num. deriv for',name,'val',deltaParmDict[name],'delta',delta)
     1665            intArrP = getPeakProfile(dataType,deltaParmDict,xdata,fixback,varyList,bakType)
     1666            deltaParmDict[name] -= 2*delta
     1667            intArrM = getPeakProfile(dataType,deltaParmDict,xdata,fixback,varyList,bakType)
     1668            dMdv[i] = 0.5 * (intArrP - intArrM) / delta
     1669        return dMdv
    16601670    if 'C' in dataType:
    16611671        shl = max(parmDict['SH/L'],0.002)
    … …  
    21122122    return True
    21132123
    2114 def getHeaderInfo(FitPgm,dataType):
     2124def getHeaderInfo(dataType):
    21152125    '''Provide parameter name, label name and formatting information for the
    21162126    contents of the peak table and where used in DoPeakFit
    … …  
    21182128    names = ['pos','int']
    21192129    lnames = ['position','intensity']
    2120     if FitPgm == 'LaueFringe':
    2121         names += ['damp','asym','sig','gam']
    2122         lnames += ['satellite\ndamping','satellite\nasym','sigma','gamma']
    2123         fmt = ["%10.5f","%10.2f","%10.3f","%10.3f","%10.3f","%10.3f"]
     2130    if 'LF' in dataType:
     2131        names = ['int','sig','gam','damp','asym','l']
     2132        lnames = ['intensity','sigma','gamma','satellite\ndamping','satellite\nasym','00l']
     2133        fmt = ["%10.2f","%10.3f","%10.3f","%10.3f","%10.3f","%4d"]
    21242134    elif 'C' in dataType:
    21252135        names += ['sig','gam']
    … …  
    21462156
    21472157    :param str FitPgm: type of fit to perform. At present this is ignored.
    2148     :param str FitPgm: type of fit to perform.
    2149       This should be 'LSQ' for a standard least-squares fit
    2150       with pseudo-Voigt peaks or 'LaueFringe' for on-specular thin-film fitting.
    21512158    :param list Peaks: a list of peaks. Each peak entry is a list with paired values:
    21522159      The number of pairs depends on the data type (see :func:`getHeaderInfo`).
    … …  
    23332340        peakDict = {}
    23342341        peakVary = []
    2335         names,_,_ = getHeaderInfo(FitPgm,dataType)
     2342        names,_,_ = getHeaderInfo(dataType)
     2343        off = 0
     2344        if 'LF' in dataType: off = 2
    23362345        for i,peak in enumerate(Peaks):
    23372346            if type(peak) is dict:
    … …  
    23402349            for j,name in enumerate(names):
    23412350                parName = name+str(i)
    2342                 peakDict[parName] = peak[2*j]
    2343                 if peak[2*j+1]:
     2351                peakDict[parName] = peak[off+2*j]
     2352                if peak[off+2*j+1]:
    23442353                    peakVary.append(parName)
    23452354        return peakDict,peakVary
    23462355               
    23472356    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
    2348         names,_,_ = getHeaderInfo(FitPgm,Inst['Type'][0])
     2357        off = 0
     2358        if 'LF' in Inst['Type'][0]:
     2359            off = 2
     2360            if 'clat' in varyList:
     2361                Peaks[-1]['clat'] = parmDict['clat']
     2362        names,_,_ = getHeaderInfo(Inst['Type'][0])
    23492363        for i,peak in enumerate(Peaks):
    23502364            if type(peak) is dict:
    23512365                continue
     2366            for j in range(len(names)):
     2367                parName = names[j]+str(i)
     2368                if parName in varyList or not peakInstPrmMode:
     2369                    peak[2*j+off] = parmDict[parName]               
    23522370            pos = parmDict['pos'+str(i)]
     2371            if 'LF' in Inst['Type'][0]: peak[0] = pos
    23532372            if 'difC' in Inst:
    23542373                dsp = pos/Inst['difC'][1]
    23552374            for j in range(len(names)):
    2356                 parName = names[j]+str(i)
    2357                 if parName in varyList or not peakInstPrmMode:
    2358                     peak[2*j] = parmDict[parName]
    2359                 elif 'alp' in parName:
     2375                if 'alp' in parName:
    23602376                    if 'T' in Inst['Type'][0]:
    23612377                        peak[2*j] = G2mth.getTOFalpha(parmDict,dsp)
    … …  
    23812397                       
    23822398    def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW):
     2399        if 'clat' in varyList:
     2400            print('c = {:.6f} esd {:.6f}'.format(parmDict['clat'],sigDict['clat']))
    23832401        print ('Peak coefficients:')
    2384         names,fmt,_ = getHeaderInfo(FitPgm,dataType)
     2402        names,fmt,_ = getHeaderInfo(dataType)
    23852403        head = 13*' '
    23862404        for name in names:
    … …  
    24472465    xBeg = np.searchsorted(x,Limits[0])
    24482466    xFin = np.searchsorted(x,Limits[1])+1
     2467    # find out what is varied
    24492468    bakType,bakDict,bakVary = SetBackgroundParms(Background)
    24502469    dataType,insDict,insVary = SetInstParms(Inst)
    … …  
    24602479    else:
    24612480        varyList = bakVary+insVary+peakVary
     2481        if 'LF' in Inst['Type'][0] and Peaks:
     2482            if Peaks[-1].get('clat-ref'): varyList += ['clat']
    24622483    fullvaryList = varyList[:]
    24632484    if not peakInstPrmMode:
    … …  
    24672488                raise Exception('Instrumental profile terms cannot be varied '+
    24682489                                    'after setPeakInstPrmMode(False) is used')
    2469     parmDict['LaueFringe'] = False
    2470     # if FitPgm == 'LaueFringe':
    2471     #     #======================================================================
    2472     #     print('Debug: reload LaueFringe')  # TODO: remove me
    2473     #     import imp
    2474     #     imp.reload(LF)
    2475     #     # TODO: remove ^^^^
    2476     #     #======================================================================
    2477     #     for v in ('U','V','W','X','Y','Z','alpha','alpha-0','alpha-1',
    2478     #         'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',):
    2479     #         if v in varyList:
    2480     #             raise Exception('Instrumental profile terms cannot be varied '+
    2481     #                                 'in Laue Fringe fits')
    2482     #     parmDict['LaueFringe'] = True
     2490    if 'LF' in Inst['Type'][0]:
     2491        warn = []
     2492        for v in ('U','V','W','X','Y','Z','alpha','alpha-0','alpha-1',
     2493            'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',):
     2494            if v in varyList:
     2495                warn.append(v)
     2496                del varyList[varyList.index(v)]
     2497        if warn:
     2498            print('Instrumental profile terms cannot be varied '+
     2499                                    'in Laue Fringe fits:',warn)
    24832500
    24842501    while not noFit:
    … …  
    24872504        Rvals = {}
    24882505        badVary = []
    2489         result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
     2506        try:
     2507            result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
    24902508               args=(x[xBeg:xFin],y[xBeg:xFin],fixback[xBeg:xFin],wtFactor*w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
     2509        except Exception as msg:
     2510            if GSASIIpath.GetConfigValue('debug'):
     2511                print('peak fit failure\n',msg)
     2512                import traceback
     2513                print (traceback.format_exc())
     2514            else:
     2515                print('peak fit failure')
     2516            return
    24912517        ncyc = int(result[2]['nfev']/2)
    24922518        runtime = time.time()-begin   
    … …  
    25222548    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
    25232549    if noFit:
     2550        GetPeaksParms(Inst,parmDict,Peaks,varyList)
    25242551        return
    25252552    sigDict = dict(zip(varyList,sig))
    … …  
    50305057    newRefs = np.array(newRefs)
    50315058    return True,newRefs
    5032    
     5059
     5060#===Laue Fringe code ===================================================================
     5061import NIST_profile as FP
     5062
     5063class profileObj(FP.FP_profile):
     5064    def conv_Lauefringe(self):
     5065        """Compute the FT of the Laue Fringe function"""
     5066           
     5067        me=self.get_function_name() #the name of this convolver,as a string
     5068        wave = self.param_dicts['conv_global']['dominant_wavelength']*1.e10 # in A       
     5069        pos = np.rad2deg(self.param_dicts["conv_global"]["twotheta0"])  # peak position as 2theta in deg
     5070        posQ = np.pi * 4 * np.sin(self.param_dicts["conv_global"]["twotheta0"]/2) / wave # peak position as Q
     5071        ttwid = self.twotheta_window_fullwidth_deg
     5072        ncell = self.param_dicts[me]['Ncells']
     5073        co2 = self.param_dicts[me]['clat'] / 2.
     5074        damp =  self.param_dicts[me]['damp']
     5075        asym =  self.param_dicts[me]['asym']
     5076        ttlist = np.linspace(pos-ttwid/2,pos+ttwid/2,len(self._epsb2))
     5077        Qs = np.pi * 4 * np.sin(np.deg2rad(ttlist/2)) / wave
     5078        w =  np.exp(-10**((damp-asym) * (Qs - posQ)**2))
     5079        w2 = np.exp(-10**((damp+asym) * (Qs - posQ)**2))
     5080        w[len(w)//2:] = w2[len(w)//2:]
     5081        weqdiv = w * np.sin(Qs * ncell * co2)**2 / (np.sin(Qs * co2)**2)
     5082        conv = FP.best_rfft(weqdiv)
     5083        conv[1::2] *= -1 #flip center
     5084        return conv
     5085       
     5086    def conv_Lorentzian(self):
     5087        """Compute the FT of a Lorentz function where gamma is the FWHM"""
     5088        ttwid = self.twotheta_window_fullwidth_deg
     5089        me=self.get_function_name() #the name of this convolver,as a string
     5090        g2gam = self.param_dicts[me]['g2gam'] # gsas-ii gamma in centidegrees
     5091        gamma = g2gam/100 # deg
     5092        ttlist = np.linspace(-ttwid/2,ttwid/2,len(self._epsb2))
     5093        eqdiv = (0.5 * gamma / np.pi) / (gamma**2/4. + ttlist**2)
     5094        conv = FP.best_rfft(eqdiv)
     5095        conv[1::2] *= -1 #flip center
     5096        return conv
     5097
     5098    def conv_Gaussian(self):
     5099        """Compute the FT of a Gaussian where sigma**2 is the variance"""
     5100        ttwid = self.twotheta_window_fullwidth_deg
     5101        me=self.get_function_name() #the name of this convolver,as a string
     5102        g2sig2 = self.param_dicts[me]['g2sig2'] # gsas-ii sigma**2 in centidegr**2
     5103        sigma = math.sqrt(g2sig2)/100.       
     5104        ttlist = np.linspace(-ttwid/2,ttwid/2,len(self._epsb2))
     5105        eqdiv = np.exp(-0.5*ttlist**2/sigma**2) / math.sqrt(2*np.pi*sigma**2)
     5106        conv = FP.best_rfft(eqdiv)
     5107        conv[1::2] *= -1 #flip center
     5108        return conv
     5109   
     5110def LaueFringePeakCalc(ttArr,intArr,lam,peakpos,intens,sigma2,gamma,shol,ncells,clat,damp,asym,calcwid,plot=False):
     5111    '''Compute the peakshape for a Laue Fringe peak convoluted with a Gaussian, Lorentzian &
     5112    an axial divergence asymmetry correction.
     5113
     5114    :param np.array ttArr: Array of two-theta values (in degrees)
     5115    :param np.array intArr: Array of intensity values (peaks are added to this)
     5116    :param float lam: wavelength in Angstrom
     5117    :param float peakpos: peak position in two-theta (deg.)
     5118    :param float intens: intensity factor for peak
     5119    :param float sigma2: Gaussian variance (in centidegrees**2) **
     5120    :param float gamma:  Lorenzian FWHM (in centidegrees) **
     5121    :param float shol: FCJ (S + H)/L where S=sample-half height, H=slit half-height, L=radius **
     5122    :param float ncells: number of unit cells in specular direction **
     5123    :param float clat: c lattice parameter **
     5124    :param float damp:
     5125    :param float asym:
     5126    :param float calcwid: two-theta (deg.) width for cutoff of peak computation.
     5127       Defaults to 5
     5128    :param bool plot: for debugging, shows contributions to peak
     5129
     5130    **  If term is <= zero, item is removed from convolution
     5131    '''
     5132    def LaueFringePeakPlot(ttArr,intArr):
     5133        import matplotlib.pyplot as plt
     5134        refColors = ['xkcd:blue','xkcd:red','xkcd:green','xkcd:cyan','xkcd:magenta','xkcd:black',
     5135                         'xkcd:pink','xkcd:brown','xkcd:teal','xkcd:orange','xkcd:grey','xkcd:violet',]
     5136        fig, ax = plt.subplots()
     5137        ax.set(title='Peak convolution functions @ 2theta={:.3f}'.format(peakpos),
     5138                   xlabel=r'$\Delta 2\theta, deg$',
     5139                   ylabel=r'Intensity (arbitrary)')
     5140        ax.set_yscale("log",nonpositive='mask')
     5141        ttmin = ttmax = 0
     5142        for i,conv in enumerate(convList):
     5143            f = NISTpk.convolver_funcs[conv]()
     5144            if f is None: continue
     5145            FFT = FP.best_irfft(f)
     5146            if f[1].real > 0: FFT = np.roll(FFT,int(len(FFT)/2.))
     5147            FFT /= FFT.max()
     5148            if i == 0:
     5149                tt = np.linspace(-NISTpk.twotheta_window_fullwidth_deg/2,
     5150                                    NISTpk.twotheta_window_fullwidth_deg/2,len(FFT))
     5151            ttmin = min(ttmin,tt[np.argmax(FFT>.005)])
     5152            ttmax = max(ttmax,tt[::-1][np.argmax(FFT[::-1]>.005)])
     5153            color = refColors[i%len(refColors)]
     5154            ax.plot(tt,FFT,color,label=conv[5:])
     5155        color = refColors[(i+1)%len(refColors)]
     5156        ax.plot(ttArr-peakpos,intArr/max(intArr),color,label='Convolution')
     5157        ax.set_xlim((ttmin,ttmax))
     5158        ax.legend(loc='best')
     5159        plt.show()
     5160    # hardcoded constants
     5161    diffRadius = 220 # diffractometer radius in mm; needed for axial divergence, etc, but should not matter
     5162    axial_factor = 1.5  # fudge factor to bring sh/l broadening to ~ agree with FPA
     5163    equatorial_divergence_deg = 0.5 # not sure exactly what this impacts
     5164    NISTparms = {
     5165        "": {
     5166            'equatorial_divergence_deg' : equatorial_divergence_deg,
     5167            'dominant_wavelength' : 1.e-10 * lam,
     5168            'diffractometer_radius' : 1e-3* diffRadius, # diffractometer radius in m
     5169            'oversampling' : 8,
     5170            },
     5171        "emission": {
     5172            'emiss_wavelengths' : 1.e-10 * np.array([lam]),
     5173            'emiss_intensities' : np.array([1.]),
     5174            'emiss_gauss_widths' : 1.e-10 * 1.e-3 * np.array([0.001]),
     5175            'emiss_lor_widths' : 1.e-10 * 1.e-3 * np.array([0.001]),
     5176            'crystallite_size_gauss' : 1.e-9 * 1e6,
     5177            'crystallite_size_lor' : 1.e-9 * 1e6
     5178            },
     5179        "axial": {
     5180            'axDiv':"full",
     5181            'slit_length_source' : 1e-3 * diffRadius * shol * axial_factor,
     5182            'slit_length_target' : 1e-3 * diffRadius * shol * 1.00001 * axial_factor, # != 'slit_length_source'
     5183            'length_sample' : 1e-3 * diffRadius * shol * axial_factor,
     5184            'n_integral_points' : 10,
     5185            'angI_deg' : 2.5,
     5186            'angD_deg': 2.5,
     5187            },
     5188        'Gaussian': {'g2sig2': sigma2},
     5189        'Lorentzian': {'g2gam': gamma},
     5190        'Lauefringe': {'Ncells': ncells, 'clat':clat, 'damp': damp, 'asym': asym},
     5191        }
     5192    NISTpk=profileObj(anglemode="twotheta",
     5193                    output_gaussian_smoother_bins_sigma=1.0,
     5194                    oversampling=NISTparms.get('oversampling',10))
     5195    NISTpk.debug_cache=False   
     5196    for key in NISTparms: #set parameters for each convolver
     5197        if key:
     5198            NISTpk.set_parameters(convolver=key,**NISTparms[key])
     5199        else:
     5200            NISTpk.set_parameters(**NISTparms[key])
     5201    # find closest point to peak location (which may be outside limits of the array)
     5202    center_bin_idx=min(ttArr.searchsorted(peakpos),len(ttArr)-1)
     5203    step = (ttArr[-1]-ttArr[0])/(len(ttArr)-1)
     5204    NISTpk.set_optimized_window(twotheta_exact_bin_spacing_deg=step,
     5205                    twotheta_window_center_deg=ttArr[center_bin_idx],
     5206                    twotheta_approx_window_fullwidth_deg=calcwid,
     5207                    )
     5208    NISTpk.set_parameters(twotheta0_deg=peakpos)
     5209    convList = ['conv_emission']
     5210    if ncells: convList += ['conv_Lauefringe']
     5211    if sigma2 > 0: convList += ['conv_Gaussian']
     5212    if gamma > 0: convList += ['conv_Lorentzian']
     5213    if shol > 0: convList += ['conv_axial']
     5214
     5215    # global deriv
     5216    # if deriv:
     5217    #     peakObj = NISTpk.compute_line_profile(convolver_names=convList,compute_derivative=True)
     5218    # else:
     5219    #     peakObj = NISTpk.compute_line_profile(convolver_names=convList)
     5220    peakObj = NISTpk.compute_line_profile(convolver_names=convList)
     5221   
     5222    pkPts = len(peakObj.peak)
     5223    pkMax = peakObj.peak.max()
     5224    startInd = center_bin_idx-(pkPts//2)
     5225    istart = None
     5226    pstart = None
     5227    iend = None
     5228    pend = None
     5229    # adjust data range if peak calc begins below data range or ends above data range
     5230    # but range of peak calc should not extend past both ends of ttArr
     5231    if startInd < 0:
     5232        iend = startInd+pkPts
     5233        pstart = -startInd
     5234    elif startInd > len(intArr):
     5235        return
     5236    elif startInd+pkPts >= len(intArr):
     5237        offset = pkPts - len( intArr[startInd:] )
     5238        istart = startInd
     5239        iend = startInd+pkPts-offset
     5240        pend = -offset
     5241    else:
     5242        istart = startInd
     5243        iend = startInd+pkPts
     5244    intArr[istart:iend] += intens * peakObj.peak[pstart:pend]/pkMax
     5245    if plot:
     5246        LaueFringePeakPlot(ttArr[istart:iend], (intens * peakObj.peak[pstart:pend]/pkMax))
     5247       
     5248def LaueSatellite(peakpos,wave,c,ncell,j=[-4,-3,-2,-1,0,1,2,3,4]):
     5249    '''Returns the locations of the Laue satellite positions relative
     5250    to the peak position
     5251   
     5252    :param float peakpos: the peak position in degrees 2theta
     5253    :param float ncell:   Laue fringe parameter, number of unit cells in layer
     5254    :param list j: the satellite order, where j=-1 is the first satellite
     5255      on the lower 2theta side and j=1 is the first satellite on the high
     5256      2theta side. j=0 gives the peak position
     5257    '''
     5258    Qpos = 4 * np.pi * np.sin(peakpos * np.pi / 360) / wave
     5259    dQvals = (2 * np.array(j) + np.sign(j)) * np.pi / (c * ncell)
     5260    return np.arcsin((Qpos+dQvals)*wave/(4*np.pi)) * (360 / np.pi)
     5261
    50335262#### testing data
    50345263NeedTestData = True

  • trunk/GSASIIpwdGUI.py

    r5265 r5271  
    834834        if not G2frame.GSASprojectfile:            #force a save of the gpx file so SaveState can write in the same directory
    835835            G2frame.OnFileSaveas(event)
    836         FitPgm = 'LSQ'
    837 #        if data.get('FitPgm',0) == 1:
    838 #            FitPgm = 'LaueFringe'
    839         wx.CallAfter(OnPeakFit,FitPgm)
     836        wx.CallAfter(OnPeakFit)
    840837       
    841838    def OnOneCycle(event):
    … …  
    844841            reflGrid.HideCellEditControl()
    845842            reflGrid.DisableCellEditControl()
    846         FitPgm = 'LSQ'
    847 #        if data.get('FitPgm',0) == 1:
    848 #            FitPgm = 'LaueFringe'
    849         wx.CallAfter(OnPeakFit,FitPgm,oneCycle=True)
     843        wx.CallAfter(OnPeakFit,oneCycle=True)
    850844       
    851845    def OnSeqPeakFit(event):
    … …  
    876870        print ('Peak Fitting with '+controls['deriv type']+' derivatives:')
    877871        oneCycle = False
    878         FitPgm = 'LSQ'
    879 #        if data.get('FitPgm',0) == 1:
    880 #            FitPgm = 'LaueFringe'
    881872        prevVaryList = []
    882873        peaks = None
    … …  
    903894                data = Pattern[1]
    904895                fixback = GetFileBackground(G2frame,data,background,scale=False)
    905                 peaks['sigDict'],result,sig,Rvals,varyList,parmDict,fullvaryList,badVary = G2pwd.DoPeakFit(FitPgm,peaks['peaks'],
     896                peaks['sigDict'],result,sig,Rvals,varyList,parmDict,fullvaryList,badVary = G2pwd.DoPeakFit(None,peaks['peaks'],
    906897                    background,limits,inst,inst2,data,fixback,prevVaryList,oneCycle,controls)   #needs wtFactor after controls?
    907898                if len(result[0]) != len(fullvaryList):
    … …  
    933924        G2plt.PlotPatterns(G2frame,plotType='PWDR')
    934925       
    935     def OnPeakFit(FitPgm,oneCycle=False,noFit=False):
     926    def OnPeakFit(oneCycle=False,noFit=False):
    936927        'Do peak fitting by least squares'
    937928        SaveState()
    … …  
    939930        if not controls:
    940931            controls = {'deriv type':'analytic','min dM/M':0.001,}     #fill in defaults if needed
    941         print ('Peak Fitting with '+controls['deriv type']+' derivatives:')
     932        #print ('Peak Fitting with '+controls['deriv type']+' derivatives:')
    942933        PatternId = G2frame.PatternId
    943934        peaks = G2frame.GPXtree.GetItemPyData(G2gd.GetGPXtreeItemId(G2frame,PatternId, 'Peak List'))
    … …  
    952943        wtFactor = Pattern[0]['wtFactor']
    953944        bxye = GetFileBackground(G2frame,data,background,scale=False)
    954         overallInfo = []
    955 #         peaks['LaueFringe'] = peaks.get('LaueFringe',{})
    956 #         peaks['LaueFringe']['satellites'] = []
    957 #         import LaueFringe as LF
    958 #         lines = peaks['LaueFringe'].get('Show')
    959 #         if 'Laue' in FitPgm:
    960 #             overallInfo = [{'ncell':peaks['LaueFringe']['ncell']}] # add overall info
    961 #             if lines:
    962 #                 for i in peaks['peaks']:
    963 #                     pks = list(range(-lines,0)) + list(range(1,lines+1))
    964 #                     peaks['LaueFringe']['satellites'].extend(
    965 #                         LF.LaueSatellite(i[0],peaks['LaueFringe']['ncell'],pks))
    966 # #======================================================================
    967 #         print('Debug: reload G2pwd')  # TODO: remove me
    968 #         import imp
    969 #         imp.reload(G2pwd)
    970 #         imp.reload(G2plt)
    971 #         # TODO: remove ^^^^
    972 #======================================================================
     945        overallInfo = {}
     946        peaks['LaueFringe'] = peaks.get('LaueFringe',{})
     947        peaks['LaueFringe']['satellites'] = []
     948        lines = peaks['LaueFringe'].get('Show')
     949        if 'LF' in inst['Type'][0]:
     950            wave = G2mth.getWave(inst)
     951            overallInfo = {'ncell':peaks['LaueFringe']['ncell'],
     952                            'clat': peaks['LaueFringe']['clat'],
     953                            'clat-ref': peaks['LaueFringe']['clat-ref']
     954                               } # add overall info
     955            if lines:
     956                for i in peaks['peaks']:
     957                    pks = list(range(-lines,0)) + list(range(1,lines+1))
     958                    peaks['LaueFringe']['satellites'].extend(
     959                        G2pwd.LaueSatellite(i[0],wave,peaks['LaueFringe']['clat'],peaks['LaueFringe']['ncell'],pks))
     960        peaksplus = peaks['peaks'] + [overallInfo]
    973961        if noFit:
    974             results = G2pwd.DoPeakFit(FitPgm,peaks['peaks']+overallInfo,background,limits,inst,inst2,data,bxye,[],oneCycle,controls,wtFactor,noFit=True)
     962            results = G2pwd.DoPeakFit(None,peaksplus,background,limits,inst,inst2,data,bxye,[],oneCycle,controls,wtFactor,noFit=True)
    975963            G2plt.PlotPatterns(G2frame,plotType='PWDR')
    976964            return
    977         dlg = wx.ProgressDialog('Residual','Peak fit Rwp = ',101.0,
    978             parent=G2frame,
    979             style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
    980         results = G2pwd.DoPeakFit(FitPgm,peaks['peaks']+overallInfo,background,limits,inst,inst2,data,bxye,[],oneCycle,controls,wtFactor,dlg)
     965        try:
     966            dlg = wx.ProgressDialog('Residual','Peak fit Rwp = ',101.0,parent=G2frame,
     967                style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
     968            results = G2pwd.DoPeakFit(None,peaksplus,background,limits,inst,inst2,data,bxye,[],oneCycle,controls,wtFactor,dlg)
     969        finally:
     970            dlg.Destroy()
     971        if results is None:
     972            return
    981973        peaks['sigDict'] = results[0]
     974        if 'LF' in inst['Type'][0] and 'clat' in overallInfo:
     975            peaks['LaueFringe']['clat'] = overallInfo['clat']
    982976        text = 'Peak fit: Rwp=%.2f%% Nobs= %d Nparm= %d Npeaks= %d'%(results[3]['Rwp'],results[1][2]['fjac'].shape[1],len(results[0]),len(peaks['peaks']))
    983977        newpeaks = copy.copy(peaks)
    … …  
    11341128            wx.CallAfter(UpdatePeakGrid,G2frame,data)
    11351129
    1136     def updateMe(*args):
    1137         'Redraw the peak listings after the mode changes'
    1138         wx.CallAfter(UpdatePeakGrid,G2frame,data)
    1139         wx.CallLater(100,RefreshPeakGrid,None)
    1140 
    11411130    def RefreshPeakGrid(event):
    11421131        'recompute & plot the peaks any time a value in the table is edited'
    1143         FitPgm = 'LSQ'
    1144 #        if data.get('FitPgm',0) == 1:
    1145 #            FitPgm = 'LaueFringe'
    1146         OnPeakFit(FitPgm,noFit=True)
     1132        if 'LF' in Inst['Type'][0]:
     1133            for i in range(len(data['LFpeaks'])):
     1134                data['peaks'][i][2:] = data['LFpeaks'][i]
     1135        OnPeakFit(noFit=True)
    11471136       
    11481137    def OnSetPeakWidMode(event):
    … …  
    11551144    #### beginning of UpdatePeakGrid
    11561145    #======================================================================
    1157     data['FitPgm'] = data.get('FitPgm',0)
    11581146    G2frame.GetStatusBar().SetStatusText('Global refine: select refine column & press Y or N',1)
    11591147    G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.PeakMenu)
    … …  
    11711159    G2frame.Bind(wx.EVT_MENU, OnResetSigGam, id=G2G.wxID_RESETSIGGAM)
    11721160    G2frame.Bind(wx.EVT_MENU, OnSetPeakWidMode, id=G2G.wxID_SETUNVARIEDWIDTHS)
     1161    OnSetPeakWidMode(None)
    11731162    if data['peaks']:
    11741163        G2frame.dataWindow.AutoSearch.Enable(False)
    … …  
    11901179    colLabels = []
    11911180    Types = []
    1192     for _,f,l in zip(*G2pwd.getHeaderInfo(
    1193         ('LSQ','LaueFringe')[data.get('FitPgm',0)],Inst['Type'][0])):
     1181    for _,f,l in zip(*G2pwd.getHeaderInfo(Inst['Type'][0])):
    11941182        colLabels.append(l)
    1195         colLabels.append('refine')
    1196         Types.append(wg.GRID_VALUE_FLOAT +
    1197                          f.replace('%',':').replace('f','').replace('.',','))
    1198         Types.append(wg.GRID_VALUE_BOOL)
     1183        if "d" in f:
     1184            Types.append(wg.GRID_VALUE_NUMBER)
     1185        else:
     1186            Types.append(wg.GRID_VALUE_FLOAT + f.replace('%',':').replace('f','').replace('.',','))
     1187            colLabels.append('refine')
     1188            Types.append(wg.GRID_VALUE_BOOL)
    11991189    T = []
    12001190    for peak in data['peaks']:
    … …  
    12071197    for key in T: X.append(D[key])
    12081198    data['peaks'] = X
    1209     # pad peak table, if needed
    1210     for j in range(len(data['peaks'])):
    1211         for i in range(len(data['peaks'][j]),len(colLabels)):
    1212             if colLabels[i] == 'refine':
    1213                 data['peaks'][j].append(False)
    1214             else:
    1215                 data['peaks'][j].append(1.0)
     1199
    12161200    G2frame.dataWindow.ClearData()
    12171201    mainSizer = wx.BoxSizer(wx.VERTICAL)
    12181202    G2frame.GPXtree.SetItemPyData(G2frame.PickId,data)
    1219     G2frame.PeakTable = G2G.Table(data['peaks'],rowLabels=rowLabels,colLabels=colLabels,types=Types)
     1203    if 'LF' in Inst['Type'][0]:
     1204        data['LFpeaks'] = []
     1205        for i in range(len(data['peaks'])):
     1206            if len(data['peaks'][i]) == 8:  # need to extend the entry
     1207                if 'Lam' in Inst:
     1208                    lam = Inst['Lam'][0]
     1209                else:
     1210                    lam = (Inst['Lam1'][0] +
     1211                            Inst['I(L2)/I(L1)'][0] * Inst['Lam2'][0]) / (
     1212                                1 + Inst['I(L2)/I(L1)'][0])
     1213                if i == 0:
     1214                    pos = data['peaks'][i][0]
     1215                    data['LaueFringe']['clat'] = 0.5 * lam / np.sin(pos*np.pi/360)
     1216                    l = 1
     1217                else:
     1218                    l = data['LaueFringe']['clat'] / (0.5 * lam / np.sin(data['peaks'][i][0]*np.pi/360))
     1219                    l = int(l + 0.5)
     1220                    data['peaks'][i][0] = 360 * np.arcsin(0.5 * lam / (data['LaueFringe']['clat']/l)) / np.pi
     1221                data['peaks'][i] += [0., False, 0., False, l, None]
     1222            data['LFpeaks'].append(data['peaks'][i][2:])
     1223        G2frame.PeakTable = G2G.Table(data['LFpeaks'],rowLabels=rowLabels,colLabels=colLabels,types=Types)
     1224    else:
     1225        G2frame.PeakTable = G2G.Table(data['peaks'],rowLabels=rowLabels,colLabels=colLabels,types=Types)
    12201226    #G2frame.SetLabel(G2frame.GetLabel().split('||')[0]+' || '+'Peak List')
    12211227    G2frame.dataWindow.currentGrids = []
    … …  
    12341240    mainSizer.Add(topSizer,0,wx.EXPAND)
    12351241    G2G.HorizontalLine(mainSizer,G2frame.dataWindow)
    1236     # if 'C' in Inst['Type'][0]:
    1237         # topSizer = wx.BoxSizer(wx.HORIZONTAL)
    1238         # topSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Fitting mode: '),0,WACV)
    1239         # pkType = G2G.G2ChoiceButton(G2frame.dataWindow,('Powder diffraction','Laue fringes'),
    1240         #             indLoc=data,indKey='FitPgm',
    1241         #             onChoice=updateMe)
    1242         # topSizer.Add(pkType,0,WACV)
    1243         # if data['FitPgm']:
    1244         #     topSizer.Add((15,-1))
    1245         #     if 'LaueFringe' not in data:
    1246         #         data['LaueFringe'] = {'ncell':20}
    1247         #     elif 'ncell' not in data['LaueFringe']:
    1248         #         data['LaueFringe']['ncell'] = 20
    1249         #     siz = G2G.G2SpinWidget(G2frame.dataWindow,data['LaueFringe'] ,'ncell',
    1250         #                                'Laue ncell',
    1251         #                                onChange=RefreshPeakGrid,onChangeArgs=[None])
    1252         #     topSizer.Add(siz,0,WACV)
    1253         #     data['LaueFringe']['Show'] =  data['LaueFringe'].get('Show',0)
    1254         #     #ch = G2G.G2CheckBox(G2frame.dataWindow,' Show',data['LaueFringe'],'Show',
    1255         #     #                        OnChange=RefreshPeakGrid)
    1256         #     topSizer.Add(wx.StaticText(G2frame.dataWindow,label='  Show '),0,WACV)
    1257         #     ch = G2G.EnumSelector(G2frame.dataWindow,data['LaueFringe'],'Show',
    1258         #                             ['None','1','2','3','4','5','6'],list(range(7)),
    1259         #                             OnChange=RefreshPeakGrid)
    1260         #     topSizer.Add(ch,0,WACV)
    1261         #     topSizer.Add(wx.StaticText(G2frame.dataWindow,label=' satellites'),0,WACV)
    1262         # mainSizer.Add(topSizer)
     1242    if 'LF' in Inst['Type'][0]:
     1243        mainSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Laue Fringe fitting'))
     1244        topSizer = wx.BoxSizer(wx.HORIZONTAL)
     1245        topSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Overall parms: c='),0,WACV)
     1246        data['LaueFringe'] = data.get('LaueFringe',{})
     1247        data['LaueFringe']['ncell'] = data['LaueFringe'].get('ncell',20)
     1248        data['LaueFringe']['clat'] =  data['LaueFringe'].get('clat',9.0)
     1249        data['LaueFringe']['clat-ref'] =  data['LaueFringe'].get('clat-ref',False)
     1250        data['LaueFringe']['Show'] =  data['LaueFringe'].get('Show',0)
     1251        cVal = G2G.ValidatedTxtCtrl(G2frame.dataWindow,data['LaueFringe'],'clat',
     1252                                    typeHint=float,nDig=(10,4),
     1253                                    OnLeave=lambda *arg,**kw:RefreshPeakGrid(None))
     1254        topSizer.Add(cVal,0,WACV)
     1255        cRef = G2G.G2CheckBox(G2frame.dataWindow,'ref',data['LaueFringe'],'clat-ref')
     1256        topSizer.Add(cRef,0,WACV)
     1257        topSizer.Add((15,-1))
     1258        siz = G2G.G2SpinWidget(G2frame.dataWindow,data['LaueFringe'] ,'ncell',
     1259                                       'Laue ncell',
     1260                                       onChange=RefreshPeakGrid,onChangeArgs=[None])
     1261        topSizer.Add(siz,0,WACV)
     1262        topSizer.Add(wx.StaticText(G2frame.dataWindow,label='  Show '),0,WACV)
     1263        ch = G2G.EnumSelector(G2frame.dataWindow,data['LaueFringe'],'Show',
     1264                                    ['None','1','2','3','4','5','6'],list(range(7)),
     1265                                    OnChange=RefreshPeakGrid)
     1266        topSizer.Add(ch,0,WACV)
     1267        topSizer.Add(wx.StaticText(G2frame.dataWindow,label=' satellites'),0,WACV)
     1268        mainSizer.Add(topSizer)
    12631269    mainSizer.Add(reflGrid,0,wx.ALL,10)
    12641270    G2frame.dataWindow.SetSizer(mainSizer)

  • trunk/imports/G2pwd_csv.py

    r5268 r5271  
    8080                    try:
    8181                        positions[j] = int(v.strip().split('=')[1])
     82                        if j == 1: positions[2] = -1
    8283                    except:
    8384                        print('Error parsing: "'+S+'"')
Note: See TracChangeset for help on using the changeset viewer.