source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/Cylinder_2D_v40.ipf @ 788

Last change on this file since 788 was 788, checked in by srkline, 12 years ago

Added patch to use the correct detector dead time based on whether VAX or ICE hardware was used. The switch is done based on the day of data acquisition using dates supplied by Jeff and Cedric for the date past which the VAX hardware was not used.

25-FEB-2010 for NG7
23-JUL-2009 for NG3

Any data collected after these dates will use the measured dead time with ICE hardware.

File size: 13.2 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
223Function Cylinder2D(cw,zw,xw,yw) : FitFunc
224        Wave cw,zw,xw,yw
225       
226        //      Variable t1=StopMSTimer(-2)
227
228#if exists("Cylinder_2DX")                      //to hide the function if XOP not installed
229
230        Make/O/D/N=12 Cyl2D_tmp                         // there seems to be no speed penalty for doing this...
231        Cyl2D_tmp[0,10] = cw
232        Cyl2D_tmp[11] = 25                                      // hard-wire the number of integration points
233        Cyl2D_tmp[5] = 0                                                // send a background of zero
234       
235        MultiThread zw = Cylinder_2DX(Cyl2D_tmp,xw,yw) + cw[5]          //add in the proper background here
236
237#endif
238
239//      Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
240       
241        return(0)
242End
243
244
245/////////////////////smeared functions //////////////////////
246
247Function SmearedCylinder2D(s)
248        Struct ResSmear_2D_AAOStruct &s
249       
250//// non-threaded, but generic calculation
251//// the last param is nord     
252//      Smear_2DModel_PP(Cylinder2D_noThread,s,10)
253
254
255//// the last param is nord
256        SmearedCylinder2D_THR(s,10)             
257
258        return(0)
259end
260
261// for the plot dependency only
262Function fSmearedCylinder2D(coefW,resultW)
263        Wave coefW,resultW
264       
265        String str = getWavesDataFolder(resultW,0)
266        String DF="root:"+str+":"
267       
268        WAVE qx = $(DF+str+"_qx")
269        WAVE qy = $(DF+str+"_qy")
270        WAVE qz = $(DF+str+"_qz")
271        WAVE sQpl = $(DF+str+"_sQpl")
272        WAVE sQpp = $(DF+str+"_sQpp")
273        WAVE shad = $(DF+str+"_fs")
274       
275        STRUCT ResSmear_2D_AAOStruct s
276        WAVE s.coefW = coefW   
277        WAVE s.zw = resultW     
278        WAVE s.xw[0] = qx
279        WAVE s.xw[1] = qy
280        WAVE s.qz = qz
281        WAVE s.sQpl = sQpl
282        WAVE s.sQpp = sQpp
283        WAVE s.fs = shad
284       
285        Variable err
286        err = SmearedCylinder2D(s)
287       
288        return (0)
289End
290
291//
292// NON-THREADED IMPLEMENTATION
293// -- same as threaded, but no MultiThread KW
294//
295ThreadSafe Function Cylinder2D_noThread(cw,zw,xw,yw)
296        Wave cw,zw,xw,yw
297       
298        //      Variable t1=StopMSTimer(-2)
299
300#if exists("Cylinder_2DX")                      //to hide the function if XOP not installed
301
302        Make/O/D/N=12 Cyl2D_tmp                         // there seems to be no speed penalty for doing this...
303        Cyl2D_tmp[0,10] = cw
304        Cyl2D_tmp[11] = 5                                       // hard-wire the number of integration points
305        Cyl2D_tmp[5] = 0                                                // send a background of zero
306       
307        zw = Cylinder_2DX(Cyl2D_tmp,xw,yw) + cw[5]              //add in the proper background here
308
309#endif
310
311//      Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
312       
313        return(0)
314End
315
316//
317// this is the threaded version, that dispatches the calculation out to threads
318//
319// must be written specific to each 2D function
320//
321Function SmearedCylinder2D_THR(s,nord)
322        Struct ResSmear_2D_AAOStruct &s
323        Variable nord
324       
325        String weightStr,zStr
326       
327// create all of the necessary quadrature waves here - rather than inside a threadsafe function
328        switch(nord)   
329                case 5:         
330                        weightStr="gauss5wt"
331                        zStr="gauss5z"
332                        if (WaveExists($weightStr) == 0)
333                                Make/O/D/N=(nord) $weightStr,$zStr
334                                Make5GaussPoints($weightStr,$zStr)     
335                        endif
336                        break                           
337                case 10:               
338                        weightStr="gauss10wt"
339                        zStr="gauss10z"
340                        if (WaveExists($weightStr) == 0)
341                                Make/O/D/N=(nord) $weightStr,$zStr
342                                Make10GaussPoints($weightStr,$zStr)     
343                        endif
344                        break                           
345                case 20:               
346                        weightStr="gauss20wt"
347                        zStr="gauss20z"
348                        if (WaveExists($weightStr) == 0)
349                                Make/O/D/N=(nord) $weightStr,$zStr
350                                Make20GaussPoints($weightStr,$zStr)     
351                        endif
352                        break
353                default:                                                       
354                        Abort "Smear_2DModel_PP_Threaded called with invalid nord value"                                       
355        endswitch
356       
357        Wave/Z wt = $weightStr
358        Wave/Z xi = $zStr               // wave references to pass
359
360        Variable npt=numpnts(s.xw[0])
361        Variable i,nthreads= ThreadProcessorCount
362        variable mt= ThreadGroupCreate(nthreads)
363
364        Variable t1=StopMSTimer(-2)
365       
366        for(i=0;i<nthreads;i+=1)
367//              Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
368                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)
369        endfor
370
371        do
372                variable tgs= ThreadGroupWait(mt,100)
373        while( tgs != 0 )
374
375        variable dummy= ThreadGroupRelease(mt)
376       
377// comment out the threading + uncomment this for testing to make sure that the single thread works
378//      nThreads=1
379//      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)
380
381        Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
382       
383        return(0)
384end
385
386//
387// - worker function for threads of Sphere2D
388//
389ThreadSafe Function SmearedCylinder2D_T(coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi,pt1,pt2,nord)
390        WAVE coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi
391        Variable pt1,pt2,nord
392       
393// now passed in....
394//      Wave wt = $weightStr
395//      Wave xi = $zStr         
396
397        Variable ii,jj,kk,num
398        Variable qx,qy,qz,qval,sx,sy,fs
399        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut
400       
401        Variable normFactor,phi,theta,maxSig,numStdDev=3
402       
403/// keep these waves local
404        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW
405       
406        // now just loop over the points as specified
407       
408        answer=0
409       
410        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt
411
412        //loop over q-values
413        for(ii=pt1;ii<(pt2+1);ii+=1)
414               
415                qx = qxw[ii]
416                qy = qyw[ii]
417                qz = qzw[ii]
418                qval = sqrt(qx^2+qy^2+qz^2)
419                spl = sxw[ii]
420                spp = syw[ii]
421                fs = fsw[ii]
422               
423                normFactor = 2*pi*spl*spp
424               
425                phi = FindPhi(qx,qy)
426               
427                apl = -numStdDev*spl + qval             //parallel = q integration limits
428                bpl = numStdDev*spl + qval
429                app = -numStdDev*spp + phi              //perpendicular = phi integration limits
430                bpp = numStdDev*spp + phi
431               
432                //make sure the limits are reasonable.
433                if(apl < 0)
434                        apl = 0
435                endif
436                // do I need to specially handle limits when phi ~ 0?
437       
438               
439                sumOut = 0
440                for(jj=0;jj<nord;jj+=1)         // call phi the "outer'
441                        phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2
442
443                        sumIn=0
444                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl
445
446                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2
447                               
448                                FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi
449                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not
450                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl
451                               
452                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) )
453                                res_tot[kk] /= normFactor
454//                              res_tot[kk] *= fs
455
456                        endfor
457                       
458                        Cylinder2D_noThread(coef,fcnRet,xptw,yptw)                      //fcn passed in is an AAO
459                       
460                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0]
461                        fcnRet *= wt[jj]*wt*res_tot
462                        //
463                        answer += (bpl-apl)/2.0*sum(fcnRet)             //
464                endfor
465
466                answer *= (bpp-app)/2.0
467                zw[ii] = answer
468        endfor
469       
470        return(0)
471end
472
473
Note: See TracBrowser for help on using the repository browser.