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

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

Changed the 2D resolution calculation to include gravity. This seems to work, but will still require some testing. It requires no modifications to the smearing calculation, still using parallel and perpendicular directions.

Added a Fake Update feature to the RealTime? update. There are specific, separate instructions for how to set this up. Usef (possibly) for summer schools or demos.

Adjusted the 2D MonteCarlo? simulation to a corrected beam center. The Gauss Peak was not symmetric around the old beam center.

Corrected the AAO resolution smearing functions to work with cursors.

File size: 10.7 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#if exists("PeakGauss2DX")                      //to hide the function if XOP not installed
136        MultiThread zw= PeakGauss2DX(cw,xw,yw)
137#else
138        MultiThread zw = I_PeakGauss2D(cw,xw,yw)
139#endif 
140       
141//      Variable npt=numpnts(yw)
142//      Variable i,nthreads= ThreadProcessorCount
143//      variable mt= ThreadGroupCreate(nthreads)
144//
145////    Variable t1=StopMSTimer(-2)
146//     
147//      for(i=0;i<nthreads;i+=1)
148//      //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
149//              ThreadStart mt,i,PeakGauss2D_T(cw,zw,xw,yw,(i*npt/nthreads),((i+1)*npt/nthreads-1))
150//      endfor
151//
152//      do
153//              variable tgs= ThreadGroupWait(mt,100)
154//      while( tgs != 0 )
155//
156//      variable dummy= ThreadGroupRelease(mt)
157//     
158////    Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
159       
160        return(0)
161End
162
163
164//non-threaded version of the function
165ThreadSafe Function PeakGauss2D_noThread(cw,zw,xw,yw)
166        WAVE cw,zw, xw,yw
167       
168#if exists("PeakGauss2DX")                      //to hide the function if XOP not installed
169        zw= PeakGauss2DX(cw,xw,yw)
170#else
171        zw = I_PeakGauss2D(cw,xw,yw)
172#endif
173
174        return 0
175End
176
177// I think I did this because when I do the quadrature loops I'm calling the AAO with 1-pt waves, so threading
178// would just be a slowdown
179Function SmearedPeakGauss2D(s)
180        Struct ResSmear_2D_AAOStruct &s
181       
182        //no threads
183//      Smear_2DModel_PP(PeakGauss2D_noThread,s,10)
184//      Smear_2DModel_PPG(PeakGauss2D_noThread,s,10)
185
186        //or threaded...
187        SmearPeakGauss2D_THR(s,10)
188       
189        return(0)
190end
191
192
193Function fSmearedPeakGauss2D(coefW,resultW)
194        Wave coefW,resultW
195       
196        String str = getWavesDataFolder(resultW,0)
197        String DF="root:"+str+":"
198       
199        WAVE qx = $(DF+str+"_qx")
200        WAVE qy = $(DF+str+"_qy")
201        WAVE qz = $(DF+str+"_qz")
202        WAVE sQpl = $(DF+str+"_sQpl")
203        WAVE sQpp = $(DF+str+"_sQpp")
204        WAVE shad = $(DF+str+"_fs")
205       
206        STRUCT ResSmear_2D_AAOStruct s
207        WAVE s.coefW = coefW   
208        WAVE s.zw = resultW     
209        WAVE s.xw[0] = qx
210        WAVE s.xw[1] = qy
211        WAVE s.qz = qz
212        WAVE s.sQpl = sQpl
213        WAVE s.sQpp = sQpp
214        WAVE s.fs = shad
215       
216        Variable err
217        err = SmearedPeakGauss2D(s)
218       
219        return (0)
220End
221
222ThreadSafe Function I_PeakGauss2D(w,x,y)
223        Wave w
224        Variable x,y
225       
226        Variable retVal,qval
227        qval = sqrt(x^2+y^2)
228               
229        retval = Peak_Gauss_modelX(w,qval)
230//      retval = fPeak_Gauss_model(w,qval)
231       
232        return(retVal)
233End
234
235///////////////
236//
237// this is the threaded version, that dispatches the calculation out to threads
238//
239// must be written specific to each 2D function
240//
241Function SmearPeakGauss2D_THR(s,nord)
242        Struct ResSmear_2D_AAOStruct &s
243        Variable nord
244       
245        String weightStr,zStr
246       
247// create all of the necessary quadrature waves here - rather than inside a threadsafe function
248        switch(nord)   
249                case 5:         
250                        weightStr="gauss5wt"
251                        zStr="gauss5z"
252                        if (WaveExists($weightStr) == 0)
253                                Make/O/D/N=(nord) $weightStr,$zStr
254                                Make5GaussPoints($weightStr,$zStr)     
255                        endif
256                        break                           
257                case 10:               
258                        weightStr="gauss10wt"
259                        zStr="gauss10z"
260                        if (WaveExists($weightStr) == 0)
261                                Make/O/D/N=(nord) $weightStr,$zStr
262                                Make10GaussPoints($weightStr,$zStr)     
263                        endif
264                        break                           
265                case 20:               
266                        weightStr="gauss20wt"
267                        zStr="gauss20z"
268                        if (WaveExists($weightStr) == 0)
269                                Make/O/D/N=(nord) $weightStr,$zStr
270                                Make20GaussPoints($weightStr,$zStr)     
271                        endif
272                        break
273                default:                                                       
274                        Abort "Smear_2DModel_PP_Threaded called with invalid nord value"                                       
275        endswitch
276       
277        Wave/Z wt = $weightStr
278        Wave/Z xi = $zStr               // wave references to pass
279
280        Variable npt=numpnts(s.xw[0])
281        Variable i,nthreads= ThreadProcessorCount
282        variable mt= ThreadGroupCreate(nthreads)
283
284        Variable t1=StopMSTimer(-2)
285       
286        for(i=0;i<nthreads;i+=1)
287        //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
288                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)
289        endfor
290
291        do
292                variable tgs= ThreadGroupWait(mt,100)
293        while( tgs != 0 )
294
295        variable dummy= ThreadGroupRelease(mt)
296       
297// comment out the threading + uncomment this for testing to make sure that the single thread works
298//      nThreads=1
299//      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)
300
301        Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
302       
303        return(0)
304end
305
306//
307// - worker function for threads of PeakGauss2D
308//
309ThreadSafe Function SmearPeakGauss2D_T(coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi,pt1,pt2,nord)
310        WAVE coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi
311        Variable pt1,pt2,nord
312       
313// now passed in....
314//      Wave wt = $weightStr
315//      Wave xi = $zStr         
316
317        Variable ii,jj,kk,num
318        Variable qx,qy,qz,qval,sx,sy,fs
319        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut
320       
321        Variable normFactor,phi,theta,maxSig,numStdDev=3
322       
323/// keep these waves local
324        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW
325       
326        // now just loop over the points as specified
327       
328        answer=0
329       
330        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt
331        Variable qperp_pt,phi_prime,q_prime
332
333        //loop over q-values
334        for(ii=pt1;ii<(pt2+1);ii+=1)
335               
336                qx = qxw[ii]
337                qy = qyw[ii]
338                qz = qzw[ii]
339                qval = sqrt(qx^2+qy^2+qz^2)
340                spl = sxw[ii]
341                spp = syw[ii]
342                fs = fsw[ii]
343               
344                normFactor = 2*pi*spl*spp
345               
346                phi = FindPhi(qx,qy)
347               
348                apl = -numStdDev*spl + qval             //parallel = q integration limits
349                bpl = numStdDev*spl + qval
350                app = -numStdDev*spp + 0                //q_perp = 0
351                bpp = numStdDev*spp + 0
352               
353                //make sure the limits are reasonable.
354                if(apl < 0)
355                        apl = 0
356                endif
357                // do I need to specially handle limits when phi ~ 0?
358       
359               
360                sumOut = 0
361                for(jj=0;jj<nord;jj+=1)         // call phi the "outer'
362                        qperp_pt = (xi[jj]*(bpp-app)+app+bpp)/2         //this is now q_perp
363
364                        sumIn=0
365                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl
366
367                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2
368                               
369                                // find QxQy given Qpl and Qperp on the grid
370                                //
371                                q_prime = sqrt(qpl_pt^2+qperp_pt^2)
372                                phi_prime = phi + qperp_pt/qpl_pt
373                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt)
374                               
375                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not
376                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl
377                               
378                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (qperp_pt)^2/spp/spp ) )
379                                res_tot[kk] /= normFactor
380//                              res_tot[kk] *= fs
381
382                        endfor
383                       
384                        PeakGauss2D_noThread(coef,fcnRet,xptw,yptw)                     //fcn passed in is an AAO
385                       
386                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0]
387                        fcnRet *= wt[jj]*wt*res_tot
388                        //
389                        answer += (bpl-apl)/2.0*sum(fcnRet)             //
390                endfor
391
392                answer *= (bpp-app)/2.0
393                zw[ii] = answer
394        endfor
395       
396        return(0)
397end
398
Note: See TracBrowser for help on using the repository browser.