source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/PeakGauss2D_v40.ipf @ 708

Last change on this file since 708 was 708, checked in by srkline, 13 years ago

SA_Includes_v410 : now include Smear_2D

PeakGauss_2D, Sphere_2D : included threaded resolution smearing calculations for testing

DataSetHandling? : Included a quick and dirty batch converter for XML->6-col. See the top
of the file for the command to run

GaussUtils? : re-define the ResSemear_2D_AAOStruct. Relocated q-value and phi calculations from
RawWindowHook? to this file so they would be available for reduction and analysis

Smear_2D : now has a generic (non-threaded) smearing routine. Threading must be done in
individual functions since FUNCREF can't be passed to threads (plus a few other issues)

PlotUtils_2D : updated loader for new QxQy? columns. Fixes to Wrapper_2D to enable smeared fits

RawWindowHook? : removed Q-value calculation functions and moved these to GaussUtils?

WriteQIS : now writes out 8-columns for QxQy? data, defining the resolution
function in terms of directions parallel and perpendicular to Q. TEMPORARILY in the data
file an error in intensity is generated that is SQRT(I), being careful to
replace any NaN, inf, or zero with an average error value

MultiScatter_MonteCarlo_2D : 4-processor aware

NCNR_Utils : 2D resolution calculation is now in terms of parallel and perpendicular
rather than x and y. Gravity is included in the y-component

