Changeset 892


Ignore:
Timestamp:
Apr 26, 2013 3:45:31 PM (9 years ago)
Author:
vondreele
Message:

develop new correction for detector penetration effects.
(1-cos2Q) dependent correction with a refinable coeff. Pretty good 1st correction. Not yet right if detector tilted.

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIimage.py

    r850 r892  
    128128    return err,makeParmsEllipse(bb)
    129129   
    130 def FitDetector(rings,p0,wave):
     130def FitDetector(rings,varyList,parmDict):
    131131       
    132132    def CalibPrint(ValSig):
     
    149149        print sigstr       
    150150       
    151     def ellipseCalcD(B,xyd,wave):
     151    def ellipseCalcD(B,xyd,varyList,parmDict):
    152152        x = xyd[0]
    153153        y = xyd[1]
    154154        dsp = xyd[2]
    155         dist,x0,y0,phi,tilt = B
     155        wave = parmDict['wave']
     156        if 'dep' in varyList:
     157            dist,x0,y0,tilt,phi,dep = B[:6]
     158        else:
     159            dist,x0,y0,tilt,phi = B[:5]
     160            dep = parmDict['dep']
     161        if 'wave' in varyList:
     162            wave = B[-1]
    156163        tth = 2.0*npasind(wave/(2.*dsp))
     164        dxy = dep*(1.-npcosd(tth))
    157165        ttth = nptand(tth)
    158         radius = dist*ttth
     166        radius = (dist+dxy)*ttth
    159167        stth = npsind(tth)
    160168        cosb = npcosd(tilt)
    161         R1 = dist*stth*npcosd(tth)*cosb/(cosb**2-stth**2)
     169        R1 = (dist+dxy)*stth*npcosd(tth)*cosb/(cosb**2-stth**2)
    162170        R0 = np.sqrt(R1*radius*cosb)
    163171        zdis = R1*ttth*nptand(tilt)
     
    168176        return (XR/R0)**2+(YR/R1)**2-1
    169177       
    170     def ellipseCalcW(C,xyd):
    171         dist,x0,y0,phi,tilt,wave = C
    172         B = dist,x0,y0,phi,tilt
    173         return ellipseCalcD(B,xyd,wave)
    174                
    175     names = ['distance','det-X0','det-Y0','rotate','tilt','wavelength']
    176     fmt = ['%12.2f','%12.2f','%12.2f','%12.2f','%12.2f','%12.5f']   
    177     result = leastsq(ellipseCalcD,p0,args=(rings.T,wave),full_output=True)
    178     result[0][3] = np.mod(result[0][3],360.0)               #remove random full circles
     178    names = ['dist','det-X','det-Y','tilt','phi','dep','wave']
     179    fmt = ['%12.2f','%12.2f','%12.2f','%12.2f','%12.2f','%12.3f','%12.5f']
     180    p0 = [parmDict[key] for key in varyList]
     181    result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True)
     182    parmDict.update(zip(varyList,result[0]))
    179183    vals = list(result[0])
    180     vals.append(wave)
    181     chi = np.sqrt(np.sum(ellipseCalcD(result[0],rings.T,wave)**2))
     184    chi = np.sqrt(np.sum(ellipseCalcD(result[0],rings.T,varyList,parmDict)**2))
    182185    sig = list(chi*np.sqrt(np.diag(result[1])))
    183     sig.append(0.)
     186    sigList = np.zeros(7)
     187    for i,name in enumerate(names):
     188        if name in varyList:
     189            sigList[i] = sig[varyList.index(name)]
    184190    ValSig = zip(names,fmt,vals,sig)
    185191    CalibPrint(ValSig)
    186     try:
    187         print 'Trial refinement of wavelength - not used for calibration'
    188         p0 = result[0]
    189         p0 = np.append(p0,wave)
    190         resultW = leastsq(ellipseCalcW,p0,args=(rings.T,),full_output=True)
    191         resultW[0][3] = np.mod(result[0][3],360.0)          #remove random full circles
    192         sig = np.sqrt(np.sum(ellipseCalcW(resultW[0],rings.T)**2))
    193         ValSig = zip(names,fmt,resultW[0],sig*np.sqrt(np.diag(resultW[1])))
    194         CalibPrint(ValSig)
    195         return result[0],resultW[0][-1]       
    196     except ValueError:
    197         print 'Bad refinement - no result'
    198         return result[0],wave
     192#    try:
     193#        print 'Trial refinement of wavelength - not used for calibration'
     194#        p0 = result[0]
     195#        print p0
     196#        print parms
     197#        p0 = np.append(p0,parms[0])
     198#        resultW = leastsq(ellipseCalcW,p0,args=(rings.T,parms[1:]),full_output=True)
     199#        resultW[0][3] = np.mod(result[0][3],360.0)          #remove random full circles
     200#        sig = np.sqrt(np.sum(ellipseCalcW(resultW[0],rings.T)**2))
     201#        ValSig = zip(names,fmt,resultW[0],sig*np.sqrt(np.diag(resultW[1])))
     202#        CalibPrint(ValSig)
     203#        return result[0],resultW[0][-1]       
     204#    except ValueError:
     205#        print 'Bad refinement - no result'
     206#        return result[0],wave
    199207                   
    200208def ImageLocalMax(image,w,Xpix,Ypix):
     
    279287    tilt = data['tilt']
    280288    phi = data['rotation']
     289    dep = data['DetDepth']
    281290    radii = [0,0]
    282291    tth = 2.0*asind(data['wavelength']/(2.*dsp))
     
    285294    ctth = cosd(tth)
    286295    cosb = cosd(tilt)
    287     radius = dist*ttth
    288     radii[1] = dist*stth*ctth*cosb/(cosb**2-stth**2)
     296    dxy = dep*(1.-npcosd(tth))
     297    radius = ttth*(dist+dxy)
     298    radii[1] = (dist+dxy)*stth*ctth*cosb/(cosb**2-stth**2)
    289299    if radii[1] > 0:
    290300        radii[0] = math.sqrt(radii[1]*radius*cosb)
     
    313323    wave = data['wavelength']
    314324    dist = data['distance']
     325    dep = data['DetDepth']
    315326    tth = 2.0*asind(wave/(2.*dsp))
     327    dxy = dep*(1.-npcosd(tth))
    316328    ttth = tand(tth)
    317     radius = dist*ttth
     329    radius = (dist+dxy)*ttth
    318330    stth = sind(tth)
    319331    cosb = cosd(tilt)
    320     R1 = dist*stth*cosd(tth)*cosb/(cosb**2-stth**2)
     332    R1 = (dist+dxy)*stth*cosd(tth)*cosb/(cosb**2-stth**2)
    321333    R0 = math.sqrt(R1*radius*cosb)
    322334    zdis = R1*ttth*tand(tilt)
     
    338350    tilt = data['tilt']
    339351    phi = data['rotation']
     352    dep = data['DetDepth']
    340353    LRazim = data['LRazimuth']
    341354    azmthoff = data['azmthOff']
     
    346359    Z = np.dot(X,makeMat(tilt,0)).T[2]
    347360    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
     361    dxy = dep*(1.-npcosd(tth))
     362    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z+dxy))
    348363    dsp = wave/(2.*npsind(tth/2.))
    349364    azm = (npatan2d(dx,-dy)+azmthoff+720.)%360.
     
    412427        HKL += hkl
    413428    HKL = G2lat.sortHKLd(HKL,True,False)
    414     wave = data['wavelength']
    415     cent = data['center']   
    416     dist = data['distance']
    417     tilt = data['tilt']
    418     phi = data['rotation']
     429    varyList = ['dist','det-X','det-Y','tilt','phi']
     430    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
     431        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
    419432    for H in HKL:
    420433        dsp = H[3]
     
    428441            continue
    429442    rings = np.concatenate((data['rings']),axis=0)
    430     p0 = [dist,cent[0],cent[1],phi,tilt]
    431     result,newWave = FitDetector(rings,p0,wave)
    432     data['distance'] = result[0]
    433     data['center'] = result[1:3]
    434     data['rotation'] = np.mod(result[3],360.0)
    435     data['tilt'] = result[4]
     443    if data['DetDepthRef']:
     444        varyList.append('dep')
     445    FitDetector(rings,varyList,parmDict)
     446    data['distance'] = parmDict['dist']
     447    data['center'] = [parmDict['det-X'],parmDict['det-Y']]
     448    data['rotation'] = np.mod(parmDict['phi'],360.0)
     449    data['tilt'] = parmDict['tilt']
     450    data['DetDepth'] = parmDict['dep']
    436451    N = len(data['ellipses'])
    437452    data['ellipses'] = []           #clear away individual ellipse fits
     
    443458    G2plt.PlotImage(self,newImage=True)       
    444459    return True
    445    
    446460           
    447461def ImageCalibrate(self,data):
     
    589603        print 'Only one ring fitted. Check your wavelength.'
    590604        return False
    591     cent = data['center'] = [xSum/Zsum,ySum/Zsum]
    592     wave = data['wavelength']
    593     dist = data['distance'] = distSum/Zsum
     605    data['center'] = [xSum/Zsum,ySum/Zsum]
     606    data['distance'] = distSum/Zsum
    594607   
    595608    #possible error if no. of rings < 3! Might need trap here
     
    608621            Zsign = -1
    609622   
    610     tilt = data['tilt'] = Zsign*tiltSum/Zsum
    611     phi = data['rotation'] = phiSum/Zsum
     623    data['tilt'] = Zsign*tiltSum/Zsum
     624    data['rotation'] = phiSum/Zsum
     625    varyList = ['dist','det-X','det-Y','tilt','phi']
     626    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
     627        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
    612628    rings = np.concatenate((data['rings']),axis=0)
    613     p0 = [dist,cent[0],cent[1],phi,tilt]
    614     result,newWave = FitDetector(rings,p0,wave)
    615     data['distance'] = result[0]
    616     data['center'] = result[1:3]
    617     data['rotation'] = np.mod(result[3],360.0)
    618     data['tilt'] = result[4]
     629    if data['DetDepthRef']:
     630        varyList.append('dep')
     631        FitDetector(rings,varyList,parmDict)
     632    data['distance'] = parmDict['dist']
     633    data['center'] = [parmDict['det-X'],parmDict['det-Y']]
     634    data['rotation'] = np.mod(parmDict['phi'],360.0)
     635    data['tilt'] = parmDict['tilt']
     636    data['DetDepth'] = parmDict['dep']
    619637    N = len(data['ellipses'])
    620638    data['ellipses'] = []           #clear away individual ellipse fits
  • trunk/GSASIIimgGUI.py

    r827 r892  
    3838    if 'GonioAngles' not in data:
    3939        data['GonioAngles'] = [0.,0.,0.]
     40    if 'DetDepth' not in data:
     41        data['DetDepth'] = 0.
     42        data['DetDepthRef'] = False
    4043#end patch
    4144
     
    360363            except ValueError:
    361364                pass
    362             waveSel.SetValue("%6.5f" % (data['wavelength']))          #reset in case of error         
     365            waveSel.SetValue("%6.5f" % (data['wavelength']))          #reset in case of error
     366           
     367        def OnDetDepthRef(event):
     368            data['DetDepthRef'] = penSel.GetValue()
     369           
     370        def OnDetDepth(event):
     371            try:
     372                data['DetDepth'] = float(penVal.GetValue())
     373            except ValueError:
     374                pass
     375            penVal.SetValue("%6.3f" % (data['DetDepth']))          #reset in case of error                     
    363376           
    364377        calibSizer = wx.FlexGridSizer(5,2,5,5)
     
    394407        rotSel.SetBackgroundColour(VERY_LIGHT_GREY)
    395408        calibSizer.Add(rotSel,0,wx.ALIGN_CENTER_VERTICAL)
     409        penSel = wx.CheckBox(parent=G2frame.dataDisplay,label='Penetration?')
     410        calibSizer.Add(penSel,0,wx.ALIGN_CENTER_VERTICAL)
     411        penSel.Bind(wx.EVT_CHECKBOX, OnDetDepthRef)
     412        penSel.SetValue(data['DetDepthRef'])
     413        penVal = wx.TextCtrl(parent=G2frame.dataDisplay,value=("%6.5f" % (data['DetDepth'])),
     414            style=wx.TE_PROCESS_ENTER)
     415        penVal.Bind(wx.EVT_TEXT_ENTER,OnDetDepth)
     416        penVal.Bind(wx.EVT_KILL_FOCUS,OnDetDepth)
     417        calibSizer.Add(penVal,0,wx.ALIGN_CENTER_VERTICAL)             
     418       
    396419        return calibSizer
    397420   
Note: See TracChangeset for help on using the changeset viewer.