source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/VC_DetectorBinning_Utils.ipf @ 954

Last change on this file since 954 was 954, checked in by srkline, 8 years ago

converted and renamed VSANS files. "VC_" prefix for VCALC related files (that's all of them right now), and moved and renamed parts of files so that the ipf names are more logical now with the contents. Deleted the "V_" prefix files. Added a lengthy routine to be able to write out a VSANS file in HDF format. This is NOT the final and approved data format, only a working version so that I can test things out...

File size: 26.2 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2
3
4// TODO: hard wired for a sphere - change this to allow minimal selections and altering of coefficients
5// TODO: add the "fake" 2D simulation to fill the panels which are then later averaged as I(Q)
6Function FillPanel_wModelData(det,qTot,type)
7        Wave det,qTot
8        String type
9
10        SetDataFolder root:Packages:NIST:VSANS:VCALC:Front
11
12        // q-values and detector arrays already allocated and calculated
13        Duplicate/O det tmpInten,tmpSig,prob_i
14               
15        Variable imon,trans,thick,sdd,pixSizeX,pixSizeY,sdd_offset
16
17        //imon = V_BeamIntensity()*CountTime
18        imon = VCALC_getImon()          //TODO: currently from the panel, not calculated
19        trans = 0.8
20        thick = 0.1
21       
22        // need SDD
23        // need pixel dimensions
24        // nominal sdd in meters, offset in mm, want result in cm !
25        sdd = VCALC_getSDD(type)*100    +  VCALC_getTopBottomSDDOffset(type) / 10               // result is sdd in [cm]
26
27        pixSizeX = VCALC_getPixSizeX(type)              // cm
28        pixSizeY = VCALC_getPixSizeY(type)
29       
30       
31        //?? pick the function from a popup on the panel? (bypass the analysis panel, or maybe it's better to
32        //  keep the panel to keep people used to using it.)
33        String funcStr = VCALC_getModelFunctionStr()
34        strswitch(funcStr)
35                case "Big Debye":
36                        tmpInten = V_Debye(10,3000,0.0001,qTot[p][q])
37                        break
38                case "Big Sphere":
39                        tmpInten = V_SphereForm(1,900,1e-6,0.01,qTot[p][q])     
40                        break
41                case "Debye":
42                        tmpInten = V_Debye(10,300,0.0001,qTot[p][q])
43                        break
44                case "Sphere":
45                        tmpInten = V_SphereForm(1,60,1e-6,0.001,qTot[p][q])     
46                        break
47                default:
48                        tmpInten = V_Debye(10,300,0.1,qTot[p][q])
49        endswitch
50
51
52///////////////
53//      // calculate the scattering cross section simply to be able to estimate the transmission
54//      Variable sig_sas=0
55//     
56//      // remember that the random deviate is the coherent portion ONLY - the incoherent background is
57//      // subtracted before the calculation.
58//      CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
59//
60//      if(sig_sas > 100)
61//              DoAlert 0,"SAS cross section > 100. Estimates of multiple scattering are unreliable. Choosing a model with a well-defined Rg may help"
62//      endif           
63//
64//      // calculate the multiple scattering fraction for display (10/2009)
65//      Variable ii,nMax=10,tau
66//      mScat=0
67//      tau = thick*sig_sas
68//      // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
69//      // neutrons out of those that were scattered
70//      for(ii=2;ii<nMax;ii+=1)
71//              mScat += tau^(ii)/factorial(ii)
72////            print tau^(ii)/factorial(ii)
73//      endfor
74//      estTrans = exp(-1*thick*sig_sas)                //thickness and sigma both in units of cm
75//      mscat *= (estTrans)/(1-estTrans)
76//
77////    if(mScat > 0.1)         //  Display warning
78//
79//      Print "Sig_sas = ",sig_sas
80////////////////////
81       
82        prob_i = trans*thick*pixSizeX*pixSizeY/(sdd)^2*tmpInten                 //probability of a neutron in q-bin(i)
83               
84        tmpInten = (imon)*prob_i                //tmpInten is not the model calculation anymore!!
85
86
87/// **** can I safely assume a Gaussian error in the count rate??
88        tmpSig = sqrt(tmpInten)         // corrected based on John's memo, from 8/9/99
89
90        tmpInten += gnoise(tmpSig)
91        tmpInten = (tmpInten[p][q] < 0) ? 0 : tmpInten[p][q]                    // MAR 2013 -- is this the right thing to do
92        tmpInten = trunc(tmpInten)
93               
94       
95        det = tmpInten
96
97// if I want "absolute" scale -- then I lose the integer nature of the detector (but keep the random)
98//      det /= trans*thick*pixSizeX*pixSizeY/(sdd)^2*imon
99
100       
101        KillWaves/Z tmpInten,tmpSig,prob_i     
102        SetDataFolder root:
103
104        return(0)
105End
106
107
108// For a given detector panel, calculate the q-values
109// -work with everything as arrays
110// Input needed:
111// detector data
112// detector type (LRTB?)
113// beam center (may be off the detector)
114// SDD
115// lambda
116//
117// pixel dimensions for detector type (global constants)
118// - data dimensions read directly from array
119//
120// --What is calculated:
121// array of Q
122// array of qx,qy,qz
123// array of error already exists
124//
125//
126// -- sdd in meters
127// -- lambda in Angstroms
128Function V_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
129        Wave data,qTot,qx,qy,qz
130        Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
131               
132        // loop over the array and calculate the values - this is done as a wave assignment
133// TODO -- be sure that it's p,q -- or maybe p+1,q+1 as used in WriteQIS.ipf   
134        qTot = V_CalcQval(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
135        qx = V_CalcQX(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
136        qy = V_CalcQY(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
137        qz = V_CalcQZ(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
138       
139        return(0)
140End
141
142
143//////////////////////
144// NOTE: The Q calculations are different than what is in GaussUtils in that they take into
145// accout the different x/y pixel sizes and the beam center not being on the detector -
146// off a different edge for each LRTB type
147/////////////////////
148
149//function to calculate the overall q-value, given all of the necesary trig inputs
150//NOTE: detector locations passed in are pixels = 0.5cm real space on the detector
151//and are in detector coordinates (1,128) rather than axis values
152//the pixel locations need not be integers, reals are ok inputs
153//sdd is in meters
154//wavelength is in Angstroms
155//
156//returned magnitude of Q is in 1/Angstroms
157//
158Function V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
159        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
160       
161        Variable dx,dy,qval,two_theta,dist
162               
163        sdd *=100               //convert to cm
164        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
165        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
166        dist = sqrt(dx^2 + dy^2)
167       
168        two_theta = atan(dist/sdd)
169
170        qval = 4*Pi/lam*sin(two_theta/2)
171       
172        return qval
173End
174
175//calculates just the q-value in the x-direction on the detector
176//input/output is the same as CalcQval()
177//ALL inputs are in detector coordinates
178//
179//NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector
180//sdd is in meters
181//wavelength is in Angstroms
182//
183// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
184// now properly accounts for qz
185//
186Function V_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
187        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
188
189        Variable qx,qval,phi,dx,dy,dist,two_theta
190       
191        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
192       
193        sdd *=100               //convert to cm
194        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
195        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
196        phi = V_FindPhi(dx,dy)
197       
198        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
199        dist = sqrt(dx^2 + dy^2)
200        two_theta = atan(dist/sdd)
201
202        qx = qval*cos(two_theta/2)*cos(phi)
203       
204        return qx
205End
206
207//calculates just the q-value in the y-direction on the detector
208//input/output is the same as CalcQval()
209//ALL inputs are in detector coordinates
210//NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector
211//sdd is in meters
212//wavelength is in Angstroms
213//
214// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
215// now properly accounts for qz
216//
217Function V_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
218        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
219       
220        Variable dy,qval,dx,phi,qy,dist,two_theta
221       
222        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
223       
224        sdd *=100               //convert to cm
225        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
226        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
227        phi = V_FindPhi(dx,dy)
228       
229        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
230        dist = sqrt(dx^2 + dy^2)
231        two_theta = atan(dist/sdd)
232       
233        qy = qval*cos(two_theta/2)*sin(phi)
234       
235        return qy
236End
237
238//calculates just the z-component of the q-vector, not measured on the detector
239//input/output is the same as CalcQval()
240//ALL inputs are in detector coordinates
241//NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector
242//sdd is in meters
243//wavelength is in Angstroms
244//
245// not actually used, but here for completeness if anyone asks
246//
247Function V_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
248        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
249       
250        Variable dy,qval,dx,phi,qz,dist,two_theta
251       
252        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
253       
254        sdd *=100               //convert to cm
255       
256        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
257        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
258        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
259        dist = sqrt(dx^2 + dy^2)
260        two_theta = atan(dist/sdd)
261       
262        qz = qval*sin(two_theta/2)
263       
264        return qz
265End
266
267//phi is defined from +x axis, proceeding CCW around [0,2Pi]
268Threadsafe Function V_FindPhi(vx,vy)
269        variable vx,vy
270       
271        variable phi
272       
273        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
274       
275        // special cases
276        if(vx==0 && vy > 0)
277                return(pi/2)
278        endif
279        if(vx==0 && vy < 0)
280                return(3*pi/2)
281        endif
282        if(vx >= 0 && vy == 0)
283                return(0)
284        endif
285        if(vx < 0 && vy == 0)
286                return(pi)
287        endif
288       
289       
290        if(vx > 0 && vy > 0)
291                return(phi)
292        endif
293        if(vx < 0 && vy > 0)
294                return(phi + pi)
295        endif
296        if(vx < 0 && vy < 0)
297                return(phi + pi)
298        endif
299        if( vx > 0 && vy < 0)
300                return(phi + 2*pi)
301        endif
302       
303        return(phi)
304end
305
306Function V_SphereForm(scale,radius,delrho,bkg,x)                               
307        Variable scale,radius,delrho,bkg
308        Variable x
309       
310        // variables are:                                                       
311        //[0] scale
312        //[1] radius (A)
313        //[2] delrho (A-2)
314        //[3] background (cm-1)
315       
316//      Variable scale,radius,delrho,bkg                               
317//      scale = w[0]
318//      radius = w[1]
319//      delrho = w[2]
320//      bkg = w[3]
321       
322       
323        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
324        // and is rescaled to give [=] cm^-1
325       
326        Variable bes,f,vol,f2
327        //
328        //handle q==0 separately
329        If(x==0)
330                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
331                return(f)
332        Endif
333       
334//      bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
335       
336        bes = 3*sqrt(pi/(2*x*radius))*BesselJ(1.5,x*radius)/(x*radius)
337       
338        vol = 4*pi/3*radius^3
339        f = vol*bes*delrho              // [=] A
340        // normalize to single particle volume, convert to 1/cm
341        f2 = f * f / vol * 1.0e8                // [=] 1/cm
342       
343        return (scale*f2+bkg)   // Scale, then add in the background
344       
345End
346
347Function V_Debye(scale,rg,bkg,x)
348        Variable scale,rg,bkg
349        Variable x
350       
351        // variables are:
352        //[0] scale factor
353        //[1] radius of gyration [A]
354        //[2] background        [cm-1]
355       
356        // calculates (scale*debye)+bkg
357        Variable Pq,qr2
358       
359        qr2=(x*rg)^2
360        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
361       
362        //scale
363        Pq *= scale
364        // then add in the background
365        return (Pq+bkg)
366End
367
368
369
370
371Function SetDeltaQ(folderStr,type)
372        String folderStr,type
373
374        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + folderStr + ":det_"+type)            // 2D detector data
375       
376        Variable xDim,yDim,delQ
377       
378        xDim=DimSize(inten,0)
379        yDim=DimSize(inten,1)
380       
381        if(xDim<yDim)
382                WAVE qx = $("root:Packages:NIST:VSANS:VCALC:" + folderStr + ":qx_"+type)
383                delQ = abs(qx[0][0] - qx[1][0])/2
384        else
385                WAVE qy = $("root:Packages:NIST:VSANS:VCALC:" + folderStr + ":qy_"+type)
386                delQ = abs(qy[0][1] - qy[0][0])/2
387        endif
388       
389        // set the global
390        Variable/G $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_"+type) = delQ
391//      Print "SET delQ = ",delQ," for ",type
392       
393        return(0)
394end
395
396
397//TODO -- need a switch here to dispatch to the averaging type
398Proc V_BinQxQy_to_1D(folderStr,type)
399        String folderStr
400        String type
401//      Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4)
402//      Prompt type,"detector identifier"
403
404
405        V_fDoBinning_QxQy2D(folderStr, type)
406
407
408/// this is for a tall, narrow slit mode       
409//      V_fBinDetector_byRows(folderStr,type)
410       
411End
412
413
414Proc V_Graph_1D_detType(folderStr,type)
415        String folderStr,type
416       
417        SetDataFolder root:Packages:NIST:VSANS:VCALC
418       
419        Display $("iBin_qxqy"+"_"+type) vs $("qBin_qxqy"+"_"+type)
420        ModifyGraph mirror=2,grid=1,log=1
421        ModifyGraph mode=4,marker=19,msize=2
422//      ErrorBars/T=0 iBin_qxqy Y,wave=(eBin2D_qxqy,eBin2D_qxqy)                // for simulations, I don't have 2D uncertainty
423        ErrorBars/T=0 $("iBin_qxqy"+"_"+type) Y,wave=($("eBin_qxqy"+"_"+type),$("eBin_qxqy"+"_"+type))
424        legend
425       
426        SetDataFolder root:
427
428End
429
430
431// see the equivalent function in PlotUtils2D_v40.ipf
432//
433//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
434//
435// this has been modified to accept different detector panels and to take arrays
436// -- type = FL or FR or...other panel identifiers
437//
438// TODO "iErr" is all messed up since it doesn't really apply here for data that is not 2D simulation
439//
440//
441Function V_fDoBinning_QxQy2D(folderStr,type)
442        String folderStr,type
443       
444        Variable nSets = 0
445        Variable xDim,yDim
446        Variable ii,jj
447        Variable qVal,nq,var,avesq,aveisq
448        Variable binIndex,val
449
450       
451        SetDataFolder root:Packages:NIST:VSANS:VCALC
452       
453// now switch on the type to determine which waves to declare and create
454// since there may be more than one panel to step through. There may be two, there may be four
455//
456
457        strswitch(type) // string switch
458                case "FL":              // execute if case matches expression
459                case "FR":
460                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_FL")
461                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+type)              // 2D detector data
462                        WAVE/Z iErr = $("iErr_"+type)                   // 2D errors -- may not exist, especially for simulation
463                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+type)                     // 2D q-values
464                        nSets = 1
465                        break   
466                                                               
467                case "FT":             
468                case "FB":
469                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_FT")
470                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+type)              // 2D detector data
471                        WAVE/Z iErr = $("iErr_"+type)                   // 2D errors -- may not exist, especially for simulation
472                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+type)                     // 2D q-values
473                        nSets = 1
474                        break
475                       
476                case "ML":             
477                case "MR":             
478                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_ML")
479                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+type)             // 2D detector data
480                        WAVE/Z iErr = $("iErr_"+type)                   // 2D errors -- may not exist, especially for simulation
481                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+type)                    // 2D q-values
482                        nSets = 1
483                        break   
484                                       
485                case "MT":             
486                case "MB":             
487                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_MT")
488                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+type)             // 2D detector data
489                        WAVE/Z iErr = $("iErr_"+type)                   // 2D errors -- may not exist, especially for simulation
490                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+type)                    // 2D q-values
491                        nSets = 1
492                        break   
493                                       
494                case "B":               
495                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_B")
496                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Back" + ":det_"+type)               // 2D detector data
497                        WAVE/Z iErr = $("iErr_"+type)                   // 2D errors -- may not exist, especially for simulation
498                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Back" +":qTot_"+type)                      // 2D q-values
499                        nSets = 1
500                        break   
501                       
502                case "FLR":
503                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_FL")
504                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FL")              // 2D detector data
505                        WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation
506                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FL")                     // 2D q-values
507                        WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FR")             // 2D detector data
508                        WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation
509                        Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FR")                    // 2D q-values
510                        nSets = 2
511                        break                   
512               
513                case "FTB":
514                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_FT")
515                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FT")              // 2D detector data
516                        WAVE/Z iErr = $("iErr_"+"FT")                   // 2D errors -- may not exist, especially for simulation
517                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FT")                     // 2D q-values
518                        WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FB")             // 2D detector data
519                        WAVE/Z iErr2 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation
520                        Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FB")                    // 2D q-values
521                        nSets = 2
522                        break           
523               
524                case "FLRTB":
525                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_FL")
526                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FL")              // 2D detector data
527                        WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation
528                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FL")                     // 2D q-values
529                        WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FR")             // 2D detector data
530                        WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation
531                        Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FR")                    // 2D q-values
532                        WAVE inten3 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FT")             // 2D detector data
533                        WAVE/Z iErr3 = $("iErr_"+"FT")                  // 2D errors -- may not exist, especially for simulation
534                        Wave qTotal3 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FT")                    // 2D q-values
535                        WAVE inten4 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FB")             // 2D detector data
536                        WAVE/Z iErr4 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation
537                        Wave qTotal4 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FB")                    // 2D q-values
538                        nSets = 4
539                        break           
540                       
541
542                case "MLR":
543                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_ML")
544                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"ML")             // 2D detector data
545                        WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation
546                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"ML")                    // 2D q-values
547                        WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MR")            // 2D detector data
548                        WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation
549                        Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MR")                   // 2D q-values
550                        nSets = 2
551                        break                   
552               
553                case "MTB":
554                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_MT")
555                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MT")             // 2D detector data
556                        WAVE/Z iErr = $("iErr_"+"MT")                   // 2D errors -- may not exist, especially for simulation
557                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MT")                    // 2D q-values
558                        WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MB")            // 2D detector data
559                        WAVE/Z iErr2 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation
560                        Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MB")                   // 2D q-values
561                        nSets = 2
562                        break                           
563               
564                case "MLRTB":
565                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_ML")
566                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"ML")             // 2D detector data
567                        WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation
568                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"ML")                    // 2D q-values
569                        WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MR")            // 2D detector data
570                        WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation
571                        Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MR")                   // 2D q-values
572                        WAVE inten3 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MT")            // 2D detector data
573                        WAVE/Z iErr3 = $("iErr_"+"MT")                  // 2D errors -- may not exist, especially for simulation
574                        Wave qTotal3 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MT")                   // 2D q-values
575                        WAVE inten4 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MB")            // 2D detector data
576                        WAVE/Z iErr4 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation
577                        Wave qTotal4 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MB")                   // 2D q-values
578                        nSets = 4
579                        break                                                                   
580                                       
581                default:
582                        nSets = 0                                                       // optional default expression executed
583                        Print "ERROR   ---- type is not recognized "
584        endswitch
585
586//      Print "delQ = ",delQ," for ",type
587
588        if(nSets == 0)
589                return(0)
590        endif
591
592
593//TODO: properly define the errors here - I'll have this if I do the simulation
594        if(WaveExists(iErr)==0)
595                Duplicate/O inten,iErr
596                Wave iErr=iErr
597//              iErr = 1+sqrt(inten+0.75)                       // can't use this -- it applies to counts, not intensity (already a count rate...)
598                iErr = sqrt(inten+0.75)                 // TODO -- here I'm just using some fictional value
599        endif
600
601        nq = 600
602
603        // note that the back panel of 320x320 (1mm res) results in 447 data points!
604        // - so I upped nq to 600
605
606//******TODO****** -- where to put the averaged data -- right now, folderStr is forced to ""   
607//      SetDataFolder $("root:"+folderStr)              //should already be here, but make sure...     
608        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"iBin_qxqy"+"_"+type)
609        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"qBin_qxqy"+"_"+type)
610        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"nBin_qxqy"+"_"+type)
611        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"iBin2_qxqy"+"_"+type)
612        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"eBin_qxqy"+"_"+type)
613        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"eBin2D_qxqy"+"_"+type)
614       
615        Wave iBin_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"iBin_qxqy_"+type)
616        Wave qBin_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"qBin_qxqy"+"_"+type)
617        Wave nBin_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"nBin_qxqy"+"_"+type)
618        Wave iBin2_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"iBin2_qxqy"+"_"+type)
619        Wave eBin_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"eBin_qxqy"+"_"+type)
620        Wave eBin2D_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"eBin2D_qxqy"+"_"+type)
621       
622       
623//      delQ = abs(sqrt(qx[2]^2+qy[2]^2+qz[2]^2) - sqrt(qx[1]^2+qy[1]^2+qz[1]^2))               //use bins of 1 pixel width
624// TODO: not sure if I want to so dQ in x or y direction...
625        // the short dimension is the 8mm tubes, use this direction as dQ?
626        // but don't use the corner of the detector, since dQ will be very different on T/B or L/R due to the location of [0,0]
627        // WRT the beam center. use qx or qy directly. Still not happy with this way...
628
629
630        qBin_qxqy[] =  p*delQ   
631        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
632
633        iBin_qxqy = 0
634        iBin2_qxqy = 0
635        eBin_qxqy = 0
636        eBin2D_qxqy = 0
637        nBin_qxqy = 0   //number of intensities added to each bin
638
639// now there are situations of:
640// 1 panel
641// 2 panels
642// 4 panels
643//
644// this needs to be a double loop now...
645
646// use set 1 (no number) only
647        if(nSets >= 1)
648                xDim=DimSize(inten,0)
649                yDim=DimSize(inten,1)
650       
651                for(ii=0;ii<xDim;ii+=1)
652                        for(jj=0;jj<yDim;jj+=1)
653                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
654                                qVal = qTotal[ii][jj]
655                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
656                                val = inten[ii][jj]
657                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
658                                        iBin_qxqy[binIndex] += val
659                                        iBin2_qxqy[binIndex] += val*val
660                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
661                                        nBin_qxqy[binIndex] += 1
662                                endif
663                        endfor
664                endfor
665               
666        endif
667
668// add in set 2 (set 1 already done)
669        if(nSets >= 2)
670                xDim=DimSize(inten2,0)
671                yDim=DimSize(inten2,1)
672       
673                for(ii=0;ii<xDim;ii+=1)
674                        for(jj=0;jj<yDim;jj+=1)
675                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
676                                qVal = qTotal2[ii][jj]
677                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
678                                val = inten2[ii][jj]
679                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
680                                        iBin_qxqy[binIndex] += val
681                                        iBin2_qxqy[binIndex] += val*val
682                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
683                                        nBin_qxqy[binIndex] += 1
684                                endif
685                        endfor
686                endfor
687               
688        endif
689
690// add in set 3 and 4 (set 1 and 2already done)
691        if(nSets == 4)
692                xDim=DimSize(inten3,0)
693                yDim=DimSize(inten3,1)
694       
695                for(ii=0;ii<xDim;ii+=1)
696                        for(jj=0;jj<yDim;jj+=1)
697                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
698                                qVal = qTotal3[ii][jj]
699                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
700                                val = inten3[ii][jj]
701                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
702                                        iBin_qxqy[binIndex] += val
703                                        iBin2_qxqy[binIndex] += val*val
704                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
705                                        nBin_qxqy[binIndex] += 1
706                                endif
707                        endfor
708                endfor
709               
710               
711                xDim=DimSize(inten4,0)
712                yDim=DimSize(inten4,1)
713       
714                for(ii=0;ii<xDim;ii+=1)
715                        for(jj=0;jj<yDim;jj+=1)
716                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
717                                qVal = qTotal4[ii][jj]
718                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
719                                val = inten4[ii][jj]
720                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
721                                        iBin_qxqy[binIndex] += val
722                                        iBin2_qxqy[binIndex] += val*val
723                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
724                                        nBin_qxqy[binIndex] += 1
725                                endif
726                        endfor
727                endfor
728               
729        endif
730
731
732// after looping through all of the data on the panels, calculate errors on I(q),
733// just like in CircSectAve.ipf
734        for(ii=0;ii<nq;ii+=1)
735                if(nBin_qxqy[ii] == 0)
736                        //no pixels in annuli, data unknown
737                        iBin_qxqy[ii] = 0
738                        eBin_qxqy[ii] = 1
739                        eBin2D_qxqy[ii] = NaN
740                else
741                        if(nBin_qxqy[ii] <= 1)
742                                //need more than one pixel to determine error
743                                iBin_qxqy[ii] /= nBin_qxqy[ii]
744                                eBin_qxqy[ii] = 1
745                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
746                        else
747                                //assume that the intensity in each pixel in annuli is normally distributed about mean...
748                                iBin_qxqy[ii] /= nBin_qxqy[ii]
749                                avesq = iBin_qxqy[ii]^2
750                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
751                                var = aveisq-avesq
752                                if(var<=0)
753                                        eBin_qxqy[ii] = 1e-6
754                                else
755                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
756                                endif
757                                // and calculate as it is propagated pixel-by-pixel
758                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
759                        endif
760                endif
761        endfor
762       
763        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo
764       
765        // find the last non-zero point, working backwards
766        val=nq
767        do
768                val -= 1
769        while((nBin_qxqy[val] == 0) && val > 0)
770       
771//      print val, nBin_qxqy[val]
772        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
773
774        if(val == 0)
775                // all the points were deleted
776                return(0)
777        endif
778       
779       
780        // since the beam center is not always on the detector, many of the low Q bins will have zero pixels
781        // find the first non-zero point, working forwards
782        val = -1
783        do
784                val += 1
785        while(nBin_qxqy[val] == 0)     
786        DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
787
788        // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
789        val = numpnts(nBin_qxqy)-1
790        do
791                if(nBin_qxqy[val] == 0)
792                        DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
793                endif
794                val -= 1
795        while(val>0)
796       
797        SetDataFolder root:
798       
799        return(0)
800End
801
802
Note: See TracBrowser for help on using the repository browser.