source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/PeakGauss2D_v40.ipf @ 794

Last change on this file since 794 was 794, checked in by srkline, 12 years ago

Lots of changes:
-2D resolution smearing
-2D error propagation

1) 2D resolution smearing has been corrected to use sigma (perp) correctly
rather than phi. This shows up in the quadrature loop in all of the 2D models
and in the Smear_2D "generic" function.

2) 1D resolutionn smearing is now corrected to account for the integration
range of +/- 3 sigma (99.73% of distribution). Now the result is divided by
0.9973 to rescale it to the correct value.

3) Smeared models are now AAO to improve speed and to allow easier use with
functions that are inherently AAO. No changes are needed, since the call is
behind the scenes, replacing Smear_Model_N() with Smear_Model_N_AAO().

4) in PlotUtils_2D added functions to re-bin the QxQy? data into a 1D format
BinQxQy_to_1D(). This also re-bins the errors in two ways, adding the per-pixel
errors in quadrature, or the deviation from the mean of the intensity. Manually
editing the intensity declaration allows 2D->1D binning of smeared models.

5) Per-pixel propagation of errors has been added through the entire data
reduction sequence. Data errors are generated on loading using Poisson
statistics (specifically tailored for accuracy at low counts) and then is
propagated through each manipulation of the data using standard error
propagation. The error matrix is linear_data_error. As a by-product, all of
the math operations on data are explicitly done on linear_data, to avoid
any potential mistakes of log/linear scaling. Results of this propagation
largely match J. Barker's /ERROR routines from the VAX, with some differences
at low data count values (as expected) and at higher count values near the
beam stop (somewhat unexpected). This per-pixel error is ouput in the QxQy_ASCII
data files. NO CHANGE has been made to the 1D data, which uses the deviation from
the mean as the error - since this is correct.

6) Added tables for the uncertainty in attenuator transmission (from JGB)

7) added two more REAL values to the VAX header to store information
necessary for error propagation. These are couting error that are part of
the transmission error and of the absolute scaling error. These add Read/Write?
functions in NCNR_DataReadWrite

The transmission error (one standard deviation) is written at byte 396 (4 bytes)

RealsRead?[41]

The Box Sum error (one standard deviation) is written at byte 400 (4 bytes)

RealsRead?[42]

