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

Last change on this file since 1100 was 1100, checked in by srkline, 4 years ago

Significant restructuring of V_ExecuteProtocol to make the logic of the flow more tractable and adaptable in the future.

Added major change to VSANS event mode reduction panel to allow F and M binned events to be saved to a copy of the correspondng data file. A new data reduction panel for multiple slice reduction is now presented for Event data. This allows reduction of a single slice (all 8 F,M panels, back ignored), or all of the slices.

File size: 14.2 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        Duplicate/O/R=[pt1,pt2] qval $(baseStr+"_qt")
152        Duplicate/O/R=[pt1,pt2] ival $(baseStr+"_it")
153        Duplicate/O/R=[pt1,pt2] sval $(baseStr+"_st")
154        WAVE qt = $(baseStr+"_qt")                      //these are trimmed based on the cursor points
155        WAVE it = $(baseStr+"_it")
156        WAVE st = $(baseStr+"_st")
157       
158        //Variable/G USANS_m = EnterSlope(baseStr)
159        Variable/G USANS_m = -4
160        Variable/G USANS_slope_numpnts = 15
161        // must reset the global for baseStr to this set, otherwise it will look for the last set loaded
162        String/G root:Packages:NIST:USANS_basestr = basestr             //this is the "current data" for slope and weights
163
164       
165        Print "---- Calculating Weighting Matrix for USANS Data ----"
166       
167        EnterSlope(basestr)
168       
169        //Deal with broken Pauseforuser
170        SetDataFolder $("root:"+basestr)
171       
172        print USANS_m
173
174        if (USANS_m == 999)
175                KillWaves/Z  $(basestr+"_res"),W1Mat,W2mat,Rmat
176                return(1)
177        endif
178       
179        Variable tref = startMSTimer
180        print "Calculating W1..."
181        W1mat = (p <= q ) && (q <= matN-2) ? CalcW1(p,q)  : 0
182        print "Calculating W2..."
183        W2mat = (p+1 <= q ) && (q <= matN-1) ?  CalcW2(p,q) : 0
184        print "Calculating Remainders..."
185        Rmat = (q == matN-1) ? CalcR(p) : 0
186//      print "Summing weights..."
187        tmpWeights = W1mat + W2mat + Rmat
188        print "Done"
189        Variable ms = stopMSTimer(tref)
190        print "Time elapsed = ", ms/1e6, "s"
191       
192        //  - put the smaller matrix into the larger one (padded w/zero)
193        // get the smaller triangular weights into the proper place in the full matrix
194        //
195        // this is necessary since the original dependency (as plotted) was set up using the
196        // full data set in the STRUCT
197        fullWeights = 0
198        fullWeights[pt1,pt2][pt1,pt2] = tmpWeights[p-pt1][q-pt1]
199       
200       
201        // put the point range of the matrix in a wave note
202        String nStr=""
203        sprintf nStr,"P1=%d;P2=%d;",pt1,pt2
204        Note/K fullWeights ,nStr
205       
206        return(0)
207End
208
209// trimStr = "" if the full (original) data set is to be used, == "t" if the set is trimmed to cursors
210Function EnterSlope(baseStr)
211        String baseStr
212       
213//      NVAR USANS_m,USANS_slope_numpnts
214        NVAR USANS_m = $("root:"+baseStr+":USANS_m")
215        NVAR USANS_slope_numpnts = $("root:"+baseStr+":USANS_slope_numpnts")
216        SVAR trimStr = $("root:"+baseStr+":trimStr")
217       
218//      Variable slope=-4
219
220//      Prompt slope "Enter a slope for the file \""+ baseStr + "\""
221//      DoPrompt "Enter Slope", slope
222//              If (V_Flag)
223//                      return (999)                    //return a bogus slope if the user canceled
224//              Endif
225//      print "slope=", slope
226//      return slope
227
228        NewPanel /W=(600,300,1000,700)/N=USANS_Slope as "USANS Slope Extrapolation"
229        SetDrawLayer UserBack
230        Button button_OK,pos={270,360},size={100,20},title="Accept Slope",font="Geneva"
231        Button button_OK,proc=USANS_Slope_ButtonProc
232        Button button_SlopeCalc,pos={270,310},size={100,20},title="Calculate",font="Geneva"
233        Button button_SlopeCalc,proc=USANS_Slope_ButtonProc
234        SetVariable setvar_numpnts, pos={20,310}, size={130,19}, title="# Points",fSize=13
235        SetVariable setvar_numpnts, value= USANS_slope_numpnts
236        SetVariable setvar_Slope,pos={160,310},size={90,19},title="Slope",fSize=13
237        SetVariable setvar_Slope,limits={-inf,inf,0},value= USANS_m
238        SetVariable setvar_slope, proc=USANS_ManualSlope
239       
240        Display/W=(9,6,402,305)/HOST=USANS_Slope $(basestr+"_i"+trimStr) vs $(basestr+"_q"+trimStr)
241        RenameWindow #,SlopePlot
242        ErrorBars/T=0/W=USANS_Slope#SlopePlot $(basestr+"_i"+trimStr), Y wave=($(basestr+"_s"+trimStr),$(basestr+"_s"+trimStr))
243        ModifyGraph log=1
244        ModifyGraph mode=3,msize=3,marker=1,rgb=(0,0,0)
245//      legend
246        SetActiveSubwindow ##
247       
248        //Print TraceNameList("USANS_Slope#SlopePlot",";",1)
249       
250        //USANS_CalculateSlope(basestr,USANS_slope_numpnts)
251       
252        PauseForUser  USANS_Slope       
253       
254       
255//      return slope
256       
257End
258
259Function USANS_CalculateSlope(basestr, nend)
260        String basestr
261        Variable nend
262       
263//      NVAR USANS_m
264        NVAR USANS_m = $("root:"+baseStr+":USANS_m")
265        SVAR trimStr = $("root:"+baseStr+":trimStr")
266       
267        Wave iw = $(basestr+"_i"+trimStr)
268        Wave qw = $(basestr+"_q"+trimStr)
269        Wave sw = $(basestr+"_s"+trimStr)
270
271        Variable num_extr=25
272        // Make extra waves for extrapolation
273        // Taken from DSM_SetExtrWaves
274               
275        Make/O/D/N=(num_extr) extr_hqq,extr_hqi
276        extr_hqi=1              //default values
277
278        //set the q-range
279        Variable qmax,num
280
281        num=numpnts(qw)
282        qmax=6*qw[num-1]
283       
284        // num-1-nend puts the extrapolation over the data set, as well as extending it
285        extr_hqq = qw[num-1-nend] + x * (qmax-qw[num-1-nend])/num_extr
286               
287       
288        // Modifed from DSM_DoExtraploate in LakeDesmearing_JB.ipf
289        // which is part of the USANS Reduction macros
290       
291               
292        //      Wave extr_lqi=extr_lqi
293        //      Wave extr_lqq=extr_lqq
294        //      Wave extr_hqi=extr_hqi
295        //      Wave extr_hqq=extr_hqq
296                Variable/G V_FitMaxIters=300
297                Variable/G V_fitOptions=4               //suppress the iteration window
298                Variable retval
299                num=numpnts(iw)
300       
301                Make/O/D P_coef={0,1,-4}                        //input
302                Make/O/T Constr={"K2<0","K2 > -8"}
303                //(set background to zero and hold fixed)
304                // initial guess
305                P_coef[1] = iw[num-1]/qw[num-1]^P_coef[2]
306               
307                CurveFit/H="100" Power kwCWave=P_coef  iw[(num-1-nend),(num-1)] /X=qw /I=1 /W=sw /C=constr
308                extr_hqi=P_coef[0]+P_coef[1]*extr_hqq^P_coef[2]
309       
310                AppendToGraph/W=USANS_Slope#SlopePlot extr_hqi vs extr_hqq
311                ModifyGraph/W=USANS_Slope#SlopePlot lsize(extr_hqi)=2
312       
313                Printf "Smeared Power law exponent = %g\r",P_coef[2]
314                Printf "**** For Desmearing, use a Power law exponent of %5.1f\r",P_coef[2]-1
315       
316                retVal = P_coef[2]-1                   
317                return(retVal)
318               
319End
320
321// adjust the slope manually on the graph (not fitted!!) just for visualization
322//
323// can be turned off by removing the action proc in the setVariable
324//
325Function USANS_ManualSlope(ctrlName,varNum,varStr,varName) : SetVariableControl
326        String ctrlName
327        Variable varNum
328        String varStr           
329        String varName
330               
331        SVAR baseStr = root:Packages:NIST:USANS_basestr
332        NVAR nFit = $("root:"+baseStr+":USANS_slope_numpnts")
333        NVAR USANS_m = $("root:"+baseStr+":USANS_m")
334        WAVE/Z extr_hqq = $("root:"+baseStr+":extr_hqq")
335        WAVE/Z extr_hqi = $("root:"+baseStr+":extr_hqi")
336        SVAR trimStr = $("root:"+baseStr+":trimStr")
337       
338        if(waveExists(extr_hqq)==0 || waveExists(extr_hqi)==0)
339                return(0)               // encourage the user to try to fit first...
340        endif
341       
342        Wave iw = $("root:"+basestr+":"+basestr+"_i"+trimStr)
343        Wave qw = $("root:"+basestr+":"+basestr+"_q"+trimStr)
344
345        Variable matchPt,num=numpnts(iw),correctedSlope
346        correctedSlope = USANS_m + 1
347       
348        matchPt = iw[num-1-nFit]/(qw[num-1-nFit]^correctedSlope)
349       
350        extr_hqi=matchPt*extr_hqq^correctedSlope
351       
352        return(0)
353End
354
355Function USANS_Slope_ButtonProc(ctrlName) : ButtonControl
356        String ctrlName
357       
358        SVAR basestr =  root:Packages:NIST:USANS_basestr
359       
360        NVAR USANS_m = $("root:"+baseStr+":USANS_m")
361        NVAR USANS_slope_numpnts = $("root:"+baseStr+":USANS_slope_numpnts")
362       
363        strswitch (ctrlName)
364                case "button_OK":
365                        ControlUpdate/W=USANS_Slope setvar_Slope
366                        DoWindow/K USANS_Slope
367                        break
368                case "button_SlopeCalc":
369                        USANS_m = USANS_CalculateSlope(basestr, USANS_slope_numpnts)
370                        ControlUpdate/W=USANS_Slope setvar_Slope
371                        break
372        endswitch
373         
374
375End
376
377// MAY 2018
378// updated to use dQv as a per-q wave, not a single global value
379Function CalcW1(i,j)
380        Variable i,j
381       
382        SVAR USANS_basestr=root:Packages:NIST:USANS_basestr
383        SVAR trimStr = $("root:"+USANS_baseStr+":trimStr")
384        SetDataFolder $("root:"+USANS_basestr)
385
386       
387// replaced with a wave
388//      NVAR dQv = USANS_dQv
389       
390        Variable UU,UL,dqj,rU,rL,wU,wL,dqw
391        Wave Qval = $(USANS_basestr+"_q"+trimStr)
392        Wave dQv = $(USANS_basestr+"_dQv"+trimStr)
393       
394        UU =sqrt(Qval[j+1]^2-Qval[i]^2)
395        UL = sqrt(Qval[j]^2-Qval[i]^2)
396        dqw = Qval[j+1]-Qval[j]
397        rU = sqrt(UU^2+Qval[i]^2)
398        rL = sqrt(UL^2+Qval[i]^2)
399       
400        wU = (1.0/dQv[i])*(Qval[j+1]*UU/dqw - 0.5*UU*rU/dqw - 0.5*Qval[i]^2*ln(UU+rU)/dqw )
401        wL = (1.0/dQv[i])*(Qval[j+1]*UL/dqw - 0.5*UL*rL/dqw - 0.5*Qval[i]^2*ln(UL+rL)/dqw )
402       
403        Return wU-wL
404
405End
406
407// MAY 2018
408// updated to use dQv as a per-q wave, not a single global value
409Function CalcW2(i,j)
410        Variable i,j
411       
412        SVAR USANS_basestr=root:Packages:NIST:USANS_basestr
413        SVAR trimStr = $("root:"+USANS_baseStr+":trimStr")
414        SetDataFolder $("root:"+USANS_basestr)
415
416
417// replaced with a wave
418//      NVAR dQv = USANS_dQv
419       
420        variable UU,UL,dqw,rU,rL,wU,wL
421       
422        Wave Qval = $(USANS_basestr+"_q"+trimStr)
423        Wave dQv = $(USANS_basestr+"_dQv"+trimStr)
424
425        UU = sqrt(Qval[j]^2-Qval[i]^2)                 
426        UL = sqrt(Qval[j-1]^2-Qval[i]^2)               
427        dqw = Qval[j]-Qval[j-1]                 
428        rU = sqrt(UU^2+Qval[i]^2)
429        rL = sqrt(UL^2+Qval[i]^2)
430        wU = (1.0/dQv[i])*( -Qval[j-1]*UU/dqw + 0.5*UU*rU/dqw + 0.5*Qval[i]^2*ln(UU+rU)/dqw )
431        wL = (1.0/dQv[i])*( -Qval[j-1]*UL/dqw + 0.5*UL*rL/dqw + 0.5*Qval[i]^2*ln(UL+rL)/dqw )
432
433        Return wU-wL
434
435End
436
437// MAY 2018
438// updated to use dQv as a per-q wave, not a single global value
439Function CalcR(i)
440        Variable i
441
442        SVAR USANS_basestr=root:Packages:NIST:USANS_basestr
443        SVAR trimStr = $("root:"+USANS_baseStr+":trimStr")
444        SetDataFolder $("root:"+USANS_basestr)
445
446
447        NVAR m = USANS_m
448        NVAR N = USANS_N
449// replaced with a wave
450//      NVAR dQv = USANS_dQv
451       
452        Variable retval
453        Wave Qval = $(USANS_basestr+"_q"+trimStr)
454        Wave Ival = $(USANS_basestr+"_i"+trimStr)
455        Variable/G USANS_intQpt = Qval[i]
456        Wave dQv = $(USANS_basestr+"_dQv"+trimStr)
457
458       
459        Variable lower = sqrt(qval[N-1]^2-qval[i]^2)
460        Variable upper
461        upper = dQv[i]
462
463        if (i == N)
464                lower = 0
465        endif
466       
467        retval = Integrate1D(Remainder,lower,upper)
468       
469        retval *= 1/dQv[i]
470       
471        Return retval
472
473End
474
475Function Remainder(i)
476       
477        Variable i
478       
479        SVAR USANS_basestr=root:Packages:NIST:USANS_basestr
480        SVAR trimStr = $("root:"+USANS_baseStr+":trimStr")
481        SetDataFolder $("root:"+USANS_basestr)
482
483        NVAR m = USANS_m
484        NVAR qi = USANS_intQpt
485        NVAR N = USANS_N
486        WAVE Qval = $(USANS_basestr+"_q"+trimStr)       
487        Variable retVal
488       
489        retVal=Qval[N-1]^(-m)*(i^2+qi^2)^(m/2)
490
491        return retval
492
493End
Note: See TracBrowser for help on using the repository browser.