source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/Cylinder_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: 13.8 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// 4.8x speedup for threading (N=4 + Nvirt=4)
18//
19
20//
21// the calculation is done as for the QxQy data set:
22// three waves XYZ, then converted to a matrix
23//
24Proc PlotCylinder2D(str)                                               
25        String str
26        Prompt str,"Pick the data folder containing the 2D data",popup,getAList(4)
27       
28        if (!exists("Cylinder_2DX") )
29                Abort "You must have the SANSAnalysis XOP installed to use 2D models"
30        endif
31        SetDataFolder $("root:"+str)
32       
33        // NOTE THAT THE COEFFICIENTS [N] ARE IN A DIFFERENT ORDER !!!
34        // Setup parameter table for model function
35        make/O/T/N=11 parameters_Cyl2D
36        Make/O/D/N=11 coef_Cyl2D
37        coef_Cyl2D[0] = 1.0
38        coef_Cyl2D[1] = 20.0
39        coef_Cyl2D[2] = 60.0
40        coef_Cyl2D[3] = 1e-6
41        coef_Cyl2D[4] = 6.3e-6
42        coef_Cyl2D[5] = 0.0
43        coef_Cyl2D[6] = 1.57
44        coef_Cyl2D[7] = 0.0
45        coef_Cyl2D[8] = 0.0
46        coef_Cyl2D[9] = 0.0
47        coef_Cyl2D[10] = 0.0
48       
49        // currently, the number of integration points is hard-wired to be 25 in Cylinder2D_T
50        //coef_Cyl2D[11] = 25
51        //
52        parameters_Cyl2D[0] = "Scale"
53        parameters_Cyl2D[1] = "Radius"
54        parameters_Cyl2D[2] = "Length"
55        parameters_Cyl2D[3] = "SLD cylinder (A^-2)"
56        parameters_Cyl2D[4] = "SLD solvent"
57        parameters_Cyl2D[5] = "Background"
58        parameters_Cyl2D[6] = "Axis Theta"
59        parameters_Cyl2D[7] = "Axis Phi"
60       
61        parameters_Cyl2D[9] = "Sigma of polydisp in Theta [rad]"                //*****
62        parameters_Cyl2D[10] = "Sigma of polydisp in Phi [rad]"                 //*****
63        parameters_Cyl2D[8] = "Sigma of polydisp in Radius [A]"         //*****
64       
65//      parameters_Cyl2D[11] = "number of integration points"
66
67        Edit parameters_Cyl2D,coef_Cyl2D                                       
68       
69        // generate the triplet representation
70        Duplicate/O $(str+"_qx") xwave_Cyl2D
71        Duplicate/O $(str+"_qy") ywave_Cyl2D,zwave_Cyl2D                       
72               
73        Variable/G g_Cyl2D=0
74        g_Cyl2D := Cylinder2D(coef_Cyl2D,zwave_Cyl2D,xwave_Cyl2D,ywave_Cyl2D)   //AAO 2D calculation
75       
76        Display ywave_Cyl2D vs xwave_Cyl2D
77        modifygraph log=0
78        ModifyGraph mode=3,marker=16,zColor(ywave_Cyl2D)={zwave_Cyl2D,*,*,YellowHot,0}
79        ModifyGraph standoff=0
80        ModifyGraph width={Aspect,1}
81        ModifyGraph lowTrip=0.001
82        Label bottom "qx (A\\S-1\\M)"
83        Label left "qy (A\\S-1\\M)"
84        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
85       
86        // generate the matrix representation
87        ConvertQxQy2Mat(xwave_Cyl2D,ywave_Cyl2D,zwave_Cyl2D,"Cyl2D_mat")
88        Duplicate/O $"Cyl2D_mat",$"Cyl2D_lin"           //keep a linear-scaled version of the data
89        // _mat is for display, _lin is the real calculation
90
91        // not a function evaluation - this simply keeps the matrix for display in sync with the triplet calculation
92        Variable/G g_Cyl2Dmat=0
93        g_Cyl2Dmat := UpdateQxQy2Mat(xwave_Cyl2D,ywave_Cyl2D,zwave_Cyl2D,Cyl2D_lin,Cyl2D_mat)
94       
95       
96        SetDataFolder root:
97        AddModelToStrings("Cylinder2D","coef_Cyl2D","parameters_Cyl2D","Cyl2D")
98End
99
100
101// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
102Proc PlotSmearedCylinder2D(str)                                         
103        String str
104        Prompt str,"Pick the data folder containing the 2D data",popup,getAList(4)
105       
106        if (!exists("Cylinder_2DX") )
107                Abort "You must have the SANSAnalysis XOP installed to use 2D models"
108        endif
109        SetDataFolder $("root:"+str)
110       
111        // NOTE THAT THE COEFFICIENTS [N] ARE IN A DIFFERENT ORDER !!!
112        // Setup parameter table for model function
113        make/O/T/N=11 smear_parameters_Cyl2D
114        Make/O/D/N=11 smear_coef_Cyl2D
115        smear_coef_Cyl2D[0] = 1.0
116        smear_coef_Cyl2D[1] = 20.0
117        smear_coef_Cyl2D[2] = 60.0
118        smear_coef_Cyl2D[3] = 1e-6
119        smear_coef_Cyl2D[4] = 6.3e-6
120        smear_coef_Cyl2D[5] = 0.0
121        smear_coef_Cyl2D[6] = 1.57
122        smear_coef_Cyl2D[7] = 0.0
123        smear_coef_Cyl2D[8] = 0.0
124        smear_coef_Cyl2D[9] = 0.0
125        smear_coef_Cyl2D[10] = 0.0
126       
127        // currently, the number of integration points is hard-wired to be 25 in Cylinder2D_T
128        //smear_coef_Cyl2D[11] = 25
129        //
130        smear_parameters_Cyl2D[0] = "Scale"
131        smear_parameters_Cyl2D[1] = "Radius"
132        smear_parameters_Cyl2D[2] = "Length"
133        smear_parameters_Cyl2D[3] = "SLD cylinder (A^-2)"
134        smear_parameters_Cyl2D[4] = "SLD solvent"
135        smear_parameters_Cyl2D[5] = "Background"
136        smear_parameters_Cyl2D[6] = "Axis Theta"
137        smear_parameters_Cyl2D[7] = "Axis Phi"
138       
139        smear_parameters_Cyl2D[9] = "Sigma of polydisp in Theta [rad]"          //*****
140        smear_parameters_Cyl2D[10] = "Sigma of polydisp in Phi [rad]"                   //*****
141        smear_parameters_Cyl2D[8] = "Sigma of polydisp in Radius [A]"           //*****
142       
143//      smear_parameters_Cyl2D[11] = "number of integration points"
144
145        Edit smear_parameters_Cyl2D,smear_coef_Cyl2D                                   
146       
147        // generate the triplet representation
148        Duplicate/O $(str+"_qx") smeared_Cyl2D
149        SetScale d,0,0,"1/cm",smeared_Cyl2D                                     
150               
151        Variable/G gs_Cyl2D=0
152        gs_Cyl2D := fSmearedCylinder2D(smear_coef_Cyl2D,smeared_Cyl2D)          //wrapper to fill the STRUCT
153       
154        Display $(str+"_qy") vs $(str+"_qx")
155        modifygraph log=0
156        ModifyGraph mode=3,marker=16,zColor($(str+"_qy"))={smeared_Cyl2D,*,*,YellowHot,0}
157        ModifyGraph standoff=0
158        ModifyGraph width={Aspect,1}
159        ModifyGraph lowTrip=0.001
160        Label bottom "qx (A\\S-1\\M)"
161        Label left "qy (A\\S-1\\M)"
162        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
163       
164        // generate the matrix representation
165        Duplicate/O $(str+"_qx"), sm_qx
166        Duplicate/O $(str+"_qy"), sm_qy         // I can't use local variables in dependencies, so I need the name (that I can't get)
167       
168        // generate the matrix representation
169        ConvertQxQy2Mat(sm_qx,sm_qy,smeared_Cyl2D,"sm_Cyl2D_mat")
170        Duplicate/O $"sm_Cyl2D_mat",$"sm_Cyl2D_lin"             //keep a linear-scaled version of the data
171        // _mat is for display, _lin is the real calculation
172
173        // not a function evaluation - this simply keeps the matrix for display in sync with the triplet calculation
174        Variable/G gs_Cyl2Dmat=0
175        gs_Cyl2Dmat := UpdateQxQy2Mat(sm_qx,sm_qy,smeared_Cyl2D,sm_Cyl2D_lin,sm_Cyl2D_mat)
176       
177       
178        SetDataFolder root:
179        AddModelToStrings("SmearedCylinder2D","smear_coef_Cyl2D","smear_parameters_Cyl2D","Cyl2D")
180End
181
182
183
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////
195//Function Cylinder2D(cw,zw,xw,yw) : FitFunc
196//      Wave cw,zw,xw,yw
197//     
198//      Variable npt=numpnts(yw)
199//      Variable ii,nthreads= ThreadProcessorCount
200//      variable mt= ThreadGroupCreate(nthreads)
201//
202////    Variable t1=StopMSTimer(-2)
203//     
204//      for(ii=0;ii<nthreads;ii+=1)
205//      //      Print (ii*npt/nthreads),((ii+1)*npt/nthreads-1)
206//              ThreadStart mt,ii,Cylinder2D_T(cw,zw,xw,yw,(ii*npt/nthreads),((ii+1)*npt/nthreads-1))
207//      endfor
208//
209//      do
210//              variable tgs= ThreadGroupWait(mt,100)
211//      while( tgs != 0 )
212//
213//      variable dummy= ThreadGroupRelease(mt)
214//     
215////    Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
216//     
217//      return(0)
218//End
219
220
221/// now using the MultiThread keyword. as of Igor 6.20, the manual threading
222// as above gives a wave read error (index out of range). Same code works fine in Igor 6.12
223//ThreadSafe Function Cylinder2D(cw,zw,xw,yw) : FitFunc
224Function Cylinder2D(cw,zw,xw,yw) : FitFunc
225        Wave cw,zw,xw,yw
226       
227        //      Variable t1=StopMSTimer(-2)
228
229#if exists("Cylinder_2DX")                      //to hide the function if XOP not installed
230
231        Make/O/D/N=12 Cyl2D_tmp                         // there seems to be no speed penalty for doing this...
232        Cyl2D_tmp[0,10] = cw
233        Cyl2D_tmp[11] = 25                                      // hard-wire the number of integration points
234        Cyl2D_tmp[5] = 0                                                // send a background of zero
235       
236        MultiThread zw = Cylinder_2DX(Cyl2D_tmp,xw,yw) + cw[5]          //add in the proper background here
237
238#endif
239
240//      Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
241       
242        return(0)
243End
244
245
246/////////////////////smeared functions //////////////////////
247
248Function SmearedCylinder2D(s)
249        Struct ResSmear_2D_AAOStruct &s
250       
251//// non-threaded, but generic calculation
252//// the last param is nord     
253//      Smear_2DModel_PP(Cylinder2D_noThread,s,10)
254
255// this is generic, but I need to declare the Cylinder2D threadsafe
256// and this calculation is significantly slower than the manually threaded calculation
257// if the function is fast to calculate. Once the function has polydispersity on 2 or more parameters
258// then this AAO calculation and the manual threading are both painfully slow, and more similar in execution time
259//
260// see the function for timing details
261//
262//      Smear_2DModel_PP_AAO(Cylinder2D,s,10)
263
264//// the last param is nord
265        SmearedCylinder2D_THR(s,10)             
266
267        return(0)
268end
269
270// for the plot dependency only
271Function fSmearedCylinder2D(coefW,resultW)
272        Wave coefW,resultW
273       
274        String str = getWavesDataFolder(resultW,0)
275        String DF="root:"+str+":"
276       
277        WAVE qx = $(DF+str+"_qx")
278        WAVE qy = $(DF+str+"_qy")
279        WAVE qz = $(DF+str+"_qz")
280        WAVE sQpl = $(DF+str+"_sQpl")
281        WAVE sQpp = $(DF+str+"_sQpp")
282        WAVE shad = $(DF+str+"_fs")
283       
284        STRUCT ResSmear_2D_AAOStruct s
285        WAVE s.coefW = coefW   
286        WAVE s.zw = resultW     
287        WAVE s.xw[0] = qx
288        WAVE s.xw[1] = qy
289        WAVE s.qz = qz
290        WAVE s.sQpl = sQpl
291        WAVE s.sQpp = sQpp
292        WAVE s.fs = shad
293       
294        Variable err
295        err = SmearedCylinder2D(s)
296       
297        return (0)
298End
299
300//
301// NON-THREADED IMPLEMENTATION
302// -- same as threaded, but no MultiThread KW
303//
304ThreadSafe Function Cylinder2D_noThread(cw,zw,xw,yw)
305        Wave cw,zw,xw,yw
306       
307        //      Variable t1=StopMSTimer(-2)
308
309#if exists("Cylinder_2DX")                      //to hide the function if XOP not installed
310
311        Make/O/D/N=12 Cyl2D_tmp                         // there seems to be no speed penalty for doing this...
312        Cyl2D_tmp[0,10] = cw
313        Cyl2D_tmp[11] = 5                                       // hard-wire the number of integration points
314        Cyl2D_tmp[5] = 0                                                // send a background of zero
315       
316        zw = Cylinder_2DX(Cyl2D_tmp,xw,yw) + cw[5]              //add in the proper background here
317
318#endif
319
320//      Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
321       
322        return(0)
323End
324
325//
326// this is the threaded version, that dispatches the calculation out to threads
327//
328// must be written specific to each 2D function
329//
330Function SmearedCylinder2D_THR(s,nord)
331        Struct ResSmear_2D_AAOStruct &s
332        Variable nord
333       
334        String weightStr,zStr
335       
336// create all of the necessary quadrature waves here - rather than inside a threadsafe function
337        switch(nord)   
338                case 5:         
339                        weightStr="gauss5wt"
340                        zStr="gauss5z"
341                        if (WaveExists($weightStr) == 0)
342                                Make/O/D/N=(nord) $weightStr,$zStr
343                                Make5GaussPoints($weightStr,$zStr)     
344                        endif
345                        break                           
346                case 10:               
347                        weightStr="gauss10wt"
348                        zStr="gauss10z"
349                        if (WaveExists($weightStr) == 0)
350                                Make/O/D/N=(nord) $weightStr,$zStr
351                                Make10GaussPoints($weightStr,$zStr)     
352                        endif
353                        break                           
354                case 20:               
355                        weightStr="gauss20wt"
356                        zStr="gauss20z"
357                        if (WaveExists($weightStr) == 0)
358                                Make/O/D/N=(nord) $weightStr,$zStr
359                                Make20GaussPoints($weightStr,$zStr)     
360                        endif
361                        break
362                default:                                                       
363                        Abort "Smear_2DModel_PP_Threaded called with invalid nord value"                                       
364        endswitch
365       
366        Wave/Z wt = $weightStr
367        Wave/Z xi = $zStr               // wave references to pass
368
369        Variable npt=numpnts(s.xw[0])
370        Variable i,nthreads= ThreadProcessorCount
371        variable mt= ThreadGroupCreate(nthreads)
372
373        Variable t1=StopMSTimer(-2)
374       
375        for(i=0;i<nthreads;i+=1)
376//              Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
377                ThreadStart mt,i,SmearedCylinder2D_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)
378        endfor
379
380        do
381                variable tgs= ThreadGroupWait(mt,100)
382        while( tgs != 0 )
383
384        variable dummy= ThreadGroupRelease(mt)
385       
386// comment out the threading + uncomment this for testing to make sure that the single thread works
387//      nThreads=1
388//      SmearedCylinder2D_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)
389
390        Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
391       
392        return(0)
393end
394
395//
396// - worker function for threads of Cylinder2D
397//
398ThreadSafe Function SmearedCylinder2D_T(coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi,pt1,pt2,nord)
399        WAVE coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi
400        Variable pt1,pt2,nord
401       
402// now passed in....
403//      Wave wt = $weightStr
404//      Wave xi = $zStr         
405
406        Variable ii,jj,kk,num
407        Variable qx,qy,qz,qval,sx,sy,fs
408        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut
409       
410        Variable normFactor,phi,theta,maxSig,numStdDev=3
411       
412/// keep these waves local
413        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW
414       
415        // now just loop over the points as specified
416       
417        answer=0
418       
419        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt
420        Variable qperp_pt,phi_prime,q_prime
421
422        //loop over q-values
423        for(ii=pt1;ii<(pt2+1);ii+=1)
424               
425                qx = qxw[ii]
426                qy = qyw[ii]
427                qz = qzw[ii]
428                qval = sqrt(qx^2+qy^2+qz^2)
429                spl = sxw[ii]
430                spp = syw[ii]
431                fs = fsw[ii]
432               
433                normFactor = 2*pi*spl*spp
434               
435                phi = FindPhi(qx,qy)
436               
437                apl = -numStdDev*spl + qval             //parallel = q integration limits
438                bpl = numStdDev*spl + qval
439                app = -numStdDev*spp + 0                //q_perp = 0
440                bpp = numStdDev*spp + 0
441               
442                //make sure the limits are reasonable.
443                if(apl < 0)
444                        apl = 0
445                endif
446                // do I need to specially handle limits when phi ~ 0?
447       
448               
449                sumOut = 0
450                for(jj=0;jj<nord;jj+=1)         // call phi the "outer'
451                        qperp_pt = (xi[jj]*(bpp-app)+app+bpp)/2         //this is now q_perp
452
453                        sumIn=0
454                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl
455
456                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2
457                               
458                                // find QxQy given Qpl and Qperp on the grid
459                                //
460                                q_prime = sqrt(qpl_pt^2+qperp_pt^2)
461                                phi_prime = phi + qperp_pt/qpl_pt
462                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt)
463                               
464                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not
465                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl
466                               
467                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (qperp_pt)^2/spp/spp ) )
468                                res_tot[kk] /= normFactor
469//                              res_tot[kk] *= fs
470
471                        endfor
472                       
473                        Cylinder2D_noThread(coef,fcnRet,xptw,yptw)                      //fcn passed in is an AAO
474                       
475                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0]
476                        fcnRet *= wt[jj]*wt*res_tot
477                        //
478                        answer += (bpl-apl)/2.0*sum(fcnRet)             //
479                endfor
480
481                answer *= (bpp-app)/2.0
482                zw[ii] = answer
483        endfor
484       
485        return(0)
486end
487
488
Note: See TracBrowser for help on using the repository browser.