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

Last change on this file since 1239 was 1239, checked in by srkline, 3 years ago

adding downstream window transmission correction -- added function to detector corrections, and added item to preference panel. Transmission value is currently set as a global since the value is not part of the VSANS header. Global value defaults to T=1= no correction.

Other changes are cleanup of TODO items that were already done and tested. inline commnents have been updated.

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