source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/Sphere_2D_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: 11.9 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 PlotSphere2D(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_sf2D = {1.,60,1e-6,6.3e-6,0.01}                                           
27        make/o/t parameters_sf2D = {"scale","Radius (A)","SLD sphere (A-2)","SLD solvent","bkgd (cm-1)"}               
28        Edit parameters_sf2D,coef_sf2D                         
29       
30        // generate the triplet representation
31        Duplicate/O $(str+"_qx") xwave_sf2D
32        Duplicate/O $(str+"_qy") ywave_sf2D,zwave_sf2D                 
33               
34        Variable/G g_sf2D=0
35        g_sf2D := Sphere2D(coef_sf2D,zwave_sf2D,xwave_sf2D,ywave_sf2D)  //AAO 2D calculation
36       
37        Display ywave_sf2D vs xwave_sf2D
38        modifygraph log=0
39        ModifyGraph mode=3,marker=16,zColor(ywave_sf2D)={zwave_sf2D,*,*,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_sf2D,ywave_sf2D,zwave_sf2D,"sf2D_mat")
49        Duplicate/O $"sf2D_mat",$"sf2D_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_sf2Dmat=0
54        g_sf2Dmat := UpdateQxQy2Mat(xwave_sf2D,ywave_sf2D,zwave_sf2D,sf2D_lin,sf2D_mat)
55       
56       
57        SetDataFolder root:
58        AddModelToStrings("Sphere2D","coef_sf2D","parameters_sf2D","sf2D")
59End
60
61// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
62Proc PlotSmearedSphere2D(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_sf2D = {1.,60,1e-6,6.3e-6,0.01}                                     
75        make/o/t smear_parameters_sf2D = {"scale","Radius (A)","SLD sphere (A-2)","SLD solvent (A-2)","bkgd (cm-1)"}
76        Edit smear_parameters_sf2D,smear_coef_sf2D                                     
77       
78        Duplicate/O $(str+"_qx") smeared_sf2D   //1d place for the smeared model
79        SetScale d,0,0,"1/cm",smeared_sf2D                                     
80               
81        Variable/G gs_sf2D=0
82        gs_sf2D := fSmearedSphere2D(smear_coef_sf2D,smeared_sf2D)       //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_sf2D,*,*,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_sf2D,"sm_sf2D_mat")
99        Duplicate/O $"sm_sf2D_mat",$"sm_sf2D_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_sf2Dmat=0
104        gs_sf2Dmat := UpdateQxQy2Mat(sm_qx,sm_qy,smeared_sf2D,sm_sf2D_lin,sm_sf2D_mat)
105       
106        SetDataFolder root:
107        AddModelToStrings("SmearedSphere2D","smear_coef_sf2D","smear_parameters_sf2D","sf2D")
108End
109
110//
111//  Fit function that is actually a wrapper to dispatch the calculation to N threads
112//
113// nthreads is 1 or an even number, typically 2
114// it doesn't matter if npt is odd. In this case, fractional point numbers are passed
115// and the wave indexing works just fine - I tested this with test waves of 7 and 8 points
116// and the points "2.5" and "3.5" evaluate correctly as 2 and 3
117//
118Function Sphere2D(cw,zw,xw,yw) : FitFunc
119        Wave cw,zw,xw,yw
120
121#if exists("Sphere_2DX")                        //to hide the function if XOP not installed
122        MultiThread     zw= Sphere_2DX(cw,xw,yw)
123#endif
124
125        return(0)
126End
127
128
129///
130//// keep this section as an example
131//
132
133//      Variable npt=numpnts(yw)
134//      Variable i,nthreads= ThreadProcessorCount
135//      variable mt= ThreadGroupCreate(nthreads)
136//
137////    Variable t1=StopMSTimer(-2)
138//     
139//      for(i=0;i<nthreads;i+=1)
140//      //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
141//              ThreadStart mt,i,Sphere2D_T(cw,zw,xw,yw,(i*npt/nthreads),((i+1)*npt/nthreads-1))
142//      endfor
143//
144//      do
145//              variable tgs= ThreadGroupWait(mt,100)
146//      while( tgs != 0 )
147//
148//      variable dummy= ThreadGroupRelease(mt)
149//     
150////    Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
151//     
152//
153////// end example of threading
154
155
156
157//threaded version of the function
158//ThreadSafe Function Sphere2D_T(cw,zw,xw,yw,p1,p2)
159//      WAVE cw,zw,xw,yw
160//      Variable p1,p2
161//     
162//#if exists("Sphere_2DX")                      //to hide the function if XOP not installed
163//     
164//      zw[p1,p2]= Sphere_2DX(cw,xw,yw)
165//
166//#endif
167//
168//      return 0
169//End
170
171
172//non-threaded version of the function, necessary for the smearing calculation
173// -- the smearing calculation can only calculate (nord) points at a time.
174//
175ThreadSafe Function Sphere2D_noThread(cw,zw,xw,yw)
176        WAVE cw,zw, xw,yw
177       
178#if exists("Sphere_2DX")                        //to hide the function if XOP not installed
179        zw= Sphere_2DX(cw,xw,yw)
180#endif
181
182        return 0
183End
184
185
186
187//// the threaded version must be specifically written, since
188//// FUNCREF can't be passed into  a threaded calc (structures can't be passed either)
189// so in this implementation, the smearing is dispatched as threads to a function that
190// can calculate the function for a range of points in the input qxqyqz. It is important
191// that the worker calls the un-threaded model function (so write one) and that in the (nord x nord)
192// loop, vectors of length (nord) are calculated rather than pointwise, since the model
193// function is AAO.
194// -- makes things rather messy to code individual functions, but I really see no other way
195// given the restrictions of what can be passed to threaded functions.
196//
197//
198// The smearing is handled this way since 1D smearing is 20 x 200 pts = 4000 evaluations
199// and the 2D is (10 x 10) x 16000 pts = 1,600,000 evaluations (if it's done like the 1D, it's 4000x slower)
200//
201//
202// - the threading gives a clean speedup of 2 for N=2, even for this simple calculation
203//  -- 4.8X speedup for N=8 (4 real cores + 4 virtual cores)
204//
205// nord = 5,10,20 allowed
206//
207Function SmearedSphere2D(s)
208        Struct ResSmear_2D_AAOStruct &s
209       
210//// non-threaded, but generic calculation
211//// the last param is nord     
212//      Smear_2DModel_PP(Sphere2D_noThread,s,10)
213
214
215//// the last param is nord
216        SmearedSphere2D_THR(s,10)               
217
218        return(0)
219end
220
221
222Function fSmearedSphere2D(coefW,resultW)
223        Wave coefW,resultW
224       
225        String str = getWavesDataFolder(resultW,0)
226        String DF="root:"+str+":"
227       
228        WAVE qx = $(DF+str+"_qx")
229        WAVE qy = $(DF+str+"_qy")
230        WAVE qz = $(DF+str+"_qz")
231        WAVE sQpl = $(DF+str+"_sQpl")
232        WAVE sQpp = $(DF+str+"_sQpp")
233        WAVE shad = $(DF+str+"_fs")
234       
235        STRUCT ResSmear_2D_AAOStruct s
236        WAVE s.coefW = coefW   
237        WAVE s.zw = resultW     
238        WAVE s.xw[0] = qx
239        WAVE s.xw[1] = qy
240        WAVE s.qz = qz
241        WAVE s.sQpl = sQpl
242        WAVE s.sQpp = sQpp
243        WAVE s.fs = shad
244       
245        Variable err
246        err = SmearedSphere2D(s)
247       
248        return (0)
249End
250
251
252
253
254//
255// this is the threaded version, that dispatches the calculation out to threads
256//
257// must be written specific to each 2D function
258//
259Function SmearedSphere2D_THR(s,nord)
260        Struct ResSmear_2D_AAOStruct &s
261        Variable nord
262       
263        String weightStr,zStr
264       
265// create all of the necessary quadrature waves here - rather than inside a threadsafe function
266        switch(nord)   
267                case 5:         
268                        weightStr="gauss5wt"
269                        zStr="gauss5z"
270                        if (WaveExists($weightStr) == 0)
271                                Make/O/D/N=(nord) $weightStr,$zStr
272                                Make5GaussPoints($weightStr,$zStr)     
273                        endif
274                        break                           
275                case 10:               
276                        weightStr="gauss10wt"
277                        zStr="gauss10z"
278                        if (WaveExists($weightStr) == 0)
279                                Make/O/D/N=(nord) $weightStr,$zStr
280                                Make10GaussPoints($weightStr,$zStr)     
281                        endif
282                        break                           
283                case 20:               
284                        weightStr="gauss20wt"
285                        zStr="gauss20z"
286                        if (WaveExists($weightStr) == 0)
287                                Make/O/D/N=(nord) $weightStr,$zStr
288                                Make20GaussPoints($weightStr,$zStr)     
289                        endif
290                        break
291                default:                                                       
292                        Abort "Smear_2DModel_PP_Threaded called with invalid nord value"                                       
293        endswitch
294       
295        Wave/Z wt = $weightStr
296        Wave/Z xi = $zStr               // wave references to pass
297
298        Variable npt=numpnts(s.xw[0])
299        Variable i,nthreads= ThreadProcessorCount
300        variable mt= ThreadGroupCreate(nthreads)
301
302        Variable t1=StopMSTimer(-2)
303       
304        for(i=0;i<nthreads;i+=1)
305//              Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
306                ThreadStart mt,i,SmearedSphere2D_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)
307        endfor
308
309        do
310                variable tgs= ThreadGroupWait(mt,100)
311        while( tgs != 0 )
312
313        variable dummy= ThreadGroupRelease(mt)
314       
315// comment out the threading + uncomment this for testing to make sure that the single thread works
316//      nThreads=1
317//      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)
318
319        Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
320       
321        return(0)
322end
323
324//
325// - worker function for threads of Sphere2D
326//
327ThreadSafe Function SmearedSphere2D_T(coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi,pt1,pt2,nord)
328        WAVE coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi
329        Variable pt1,pt2,nord
330       
331// now passed in....
332//      Wave wt = $weightStr
333//      Wave xi = $zStr         
334
335        Variable ii,jj,kk,num
336        Variable qx,qy,qz,qval,sx,sy,fs
337        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut
338       
339        Variable normFactor,phi,theta,maxSig,numStdDev=3
340       
341/// keep these waves local
342        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW
343       
344        // now just loop over the points as specified
345       
346        answer=0
347       
348        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt
349        Variable qperp_pt,phi_prime,q_prime
350
351        //loop over q-values
352        for(ii=pt1;ii<(pt2+1);ii+=1)
353               
354                qx = qxw[ii]
355                qy = qyw[ii]
356                qz = qzw[ii]
357                qval = sqrt(qx^2+qy^2+qz^2)
358                spl = sxw[ii]
359                spp = syw[ii]
360                fs = fsw[ii]
361
362               
363                normFactor = 2*pi*spl*spp
364               
365                phi = FindPhi(qx,qy)
366               
367                apl = -numStdDev*spl + qval             //parallel = q integration limits
368                bpl = numStdDev*spl + qval
369///             app = -numStdDev*spp + phi              //perpendicular = phi integration limits (WRONG)
370///             bpp = numStdDev*spp + phi
371                app = -numStdDev*spp + 0                //q_perp = 0
372                bpp = numStdDev*spp + 0
373               
374                //make sure the limits are reasonable.
375                if(apl < 0)
376                        apl = 0
377                endif
378                // do I need to specially handle limits when phi ~ 0?
379       
380               
381                sumOut = 0
382                for(jj=0;jj<nord;jj+=1)         // call phi the "outer'
383///                     phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2
384                        qperp_pt = (xi[jj]*(bpp-app)+app+bpp)/2         //this is now q_perp
385                       
386                        sumIn=0
387                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl
388
389                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2
390                               
391///                             FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi
392
393                                // find QxQy given Qpl and Qperp on the grid
394                                //
395                                q_prime = sqrt(qpl_pt^2+qperp_pt^2)
396                                phi_prime = phi + qperp_pt/q_prime
397                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt)
398                               
399                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not
400                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl
401
402                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (qperp_pt)^2/spp/spp ) )
403///                             res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) )
404                                res_tot[kk] /= normFactor
405//                              res_tot[kk] *= fs
406
407                        endfor
408                       
409                        Sphere2D_noThread(coef,fcnRet,xptw,yptw)                        //fcn passed in is an AAO
410                       
411                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0]
412                        fcnRet *= wt[jj]*wt*res_tot
413                        //
414                        answer += (bpl-apl)/2.0*sum(fcnRet)             //
415                endfor
416
417                answer *= (bpp-app)/2.0
418                zw[ii] = answer
419        endfor
420       
421        return(0)
422end
423
424
Note: See TracBrowser for help on using the repository browser.