File size: 10.7 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4//
5// The plotting macro sets up TWO dependencies
6// - one for the triplet calculation
7// - one for a matrix to display, a copy of the triplet
8//
9// For display, there are two copies of the matrix. One matrix is linear, and is a copy of the
10// triplet (which is ALWAYS linear). The other matrix is toggled log/lin for display
11// in the same way the 2D SANS data matrix is handled.
12//
13
14///  REQUIRES XOP for 2D FUNCTIONS
15
16//
17// the calculation is done as for the QxQy data set:
18// three waves XYZ, then converted to a matrix
19//
20Proc PlotPeakGauss2D(str)                                               
21        String str
22        Prompt str,"Pick the data folder containing the 2D data",popup,getAList(4)
23       
24        SetDataFolder $("root:"+str)
25       
26        Make/O/D coef_pkGau2D = {100.0, 0.25,0.005, 1.0}                                               
27        make/o/t parameters_pkGau2D = {"Scale Factor, I0 ", "Peak position (A^-1)", "Std Dev (A^-1)","Incoherent Bgd (cm-1)"}           
28        Edit parameters_pkGau2D,coef_pkGau2D                           
29       
30        // generate the triplet representation
31        Duplicate/O $(str+"_qx") xwave_pkGau2D
32        Duplicate/O $(str+"_qy") ywave_pkGau2D,zwave_pkGau2D                   
33               
34        Variable/G g_pkGau2D=0
35        g_pkGau2D := PeakGauss2D(coef_pkGau2D,zwave_pkGau2D,xwave_pkGau2D,ywave_pkGau2D)        //AAO 2D calculation
36       
37        Display ywave_pkGau2D vs xwave_pkGau2D
38        modifygraph log=0
39        ModifyGraph mode=3,marker=16,zColor(ywave_pkGau2D)={zwave_pkGau2D,*,*,YellowHot,0}
40        ModifyGraph standoff=0
41        ModifyGraph width={Aspect,1}
42        ModifyGraph lowTrip=0.001
43        Label bottom "qx (A\\S-1\\M)"
44        Label left "qy (A\\S-1\\M)"
45        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
46       
47        // generate the matrix representation
48        ConvertQxQy2Mat(xwave_pkGau2D,ywave_pkGau2D,zwave_pkGau2D,"pkGau2D_mat")
49        Duplicate/O $"pkGau2D_mat",$"pkGau2D_lin"               //keep a linear-scaled version of the data
50        // _mat is for display, _lin is the real calculation
51
52        // not a function evaluation - this simply keeps the matrix for display in sync with the triplet calculation
53        Variable/G g_pkGau2Dmat=0
54        g_pkGau2Dmat := UpdateQxQy2Mat(xwave_pkGau2D,ywave_pkGau2D,zwave_pkGau2D,pkGau2D_lin,pkGau2D_mat)
55       
56       
57        SetDataFolder root:
58        AddModelToStrings("PeakGauss2D","coef_pkGau2D","parameters_pkGau2D","pkGau2D")
59End
60
61// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
62Proc PlotSmearedPeakGauss2D(str)                                                               
63        String str
64        Prompt str,"Pick the data folder containing the 2D data",popup,getAList(4)
65       
66        // if any of the resolution waves are missing => abort
67//      if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
68//              Abort
69//      endif
70       
71        SetDataFolder $("root:"+str)
72       
73        // Setup parameter table for model function
74        Make/O/D smear_coef_pkGau2D = {100.0, 0.25,0.005, 1.0}                                 
75        make/o/t smear_parameters_pkGau2D = {"Scale Factor, I0 ", "Peak position (A^-1)", "Std Dev (A^-1)","Incoherent Bgd (cm-1)"}
76        Edit smear_parameters_pkGau2D,smear_coef_pkGau2D                                       
77       
78        Duplicate/O $(str+"_qx") smeared_pkGau2D        //1d place for the smeared model
79        SetScale d,0,0,"1/cm",smeared_pkGau2D                                   
80               
81        Variable/G gs_pkGau2D=0
82        gs_pkGau2D := fSmearedPeakGauss2D(smear_coef_pkGau2D,smeared_pkGau2D)   //this wrapper fills the STRUCT
83
84        Display $(str+"_qy") vs $(str+"_qx")
85        modifygraph log=0
86        ModifyGraph mode=3,marker=16,zColor($(str+"_qy"))={smeared_pkGau2D,*,*,YellowHot,0}
87        ModifyGraph standoff=0
88        ModifyGraph width={Aspect,1}
89        ModifyGraph lowTrip=0.001
90        Label bottom "qx (A\\S-1\\M)"
91        Label left "qy (A\\S-1\\M)"
92        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
93       
94        // generate the matrix representation
95        Duplicate/O $(str+"_qx"), sm_qx
96        Duplicate/O $(str+"_qy"), sm_qy         // I can't use local variables in dependencies, so I need the name (that I can't get)
97       
98        ConvertQxQy2Mat(sm_qx,sm_qy,smeared_pkGau2D,"sm_pkGau2D_mat")
99        Duplicate/O $"sm_pkGau2D_mat",$"sm_pkGau2D_lin"                 //keep a linear-scaled version of the data
100        // _mat is for display, _lin is the real calculation
101
102        // not a function evaluation - this simply keeps the matrix for display in sync with the triplet calculation
103        Variable/G gs_pkGau2Dmat=0
104        gs_pkGau2Dmat := UpdateQxQy2Mat(sm_qx,sm_qy,smeared_pkGau2D,sm_pkGau2D_lin,sm_pkGau2D_mat)
105       
106        SetDataFolder root:
107        AddModelToStrings("SmearedPeakGauss2D","smear_coef_pkGau2D","smear_parameters_pkGau2D","pkGau2D")
108End
109
110
111//threaded version of the function
112ThreadSafe Function PeakGauss2D_T(cw,zw,xw,yw,p1,p2)
113        WAVE cw,zw,xw,yw
114        Variable p1,p2
115       
116#if exists("PeakGauss2DX")                      //to hide the function if XOP not installed
117       
118        zw[p1,p2]= PeakGauss2DX(cw,xw,yw)
119
120#else
121        zw[p1,p2] =  I_PeakGauss2D(cw,xw,yw)
122#endif
123
124        return 0
125End
126
127//
128//  Fit function that is actually a wrapper to dispatch the calculation to N threads
129//
130// not used
131//
132Function PeakGauss2D(cw,zw,xw,yw) : FitFunc
133        Wave cw,zw,xw,yw
134       
135#if exists("PeakGauss2DX")                      //to hide the function if XOP not installed
136        MultiThread zw= PeakGauss2DX(cw,xw,yw)
137#else
138        MultiThread zw = I_PeakGauss2D(cw,xw,yw)
139#endif 
140       
141//      Variable npt=numpnts(yw)
142//      Variable i,nthreads= ThreadProcessorCount
143//      variable mt= ThreadGroupCreate(nthreads)
144//
145////    Variable t1=StopMSTimer(-2)
146//     
147//      for(i=0;i<nthreads;i+=1)
148//      //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
149//              ThreadStart mt,i,PeakGauss2D_T(cw,zw,xw,yw,(i*npt/nthreads),((i+1)*npt/nthreads-1))
150//      endfor
151//
152//      do
153//              variable tgs= ThreadGroupWait(mt,100)
154//      while( tgs != 0 )
155//
156//      variable dummy= ThreadGroupRelease(mt)
157//     
158////    Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
159       
160        return(0)
161End
162
163
164//non-threaded version of the function
165ThreadSafe Function PeakGauss2D_noThread(cw,zw,xw,yw)
166        WAVE cw,zw, xw,yw
167       
168#if exists("PeakGauss2DX")                      //to hide the function if XOP not installed
169        zw= PeakGauss2DX(cw,xw,yw)
170#else
171        zw = I_PeakGauss2D(cw,xw,yw)
172#endif
173
174        return 0
175End
176
177// I think I did this because when I do the quadrature loops I'm calling the AAO with 1-pt waves, so threading
178// would just be a slowdown
179Function SmearedPeakGauss2D(s)
180        Struct ResSmear_2D_AAOStruct &s
181       
182        //no threads
183//      Smear_2DModel_PP(PeakGauss2D_noThread,s,10)
184
185        //or threaded...
186        SmearPeakGauss2D_THR(s,10)
187       
188        return(0)
189end
190
191
192Function fSmearedPeakGauss2D(coefW,resultW)
193        Wave coefW,resultW
194       
195        String str = getWavesDataFolder(resultW,0)
196        String DF="root:"+str+":"
197       
198        WAVE qx = $(DF+str+"_qx")
199        WAVE qy = $(DF+str+"_qy")
200        WAVE qz = $(DF+str+"_qz")
201        WAVE sQpl = $(DF+str+"_sQpl")
202        WAVE sQpp = $(DF+str+"_sQpp")
203        WAVE shad = $(DF+str+"_fs")
204       
205        STRUCT ResSmear_2D_AAOStruct s
206        WAVE s.coefW = coefW   
207        WAVE s.zw = resultW     
208        WAVE s.xw[0] = qx
209        WAVE s.xw[1] = qy
210        WAVE s.qz = qz
211        WAVE s.sQpl = sQpl
212        WAVE s.sQpp = sQpp
213        WAVE s.fs = shad
214       
215        Variable err
216        err = SmearedPeakGauss2D(s)
217       
218        return (0)
219End
220
221ThreadSafe Function I_PeakGauss2D(w,x,y)
222        Wave w
223        Variable x,y
224       
225        Variable retVal,qval
226        qval = sqrt(x^2+y^2)
227               
228        retval = Peak_Gauss_modelX(w,qval)
229//      retval = fPeak_Gauss_model(w,qval)
230       
231        return(retVal)
232End
233
234///////////////
235//
236// this is the threaded version, that dispatches the calculation out to threads
237//
238// must be written specific to each 2D function
239//
240Function SmearPeakGauss2D_THR(s,nord)
241        Struct ResSmear_2D_AAOStruct &s
242        Variable nord
243       
244        String weightStr,zStr
245       
246// create all of the necessary quadrature waves here - rather than inside a threadsafe function
247        switch(nord)   
248                case 5:         
249                        weightStr="gauss5wt"
250                        zStr="gauss5z"
251                        if (WaveExists($weightStr) == 0)
252                                Make/O/D/N=(nord) $weightStr,$zStr
253                                Make5GaussPoints($weightStr,$zStr)     
254                        endif
255                        break                           
256                case 10:               
257                        weightStr="gauss10wt"
258                        zStr="gauss10z"
259                        if (WaveExists($weightStr) == 0)
260                                Make/O/D/N=(nord) $weightStr,$zStr
261                                Make10GaussPoints($weightStr,$zStr)     
262                        endif
263                        break                           
264                case 20:               
265                        weightStr="gauss20wt"
266                        zStr="gauss20z"
267                        if (WaveExists($weightStr) == 0)
268                                Make/O/D/N=(nord) $weightStr,$zStr
269                                Make20GaussPoints($weightStr,$zStr)     
270                        endif
271                        break
272                default:                                                       
273                        Abort "Smear_2DModel_PP_Threaded called with invalid nord value"                                       
274        endswitch
275       
276        Wave/Z wt = $weightStr
277        Wave/Z xi = $zStr               // wave references to pass
278
279        Variable npt=numpnts(s.xw[0])
280        Variable i,nthreads= ThreadProcessorCount
281        variable mt= ThreadGroupCreate(nthreads)
282
283        Variable t1=StopMSTimer(-2)
284       
285        for(i=0;i<nthreads;i+=1)
286        //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
287                ThreadStart mt,i,SmearPeakGauss2D_T(s.coefW,s.xw[0],s.xw[1],s.qz,s.sQpl,s.sQpp,s.fs,s.zw,wt,xi,(i*npt/nthreads),((i+1)*npt/nthreads-1),nord)
288        endfor
289
290        do
291                variable tgs= ThreadGroupWait(mt,100)
292        while( tgs != 0 )
293
294        variable dummy= ThreadGroupRelease(mt)
295       
296// comment out the threading + uncomment this for testing to make sure that the single thread works
297//      nThreads=1
298//      SmearSphere2D_T(s.coefW,s.xw[0],s.xw[1],s.qz,s.sQpl,s.sQpp,s.fs,s.zw,wt,xi,(i*npt/nthreads),((i+1)*npt/nthreads-1),nord)
299
300        Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
301       
302        return(0)
303end
304
305//
306// - worker function for threads of PeakGauss2D
307//
308ThreadSafe Function SmearPeakGauss2D_T(coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi,pt1,pt2,nord)
309        WAVE coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi
310        Variable pt1,pt2,nord
311       
312// now passed in....
313//      Wave wt = $weightStr
314//      Wave xi = $zStr         
315
316        Variable ii,jj,kk,num
317        Variable qx,qy,qz,qval,sx,sy,fs
318        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut
319       
320        Variable normFactor,phi,theta,maxSig,numStdDev=3
321       
322/// keep these waves local
323        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW
324       
325        // now just loop over the points as specified
326       
327        answer=0
328       
329        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt
330        Variable qperp_pt,phi_prime,q_prime
331
332        //loop over q-values
333        for(ii=pt1;ii<(pt2+1);ii+=1)
334               
335                qx = qxw[ii]
336                qy = qyw[ii]
337                qz = qzw[ii]
338                qval = sqrt(qx^2+qy^2+qz^2)
339                spl = sxw[ii]
340                spp = syw[ii]
341                fs = fsw[ii]
342               
343                normFactor = 2*pi*spl*spp
344               
345                phi = FindPhi(qx,qy)
346               
347                apl = -numStdDev*spl + qval             //parallel = q integration limits
348                bpl = numStdDev*spl + qval
349                app = -numStdDev*spp + 0                //q_perp = 0
350                bpp = numStdDev*spp + 0
351               
352                //make sure the limits are reasonable.
353                if(apl < 0)
354                        apl = 0
355                endif
356                // do I need to specially handle limits when phi ~ 0?
357       
358               
359                sumOut = 0
360                for(jj=0;jj<nord;jj+=1)         // call phi the "outer'
361                        qperp_pt = (xi[jj]*(bpp-app)+app+bpp)/2         //this is now q_perp
362
363                        sumIn=0
364                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl
365
366                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2
367                               
368                                // find QxQy given Qpl and Qperp on the grid
369                                //
370                                q_prime = sqrt(qpl_pt^2+qperp_pt^2)
371                                phi_prime = phi + qperp_pt/qpl_pt
372                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt)
373                               
374                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not
375                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl
376                               
377                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (qperp_pt)^2/spp/spp ) )
378                                res_tot[kk] /= normFactor
379//                              res_tot[kk] *= fs
380
381                        endfor
382                       
383                        PeakGauss2D_noThread(coef,fcnRet,xptw,yptw)                     //fcn passed in is an AAO
384                       
385                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0]
386                        fcnRet *= wt[jj]*wt*res_tot
387                        //
388                        answer += (bpl-apl)/2.0*sum(fcnRet)             //
389                endfor
390
391                answer *= (bpp-app)/2.0
392                zw[ii] = answer
393        endfor
394       
395        return(0)
396end
397
Note: See TracBrowser for help on using the repository browser.