source: sans/Dev/trunk/NCNR_User_Procedures/Common/Packages/PlotManager/USANS_SlitSmearing_v40.ipf @ 1247

Last change on this file since 1247 was 1111, checked in by srkline, 5 years ago

two significant changes:

1- in the averaging, before the data is finally written out, duplicate q-values (within 0.1%) are averaged so that duplicate q-values do not appear in the data file. This has a very bad effect on the calculation of the smearing matrix (dq=0).

2-for the data from the back (high res) detector, two steps have been added to the processing at the point where the data is converted to WORK. First, a constant read noise value of 200 cts/pixel is subtracted (found from average of multiple runs with beam off) Second, a 3x3 median filter is applied to the whole image to eliminate the stray bright pixels. Some are saturatd, some are simply outliers. Very effective.
File size: 14.3 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=3.00
3#pragma IgorVersion=6.1
4
5///////
6// JAN-MAY 2018
7//
8//
9// This is a first approximation to model slit smeared VSANS data
10// - the dQv value is now a "per-Q" quantity, rather than a single global value
11// - the wavelength smearing must be accounted for by a separate numerical integral
12//   BEFORE this step is applied.
13//  -- this is significant issue for White Beam conditions
14//
15
16
17
18
19
20//Functions for doing USANS Slit smearing by method of weight matrix
21//Routines originally from J Barker fortran code
22//Translated to IGOR by M-H Kim
23//Updated to use IGOR features and integrated into SANS Macros by A J Jackson
24//
25// AJJ
26// July 26 2007 - Modified functions to work with new SANS Analysis Macros
27//     - pass basestr to functions to determine wave names and avoid globals
28//     - pass dQv to functions to  avoid globals
29//     - pass N to CalcR to avoid globals
30
31//////////////
32//
33//// July 2008 SRK
34//
35// For fitting withe cursors, the matrix must be recalculated for the exact range of the data
36// - but the inital dependency of the model (any number of them) was (were) set up to use the full matrix
37// -- so:
38//      $_res is the resolution for the dependencies. It is always the full size of the data set. If cursors are used,
39//                      it is padded with zeroes at the edges to fill.
40//
41// weights_save: is a pristine copy of the resolution matrix for the full data set, as when loaded
42// $_qt, $_qi, $_qs, $_rest: are "trimmed" sets that use the range specified by the cursors. They are created
43//                      at load time and are initially the full data range.
44//
45// there is a wave note attached to $_res that has the current point range for the matrix. This is what the
46//              $_rest matrix is expected to be too, so keep these two matrices in sync.
47//
48//      during fitting, $_rest is used in the structure if cursors are used, since the matrix must be the same dimension (N)
49//                      as the (trimmed) data range. It may be the full data range if the cursors are at the ends of the data set.
50//                      If no cursors are used, then the $_res wave is used in the structure
51//////////
52
53
54
55
56// called only by the main file loader
57//
58Function USANS_CalcWeights(basestr, dQv)
59        String basestr
60        Variable dQv
61
62        Variable/G USANS_N=numpnts($(basestr+"_q"))
63        Variable/G USANS_dQv = dQv
64        String/G root:Packages:NIST:USANS_basestr = basestr             //this is the "current data" for slope and weights
65
66        Make/D/O/N=(USANS_N,USANS_N) $(basestr+"_res")
67        Make/D/O/N=(USANS_N,USANS_N) W1mat
68        Make/D/O/N=(USANS_N,USANS_N) W2mat
69        Make/D/O/N=(USANS_N,USANS_N) Rmat
70        Wave weights = $(basestr+"_res")
71
72        //Variable/G USANS_m = EnterSlope(baseStr)
73        Variable/G USANS_m = -4
74        Variable/G USANS_slope_numpnts = 15
75        String/G trimStr=""                     //null string if NOT using cursors, "t" if yes (and re-calculating)
76       
77        Print "---- Calculating Weighting Matrix for USANS Data ----"
78       
79        EnterSlope(basestr)
80       
81        //Deal with broken Pauseforuser
82        SetDataFolder $("root:"+basestr)
83       
84        print USANS_m
85
86        if (USANS_m == 999)
87                KillWaves/Z  $(basestr+"_res"),W1Mat,W2mat,Rmat
88                return(1)
89        endif
90       
91        Variable tref = startMSTimer
92        print "Calculating W1..."
93        W1mat = (p <= q ) && (q <= USANS_N-2) ? CalcW1(p,q)  : 0
94        print "Calculating W2..."
95        W2mat = (p+1 <= q ) && (q <= USANS_N-1) ?  CalcW2(p,q) : 0
96        print "Calculating Remainders..."
97        Rmat = (q == USANS_N-1) ? CalcR(p) : 0
98//      print "Summing weights..."
99        weights = W1mat + W2mat + Rmat
100        print "Done"
101        Variable ms = stopMSTimer(tref)
102        print "Time elapsed = ", ms/1e6, "s"
103       
104        // put the point range of the matrix in a wave note (this is the full range)
105        String nStr=""
106        sprintf nStr,"P1=%d;P2=%d;",0,USANS_N-1
107        Note weights ,nStr
108       
109        // make a set of "trimmed" data that is currently the full data set
110        // but will be trimmed as needed for use with cursors
111        // this is necessary to handle the base case of cursors that use the full range of the data!
112        Duplicate/O $(baseStr+"_q") $(baseStr+"_qt")
113        Duplicate/O $(baseStr+"_i") $(baseStr+"_it")
114        Duplicate/O $(baseStr+"_s") $(baseStr+"_st")
115        Duplicate/O weights, $(basestr+"_res"+"t")
116
117        // save a copy of the untainted matrix to return to...
118        Duplicate/O weights, weights_save
119       
120        return(0)
121End
122
123// re-calculate the weighting matrix, maybe easier to do in a separate function
124//
125// supposedly I'm in the data folder, not root:, double check to make sure
126//
127Function USANS_RE_CalcWeights(baseStr,pt1,pt2)
128        String baseStr
129        Variable pt1,pt2
130
131        SetDataFolder $("root:"+baseStr)
132        NVAR USANS_dQv = USANS_dQv
133        SVAR trimStr = trimStr
134        trimStr="t"             //yes, we're re-calculating, using trimmed data sets
135       
136        Variable matN
137        matN=pt2-pt1+1          //new size of matrix
138
139        Make/D/O/N=(matN,matN) $(basestr+"_res"+trimStr)                //this is a temporary matrix
140        Make/D/O/N=(matN,matN) W1mat
141        Make/D/O/N=(matN,matN) W2mat
142        Make/D/O/N=(matN,matN) Rmat
143        Wave tmpWeights = $(basestr+"_res"+trimStr)
144        Wave fullWeights = $(basestr+"_res")
145
146        // make a trimmed set of QIS waves to match the size selected by the cursors
147        // == "t" added to the suffix
148        WAVE qval = $(baseStr+"_q")
149        WAVE ival = $(baseStr+"_i")
150        WAVE sval = $(baseStr+"_s")
151        WAVE dQvW = $(baseStr+"_dQv")
152        Duplicate/O/R=[pt1,pt2] qval $(baseStr+"_qt")
153        Duplicate/O/R=[pt1,pt2] ival $(baseStr+"_it")
154        Duplicate/O/R=[pt1,pt2] sval $(baseStr+"_st")
155        Duplicate/O/R=[pt1,pt2] dQvW $(baseStr+"_dQvt")
156        WAVE qt = $(baseStr+"_qt")                      //these are trimmed based on the cursor points
157        WAVE it = $(baseStr+"_it")
158        WAVE st = $(baseStr+"_st")
159        WAVE dQvt = $(baseStr+"_dQvt")
160       
161        //Variable/G USANS_m = EnterSlope(baseStr)
162        Variable/G USANS_m = -4
163        Variable/G USANS_slope_numpnts = 15
164        // must reset the global for baseStr to this set, otherwise it will look for the last set loaded
165        String/G root:Packages:NIST:USANS_basestr = basestr             //this is the "current data" for slope and weights
166
167       
168        Print "---- Calculating Weighting Matrix for USANS Data ----"
169       
170        EnterSlope(basestr)
171       
172        //Deal with broken Pauseforuser
173        SetDataFolder $("root:"+basestr)
174       
175        print USANS_m
176
177        if (USANS_m == 999)
178                KillWaves/Z  $(basestr+"_res"),W1Mat,W2mat,Rmat
179                return(1)
180        endif
181       
182        Variable tref = startMSTimer
183        print "Calculating W1..."
184        W1mat = (p <= q ) && (q <= matN-2) ? CalcW1(p,q)  : 0
185        print "Calculating W2..."
186        W2mat = (p+1 <= q ) && (q <= matN-1) ?  CalcW2(p,q) : 0
187        print "Calculating Remainders..."
188        Rmat = (q == matN-1) ? CalcR(p) : 0
189//      print "Summing weights..."
190        tmpWeights = W1mat + W2mat + Rmat
191        print "Done"
192        Variable ms = stopMSTimer(tref)
193        print "Time elapsed = ", ms/1e6, "s"
194       
195        //  - put the smaller matrix into the larger one (padded w/zero)
196        // get the smaller triangular weights into the proper place in the full matrix
197        //
198        // this is necessary since the original dependency (as plotted) was set up using the
199        // full data set in the STRUCT
200        fullWeights = 0
201        fullWeights[pt1,pt2][pt1,pt2] = tmpWeights[p-pt1][q-pt1]
202       
203       
204        // put the point range of the matrix in a wave note
205        String nStr=""
206        sprintf nStr,"P1=%d;P2=%d;",pt1,pt2
207        Note/K fullWeights ,nStr
208       
209        return(0)
210End
211
212// trimStr = "" if the full (original) data set is to be used, == "t" if the set is trimmed to cursors
213Function EnterSlope(baseStr)
214        String baseStr
215       
216//      NVAR USANS_m,USANS_slope_numpnts
217        NVAR USANS_m = $("root:"+baseStr+":USANS_m")
218        NVAR USANS_slope_numpnts = $("root:"+baseStr+":USANS_slope_numpnts")
219        SVAR trimStr = $("root:"+baseStr+":trimStr")
220       
221//      Variable slope=-4
222
223//      Prompt slope "Enter a slope for the file \""+ baseStr + "\""
224//      DoPrompt "Enter Slope", slope
225//              If (V_Flag)
226//                      return (999)                    //return a bogus slope if the user canceled
227//              Endif
228//      print "slope=", slope
229//      return slope
230
231        NewPanel /W=(600,300,1000,700)/N=USANS_Slope as "USANS Slope Extrapolation"
232        SetDrawLayer UserBack
233        Button button_OK,pos={270,360},size={100,20},title="Accept Slope",font="Geneva"
234        Button button_OK,proc=USANS_Slope_ButtonProc
235        Button button_SlopeCalc,pos={270,310},size={100,20},title="Calculate",font="Geneva"
236        Button button_SlopeCalc,proc=USANS_Slope_ButtonProc
237        SetVariable setvar_numpnts, pos={20,310}, size={130,19}, title="# Points",fSize=13
238        SetVariable setvar_numpnts, value= USANS_slope_numpnts
239        SetVariable setvar_Slope,pos={160,310},size={90,19},title="Slope",fSize=13
240        SetVariable setvar_Slope,limits={-inf,inf,0},value= USANS_m
241        SetVariable setvar_slope, proc=USANS_ManualSlope
242       
243        Display/W=(9,6,402,305)/HOST=USANS_Slope $(basestr+"_i"+trimStr) vs $(basestr+"_q"+trimStr)
244        RenameWindow #,SlopePlot
245        ErrorBars/T=0/W=USANS_Slope#SlopePlot $(basestr+"_i"+trimStr), Y wave=($(basestr+"_s"+trimStr),$(basestr+"_s"+trimStr))
246        ModifyGraph log=1
247        ModifyGraph mode=3,msize=3,marker=1,rgb=(0,0,0)
248//      legend
249        SetActiveSubwindow ##
250       
251        //Print TraceNameList("USANS_Slope#SlopePlot",";",1)
252       
253        //USANS_CalculateSlope(basestr,USANS_slope_numpnts)
254       
255        PauseForUser  USANS_Slope       
256       
257       
258//      return slope
259       
260End
261
262Function USANS_CalculateSlope(basestr, nend)
263        String basestr
264        Variable nend
265       
266//      NVAR USANS_m
267        NVAR USANS_m = $("root:"+baseStr+":USANS_m")
268        SVAR trimStr = $("root:"+baseStr+":trimStr")
269       
270        Wave iw = $(basestr+"_i"+trimStr)
271        Wave qw = $(basestr+"_q"+trimStr)
272        Wave sw = $(basestr+"_s"+trimStr)
273
274        Variable num_extr=25
275        // Make extra waves for extrapolation
276        // Taken from DSM_SetExtrWaves
277               
278        Make/O/D/N=(num_extr) extr_hqq,extr_hqi
279        extr_hqi=1              //default values
280
281        //set the q-range
282        Variable qmax,num
283
284        num=numpnts(qw)
285        qmax=6*qw[num-1]
286       
287        // num-1-nend puts the extrapolation over the data set, as well as extending it
288        extr_hqq = qw[num-1-nend] + x * (qmax-qw[num-1-nend])/num_extr
289               
290       
291        // Modifed from DSM_DoExtraploate in LakeDesmearing_JB.ipf
292        // which is part of the USANS Reduction macros
293       
294               
295        //      Wave extr_lqi=extr_lqi
296        //      Wave extr_lqq=extr_lqq
297        //      Wave extr_hqi=extr_hqi
298        //      Wave extr_hqq=extr_hqq
299                Variable/G V_FitMaxIters=300
300                Variable/G V_fitOptions=4               //suppress the iteration window
301                Variable retval
302                num=numpnts(iw)
303       
304                Make/O/D P_coef={0,1,-4}                        //input
305                Make/O/T Constr={"K2<0","K2 > -8"}
306                //(set background to zero and hold fixed)
307                // initial guess
308                P_coef[1] = iw[num-1]/qw[num-1]^P_coef[2]
309               
310                CurveFit/H="100" Power kwCWave=P_coef  iw[(num-1-nend),(num-1)] /X=qw /I=1 /W=sw /C=constr
311                extr_hqi=P_coef[0]+P_coef[1]*extr_hqq^P_coef[2]
312       
313                AppendToGraph/W=USANS_Slope#SlopePlot extr_hqi vs extr_hqq
314                ModifyGraph/W=USANS_Slope#SlopePlot lsize(extr_hqi)=2
315       
316                Printf "Smeared Power law exponent = %g\r",P_coef[2]
317                Printf "**** For Desmearing, use a Power law exponent of %5.1f\r",P_coef[2]-1
318       
319                retVal = P_coef[2]-1                   
320                return(retVal)
321               
322End
323
324// adjust the slope manually on the graph (not fitted!!) just for visualization
325//
326// can be turned off by removing the action proc in the setVariable
327//
328Function USANS_ManualSlope(ctrlName,varNum,varStr,varName) : SetVariableControl
329        String ctrlName
330        Variable varNum
331        String varStr           
332        String varName
333               
334        SVAR baseStr = root:Packages:NIST:USANS_basestr
335        NVAR nFit = $("root:"+baseStr+":USANS_slope_numpnts")
336        NVAR USANS_m = $("root:"+baseStr+":USANS_m")
337        WAVE/Z extr_hqq = $("root:"+baseStr+":extr_hqq")
338        WAVE/Z extr_hqi = $("root:"+baseStr+":extr_hqi")
339        SVAR trimStr = $("root:"+baseStr+":trimStr")
340       
341        if(waveExists(extr_hqq)==0 || waveExists(extr_hqi)==0)
342                return(0)               // encourage the user to try to fit first...
343        endif
344       
345        Wave iw = $("root:"+basestr+":"+basestr+"_i"+trimStr)
346        Wave qw = $("root:"+basestr+":"+basestr+"_q"+trimStr)
347
348        Variable matchPt,num=numpnts(iw),correctedSlope
349        correctedSlope = USANS_m + 1
350       
351        matchPt = iw[num-1-nFit]/(qw[num-1-nFit]^correctedSlope)
352       
353        extr_hqi=matchPt*extr_hqq^correctedSlope
354       
355        return(0)
356End
357
358Function USANS_Slope_ButtonProc(ctrlName) : ButtonControl
359        String ctrlName
360       
361        SVAR basestr =  root:Packages:NIST:USANS_basestr
362       
363        NVAR USANS_m = $("root:"+baseStr+":USANS_m")
364        NVAR USANS_slope_numpnts = $("root:"+baseStr+":USANS_slope_numpnts")
365       
366        strswitch (ctrlName)
367                case "button_OK":
368                        ControlUpdate/W=USANS_Slope setvar_Slope
369                        DoWindow/K USANS_Slope
370                        break
371                case "button_SlopeCalc":
372                        USANS_m = USANS_CalculateSlope(basestr, USANS_slope_numpnts)
373                        ControlUpdate/W=USANS_Slope setvar_Slope
374                        break
375        endswitch
376         
377
378End
379
380// MAY 2018
381// updated to use dQv as a per-q wave, not a single global value
382Function CalcW1(i,j)
383        Variable i,j
384       
385        SVAR USANS_basestr=root:Packages:NIST:USANS_basestr
386        SVAR trimStr = $("root:"+USANS_baseStr+":trimStr")
387        SetDataFolder $("root:"+USANS_basestr)
388
389       
390// replaced with a wave
391//      NVAR dQv = USANS_dQv
392       
393        Variable UU,UL,dqj,rU,rL,wU,wL,dqw
394        Wave Qval = $(USANS_basestr+"_q"+trimStr)
395        Wave dQv = $(USANS_basestr+"_dQv"+trimStr)
396       
397        UU =sqrt(Qval[j+1]^2-Qval[i]^2)
398        UL = sqrt(Qval[j]^2-Qval[i]^2)
399        dqw = Qval[j+1]-Qval[j]
400        rU = sqrt(UU^2+Qval[i]^2)
401        rL = sqrt(UL^2+Qval[i]^2)
402       
403        wU = (1.0/dQv[i])*(Qval[j+1]*UU/dqw - 0.5*UU*rU/dqw - 0.5*Qval[i]^2*ln(UU+rU)/dqw )
404        wL = (1.0/dQv[i])*(Qval[j+1]*UL/dqw - 0.5*UL*rL/dqw - 0.5*Qval[i]^2*ln(UL+rL)/dqw )
405       
406        Return wU-wL
407
408End
409
410// MAY 2018
411// updated to use dQv as a per-q wave, not a single global value
412Function CalcW2(i,j)
413        Variable i,j
414       
415        SVAR USANS_basestr=root:Packages:NIST:USANS_basestr
416        SVAR trimStr = $("root:"+USANS_baseStr+":trimStr")
417        SetDataFolder $("root:"+USANS_basestr)
418
419
420// replaced with a wave
421//      NVAR dQv = USANS_dQv
422       
423        variable UU,UL,dqw,rU,rL,wU,wL
424       
425        Wave Qval = $(USANS_basestr+"_q"+trimStr)
426        Wave dQv = $(USANS_basestr+"_dQv"+trimStr)
427
428        UU = sqrt(Qval[j]^2-Qval[i]^2)                 
429        UL = sqrt(Qval[j-1]^2-Qval[i]^2)               
430        dqw = Qval[j]-Qval[j-1]                 
431        rU = sqrt(UU^2+Qval[i]^2)
432        rL = sqrt(UL^2+Qval[i]^2)
433        wU = (1.0/dQv[i])*( -Qval[j-1]*UU/dqw + 0.5*UU*rU/dqw + 0.5*Qval[i]^2*ln(UU+rU)/dqw )
434        wL = (1.0/dQv[i])*( -Qval[j-1]*UL/dqw + 0.5*UL*rL/dqw + 0.5*Qval[i]^2*ln(UL+rL)/dqw )
435
436        Return wU-wL
437
438End
439
440// MAY 2018
441// updated to use dQv as a per-q wave, not a single global value
442Function CalcR(i)
443        Variable i
444
445        SVAR USANS_basestr=root:Packages:NIST:USANS_basestr
446        SVAR trimStr = $("root:"+USANS_baseStr+":trimStr")
447        SetDataFolder $("root:"+USANS_basestr)
448
449
450        NVAR m = USANS_m
451        NVAR N = USANS_N
452// replaced with a wave
453//      NVAR dQv = USANS_dQv
454       
455        Variable retval
456        Wave Qval = $(USANS_basestr+"_q"+trimStr)
457        Wave Ival = $(USANS_basestr+"_i"+trimStr)
458        Variable/G USANS_intQpt = Qval[i]
459        Wave dQv = $(USANS_basestr+"_dQv"+trimStr)
460
461       
462        Variable lower = sqrt(qval[N-1]^2-qval[i]^2)
463        Variable upper
464        upper = dQv[i]
465
466        if (i == N)
467                lower = 0
468        endif
469       
470        retval = Integrate1D(Remainder,lower,upper)
471       
472        retval *= 1/dQv[i]
473       
474        Return retval
475
476End
477
478Function Remainder(i)
479       
480        Variable i
481       
482        SVAR USANS_basestr=root:Packages:NIST:USANS_basestr
483        SVAR trimStr = $("root:"+USANS_baseStr+":trimStr")
484        SetDataFolder $("root:"+USANS_basestr)
485
486        NVAR m = USANS_m
487        NVAR qi = USANS_intQpt
488        NVAR N = USANS_N
489        WAVE Qval = $(USANS_basestr+"_q"+trimStr)       
490        Variable retVal
491       
492        retVal=Qval[N-1]^(-m)*(i^2+qi^2)^(m/2)
493
494        return retval
495
496End
Note: See TracBrowser for help on using the repository browser.