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