source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_WhiteBeamSmear.ipf @ 1130

Last change on this file since 1130 was 1100, checked in by srkline, 5 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: 17.5 KB
Line 
1#pragma TextEncoding = "MacRoman"
2#pragma rtGlobals=3             // Use modern global access method and strict wave access.
3
4//
5//
6// testing routines to compare various integration methods and approximations
7// for calculating the resolution smearing from the white beam wavelength distribution
8//
9//
10//
11// IntegrateFn_N is something that I wrote (in GaussUtils) for quadrature with any number of
12//  points (user-selected)
13//
14// 2018:
15// my quadrature and the built-in function are equivalent. Romberg may be useful in some cases
16// especially for multiple integrals. then number of points and timing can be optimized. But either
17// method can be used. 
18//     
19//              answer = IntegrateFn_N(V_WB_testKernel,loLim,upLim,cw,qVals,nord)
20//      answer_Rom_WB = Integrate_BuiltIn(cw,loLim,upLim,qVals)
21
22// using a matrix multiplication for this calculation of the white beam wavelength smearing is NOT
23// recommended -- the calculation is not nearly accurate enough.
24//
25//
26// Using my built-in quadrature routines (see V_TestWavelengthIntegral) may be of use when
27// writing fitting functions for all of these cases. The built-in Integrate may be limited
28//
29// TODO -- beware what might happen to the calculations since there is a single global string
30//   containing the function name.
31//
32// TODO:
33// -- a significant problem with using the coef waves that are used in the wrapper are that
34//   they are set up with a dependency, so doing the WB calculation also does the "regular"
35//   smeared calculation, doubling the time required...
36
37
38//
39// needs V_DummyFunctions for the FUNCREF to work - since it fails if I simply call the XOP
40//
41//
42// SANSModel_proto(w,x)         is in GaussUtils_v40.ipf
43//
44// FUNCREF SANSModel_proto fcn
45//
46
47
48
49// call the calculation
50// see DoTheFitButton in Wrapper_v40.ipf
51//
52//
53Proc V_Calc_WB_Smearing_top()
54
55        String folderStr,funcStr,coefStr
56       
57        ControlInfo/W=WrapperPanel popup_0
58        folderStr=S_Value
59       
60        ControlInfo/W=WrapperPanel popup_1
61        funcStr=S_Value
62       
63        ControlInfo/W=WrapperPanel popup_2
64        coefStr=S_Value
65
66        V_DoWavelengthIntegral_top(folderStr,funcStr,coefStr)
67
68        SetDataFolder root:
69End
70
71
72Proc V_Calc_WB_Smearing_mid()
73
74        String folderStr,funcStr,coefStr
75       
76        ControlInfo/W=WrapperPanel popup_0
77        folderStr=S_Value
78       
79        ControlInfo/W=WrapperPanel popup_1
80        funcStr=S_Value
81       
82        ControlInfo/W=WrapperPanel popup_2
83        coefStr=S_Value
84
85        V_DoWavelengthIntegral_mid(folderStr,funcStr,coefStr)
86
87        SetDataFolder root:
88End
89
90Proc V_Calc_WB_Smearing_interp()
91
92        String folderStr,funcStr,coefStr
93       
94        ControlInfo/W=WrapperPanel popup_0
95        folderStr=S_Value
96       
97        ControlInfo/W=WrapperPanel popup_1
98        funcStr=S_Value
99       
100        ControlInfo/W=WrapperPanel popup_2
101        coefStr=S_Value
102
103        V_DoWavelengthIntegral_interp(folderStr,funcStr,coefStr)
104
105        SetDataFolder root:
106End
107
108Proc V_Calc_WB_Smearing_triang()
109
110        String folderStr,funcStr,coefStr
111       
112        ControlInfo/W=WrapperPanel popup_0
113        folderStr=S_Value
114       
115        ControlInfo/W=WrapperPanel popup_1
116        funcStr=S_Value
117       
118        ControlInfo/W=WrapperPanel popup_2
119        coefStr=S_Value
120
121        V_DoWavelengthIntegral_triang(folderStr,funcStr,coefStr)
122
123        SetDataFolder root:
124End
125
126
127
128// uses built-in Integrate1d()
129//
130Function V_DoWavelengthIntegral_top(folderStr,funcStr,coefStr)
131        String folderStr,funcStr,coefStr
132
133        SetDataFolder $("root:"+folderStr)
134               
135        // gather the input waves
136        WAVE qVals = $(folderStr+"_q")
137//      WAVE cw = smear_coef_BroadPeak
138        WAVE cw = $coefStr
139       
140        funcStr = V_getXFuncStrFromCoef(cw)+"_"         //get the modelX name, tag on "_"
141        String/G root:gFunctionString = funcStr         // need a global reference to pass to Integrate1D
142       
143        // make a wave for the answer
144        Duplicate/O qvals answer_Rom_WB_top
145       
146        // do the integration
147        Variable loLim,upLim
148               
149        // define limits based on lo/mean, hi/mean of the wavelength distribution
150        // using the empirical definition, "top" of the peaks
151        loLim = 3.37/5.3
152        upLim = 8.25/5.3
153       
154//      // using the "middle"
155//      loLim = 3.37/5.3
156//      upLim = 8.37/5.3       
157//     
158//      // using the interpolated distribution (must change the function call)
159//      lolim = 3/5.3
160//      uplim = 9/5.3
161
162// using the "triangular" distribution (must change the function call)
163//      loLim = 4/5.3
164//      upLim = 8/5.3
165       
166        answer_Rom_WB_top = V_Integrate_BuiltIn_top(cw,loLim,upLim,qVals)
167
168// why do I need this? Is this because this is defined as the mean of the distribution
169//  and is needed to normalize the integral? verify this on paper.     
170        answer_Rom_WB_top *= 5.3
171
172// normalize the integral       
173        answer_Rom_WB_top /= 20926              // "top"  of peaks
174//      answer_Rom_WB /= 19933          // "middle"  of peaks
175//      answer_Rom_WB /= 20051          // interpolated distribution
176//      answer_Rom_WB /= 1              // triangular distribution (it's already normalized)
177
178// additional normalization???
179        answer_Rom_WB_top /= 1.05               //
180       
181        SetDataFolder root:
182       
183        return 0
184End
185
186
187// uses built-in Integrate1d()
188//
189Function V_DoWavelengthIntegral_mid(folderStr,funcStr,coefStr)
190        String folderStr,funcStr,coefStr
191
192        SetDataFolder $("root:"+folderStr)
193               
194        // gather the input waves
195        WAVE qVals = $(folderStr+"_q")
196//      WAVE cw = smear_coef_BroadPeak
197        WAVE cw = $coefStr
198       
199        funcStr = V_getXFuncStrFromCoef(cw)+"_"         //get the modelX name, tag on "_"
200        String/G root:gFunctionString = funcStr         // need a global reference to pass to Integrate1D
201       
202        // make a wave for the answer
203        Duplicate/O qvals answer_Rom_WB_mid
204       
205        // do the integration
206        Variable loLim,upLim
207               
208        // define limits based on lo/mean, hi/mean of the wavelength distribution
209        // using the empirical definition, "top" of the peaks
210//      loLim = 3.37/5.3
211//      upLim = 8.25/5.3
212       
213//      // using the "middle"
214        loLim = 3.37/5.3
215        upLim = 8.37/5.3       
216//     
217//      // using the interpolated distribution (must change the function call)
218//      lolim = 3/5.3
219//      uplim = 9/5.3
220
221// using the "triangular" distribution (must change the function call)
222//      loLim = 4/5.3
223//      upLim = 8/5.3
224       
225        answer_Rom_WB_mid = V_Integrate_BuiltIn_mid(cw,loLim,upLim,qVals)
226
227// why do I need this? Is this because this is defined as the mean of the distribution
228//  and is needed to normalize the integral? verify this on paper.     
229        answer_Rom_WB_mid *= 5.3
230
231// normalize the integral       
232//      answer_Rom_WB /= 20926          // "top"  of peaks
233        answer_Rom_WB_mid /= 19933              // "middle"  of peaks
234//      answer_Rom_WB /= 20051          // interpolated distribution
235//      answer_Rom_WB /= 1              // triangular distribution (it's already normalized)
236
237// additional normalization???
238        answer_Rom_WB_mid /= 1.05               // "middle"  of peaks
239       
240        SetDataFolder root:
241       
242        return 0
243End
244
245// uses built-in Integrate1d()
246//
247Function V_DoWavelengthIntegral_interp(folderStr,funcStr,coefStr)
248        String folderStr,funcStr,coefStr
249
250        SetDataFolder $("root:"+folderStr)
251               
252        // gather the input waves
253        WAVE qVals = $(folderStr+"_q")
254//      WAVE cw = smear_coef_BroadPeak
255        WAVE cw = $coefStr
256       
257        funcStr = V_getXFuncStrFromCoef(cw)+"_"         //get the modelX name, tag on "_"
258        String/G root:gFunctionString = funcStr         // need a global reference to pass to Integrate1D
259       
260        // make a wave for the answer
261        Duplicate/O qvals answer_Rom_WB_interp
262       
263        // do the integration
264        Variable loLim,upLim
265               
266        // define limits based on lo/mean, hi/mean of the wavelength distribution
267        // using the empirical definition, "top" of the peaks
268//      loLim = 3.37/5.3
269//      upLim = 8.25/5.3
270       
271//      // using the "middle"
272//      loLim = 3.37/5.3
273//      upLim = 8.37/5.3       
274//     
275//      // using the interpolated distribution (must change the function call)
276        lolim = 3/5.3
277        uplim = 9/5.3
278
279// using the "triangular" distribution (must change the function call)
280//      loLim = 4/5.3
281//      upLim = 8/5.3
282       
283        answer_Rom_WB_interp = V_Integrate_BuiltIn_interp(cw,loLim,upLim,qVals)
284
285// why do I need this? Is this because this is defined as the mean of the distribution
286//  and is needed to normalize the integral? verify this on paper.     
287        answer_Rom_WB_interp *= 5.3
288
289// normalize the integral       
290//      answer_Rom_WB /= 20926          // "top"  of peaks
291//      answer_Rom_WB /= 19933          // "middle"  of peaks
292        answer_Rom_WB_interp /= 20051           // interpolated distribution
293//      answer_Rom_WB /= 1              // triangular distribution (it's already normalized)
294
295// additional normalization???
296        answer_Rom_WB_interp /= 1.05            // "middle"  of peaks
297       
298        SetDataFolder root:
299       
300        return 0
301End
302
303// uses built-in Integrate1d()
304//
305Function V_DoWavelengthIntegral_triang(folderStr,funcStr,coefStr)
306        String folderStr,funcStr,coefStr
307
308        SetDataFolder $("root:"+folderStr)
309               
310        // gather the input waves
311        WAVE qVals = $(folderStr+"_q")
312//      WAVE cw = smear_coef_BroadPeak
313        WAVE cw = $coefStr
314       
315        funcStr = V_getXFuncStrFromCoef(cw)+"_"         //get the modelX name, tag on "_"
316        String/G root:gFunctionString = funcStr         // need a global reference to pass to Integrate1D
317       
318        // make a wave for the answer
319        Duplicate/O qvals answer_Rom_WB_triang
320       
321        // do the integration
322        Variable loLim,upLim
323               
324        // define limits based on lo/mean, hi/mean of the wavelength distribution
325        // using the empirical definition, "top" of the peaks
326//      loLim = 3.37/5.3
327//      upLim = 8.25/5.3
328       
329//      // using the "middle"
330//      loLim = 3.37/5.3
331//      upLim = 8.37/5.3       
332//     
333//      // using the interpolated distribution (must change the function call)
334//      lolim = 3/5.3
335//      uplim = 9/5.3
336
337// using the "triangular" distribution (must change the function call)
338        loLim = 4/5.3
339        upLim = 8/5.3
340       
341        answer_Rom_WB_triang = V_Integrate_BuiltIn_triangle(cw,loLim,upLim,qVals)
342
343// why do I need this? Is this because this is defined as the mean of the distribution
344//  and is needed to normalize the integral? verify this on paper.     
345        answer_Rom_WB_triang *= 5.3
346
347// normalize the integral       
348//      answer_Rom_WB /= 20926          // "top"  of peaks
349//      answer_Rom_WB /= 19933          // "middle"  of peaks
350//      answer_Rom_WB /= 20051          // interpolated distribution
351        answer_Rom_WB_triang /= 1               // triangular distribution (it's already normalized)
352
353// additional normalization???
354        answer_Rom_WB_triang /= 1.1             //
355       
356        SetDataFolder root:
357       
358        return 0
359End
360
361
362//
363// not used anymore - the built-in works fine, but this
364// may be of use if I convert all of these to fitting functions.
365//
366Function V_TestWavelengthIntegral(folderStr)
367        String folderStr
368
369        SetDataFolder $("root:"+folderStr)
370       
371       
372        // gather the input waves
373        WAVE qVals = $(folderStr+"_q")
374//      WAVE cw = smear_coef_sf
375//      WAVE cw = smear_coef_pgs
376        WAVE cw = smear_coef_BroadPeak
377       
378        // make a wave for the answer
379//      Duplicate/O qvals answer, answer_builtIn
380        Duplicate/O qvals answer_Quad
381       
382        // do the integration
383        // Function IntegrateFn_N(fcn,loLim,upLim,w,x,nord)                             
384
385        Variable loLim,upLim,nord
386       
387        nord = 76       // 20 quadrature points not enough for white beam (especially AgBeh test)
388       
389        loLim = 4/5.3
390        upLim = 8/5.3
391
392// 2018:
393// my quadrature and the built-in function are equivalent. Romberg may be useful in some cases
394// especially for multiple integrals. then number of points and timing can be optimized. But either
395// method can be used. 
396        answer_Quad = IntegrateFn_N(V_WB_testKernel,loLim,upLim,cw,qVals,nord)
397       
398
399// why do I need this? Is this because this is defined as the mean of the distribution
400//  and is needed to normalize the integral? verify this on paper.     
401        answer_Quad *= 5.3
402//      answer_builtIn *= 5.3
403       
404        SetDataFolder root:
405       
406        return 0
407End
408
409
410
411Function V_WB_testKernel(cw,x,dum)
412        Wave cw
413        Variable x              // the q-value for the calculation
414        Variable dum    // the dummy integration variable
415
416        Variable val
417        SVAR funcStr = root:gFunctionString
418        FUNCREF SANSModel_proto func = $funcStr
419               
420//      val = (1-dum*5.3/8)*BroadPeakX(cw,x/dum)
421        val = (1-dum*5.3/8)*func(cw,x/dum)
422       
423//      val = V_WhiteBeamDist(dum*5.3)*BroadPeakX(cw,x/dum)
424//      val = V_WhiteBeamDist(dum*5.3)*func(cw,x/dum)
425
426        return (val)
427End
428
429Proc WBDistr()
430
431        make/O/D distr
432        SetScale/I x 0.755,1.509,"", distr
433        distr = (1-x*5.3/8)
434        display distr
435
436end
437
438// the trick here is that declaring the last qVal wave as a variable
439// since this is implicitly called N times in the wave assignment of the answer wave
440Function V_Integrate_BuiltIn_top(cw,loLim,upLim,qVal)
441        Wave cw
442        Variable loLim,upLim
443        Variable qVal
444       
445        Variable/G root:qq = qval
446        Variable ans
447       
448//      ans = Integrate1D(V_intgrnd_top,lolim,uplim,2,0,cw)             //adaptive quadrature
449        ans = Integrate1D(V_intgrnd_top,lolim,uplim,1,0,cw)             // Romberg integration
450       
451        return ans
452end
453
454// the trick here is that declaring the last qVal wave as a variable
455// since this is implicitly called N times in the wave assignment of the answer wave
456Function V_Integrate_BuiltIn_mid(cw,loLim,upLim,qVal)
457        Wave cw
458        Variable loLim,upLim
459        Variable qVal
460       
461        Variable/G root:qq = qval
462        Variable ans
463       
464//      ans = Integrate1D(V_intgrnd_mid,lolim,uplim,2,0,cw)             //adaptive quadrature
465        ans = Integrate1D(V_intgrnd_mid,lolim,uplim,1,0,cw)             // Romberg integration
466       
467        return ans
468end
469
470// the trick here is that declaring the last qVal wave as a variable
471// since this is implicitly called N times in the wave assignment of the answer wave
472Function V_Integrate_BuiltIn_triangle(cw,loLim,upLim,qVal)
473        Wave cw
474        Variable loLim,upLim
475        Variable qVal
476       
477        Variable/G root:qq = qval
478        Variable ans
479       
480//      ans = Integrate1D(V_intgrnd_triangle,lolim,uplim,2,0,cw)                //adaptive quadrature
481        ans = Integrate1D(V_intgrnd_triangle,lolim,uplim,1,0,cw)                // Romberg integration
482       
483        return ans
484end
485
486// the trick here is that declaring the last qVal wave as a variable
487// since this is implicitly called N times in the wave assignment of the answer wave
488Function V_Integrate_BuiltIn_interp(cw,loLim,upLim,qVal)
489        Wave cw
490        Variable loLim,upLim
491        Variable qVal
492       
493        Variable/G root:qq = qval
494        Variable ans
495       
496//      ans = Integrate1D(V_intgrnd_interp,lolim,uplim,2,0,cw)          //adaptive quadrature
497        ans = Integrate1D(V_intgrnd_interp,lolim,uplim,1,0,cw)          // Romberg integration
498       
499        return ans
500end
501
502//
503// See V_DummyFunctions.ipf for the full list
504//
505//Function BroadPeakX_(cw,x)
506//      Wave cw
507//      Variable x
508//     
509//      return(BroadPeakX(cw,x))
510//end
511
512Function V_intgrnd_top(cw,dum)
513        Wave cw
514        Variable dum            // the dummy of the integration
515
516        Variable val
517        NVAR qq = root:qq               //the q-value of the integration, not part of cw, so pass global
518        SVAR funcStr = root:gFunctionString
519        FUNCREF SANSModel_proto func = $funcStr
520
521//      val = (1-dum*5.3/8)*BroadPeakX(cw,qq/dum)       
522//      val = (1-dum*5.3/8)*func(cw,qq/dum)
523
524//      val = V_WhiteBeamDist(dum*5.3)*BroadPeakX(cw,qq/dum)
525        val = V_WhiteBeamDist_top(dum*5.3)*func(cw,qq/dum)
526       
527//      val = V_WhiteBeamInterp(dum*5.3)*func(cw,qq/dum)
528
529        return (val)
530End
531
532Function V_intgrnd_mid(cw,dum)
533        Wave cw
534        Variable dum            // the dummy of the integration
535
536        Variable val
537        NVAR qq = root:qq               //the q-value of the integration, not part of cw, so pass global
538        SVAR funcStr = root:gFunctionString
539        FUNCREF SANSModel_proto func = $funcStr
540
541//      val = (1-dum*5.3/8)*BroadPeakX(cw,qq/dum)       
542//      val = (1-dum*5.3/8)*func(cw,qq/dum)
543
544//      val = V_WhiteBeamDist(dum*5.3)*BroadPeakX(cw,qq/dum)
545        val = V_WhiteBeamDist_mid(dum*5.3)*func(cw,qq/dum)
546       
547//      val = V_WhiteBeamInterp(dum*5.3)*func(cw,qq/dum)
548
549        return (val)
550End
551
552Function V_intgrnd_triangle(cw,dum)
553        Wave cw
554        Variable dum            // the dummy of the integration
555
556        Variable val
557        NVAR qq = root:qq               //the q-value of the integration, not part of cw, so pass global
558        SVAR funcStr = root:gFunctionString
559        FUNCREF SANSModel_proto func = $funcStr
560
561//      val = (1-dum*5.3/8)*BroadPeakX(cw,qq/dum)       
562        val = (1-dum*5.3/8)*func(cw,qq/dum)
563
564//      val = V_WhiteBeamDist(dum*5.3)*BroadPeakX(cw,qq/dum)
565//      val = V_WhiteBeamDist(dum*5.3)*func(cw,qq/dum)
566       
567//      val = V_WhiteBeamInterp(dum*5.3)*func(cw,qq/dum)
568
569        return (val)
570End
571
572Function V_intgrnd_interp(cw,dum)
573        Wave cw
574        Variable dum            // the dummy of the integration
575
576        Variable val
577        NVAR qq = root:qq               //the q-value of the integration, not part of cw, so pass global
578        SVAR funcStr = root:gFunctionString
579        FUNCREF SANSModel_proto func = $funcStr
580
581//      val = (1-dum*5.3/8)*BroadPeakX(cw,qq/dum)       
582//      val = (1-dum*5.3/8)*func(cw,qq/dum)
583
584//      val = V_WhiteBeamDist(dum*5.3)*BroadPeakX(cw,qq/dum)
585//      val = V_WhiteBeamDist(dum*5.3)*func(cw,qq/dum)
586       
587        val = V_WhiteBeamInterp(dum*5.3)*func(cw,qq/dum)
588
589        return (val)
590End
591
592
593////////////////////////////
594
595// need a function to return the model function name
596// given the coefficient wave
597//
598// want the function NameX for use in the integration, not the AAO function
599//
600
601// from the name of the coefficient wave, get the function name
602// be sure that there is no "Smeared" at the beginning of the name
603// tag X to the end of the name string
604//
605// then the funcString must be passed in as a global to the built-in integration function.
606//
607Function/S V_getXFuncStrFromCoef(cw)
608        Wave cw
609       
610        String cwStr = NameOfWave(cw)
611        String outStr = "",extStr=""
612       
613//      String convStr = ReplaceString("_",cwStr,".")           // change the _ to .   
614//      extStr = ParseFilePath(4, convStr, ":", 0, 0)           // extracts the last .nnn, without the .
615
616// go through the list of coefKWStr pairs
617// look for the cwStr
618// take up to the = (that is the funcStr)
619// remove "Smeared" if needed   
620        SVAR coefList=root:Packages:NIST:coefKWStr
621
622        Variable ii,num
623        String item
624       
625        num=ItemsInList(coefList,";")
626        ii=0
627        do
628                item = StringFromList(ii, coefList, ";")
629               
630                if(strsearch(item,cwStr,0) != -1)               //match
631                        item = ReplaceString("=",item,".")              //replace the = with .
632                        outStr = ParseFilePath(3, item, ":", 0, 0)              // extract file name without extension
633                        outStr = ReplaceString("Smeared",outStr,"")             // replace "Smeared" with null, if it's there
634                        ii = num + 1
635                endif
636               
637                ii+=1
638        while(ii<num)
639       
640        return(outStr+"X")
641end
642
643//////////////////////////////////////////
644// generates dummy functions of the form:
645//
646//Function BroadPeakX_(cw,x)
647//      Wave cw
648//      Variable x
649//      return(BroadPeakX(cw,x))
650//End
651//
652// so that I can use the FUNCREF
653// which fails for some reason when I just use the XOP name?
654//
655//
656// not everything ending in X is a model function - trimmed list is in V_DummyFunctions.ipf
657//
658Function V_generateDummyFuncs()
659
660        String list = FunctionList("*X",";","KIND:4")
661        Variable ii,num
662        String item,str
663       
664        num=ItemsInList(list,";")
665
666        NewNotebook/N=Notebook1/F=0
667       
668       
669        for(ii=0;ii<num;ii+=1)
670                item = StringFromList(ii,list,";")
671                str = "\r"
672                str = "Function "+item+"_(cw,x)\r"
673                str += "\tWave cw\r"
674                str += "\tVariable x\r"
675                str += "\treturn("+item+"(cw,x))\r"
676                str += "End\r\r"
677               
678                //print str
679       
680                Notebook $"", text=str
681               
682        endfor
683        return(0)
684       
685End
Note: See TracBrowser for help on using the repository browser.