source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/Ellipsoid2D_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: 12.5 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4//
5// !!! FOR THE ELLIPSOID, THE ANGLE THETA IS DEFINED FROM ????
6//
7// The plotting macro sets up TWO dependencies
8// - one for the triplet calculation
9// - one for a matrix to display, a copy of the triplet
10//
11// For display, there are two copies of the matrix. One matrix is linear, and is a copy of the
12// triplet (which is ALWAYS linear). The other matrix is toggled log/lin for display
13// in the same way the 2D SANS data matrix is handled.
14//
15
16///  REQUIRES DANSE XOP for 2D FUNCTIONS
17
18//
19// the calculation is done as for the QxQy data set:
20// three waves XYZ, then converted to a matrix
21//
22Proc PlotEllipsoid2D(str)                                               
23        String str
24        Prompt str,"Pick the data folder containing the 2D data",popup,getAList(4)
25
26        if (!exists("Ellipsoid_2DX"))
27                Abort "You must have the SANSAnalysis XOP installed to use 2D models"
28        endif
29       
30        SetDataFolder $("root:"+str)
31       
32        // Setup parameter table for model function
33//      make/O/T/N=12 parameters_Ellip2D
34//      Make/O/D/N=12 coef_Ellip2D
35        make/O/T/N=12 parameters_Ellip2D
36        Make/O/D/N=12 coef_Ellip2D
37       
38        coef_Ellip2D[0] = 1.0
39        coef_Ellip2D[1] = 20.0
40        coef_Ellip2D[2] = 60.0
41        coef_Ellip2D[3] = 1.0e-6
42        coef_Ellip2D[4] = 6.3e-6
43        coef_Ellip2D[5] = 0.0
44        coef_Ellip2D[6] = 1.57
45        coef_Ellip2D[7] = 0.0
46        coef_Ellip2D[8] = 0.0
47        coef_Ellip2D[9] = 0.0
48        coef_Ellip2D[10] = 0.0
49        coef_Ellip2D[11] = 0.0
50        // hard-wire the number of integration points
51//      coef_Ellip2D[12] = 10
52       
53        parameters_Ellip2D[0] = "Scale"
54        parameters_Ellip2D[1] = "Radius_a (rotation axis)"
55        parameters_Ellip2D[2] = "Radius_b"
56        parameters_Ellip2D[3] = "SLD cylinder (A^-2)"
57        parameters_Ellip2D[4] = "SLD solvent"
58        parameters_Ellip2D[5] = "Background"
59        parameters_Ellip2D[6] = "Axis Theta"
60        parameters_Ellip2D[7] = "Axis Phi"
61        parameters_Ellip2D[8] = "Sigma of polydisp in R_a [Angstrom]"
62        parameters_Ellip2D[9] = "Sigma of polydisp in R_b [Angstrom]"
63        parameters_Ellip2D[10] = "Sigma of polydisp in Theta [rad]"
64        parameters_Ellip2D[11] = "Sigma of polydisp in Phi [rad]"
65       
66//      parameters_Ellip2D[12] = "Num of polydisp points"
67
68       
69        Edit parameters_Ellip2D,coef_Ellip2D                                   
70       
71        // generate the triplet representation
72        Duplicate/O $(str+"_qx") xwave_Ellip2D
73        Duplicate/O $(str+"_qy") ywave_Ellip2D,zwave_Ellip2D                   
74               
75        Variable/G g_Ellip2D=0
76        g_Ellip2D := Ellipsoid2D(coef_Ellip2D,zwave_Ellip2D,xwave_Ellip2D,ywave_Ellip2D)        //AAO 2D calculation
77       
78        Display ywave_Ellip2D vs xwave_Ellip2D
79        modifygraph log=0
80        ModifyGraph mode=3,marker=16,zColor(ywave_Ellip2D)={zwave_Ellip2D,*,*,YellowHot,0}
81        ModifyGraph standoff=0
82        ModifyGraph width={Aspect,1}
83        ModifyGraph lowTrip=0.001
84        Label bottom "qx (A\\S-1\\M)"
85        Label left "qy (A\\S-1\\M)"
86        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
87       
88        // generate the matrix representation
89        ConvertQxQy2Mat(xwave_Ellip2D,ywave_Ellip2D,zwave_Ellip2D,"Ellip2D_mat")
90        Duplicate/O $"Ellip2D_mat",$"Ellip2D_lin"               //keep a linear-scaled version of the data
91        // _mat is for display, _lin is the real calculation
92
93        // not a function evaluation - this simply keeps the matrix for display in sync with the triplet calculation
94        Variable/G g_Ellip2Dmat=0
95        g_Ellip2Dmat := UpdateQxQy2Mat(xwave_Ellip2D,ywave_Ellip2D,zwave_Ellip2D,Ellip2D_lin,Ellip2D_mat)
96       
97       
98        SetDataFolder root:
99        AddModelToStrings("Ellipsoid2D","coef_Ellip2D","parameters_Ellip2D","Ellip2D")
100End
101
102
103// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
104Proc PlotSmearedEllipsoid2D(str)                                               
105        String str
106        Prompt str,"Pick the data folder containing the 2D data",popup,getAList(4)
107
108        if (!exists("Ellipsoid_2DX"))
109                Abort "You must have the SANSAnalysis XOP installed to use 2D models"
110        endif
111       
112        SetDataFolder $("root:"+str)
113       
114        // Setup parameter table for model function
115        make/O/T/N=12 smear_parameters_Ellip2D
116        Make/O/D/N=12 smear_coef_Ellip2D
117       
118        smear_coef_Ellip2D[0] = 1.0
119        smear_coef_Ellip2D[1] = 20.0
120        smear_coef_Ellip2D[2] = 60.0
121        smear_coef_Ellip2D[3] = 1.0e-6
122        smear_coef_Ellip2D[4] = 6.3e-6
123        smear_coef_Ellip2D[5] = 0.0
124        smear_coef_Ellip2D[6] = 1.57
125        smear_coef_Ellip2D[7] = 0.0
126        smear_coef_Ellip2D[8] = 0.0
127        smear_coef_Ellip2D[9] = 0.0
128        smear_coef_Ellip2D[10] = 0.0
129        smear_coef_Ellip2D[11] = 0.0
130        // hard-wire the number of integration points
131//      smear_coef_Ellip2D[12] = 10
132       
133        smear_parameters_Ellip2D[0] = "Scale"
134        smear_parameters_Ellip2D[1] = "Radius_a (rotation axis)"
135        smear_parameters_Ellip2D[2] = "Radius_b"
136        smear_parameters_Ellip2D[3] = "SLD cylinder (A^-2)"
137        smear_parameters_Ellip2D[4] = "SLD solvent"
138        smear_parameters_Ellip2D[5] = "Background"
139        smear_parameters_Ellip2D[6] = "Axis Theta"
140        smear_parameters_Ellip2D[7] = "Axis Phi"
141        smear_parameters_Ellip2D[8] = "Sigma of polydisp in R_a [Angstrom]"
142        smear_parameters_Ellip2D[9] = "Sigma of polydisp in R_b [Angstrom]"
143        smear_parameters_Ellip2D[10] = "Sigma of polydisp in Theta [rad]"
144        smear_parameters_Ellip2D[11] = "Sigma of polydisp in Phi [rad]"
145       
146//      smear_parameters_Ellip2D[12] = "Num of polydisp points"
147
148       
149        Edit smear_parameters_Ellip2D,smear_coef_Ellip2D                                       
150       
151        // generate the triplet representation
152        Duplicate/O $(str+"_qx") smeared_Ellip2D
153        SetScale d,0,0,"1/cm",smeared_Ellip2D                                   
154               
155        Variable/G gs_Ellip2D=0
156        gs_Ellip2D := fSmearedEllipsoid2D(smear_coef_Ellip2D,smeared_Ellip2D)           //wrapper to fill the STRUCT
157       
158        Display $(str+"_qy") vs $(str+"_qx")
159        modifygraph log=0
160        ModifyGraph mode=3,marker=16,zColor($(str+"_qy"))={smeared_Ellip2D,*,*,YellowHot,0}
161        ModifyGraph standoff=0
162        ModifyGraph width={Aspect,1}
163        ModifyGraph lowTrip=0.001
164        Label bottom "qx (A\\S-1\\M)"
165        Label left "qy (A\\S-1\\M)"
166        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
167       
168        // generate the matrix representation
169        Duplicate/O $(str+"_qx"), sm_qx
170        Duplicate/O $(str+"_qy"), sm_qy         // I can't use local variables in dependencies, so I need the name (that I can't get)
171       
172        // generate the matrix representation
173        ConvertQxQy2Mat(sm_qx,sm_qy,smeared_Ellip2D,"sm_Ellip2D_mat")
174        Duplicate/O $"sm_Ellip2D_mat",$"sm_Ellip2D_lin"                 //keep a linear-scaled version of the data
175        // _mat is for display, _lin is the real calculation
176
177        // not a function evaluation - this simply keeps the matrix for display in sync with the triplet calculation
178        Variable/G gs_Ellip2Dmat=0
179        gs_Ellip2Dmat := UpdateQxQy2Mat(sm_qx,sm_qy,smeared_Ellip2D,sm_Ellip2D_lin,sm_Ellip2D_mat)
180       
181        SetDataFolder root:
182        AddModelToStrings("SmearedEllipsoid2D","smear_coef_Ellip2D","smear_parameters_Ellip2D","Ellip2D")
183End
184
185
186
187//
188//  Fit function that is actually a wrapper to dispatch the calculation to N threads
189//
190// nthreads is 1 or an even number, typically 2
191// it doesn't matter if npt is odd. In this case, fractional point numbers are passed
192// and the wave indexing works just fine - I tested this with test waves of 7 and 8 points
193// and the points "2.5" and "3.5" evaluate correctly as 2 and 3
194//
195Function Ellipsoid2D(cw,zw,xw,yw) : FitFunc
196        Wave cw,zw,xw,yw
197       
198#if exists("Ellipsoid_2DX")                     //to hide the function if XOP not installed
199
200        Make/O/D/N=13 Ellip2D_tmp
201        Ellip2D_tmp[0,11] = cw
202        Ellip2D_tmp[12] = 25
203        Ellip2D_tmp[5] = 0              //pass in a zero background and add it in later
204       
205        MultiThread zw= Ellipsoid_2DX(Ellip2D_tmp,xw,yw) + cw[5]
206       
207#endif
208       
209        return(0)
210End
211
212
213/////////////////////smeared functions //////////////////////
214
215Function SmearedEllipsoid2D(s)
216        Struct ResSmear_2D_AAOStruct &s
217       
218//// non-threaded, but generic calculation
219//// the last param is nord     
220//      Smear_2DModel_PP(Ellipsoid2D_noThread,s,10)
221
222
223//// the last param is nord
224        SmearedEllipsoid2D_THR(s,10)           
225
226        return(0)
227end
228
229// for the plot dependency only
230Function fSmearedEllipsoid2D(coefW,resultW)
231        Wave coefW,resultW
232       
233        String str = getWavesDataFolder(resultW,0)
234        String DF="root:"+str+":"
235       
236        WAVE qx = $(DF+str+"_qx")
237        WAVE qy = $(DF+str+"_qy")
238        WAVE qz = $(DF+str+"_qz")
239        WAVE sQpl = $(DF+str+"_sQpl")
240        WAVE sQpp = $(DF+str+"_sQpp")
241        WAVE shad = $(DF+str+"_fs")
242       
243        STRUCT ResSmear_2D_AAOStruct s
244        WAVE s.coefW = coefW   
245        WAVE s.zw = resultW     
246        WAVE s.xw[0] = qx
247        WAVE s.xw[1] = qy
248        WAVE s.qz = qz
249        WAVE s.sQpl = sQpl
250        WAVE s.sQpp = sQpp
251        WAVE s.fs = shad
252       
253        Variable err
254        err = SmearedEllipsoid2D(s)
255       
256        return (0)
257End
258
259//
260// NON-THREADED IMPLEMENTATION
261// -- same as threaded, but no MultiThread KW
262//
263ThreadSafe Function Ellipsoid2D_noThread(cw,zw,xw,yw)
264        Wave cw,zw,xw,yw
265       
266        //      Variable t1=StopMSTimer(-2)
267
268#if exists("Ellipsoid_2DX")                     //to hide the function if XOP not installed
269
270        Make/O/D/N=13 Ellip2D_tmp
271        Ellip2D_tmp[0,11] = cw
272        Ellip2D_tmp[12] = 5                     //use a small number of integration points since smearing is used
273        Ellip2D_tmp[5] = 0              //pass in a zero background and add it in later
274       
275        zw= Ellipsoid_2DX(Ellip2D_tmp,xw,yw) + cw[5]
276       
277#endif
278
279//      Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
280       
281        return(0)
282End
283
284//
285// this is the threaded version, that dispatches the calculation out to threads
286//
287// must be written specific to each 2D function
288//
289Function SmearedEllipsoid2D_THR(s,nord)
290        Struct ResSmear_2D_AAOStruct &s
291        Variable nord
292       
293        String weightStr,zStr
294       
295// create all of the necessary quadrature waves here - rather than inside a threadsafe function
296        switch(nord)   
297                case 5:         
298                        weightStr="gauss5wt"
299                        zStr="gauss5z"
300                        if (WaveExists($weightStr) == 0)
301                                Make/O/D/N=(nord) $weightStr,$zStr
302                                Make5GaussPoints($weightStr,$zStr)     
303                        endif
304                        break                           
305                case 10:               
306                        weightStr="gauss10wt"
307                        zStr="gauss10z"
308                        if (WaveExists($weightStr) == 0)
309                                Make/O/D/N=(nord) $weightStr,$zStr
310                                Make10GaussPoints($weightStr,$zStr)     
311                        endif
312                        break                           
313                case 20:               
314                        weightStr="gauss20wt"
315                        zStr="gauss20z"
316                        if (WaveExists($weightStr) == 0)
317                                Make/O/D/N=(nord) $weightStr,$zStr
318                                Make20GaussPoints($weightStr,$zStr)     
319                        endif
320                        break
321                default:                                                       
322                        Abort "Smear_2DModel_PP_Threaded called with invalid nord value"                                       
323        endswitch
324       
325        Wave/Z wt = $weightStr
326        Wave/Z xi = $zStr               // wave references to pass
327
328        Variable npt=numpnts(s.xw[0])
329        Variable i,nthreads= ThreadProcessorCount
330        variable mt= ThreadGroupCreate(nthreads)
331
332        Variable t1=StopMSTimer(-2)
333       
334        for(i=0;i<nthreads;i+=1)
335//              Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
336                ThreadStart mt,i,SmearedEllipsoid2D_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)
337        endfor
338
339        do
340                variable tgs= ThreadGroupWait(mt,100)
341        while( tgs != 0 )
342
343        variable dummy= ThreadGroupRelease(mt)
344       
345// comment out the threading + uncomment this for testing to make sure that the single thread works
346//      nThreads=1
347//      SmearedEllipsoid2D_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)
348
349        Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
350       
351        return(0)
352end
353
354//
355// - worker function for threads of Sphere2D
356//
357ThreadSafe Function SmearedEllipsoid2D_T(coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi,pt1,pt2,nord)
358        WAVE coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi
359        Variable pt1,pt2,nord
360       
361// now passed in....
362//      Wave wt = $weightStr
363//      Wave xi = $zStr         
364
365        Variable ii,jj,kk,num
366        Variable qx,qy,qz,qval,sx,sy,fs
367        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut
368       
369        Variable normFactor,phi,theta,maxSig,numStdDev=3
370       
371/// keep these waves local
372        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW
373       
374        // now just loop over the points as specified
375       
376        answer=0
377       
378        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt
379        Variable qperp_pt,phi_prime,q_prime
380
381        //loop over q-values
382        for(ii=pt1;ii<(pt2+1);ii+=1)
383               
384                qx = qxw[ii]
385                qy = qyw[ii]
386                qz = qzw[ii]
387                qval = sqrt(qx^2+qy^2+qz^2)
388                spl = sxw[ii]
389                spp = syw[ii]
390                fs = fsw[ii]
391               
392                normFactor = 2*pi*spl*spp
393               
394                phi = FindPhi(qx,qy)
395               
396                apl = -numStdDev*spl + qval             //parallel = q integration limits
397                bpl = numStdDev*spl + qval
398                app = -numStdDev*spp + 0                //q_perp = 0
399                bpp = numStdDev*spp + 0
400               
401                //make sure the limits are reasonable.
402                if(apl < 0)
403                        apl = 0
404                endif
405                // do I need to specially handle limits when phi ~ 0?
406       
407               
408                sumOut = 0
409                for(jj=0;jj<nord;jj+=1)         // call phi the "outer'
410                        qperp_pt = (xi[jj]*(bpp-app)+app+bpp)/2         //this is now q_perp
411
412                        sumIn=0
413                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl
414
415                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2
416                               
417                                // find QxQy given Qpl and Qperp on the grid
418                                //
419                                q_prime = sqrt(qpl_pt^2+qperp_pt^2)
420                                phi_prime = phi + qperp_pt/qpl_pt
421                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt)
422                               
423                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not
424                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl
425                               
426                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (qperp_pt)^2/spp/spp ) )
427                                res_tot[kk] /= normFactor
428//                              res_tot[kk] *= fs
429
430                        endfor
431                       
432                        Ellipsoid2D_noThread(coef,fcnRet,xptw,yptw)                     //fcn passed in is an AAO
433                       
434                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0]
435                        fcnRet *= wt[jj]*wt*res_tot
436                        //
437                        answer += (bpl-apl)/2.0*sum(fcnRet)             //
438                endfor
439
440                answer *= (bpp-app)/2.0
441                zw[ii] = answer
442        endfor
443       
444        return(0)
445end
446
447
Note: See TracBrowser for help on using the repository browser.