source: sans/Dev/trunk/NCNR_User_Procedures/Common/Smear_2D.ipf @ 708

Last change on this file since 708 was 708, checked in by srkline, 12 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: 11.5 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4
5
6// there are utlity functions here (commented out) for plotting the 2D resolution function
7// at any point on the detector
8
9
10// I can only calculate (nord) points in an AAO fashion, so threading the function calculation
11// doesn't help. Even if I can rewrite this to calculate nord*nord AAO, that will typically be 10*10=100
12// and that is not enough points to thread with any benefit
13
14// but the calculation of each q value is independent, so I can split the total number of q-values between processors
15//
16// -- but I must be careful to pass this a function that is not already threaded!
17// --- BUT - I can't pass either a function reference, OR a structure to a thread!
18// ---- So sadly, each 2D resolution calculation must be threaded by hand.
19//
20// See the Sphere_2D function for an example of how to thread the smearing
21//
22//
23// SRK May 2010
24
25//
26// NOTE: there is a definition of FindTheta = -1* FindPhi() that is a duplicate of what is in RawWindowHook.ipf
27// and should eventually be in a common location for analysis and reduction packages
28//
29
30
31//// this is the completely generic 2D smearing, not threaded, but takes FUNCREF and STRUCT parameters
32//
33// uses resolution ellipse defined perpendicular "Y" and parallel "X",
34// rotating the ellipse into its proper orientaiton based on qx,qy
35//
36// 5 gauss points is not enough - it gives artifacts as a function of phi
37// 10 gauss points is minimally sufficient
38// 20 gauss points are needed if lots of oscillations (just like in 1D)
39// even more may be necessary for highly peaked functions
40//
41//
42Function Smear_2DModel_PP(fcn,s,nord)
43        FUNCREF SANS_2D_ModelAAO_proto fcn
44        Struct ResSmear_2D_AAOStruct &s
45        Variable nord
46       
47        String weightStr,zStr
48
49        Variable ii,jj,kk,num
50        Variable qx,qy,qz,qval,fs
51        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut
52       
53        Variable a,b,c,normFactor,phi,theta,maxSig,numStdDev=3
54       
55        switch(nord)   
56                case 5:         
57                        weightStr="gauss5wt"
58                        zStr="gauss5z"
59                        if (WaveExists($weightStr) == 0)
60                                Make/O/D/N=(nord) $weightStr,$zStr
61                                Make5GaussPoints($weightStr,$zStr)     
62                        endif
63                        break                           
64                case 10:               
65                        weightStr="gauss10wt"
66                        zStr="gauss10z"
67                        if (WaveExists($weightStr) == 0)
68                                Make/O/D/N=(nord) $weightStr,$zStr
69                                Make10GaussPoints($weightStr,$zStr)     
70                        endif
71                        break                           
72                case 20:               
73                        weightStr="gauss20wt"
74                        zStr="gauss20z"
75                        if (WaveExists($weightStr) == 0)
76                                Make/O/D/N=(nord) $weightStr,$zStr
77                                Make20GaussPoints($weightStr,$zStr)     
78                        endif
79                        break
80                default:                                                       
81                        Abort "Smear_2DModel_PP called with invalid nord value"                                 
82        endswitch
83       
84        Wave/Z wt = $weightStr
85        Wave/Z xi = $zStr
86       
87/// keep these waves local
88//      Make/O/D/N=1 yPtW
89        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW
90       
91        // now just loop over the points as specified
92        num=numpnts(s.xw[0])
93       
94        answer=0
95       
96        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt
97        Variable t1=StopMSTimer(-2)
98
99        //loop over q-values
100        for(ii=0;ii<num;ii+=1)
101       
102//              if(mod(ii, 1000 ) == 0)
103//                      Print "ii= ",ii
104//              endif
105               
106                qx = s.xw[0][ii]
107                qy = s.xw[1][ii]
108                qz = s.qz[ii]
109                qval = sqrt(qx^2+qy^2+qz^2)
110                spl = s.sQpl[ii]
111                spp = s.sQpp[ii]
112                fs = s.fs[ii]
113               
114                normFactor = 2*pi*spl*spp
115               
116                phi = -1*FindTheta(qx,qy)               //Findtheta is an exact duplicate of FindPhi() * -1
117               
118                apl = -numStdDev*spl + qval             //parallel = q integration limits
119                bpl = numStdDev*spl + qval
120                app = -numStdDev*spp + phi              //perpendicular = phi integration limits
121                bpp = numStdDev*spp + phi
122               
123                //make sure the limits are reasonable.
124                if(apl < 0)
125                        apl = 0
126                endif
127                // do I need to specially handle limits when phi ~ 0?
128       
129               
130                sumOut = 0
131                for(jj=0;jj<nord;jj+=1)         // call phi the "outer'
132                        phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2
133
134                        sumIn=0
135                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl
136
137                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2
138                               
139                                FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi
140                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not
141                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl
142                               
143                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) )
144                                res_tot[kk] /= normFactor
145//                              res_tot[kk] *= fs
146
147                        endfor
148                       
149                        fcn(s.coefW,fcnRet,xptw,yptw)                   //calculate nord pts at a time
150                       
151                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0]
152                        fcnRet *= wt[jj]*wt*res_tot
153                        //
154                        answer += (bpl-apl)/2.0*sum(fcnRet)             //
155                endfor
156
157                answer *= (bpp-app)/2.0
158                s.zw[ii] = answer
159        endfor
160       
161        Variable elap = (StopMSTimer(-2) - t1)/1e6
162        Print "elapsed time = ",elap
163       
164        return(0)
165end
166
167
168
169
170
171//phi is defined from +x axis, proceeding CCW around [0,2Pi]
172//rotate the resolution function by theta,  = -phi
173//
174// this is only different by (-1) from FindPhi
175// I'd just call FindPhi, but it's awkward to include
176//
177Threadsafe Function FindTheta(vx,vy)
178        variable vx,vy
179
180       
181        return(-1 * FindPhi(vx,vy))
182end
183
184
185
186       
187//////////////////////////////////////////////////////////
188//////////////////////////////////////////////////////////
189/// utility functions to test the calculation of the resolution function and
190//       plotting of it for visualization
191//
192// -- these need to have the reduction package included for it to compile
193//
194// this works with the raw 2D data, calculating resolution in place
195// type is "RAW"
196// xx,yy are in detector coordinates
197//
198// call as:
199// PlotResolution_atPixel("RAW",pcsr(A),qcsr(A))
200//
201// then:
202// Display;AppendImage res;AppendMatrixContour res;ModifyContour res labels=0,autoLevels={*,*,3}
203//
204
205//Function PlotResolution_atPixel(type,xx,yy)
206//      String type
207//      Variable xx,yy
208//     
209/////from QxQyExport
210//      String destStr="",typeStr=""
211//      Variable step=1,refnum
212//      destStr = "root:Packages:NIST:"+type
213//     
214//      //must select the linear_data to export
215//      NVAR isLog = $(destStr+":gIsLogScale")
216//      if(isLog==1)
217//              typeStr = ":linear_data"
218//      else
219//              typeStr = ":data"
220//      endif
221//     
222//      NVAR pixelsX = root:myGlobals:gNPixelsX
223//      NVAR pixelsY = root:myGlobals:gNPixelsY
224//     
225//      Wave data=$(destStr+typeStr)
226//      WAVE intw=$(destStr + ":integersRead")
227//      WAVE rw=$(destStr + ":realsRead")
228//      WAVE/T textw=$(destStr + ":textRead")
229//     
230////    Duplicate/O data,qx_val,qy_val,z_val,qval,qz_val,phi,r_dist
231//      Variable qx_val,qy_val,z_val,qval,qz_val,phi,r_dist
232//
233//      Variable xctr,yctr,sdd,lambda,pixSize
234//      xctr = rw[16]
235//      yctr = rw[17]
236//      sdd = rw[18]
237//      lambda = rw[26]
238//      pixSize = rw[13]/10             //convert mm to cm (x and y are the same size pixels)
239//     
240////    qx_val = CalcQx(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)          //+1 converts to detector coordinate system
241////    qy_val = CalcQy(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)
242//     
243//      qx_val = CalcQx(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)                //+1 converts to detector coordinate system
244//      qy_val = CalcQy(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)
245//     
246////    Redimension/N=(pixelsX*pixelsY) qx_val,qy_val,z_val
247//
248/////************
249//// do everything to write out the resolution too
250//      // un-comment these if you want to write out qz_val and qval too, then use the proper save command
251//      qval = CalcQval(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)
252//      qz_val = CalcQz(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)
253//      phi = FindPhi( pixSize*((xx+1)-xctr) , pixSize*((yy+1)-yctr))           //(dx,dy)
254//      r_dist = sqrt(  (pixSize*((xx+1)-xctr))^2 +  (pixSize*((yy+1)-yctr))^2 )                //radial distance from ctr to pt
255////    Redimension/N=(pixelsX*pixelsY) qz_val,qval,phi,r_dist
256//      //everything in 1D now
257////    Duplicate/O qval SigmaQX,SigmaQY,fsubS
258//      Variable SigmaQX,SigmaQY,fsubS
259//
260//      Variable L2 = rw[18]
261//      Variable BS = rw[21]
262//      Variable S1 = rw[23]
263//      Variable S2 = rw[24]
264//      Variable L1 = rw[25]
265//      Variable lambdaWidth = rw[27]   
266//      Variable usingLenses = rw[28]           //new 2007
267//
268//      //Two parameters DDET and APOFF are instrument dependent.  Determine
269//      //these from the instrument name in the header.
270//      //From conversation with JB on 01.06.99 these are the current good values
271//      Variable DDet
272//      NVAR apOff = root:myGlobals:apOff               //in cm
273//      DDet = rw[10]/10                        // header value (X) is in mm, want cm here
274//
275//      Variable ret1,ret2,ret3,del_r
276//      del_r = rw[10]
277//     
278//      get2DResolution(qval,phi,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,r_dist,ret1,ret2,ret3)
279//      SigmaQX = ret1 
280//      SigmaQY = ret2 
281//      fsubs = ret3   
282///////
283//
284////    Variable theta,phi,qx,qy,sx,sy,a,b,c,val,ii,jj,num=15,x0,y0
285//      Variable theta,qx,qy,sx,sy,a,b,c,val,ii,jj,num=15,x0,y0,maxSig,nStdDev=3,normFactor
286//      Variable qx_ret,qy_ret
287//     
288////    theta = FindPhi(qx_val,qy_val)
289//// need to rotate properly - theta is defined a starting from +y axis, moving CW
290//// we define phi starting from +x and moving CCW
291//      theta = -phi                    //seems to give the right behavior...
292//     
293//     
294//      Print qx_val,qy_val,qval
295//      Print "phi, theta",phi,theta
296//     
297//      FindQxQy(qval,phi,qx_ret,qy_ret)
298//     
299//      sx = SigmaQx
300//      sy = sigmaQy
301//      x0 = qx_val
302//      y0 = qy_val
303//     
304//      a = cos(theta)^2/(2*sx*sx) + sin(theta)^2/(2*sy*sy)
305//      b = -1*sin(2*theta)/(4*sx*sx) + sin(2*theta)/(4*sy*sy)
306//      c = sin(theta)^2/(2*sx*sx) + cos(theta)^2/(2*sy*sy)
307//     
308//      normFactor = pi/sqrt(a*c-b*b)
309//
310//      Make/O/D/N=(num,num) res
311//      // so the resolution function 'looks right' on a 2D plot - otherwise it always looks like a circle
312//      maxSig = max(sx,sy)
313//      Setscale/I x -nStdDev*maxSig+x0,nStdDev*maxSig+x0,res
314//      Setscale/I y -nStdDev*maxSig+y0,nStdDev*maxSig+y0,res
315////    Setscale/I x -nStdDev*sx+x0,nStdDev*sx+x0,res
316////    Setscale/I y -nStdDev*sy+y0,nStdDev*sy+y0,res
317//     
318//      Variable xPt,yPt,delx,dely,offx,offy
319//      delx = DimDelta(res,0)
320//      dely = DimDelta(res,1)
321//      offx = DimOffset(res,0)
322//      offy = DimOffset(res,1)
323//
324//      Print "sx,sy = ",sx,sy
325//      for(ii=0;ii<num;ii+=1)
326//              xPt = offx + ii*delx
327//              for(jj=0;jj<num;jj+=1)
328//                      yPt = offy + jj*dely
329//                      res[ii][jj] = exp(-1*(a*(xPt-x0)^2 + 2*b*(xPt-x0)*(yPt-y0) + c*(yPt-y0)^2))
330//              endfor
331//      endfor 
332//      res /= normFactor
333//     
334//      //Print sum(res,-inf,inf)*delx*dely
335//      if(WaveExists($"coef")==0)
336//              Make/O/D/N=6 coef
337//      endif
338//      Wave coef=coef
339//      coef[0] = 1
340//      coef[1] = qx_val
341//      coef[2] = qy_val
342//      coef[3] = sx
343//      coef[4] = sy
344//      coef[5] = theta
345//
346////    Variable t1=StopMSTimer(-2)
347//
348////
349////    do2dIntegrationGauss(-nStdDev*maxSig+x0,nStdDev*maxSig+x0,-nStdDev*maxSig+y0,nStdDev*maxSig+y0)
350////
351//
352////    Variable elap = (StopMSTimer(-2) - t1)/1e6
353////    Print "elapsed time = ",elap
354////    Print "time for 16384 = (minutes)",16384*elap/60
355//      return(0)
356//End
357
358
359//// this is called each time to integrate the gaussian
360//Function do2dIntegrationGauss(xMin,xMax,yMin,yMax)
361//      Variable xMin,xMax,yMin,yMax
362//     
363//      Variable/G globalXmin=xMin
364//      Variable/G globalXmax=xMax
365//      Variable/G globalY
366//                     
367//      Variable result=Integrate1d(Gauss2DFuncOuter,yMin,yMax,2,5)       
368//      KillVariables/z globalXmax,globalXmin,globalY
369//      print "integration of 2D = ",result
370//End
371//
372//Function Gauss2DFuncOuter(inY)
373//      Variable inY
374//     
375//      NVAR globalXmin,globalXmax,globalY
376//      globalY=inY
377//     
378//      return integrate1D(Gauss2DFuncInner,globalXmin,globalXmax,2,5)         
379//End
380//
381//Function Gauss2DFuncInner(inX)
382//      Variable inX
383//     
384//      NVAR globalY
385//      Wave coef=coef
386//     
387//      return Gauss2D_theta(coef,inX,GlobalY)
388//End
389//
390//Function Gauss2D_theta(w,x,y)
391//      Wave w
392//      Variable x,y
393//     
394//      Variable val,a,b,c
395//      Variable scale,x0,y0,sx,sy,theta,normFactor
396//     
397//      scale = w[0]
398//      x0 = w[1]
399//      y0 = w[2]
400//      sx = w[3]
401//      sy = w[4]
402//      theta = w[5]
403//     
404//      a = cos(theta)^2/(2*sx*sx) + sin(theta)^2/(2*sy*sy)
405//      b = -1*sin(2*theta)/(4*sx*sx) + sin(2*theta)/(4*sy*sy)
406//      c = sin(theta)^2/(2*sx*sx) + cos(theta)^2/(2*sy*sy)
407//     
408//      val = exp(-1*(a*(x-x0)^2 + 2*b*(x-x0)*(y-y0) + c*(y-y0)^2))
409//     
410//      normFactor = pi/sqrt(a*c-b*b)
411//     
412//      return(scale*val/normFactor)
413//end
414//
415////////////////////////////////////////////////////////////
416////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.