File size: 10.4 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 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 PlotPeakGauss2D(str)                                               
21        String str
22        Prompt str,"Pick the data folder containing the 2D data",popup,getAList(4)
23       
24        SetDataFolder $("root:"+str)
25       
26        Make/O/D coef_pkGau2D = {100.0, 0.25,0.005, 1.0}                                               
27        make/o/t parameters_pkGau2D = {"Scale Factor, I0 ", "Peak position (A^-1)", "Std Dev (A^-1)","Incoherent Bgd (cm-1)"}           
28        Edit parameters_pkGau2D,coef_pkGau2D                           
29       
30        // generate the triplet representation
31        Duplicate/O $(str+"_qx") xwave_pkGau2D
32        Duplicate/O $(str+"_qy") ywave_pkGau2D,zwave_pkGau2D                   
33               
34        Variable/G g_pkGau2D=0
35        g_pkGau2D := PeakGauss2D(coef_pkGau2D,zwave_pkGau2D,xwave_pkGau2D,ywave_pkGau2D)        //AAO 2D calculation
36       
37        Display ywave_pkGau2D vs xwave_pkGau2D
38        modifygraph log=0
39        ModifyGraph mode=3,marker=16,zColor(ywave_pkGau2D)={zwave_pkGau2D,*,*,YellowHot,0}
40        ModifyGraph standoff=0
41        ModifyGraph width={Aspect,1}
42        ModifyGraph lowTrip=0.001
43        Label bottom "qx (A\\S-1\\M)"
44        Label left "qy (A\\S-1\\M)"
45        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
46       
47        // generate the matrix representation
48        ConvertQxQy2Mat(xwave_pkGau2D,ywave_pkGau2D,zwave_pkGau2D,"pkGau2D_mat")
49        Duplicate/O $"pkGau2D_mat",$"pkGau2D_lin"               //keep a linear-scaled version of the data
50        // _mat is for display, _lin is the real calculation
51
52        // not a function evaluation - this simply keeps the matrix for display in sync with the triplet calculation
53        Variable/G g_pkGau2Dmat=0
54        g_pkGau2Dmat := UpdateQxQy2Mat(xwave_pkGau2D,ywave_pkGau2D,zwave_pkGau2D,pkGau2D_lin,pkGau2D_mat)
55       
56       
57        SetDataFolder root:
58        AddModelToStrings("PeakGauss2D","coef_pkGau2D","parameters_pkGau2D","pkGau2D")
59End
60
61// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
62Proc PlotSmearedPeakGauss2D(str)                                                               
63        String str
64        Prompt str,"Pick the data folder containing the 2D data",popup,getAList(4)
65       
66        // if any of the resolution waves are missing => abort
67//      if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
68//              Abort
69//      endif
70       
71        SetDataFolder $("root:"+str)
72       
73        // Setup parameter table for model function
74        Make/O/D smear_coef_pkGau2D = {100.0, 0.25,0.005, 1.0}                                 
75        make/o/t smear_parameters_pkGau2D = {"Scale Factor, I0 ", "Peak position (A^-1)", "Std Dev (A^-1)","Incoherent Bgd (cm-1)"}
76        Edit smear_parameters_pkGau2D,smear_coef_pkGau2D                                       
77       
78        Duplicate/O $(str+"_qx") smeared_pkGau2D        //1d place for the smeared model
79        SetScale d,0,0,"1/cm",smeared_pkGau2D                                   
80               
81        Variable/G gs_pkGau2D=0
82        gs_pkGau2D := fSmearedPeakGauss2D(smear_coef_pkGau2D,smeared_pkGau2D)   //this wrapper fills the STRUCT
83
84        Display $(str+"_qy") vs $(str+"_qx")
85        modifygraph log=0
86        ModifyGraph mode=3,marker=16,zColor($(str+"_qy"))={smeared_pkGau2D,*,*,YellowHot,0}
87        ModifyGraph standoff=0
88        ModifyGraph width={Aspect,1}
89        ModifyGraph lowTrip=0.001
90        Label bottom "qx (A\\S-1\\M)"
91        Label left "qy (A\\S-1\\M)"
92        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
93       
94        // generate the matrix representation
95        Duplicate/O $(str+"_qx"), sm_qx
96        Duplicate/O $(str+"_qy"), sm_qy         // I can't use local variables in dependencies, so I need the name (that I can't get)
97       
98        ConvertQxQy2Mat(sm_qx,sm_qy,smeared_pkGau2D,"sm_pkGau2D_mat")
99        Duplicate/O $"sm_pkGau2D_mat",$"sm_pkGau2D_lin"                 //keep a linear-scaled version of the data
100        // _mat is for display, _lin is the real calculation
101
102        // not a function evaluation - this simply keeps the matrix for display in sync with the triplet calculation
103        Variable/G gs_pkGau2Dmat=0
104        gs_pkGau2Dmat := UpdateQxQy2Mat(sm_qx,sm_qy,smeared_pkGau2D,sm_pkGau2D_lin,sm_pkGau2D_mat)
105       
106        SetDataFolder root:
107        AddModelToStrings("SmearedPeakGauss2D","smear_coef_pkGau2D","smear_parameters_pkGau2D","pkGau2D")
108End
109
110
111//threaded version of the function
112ThreadSafe Function PeakGauss2D_T(cw,zw,xw,yw,p1,p2)
113        WAVE cw,zw,xw,yw
114        Variable p1,p2
115       
116#if exists("PeakGauss2DX")                      //to hide the function if XOP not installed
117       
118        zw[p1,p2]= PeakGauss2DX(cw,xw,yw)
119
120#else
121        zw[p1,p2] =  I_PeakGauss2D(cw,xw,yw)
122#endif
123
124        return 0
125End
126
127//
128//  Fit function that is actually a wrapper to dispatch the calculation to N threads
129//
130// not used
131//
132Function PeakGauss2D(cw,zw,xw,yw) : FitFunc
133        Wave cw,zw,xw,yw
134       
135        Variable npt=numpnts(yw)
136        Variable i,nthreads= ThreadProcessorCount
137        variable mt= ThreadGroupCreate(nthreads)
138
139//      Variable t1=StopMSTimer(-2)
140       
141        for(i=0;i<nthreads;i+=1)
142        //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
143                ThreadStart mt,i,PeakGauss2D_T(cw,zw,xw,yw,(i*npt/nthreads),((i+1)*npt/nthreads-1))
144        endfor
145
146        do
147                variable tgs= ThreadGroupWait(mt,100)
148        while( tgs != 0 )
149
150        variable dummy= ThreadGroupRelease(mt)
151       
152//      Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
153       
154        return(0)
155End
156
157
158//non-threaded version of the function
159ThreadSafe Function PeakGauss2D_noThread(cw,zw,xw,yw)
160        WAVE cw,zw, xw,yw
161       
162#if exists("PeakGauss2DX")                      //to hide the function if XOP not installed
163        zw= PeakGauss2DX(cw,xw,yw)
164#else
165        zw = I_PeakGauss2D(cw,xw,yw)
166#endif
167
168        return 0
169End
170
171// I think I did this because when I do the quadrature loops I'm calling the AAO with 1-pt waves, so threading
172// would just be a slowdown
173Function SmearedPeakGauss2D(s)
174        Struct ResSmear_2D_AAOStruct &s
175       
176        //no threads
177//      Smear_2DModel_PP(PeakGauss2D_noThread,s,10)
178
179        //or threaded...
180        SmearPeakGauss2D_THR(s,10)
181       
182        return(0)
183end
184
185
186Function fSmearedPeakGauss2D(coefW,resultW)
187        Wave coefW,resultW
188       
189        String str = getWavesDataFolder(resultW,0)
190        String DF="root:"+str+":"
191       
192        WAVE qx = $(DF+str+"_qx")
193        WAVE qy = $(DF+str+"_qy")
194        WAVE qz = $(DF+str+"_qz")
195        WAVE sQpl = $(DF+str+"_sQpl")
196        WAVE sQpp = $(DF+str+"_sQpp")
197        WAVE shad = $(DF+str+"_fs")
198       
199        STRUCT ResSmear_2D_AAOStruct s
200        WAVE s.coefW = coefW   
201        WAVE s.zw = resultW     
202        WAVE s.xw[0] = qx
203        WAVE s.xw[1] = qy
204        WAVE s.qz = qz
205        WAVE s.sQpl = sQpl
206        WAVE s.sQpp = sQpp
207        WAVE s.fs = shad
208       
209        Variable err
210        err = SmearedPeakGauss2D(s)
211       
212        return (0)
213End
214
215ThreadSafe Function I_PeakGauss2D(w,x,y)
216        Wave w
217        Variable x,y
218       
219        Variable retVal,qval
220        qval = sqrt(x^2+y^2)
221               
222        retval = Peak_Gauss_modelX(w,qval)
223//      retval = fPeak_Gauss_model(w,qval)
224       
225        return(retVal)
226End
227
228///////////////
229//
230// this is the threaded version, that dispatches the calculation out to threads
231//
232// must be written specific to each 2D function
233//
234Function SmearPeakGauss2D_THR(s,nord)
235        Struct ResSmear_2D_AAOStruct &s
236        Variable nord
237       
238        String weightStr,zStr
239       
240// create all of the necessary quadrature waves here - rather than inside a threadsafe function
241        switch(nord)   
242                case 5:         
243                        weightStr="gauss5wt"
244                        zStr="gauss5z"
245                        if (WaveExists($weightStr) == 0)
246                                Make/O/D/N=(nord) $weightStr,$zStr
247                                Make5GaussPoints($weightStr,$zStr)     
248                        endif
249                        break                           
250                case 10:               
251                        weightStr="gauss10wt"
252                        zStr="gauss10z"
253                        if (WaveExists($weightStr) == 0)
254                                Make/O/D/N=(nord) $weightStr,$zStr
255                                Make10GaussPoints($weightStr,$zStr)     
256                        endif
257                        break                           
258                case 20:               
259                        weightStr="gauss20wt"
260                        zStr="gauss20z"
261                        if (WaveExists($weightStr) == 0)
262                                Make/O/D/N=(nord) $weightStr,$zStr
263                                Make20GaussPoints($weightStr,$zStr)     
264                        endif
265                        break
266                default:                                                       
267                        Abort "Smear_2DModel_PP_Threaded called with invalid nord value"                                       
268        endswitch
269       
270        Wave/Z wt = $weightStr
271        Wave/Z xi = $zStr               // wave references to pass
272
273        Variable npt=numpnts(s.xw[0])
274        Variable i,nthreads= ThreadProcessorCount
275        variable mt= ThreadGroupCreate(nthreads)
276
277        Variable t1=StopMSTimer(-2)
278       
279        for(i=0;i<nthreads;i+=1)
280        //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
281                ThreadStart mt,i,SmearPeakGauss2D_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)
282        endfor
283
284        do
285                variable tgs= ThreadGroupWait(mt,100)
286        while( tgs != 0 )
287
288        variable dummy= ThreadGroupRelease(mt)
289       
290// comment out the threading + uncomment this for testing to make sure that the single thread works
291//      nThreads=1
292//      SmearSphere2D_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)
293
294        Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
295       
296        return(0)
297end
298
299//
300// - worker function for threads of PeakGauss2D
301//
302ThreadSafe Function SmearPeakGauss2D_T(coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi,pt1,pt2,nord)
303        WAVE coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi
304        Variable pt1,pt2,nord
305       
306// now passed in....
307//      Wave wt = $weightStr
308//      Wave xi = $zStr         
309
310        Variable ii,jj,kk,num
311        Variable qx,qy,qz,qval,sx,sy,fs
312        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut
313       
314        Variable normFactor,phi,theta,maxSig,numStdDev=3
315       
316/// keep these waves local
317        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW
318       
319        // now just loop over the points as specified
320       
321        answer=0
322       
323        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt
324
325        //loop over q-values
326        for(ii=pt1;ii<(pt2+1);ii+=1)
327               
328                qx = qxw[ii]
329                qy = qyw[ii]
330                qz = qzw[ii]
331                qval = sqrt(qx^2+qy^2+qz^2)
332                spl = sxw[ii]
333                spp = syw[ii]
334                fs = fsw[ii]
335               
336                normFactor = 2*pi*spl*spp
337               
338                phi = FindPhi(qx,qy)
339               
340                apl = -numStdDev*spl + qval             //parallel = q integration limits
341                bpl = numStdDev*spl + qval
342                app = -numStdDev*spp + phi              //perpendicular = phi integration limits
343                bpp = numStdDev*spp + phi
344               
345                //make sure the limits are reasonable.
346                if(apl < 0)
347                        apl = 0
348                endif
349                // do I need to specially handle limits when phi ~ 0?
350       
351               
352                sumOut = 0
353                for(jj=0;jj<nord;jj+=1)         // call phi the "outer'
354                        phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2
355
356                        sumIn=0
357                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl
358
359                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2
360                               
361                                FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi
362                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not
363                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl
364                               
365                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) )
366                                res_tot[kk] /= normFactor
367//                              res_tot[kk] *= fs
368
369                        endfor
370                       
371                        PeakGauss2D_noThread(coef,fcnRet,xptw,yptw)                     //fcn passed in is an AAO
372                       
373                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0]
374                        fcnRet *= wt[jj]*wt*res_tot
375                        //
376                        answer += (bpl-apl)/2.0*sum(fcnRet)             //
377                endfor
378
379                answer *= (bpp-app)/2.0
380                zw[ii] = answer
381        endfor
382       
383        return(0)
384end
385
Note: See TracBrowser for help on using the repository browser.