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