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

Last change on this file since 1019 was 1019, checked in by srkline, 6 years ago

adding procedures for:

simple save of a DIV file. no functionality to generate a DIV file yet, since I don't know how this will happen.

Simple dump of the file structure in a data "tree"

Verified that the error bars on the I(q) data are correctly calculated as standard error of the mean. There was never an issue with this, or with SANS calculations.

Started filling in "Correct" routines based on the SANS version. Only stubs present currently.

File size: 39.6 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2
3/////////////////////////
4//
5// Utility functions to:
6//              calculate Q, Qx, Qy, Qz
7//              fill the detector panels with simulated data (the model functions are here)
8//              bin the 2D detector to 1D I(Q) based on Q and deltaQ (bin width)
9//
10/////////////////////////
11
12
13
14
15// TODO: hard wired for a sphere - change this to allow minimal selections and altering of coefficients
16// TODO: add the "fake" 2D simulation to fill the panels which are then later averaged as I(Q)
17//
18// NOTE - this is a VCALC only routine, so it's not been made completely generic
19//
20Function FillPanel_wModelData(det,qTot,type)
21        Wave det,qTot
22        String type
23
24//      SetDataFolder root:Packages:NIST:VSANS:VCALC:Front
25
26        // q-values and detector arrays already allocated and calculated
27        Duplicate/O det tmpInten,tmpSig,prob_i
28               
29        Variable imon,trans,thick,sdd,pixSizeX,pixSizeY,sdd_offset
30
31        //imon = VC_BeamIntensity()*CountTime
32        imon = VCALC_getImon()          //TODO: currently from the panel, not calculated
33        trans = 0.8
34        thick = 0.1
35       
36        // need SDD
37        // need pixel dimensions
38        // nominal sdd in meters, offset in mm, want result in cm !
39        sdd = VCALC_getSDD(type)*100    +  VCALC_getTopBottomSDDOffset(type) / 10               // result is sdd in [cm]
40
41        pixSizeX = VCALC_getPixSizeX(type)              // cm
42        pixSizeY = VCALC_getPixSizeY(type)
43       
44       
45        //?? pick the function from a popup on the panel? (bypass the analysis panel, or maybe it's better to
46        //  keep the panel to keep people used to using it.)
47        // peak @ 0.1 ~ AgBeh
48        //      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}               
49        //
50        // peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
51        //      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}             
52        String funcStr = VCALC_getModelFunctionStr()
53        strswitch(funcStr)
54                case "Big Debye":
55                        tmpInten = VC_Debye(100,3000,0.0001,qTot[p][q])
56                        break
57                case "Big Sphere":
58                        tmpInten = VC_SphereForm(1,900,1e-6,0.01,qTot[p][q])   
59                        break
60                case "Debye":
61                        tmpInten = VC_Debye(10,300,0.0001,qTot[p][q])
62                        break
63                case "Sphere":
64                        tmpInten = VC_SphereForm(1,60,1e-6,0.001,qTot[p][q])   
65                        break
66                case "AgBeh":
67                        tmpInten = VC_BroadPeak(1e-9,3,20,100.0,0.1,3,0.1,qTot[p][q])
68                        break
69                case "Vycor":
70                        tmpInten = VC_BroadPeak(1e-9,3,20,500.0,0.015,3,0.1,qTot[p][q])
71                        break   
72                case "Empty Cell":
73                        tmpInten = VC_EC_Empirical(2.2e-8,3.346,0.0065,9.0,0.016,qTot[p][q])
74                        break
75                case "Blocked Beam":
76                        tmpInten = VC_BlockedBeam(1,qTot[p][q])
77                        break
78                default:
79                        tmpInten = VC_Debye(10,300,0.1,qTot[p][q])
80        endswitch
81
82// TODO: this is faked to get around the singularity at the center of the back detector
83//
84//
85        if(cmpstr(type,"B") == 0)
86                Variable nx,ny,px,py
87                nx = VCALC_get_nPix_X(type)
88                ny = VCALC_get_nPix_Y(type)
89                px = trunc(nx/2)
90                py = trunc(ny/2)
91               
92                tmpInten[px][py] = (tmpInten[px][py+1] + tmpInten[px][py-1])/2
93        endif
94
95
96
97///////////////
98//      // calculate the scattering cross section simply to be able to estimate the transmission
99//      Variable sig_sas=0
100//     
101//      // remember that the random deviate is the coherent portion ONLY - the incoherent background is
102//      // subtracted before the calculation.
103//      CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
104//
105//      if(sig_sas > 100)
106//              DoAlert 0,"SAS cross section > 100. Estimates of multiple scattering are unreliable. Choosing a model with a well-defined Rg may help"
107//      endif           
108//
109//      // calculate the multiple scattering fraction for display (10/2009)
110//      Variable ii,nMax=10,tau
111//      mScat=0
112//      tau = thick*sig_sas
113//      // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
114//      // neutrons out of those that were scattered
115//      for(ii=2;ii<nMax;ii+=1)
116//              mScat += tau^(ii)/factorial(ii)
117////            print tau^(ii)/factorial(ii)
118//      endfor
119//      estTrans = exp(-1*thick*sig_sas)                //thickness and sigma both in units of cm
120//      mscat *= (estTrans)/(1-estTrans)
121//
122////    if(mScat > 0.1)         //  Display warning
123//
124//      Print "Sig_sas = ",sig_sas
125////////////////////
126       
127        prob_i = trans*thick*pixSizeX*pixSizeY/(sdd)^2*tmpInten                 //probability of a neutron in q-bin(i)
128               
129        tmpInten = (imon)*prob_i                //tmpInten is not the model calculation anymore!!
130
131
132/// **** can I safely assume a Gaussian error in the count rate??
133        tmpSig = sqrt(tmpInten)         // corrected based on John's memo, from 8/9/99
134
135        tmpInten += gnoise(tmpSig)
136        tmpInten = (tmpInten[p][q] < 0) ? 0 : tmpInten[p][q]                    // MAR 2013 -- is this the right thing to do
137        tmpInten = trunc(tmpInten)
138               
139       
140        det = tmpInten
141
142// if I want "absolute" scale -- then I lose the integer nature of the detector (but keep the random)
143//      det /= trans*thick*pixSizeX*pixSizeY/(sdd)^2*imon
144
145       
146        KillWaves/Z tmpInten,tmpSig,prob_i     
147        SetDataFolder root:
148
149        return(0)
150End
151
152
153// For a given detector panel, calculate the q-values
154// -work with everything as arrays
155// Input needed:
156// detector data
157// detector type (LRTB?)
158// beam center (may be off the detector)
159// SDD
160// lambda
161//
162// pixel dimensions for detector type (global constants)
163// - data dimensions read directly from array
164//
165// --What is calculated:
166// array of Q
167// array of qx,qy,qz
168// array of error already exists
169//
170//
171// -- sdd in meters
172// -- lambda in Angstroms
173Function VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
174        Wave data,qTot,qx,qy,qz
175        Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
176               
177        // loop over the array and calculate the values - this is done as a wave assignment
178// TODO -- be sure that it's p,q -- or maybe p+1,q+1 as used in WriteQIS.ipf   
179        qTot = VC_CalcQval(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
180        qx = VC_CalcQX(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
181        qy = VC_CalcQY(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
182        qz = VC_CalcQZ(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
183       
184        return(0)
185End
186
187
188// for testing, a version that will calculate the q-arrays for VCALC based on whatever nonlinear coefficients
189// exist in the RAW data folder
190//
191// reverts to the "regular" linear detector if waves not found or a flag is set
192//
193// need to convert the beam center to mm
194// need to call the VSANS V_CalcQval routines (these use the real-space distance, not pixel dims)
195//
196// TODO:
197// -- tube width is hard-wired in
198//
199//
200Function VC_Detector_2Q_NonLin(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY,detStr)
201        Wave data,qTot,qx,qy,qz
202        Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
203        String detStr
204       
205        String destPath = "root:Packages:NIST:VSANS:RAW"
206                       
207        Wave/Z data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
208        Wave/Z data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
209        NVAR gUseNonLinearDet = root:Packages:NIST:VSANS:VCALC:gUseNonLinearDet
210       
211        if(gUseNonLinearDet && WaveExists(data_realDistX) && WaveExists(data_realDistY))
212                // convert the beam centers to mm
213                String orientation
214                Variable dimX,dimY,newX,newY
215                dimX = DimSize(data_realDistX,0)
216                dimY = DimSize(data_realDistX,1)
217                if(dimX > dimY)
218                        orientation = "horizontal"
219                else
220                        orientation = "vertical"
221                endif
222               
223                Variable tube_width = 8.4               //mm
224       
225        //
226                if(cmpstr(orientation,"vertical")==0)
227                        //      this is data dimensioned as (Ntubes,Npix)
228                        newX = tube_width*xCtr
229                        newY = data_realDistY[0][yCtr]
230                else
231                        //      this is data (horizontal) dimensioned as (Npix,Ntubes)
232                        newX = data_realDistX[xCtr][0]
233                        newY = tube_width*yCtr
234                endif   
235
236                //if detector "B", different calculation for the centers (not tubes)
237                if(cmpstr(detStr,"B")==0)
238                        newX = data_realDistX[xCtr][0]
239                        newY = data_realDistY[0][yCtr]
240                        //newX = xCtr
241                        //newY = yCtr
242                endif           
243                               
244                // calculate all of the q-values
245                qTot = V_CalcQval(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
246                qx = V_CalcQX(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
247                qy = V_CalcQY(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
248                qz = V_CalcQZ(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
249       
250//              Print "det, x_mm, y_mm ",detStr,num2str(newX),num2str(newY)
251//              Print "det, x_pix, y_pix ",detStr,num2str(xCtr),num2str(yCtr)
252        else
253                // do the q-calculation using linear detector
254                VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
255        endif
256       
257        return(0)
258End
259
260
261//////////////////////
262// NOTE: The Q calculations are different than what is in GaussUtils in that they take into
263// accout the different x/y pixel sizes and the beam center not being on the detector -
264// off a different edge for each LRTB type
265/////////////////////
266
267//function to calculate the overall q-value, given all of the necesary trig inputs
268//and are in detector coordinates (1,128) rather than axis values
269//the pixel locations need not be integers, reals are ok inputs
270//sdd is in meters
271//wavelength is in Angstroms
272//
273//returned magnitude of Q is in 1/Angstroms
274//
275Function VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
276        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
277       
278        Variable dx,dy,qval,two_theta,dist
279               
280        sdd *=100               //convert to cm
281        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
282        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
283        dist = sqrt(dx^2 + dy^2)
284       
285        two_theta = atan(dist/sdd)
286
287        qval = 4*Pi/lam*sin(two_theta/2)
288       
289        return qval
290End
291
292//calculates just the q-value in the x-direction on the detector
293//input/output is the same as CalcQval()
294//ALL inputs are in detector coordinates
295//
296//sdd is in meters
297//wavelength is in Angstroms
298//
299// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
300// now properly accounts for qz
301//
302Function VC_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
303        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
304
305        Variable qx,qval,phi,dx,dy,dist,two_theta
306       
307        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
308       
309        sdd *=100               //convert to cm
310        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
311        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
312        phi = V_FindPhi(dx,dy)
313       
314        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
315        dist = sqrt(dx^2 + dy^2)
316        two_theta = atan(dist/sdd)
317
318        qx = qval*cos(two_theta/2)*cos(phi)
319       
320        return qx
321End
322
323//calculates just the q-value in the y-direction on the detector
324//input/output is the same as CalcQval()
325//ALL inputs are in detector coordinates
326//sdd is in meters
327//wavelength is in Angstroms
328//
329// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
330// now properly accounts for qz
331//
332Function VC_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
333        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
334       
335        Variable dy,qval,dx,phi,qy,dist,two_theta
336       
337        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
338       
339        sdd *=100               //convert to cm
340        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
341        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
342        phi = V_FindPhi(dx,dy)
343       
344        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
345        dist = sqrt(dx^2 + dy^2)
346        two_theta = atan(dist/sdd)
347       
348        qy = qval*cos(two_theta/2)*sin(phi)
349       
350        return qy
351End
352
353//calculates just the z-component of the q-vector, not measured on the detector
354//input/output is the same as CalcQval()
355//ALL inputs are in detector coordinates
356//sdd is in meters
357//wavelength is in Angstroms
358//
359// not actually used, but here for completeness if anyone asks
360//
361Function VC_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
362        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
363       
364        Variable dy,qval,dx,phi,qz,dist,two_theta
365       
366        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
367       
368        sdd *=100               //convert to cm
369       
370        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
371        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
372        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
373        dist = sqrt(dx^2 + dy^2)
374        two_theta = atan(dist/sdd)
375       
376        qz = qval*sin(two_theta/2)
377       
378        return qz
379End
380
381//phi is defined from +x axis, proceeding CCW around [0,2Pi]
382Threadsafe Function V_FindPhi(vx,vy)
383        variable vx,vy
384       
385        variable phi
386       
387        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
388       
389        // special cases
390        if(vx==0 && vy > 0)
391                return(pi/2)
392        endif
393        if(vx==0 && vy < 0)
394                return(3*pi/2)
395        endif
396        if(vx >= 0 && vy == 0)
397                return(0)
398        endif
399        if(vx < 0 && vy == 0)
400                return(pi)
401        endif
402       
403       
404        if(vx > 0 && vy > 0)
405                return(phi)
406        endif
407        if(vx < 0 && vy > 0)
408                return(phi + pi)
409        endif
410        if(vx < 0 && vy < 0)
411                return(phi + pi)
412        endif
413        if( vx > 0 && vy < 0)
414                return(phi + 2*pi)
415        endif
416       
417        return(phi)
418end
419
420Function VC_SphereForm(scale,radius,delrho,bkg,x)                               
421        Variable scale,radius,delrho,bkg
422        Variable x
423       
424        // variables are:                                                       
425        //[0] scale
426        //[1] radius (A)
427        //[2] delrho (A-2)
428        //[3] background (cm-1)
429       
430//      Variable scale,radius,delrho,bkg                               
431//      scale = w[0]
432//      radius = w[1]
433//      delrho = w[2]
434//      bkg = w[3]
435       
436       
437        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
438        // and is rescaled to give [=] cm^-1
439       
440        Variable bes,f,vol,f2
441        //
442        //handle q==0 separately
443        If(x==0)
444                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
445                return(f)
446        Endif
447       
448//      bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
449       
450        bes = 3*sqrt(pi/(2*x*radius))*BesselJ(1.5,x*radius)/(x*radius)
451       
452        vol = 4*pi/3*radius^3
453        f = vol*bes*delrho              // [=] A
454        // normalize to single particle volume, convert to 1/cm
455        f2 = f * f / vol * 1.0e8                // [=] 1/cm
456       
457        return (scale*f2+bkg)   // Scale, then add in the background
458       
459End
460
461Function VC_Debye(scale,rg,bkg,x)
462        Variable scale,rg,bkg
463        Variable x
464       
465        // variables are:
466        //[0] scale factor
467        //[1] radius of gyration [A]
468        //[2] background        [cm-1]
469       
470        // calculates (scale*debye)+bkg
471        Variable Pq,qr2
472       
473        qr2=(x*rg)^2
474        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
475       
476        //scale
477        Pq *= scale
478        // then add in the background
479        return (Pq+bkg)
480End
481
482// a sum of a power law and debye to approximate the scattering from a real empty cell
483//
484//      make/O/D coef_ECEmp = {2.2e-8,3.346,0.0065,9.0,0.016}
485//
486Function VC_EC_Empirical(aa,mm,scale,rg,bkg,x)
487        Variable aa,mm,scale,rg,bkg
488        Variable x
489       
490        // variables are:
491        //[0] = A
492        //[1] = power m
493        //[2] scale factor
494        //[3] radius of gyration [A]
495        //[4] background        [cm-1]
496       
497        Variable Iq
498       
499        // calculates (scale*debye)+bkg
500        Variable Pq,qr2
501       
502//      if(x*Rg < 1e-3)         //added Oct 2008 to avoid numerical errors at low arg values
503//              return(scale+bkg)
504//      endif
505       
506        Iq = aa*x^-mm
507       
508        qr2=(x*rg)^2
509        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
510       
511        //scale
512        Pq *= scale
513        // then add the terms up
514        return (Iq + Pq + bkg)
515End
516
517// blocked beam
518//
519Function VC_BlockedBeam(bkg,x)
520        Variable bkg
521        Variable x
522       
523        return (bkg)
524End
525
526
527//
528// a broad peak to simulate silver behenate or vycor
529//
530// peak @ 0.1 ~ AgBeh
531//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}               
532//
533//
534// peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
535//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}             
536//
537//
538Function VC_BroadPeak(aa,nn,cc,LL,Qzero,mm,bgd,x)
539        Variable aa,nn,cc,LL,Qzero,mm,bgd
540        Variable x
541       
542        // variables are:                                                       
543        //[0] Porod term scaling
544        //[1] Porod exponent
545        //[2] Lorentzian term scaling
546        //[3] Lorentzian screening length [A]
547        //[4] peak location [1/A]
548        //[5] Lorentzian exponent
549        //[6] background
550       
551//      local variables
552        Variable inten, qval
553//      x is the q-value for the calculation
554        qval = x
555//      do the calculation and return the function value
556       
557        inten = aa/(qval)^nn + cc/(1 + (abs(qval-Qzero)*LL)^mm) + bgd
558
559        Return (inten)
560       
561End
562
563//
564// updated to new folder structure Feb 2016
565// folderStr = RAW,SAM, VCALC or other
566// detStr is the panel identifer "ML", etc.
567//
568Function SetDeltaQ(folderStr,detStr)
569        String folderStr,detStr
570
571        Variable isVCALC
572        if(cmpstr(folderStr,"VCALC") == 0)
573                isVCALC = 1
574        endif
575       
576        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
577        String instPath = ":entry:instrument:detector_"
578               
579        if(isVCALC)
580                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)               // 2D detector data
581        else
582                Wave inten = V_getDetectorDataW(folderStr,detStr)
583        endif
584
585        Wave qx = $(folderPath+instPath+detStr+":qx_"+detStr)
586        Wave qy = $(folderPath+instPath+detStr+":qy_"+detStr)
587       
588        Variable xDim,yDim,delQ
589       
590        xDim=DimSize(inten,0)
591        yDim=DimSize(inten,1)
592       
593        if(xDim<yDim)
594                delQ = abs(qx[0][0] - qx[1][0])/2
595        else
596                delQ = abs(qy[0][1] - qy[0][0])/2
597        endif
598       
599        // set the global
600        Variable/G $(folderPath+instPath+detStr+":gDelQ_"+detStr) = delQ
601//      Print "SET delQ = ",delQ," for ",type
602       
603        return(delQ)
604end
605
606
607//TODO -- need a switch here to dispatch to the averaging type
608Proc VC_BinQxQy_to_1D(folderStr,type)
609        String folderStr
610        String type
611//      Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4)
612//      Prompt type,"detector identifier"
613
614
615        VC_fDoBinning_QxQy2D(folderStr, type)
616
617
618/// this is for a tall, narrow slit mode       
619//      VC_fBinDetector_byRows(folderStr,type)
620       
621End
622
623
624// folderStr is RAW, VCALC, SAM, etc.
625// type is "B", "FL" for single binning, "FLR", or "MLRTB" or similar if multiple panels are combined
626//
627Proc VC_Graph_1D_detType(folderStr,type)
628        String folderStr,type
629       
630        SetDataFolder $("root:Packages:NIST:VSANS:"+folderStr)
631       
632        Display $("iBin_qxqy"+"_"+type) vs $("qBin_qxqy"+"_"+type)
633        ModifyGraph mirror=2,grid=1,log=1
634        ModifyGraph mode=4,marker=19,msize=2
635//      ErrorBars/T=0 iBin_qxqy Y,wave=(eBin2D_qxqy,eBin2D_qxqy)                // for simulations, I don't have 2D uncertainty
636        ErrorBars/T=0 $("iBin_qxqy"+"_"+type) Y,wave=($("eBin_qxqy"+"_"+type),$("eBin_qxqy"+"_"+type))
637        legend
638       
639        SetDataFolder root:
640
641End
642
643
644
645//////////
646//
647//              Function that bins a 2D detctor panel into I(q) based on the q-value of the pixel
648//              - each pixel QxQyQz has been calculated beforehand
649//              - if multiple panels are selected to be combined, it is done here during the binning
650//              - the setting of deltaQ step is still a little suspect (TODO)
651//
652//
653// see the equivalent function in PlotUtils2D_v40.ipf
654//
655//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
656//
657// this has been modified to accept different detector panels and to take arrays
658// -- type = FL or FR or...other panel identifiers
659//
660// TODO "iErr" is not always defined correctly since it doesn't really apply here for data that is not 2D simulation
661//
662//
663// updated Feb2016 to take new folder structure
664// TODO
665// -- VERIFY
666// -- figure out what the best location is to put the averaged data? currently @ top level of WORK folder
667//    but this is a lousy choice.
668// x- binning is now Mask-aware. If mask is not present, all data is used. If data is from VCALC, all data is used
669// -- Where do I put the solid angle correction? In here as a weight for each point, or later on as
670//    a blanket correction (matrix multiply) for an entire panel?
671//
672//
673//
674// folderStr = WORK folder, type = the binning type (may include multiple detectors)
675Function VC_fDoBinning_QxQy2D(folderStr,type)
676        String folderStr,type
677       
678        Variable nSets = 0
679        Variable xDim,yDim
680        Variable ii,jj
681        Variable qVal,nq,var,avesq,aveisq
682        Variable binIndex,val,isVCALC=0,maskMissing
683
684        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
685        String instPath = ":entry:instrument:detector_"
686        String detStr
687               
688        if(cmpstr(folderStr,"VCALC") == 0)
689                isVCALC = 1
690        endif
691       
692// now switch on the type to determine which waves to declare and create
693// since there may be more than one panel to step through. There may be two, there may be four
694//
695
696// TODO:
697// -- Solid_Angle -- waves will be present for WORK data other than RAW, but not for RAW
698//
699// assume that the mask files are missing unless we can find them. If VCALC data,
700//  then the Mask is missing by definition
701        maskMissing = 1
702
703        strswitch(type) // string switch
704                case "FL":              // execute if case matches expression
705                case "FR":
706                        detStr = type
707                        if(isVCALC)
708                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
709                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation
710                        else
711                                Wave inten = V_getDetectorDataW(folderStr,detStr)
712                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
713                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
714                                if(WaveExists(mask) == 1)
715                                        maskMissing = 0
716                                endif
717                               
718                        endif
719                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
720                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
721                        nSets = 1
722                        break   
723                                                               
724                case "FT":             
725                case "FB":
726                        detStr = type
727                        if(isVCALC)
728                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
729                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
730                        else
731                                Wave inten = V_getDetectorDataW(folderStr,detStr)
732                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
733                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
734                                if(WaveExists(mask) == 1)
735                                        maskMissing = 0
736                                endif
737                        endif
738                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
739                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
740                        nSets = 1
741                        break
742                       
743                case "ML":             
744                case "MR":
745                        detStr = type
746                        if(isVCALC)
747                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
748                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
749                        else
750                                Wave inten = V_getDetectorDataW(folderStr,detStr)
751                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
752                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
753                                if(WaveExists(mask) == 1)
754                                        maskMissing = 0
755                                endif
756                        endif   
757                        //TODO:
758                        // -- decide on the proper deltaQ for binning. either nominal value for LR, or one
759                        //    determined specifically for that panel (currently using one tube width as deltaQ)
760                        // -- this is repeated multiple times in this switch
761                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
762//                      NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
763                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
764                        nSets = 1
765                        break   
766                                       
767                case "MT":             
768                case "MB":
769                        detStr = type
770                        if(isVCALC)
771                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
772                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
773                        else
774                                Wave inten = V_getDetectorDataW(folderStr,detStr)
775                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
776                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
777                                if(WaveExists(mask) == 1)
778                                        maskMissing = 0
779                                endif
780                        endif   
781                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
782                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
783                        nSets = 1
784                        break   
785                                       
786                case "B":       
787                        detStr = type
788                        if(isVCALC)
789                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
790                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
791                        else
792                                Wave inten = V_getDetectorDataW(folderStr,detStr)
793                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
794                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
795                                if(WaveExists(mask) == 1)
796                                        maskMissing = 0
797                                endif
798                        endif   
799                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_B")
800                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values 
801                        nSets = 1
802                        break   
803                       
804                case "FLR":
805                // detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
806                // TODO
807                // -- see if I can un-hard-wire some of this below when more than one panel is combined
808                        if(isVCALC)
809                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
810                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
811                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
812                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
813                        else
814                                Wave inten = V_getDetectorDataW(folderStr,"FL")
815                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
816                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
817                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
818                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
819                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
820                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
821                                        maskMissing = 0
822                                endif
823                        endif   
824                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
825                       
826                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
827                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
828               
829                        nSets = 2
830                        break                   
831               
832                case "FTB":
833                        if(isVCALC)
834                                WAVE inten = $(folderPath+instPath+"FT"+":det_"+"FT")
835                                WAVE/Z iErr = $("iErr_"+"FT")                   // 2D errors -- may not exist, especially for simulation               
836                                WAVE inten2 = $(folderPath+instPath+"FB"+":det_"+"FB")
837                                WAVE/Z iErr2 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
838                        else
839                                Wave inten = V_getDetectorDataW(folderStr,"FT")
840                                Wave iErr = V_getDetectorDataErrW(folderStr,"FT")
841                                Wave inten2 = V_getDetectorDataW(folderStr,"FB")
842                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FB")
843                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
844                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
845                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
846                                        maskMissing = 0
847                                endif
848                        endif   
849                        NVAR delQ = $(folderPath+instPath+"FT"+":gDelQ_FT")
850                       
851                        Wave qTotal = $(folderPath+instPath+"FT"+":qTot_"+"FT")                 // 2D q-values 
852                        Wave qTotal2 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
853       
854                        nSets = 2
855                        break           
856               
857                case "FLRTB":
858                        if(isVCALC)
859                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
860                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
861                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
862                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
863                                WAVE inten3 = $(folderPath+instPath+"FT"+":det_"+"FT")
864                                WAVE/Z iErr3 = $("iErr_"+"FT")                  // 2D errors -- may not exist, especially for simulation               
865                                WAVE inten4 = $(folderPath+instPath+"FB"+":det_"+"FB")
866                                WAVE/Z iErr4 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
867                        else
868                                Wave inten = V_getDetectorDataW(folderStr,"FL")
869                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
870                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
871                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
872                                Wave inten3 = V_getDetectorDataW(folderStr,"FT")
873                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"FT")
874                                Wave inten4 = V_getDetectorDataW(folderStr,"FB")
875                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"FB")
876                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
877                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
878                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
879                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
880                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
881                                        maskMissing = 0
882                                endif
883                        endif   
884                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
885                       
886                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
887                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
888                        Wave qTotal3 = $(folderPath+instPath+"FT"+":qTot_"+"FT")                        // 2D q-values 
889                        Wave qTotal4 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
890               
891                        nSets = 4
892                        break           
893                       
894                case "MLR":
895                        if(isVCALC)
896                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
897                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
898                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
899                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
900                        else
901                                Wave inten = V_getDetectorDataW(folderStr,"ML")
902                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
903                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
904                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
905                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
906                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
907                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
908                                        maskMissing = 0
909                                endif
910                        endif   
911                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
912                       
913                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
914                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
915               
916                        nSets = 2
917                        break                   
918               
919                case "MTB":
920                        if(isVCALC)
921                                WAVE inten = $(folderPath+instPath+"MT"+":det_"+"MT")
922                                WAVE/Z iErr = $("iErr_"+"MT")                   // 2D errors -- may not exist, especially for simulation               
923                                WAVE inten2 = $(folderPath+instPath+"MB"+":det_"+"MB")
924                                WAVE/Z iErr2 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
925                        else
926                                Wave inten = V_getDetectorDataW(folderStr,"MT")
927                                Wave iErr = V_getDetectorDataErrW(folderStr,"MT")
928                                Wave inten2 = V_getDetectorDataW(folderStr,"MB")
929                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MB")
930                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
931                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
932                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
933                                        maskMissing = 0
934                                endif
935                        endif   
936                        NVAR delQ = $(folderPath+instPath+"MT"+":gDelQ_MT")
937                       
938                        Wave qTotal = $(folderPath+instPath+"MT"+":qTot_"+"MT")                 // 2D q-values 
939                        Wave qTotal2 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
940               
941                        nSets = 2
942                        break                           
943               
944                case "MLRTB":
945                        if(isVCALC)
946                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
947                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
948                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
949                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
950                                WAVE inten3 = $(folderPath+instPath+"MT"+":det_"+"MT")
951                                WAVE/Z iErr3 = $("iErr_"+"MT")                  // 2D errors -- may not exist, especially for simulation               
952                                WAVE inten4 = $(folderPath+instPath+"MB"+":det_"+"MB")
953                                WAVE/Z iErr4 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
954                        else
955                                Wave inten = V_getDetectorDataW(folderStr,"ML")
956                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
957                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
958                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
959                                Wave inten3 = V_getDetectorDataW(folderStr,"MT")
960                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"MT")
961                                Wave inten4 = V_getDetectorDataW(folderStr,"MB")
962                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"MB")
963                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
964                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
965                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
966                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
967                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
968                                        maskMissing = 0
969                                endif
970                        endif   
971                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
972                       
973                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
974                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
975                        Wave qTotal3 = $(folderPath+instPath+"MT"+":qTot_"+"MT")                        // 2D q-values 
976                        Wave qTotal4 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
977               
978                        nSets = 4
979                        break                                                                   
980                                       
981                default:
982                        nSets = 0                                                       // optional default expression executed
983                        Print "ERROR   ---- type is not recognized "
984        endswitch
985
986//      Print "delQ = ",delQ," for ",type
987
988        if(nSets == 0)
989                SetDataFolder root:
990                return(0)
991        endif
992
993
994//TODO: properly define the errors here - I'll have this if I do the simulation
995// -- need to propagate the errors up to this point
996//
997        if(WaveExists(iErr)==0  && WaveExists(inten) != 0)
998                Duplicate/O inten,iErr
999                Wave iErr=iErr
1000//              iErr = 1+sqrt(inten+0.75)                       // can't use this -- it applies to counts, not intensity (already a count rate...)
1001                iErr = sqrt(inten+0.75)                 // TODO -- here I'm just using some fictional value
1002        endif
1003        if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0)
1004                Duplicate/O inten2,iErr2
1005                Wave iErr2=iErr2
1006//              iErr2 = 1+sqrt(inten2+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1007                iErr2 = sqrt(inten2+0.75)                       // TODO -- here I'm just using some fictional value
1008        endif
1009        if(WaveExists(iErr3)==0  && WaveExists(inten3) != 0)
1010                Duplicate/O inten3,iErr3
1011                Wave iErr3=iErr3
1012//              iErr3 = 1+sqrt(inten3+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1013                iErr3 = sqrt(inten3+0.75)                       // TODO -- here I'm just using some fictional value
1014        endif
1015        if(WaveExists(iErr4)==0  && WaveExists(inten4) != 0)
1016                Duplicate/O inten4,iErr4
1017                Wave iErr4=iErr4
1018//              iErr4 = 1+sqrt(inten4+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1019                iErr4 = sqrt(inten4+0.75)                       // TODO -- here I'm just using some fictional value
1020        endif
1021
1022        nq = 600
1023
1024        // note that the back panel of 320x320 (1mm res) results in 447 data points!
1025        // - so I upped nq to 600
1026
1027//******TODO****** -- where to put the averaged data -- right now, folderStr is forced to ""   
1028//      SetDataFolder $("root:"+folderStr)              //should already be here, but make sure...     
1029        Make/O/D/N=(nq)  $(folderPath+":"+"iBin_qxqy"+"_"+type)
1030        Make/O/D/N=(nq)  $(folderPath+":"+"qBin_qxqy"+"_"+type)
1031        Make/O/D/N=(nq)  $(folderPath+":"+"nBin_qxqy"+"_"+type)
1032        Make/O/D/N=(nq)  $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1033        Make/O/D/N=(nq)  $(folderPath+":"+"eBin_qxqy"+"_"+type)
1034        Make/O/D/N=(nq)  $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1035       
1036        Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+type)
1037        Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type)
1038        Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy"+"_"+type)
1039        Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1040        Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy"+"_"+type)
1041        Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1042       
1043       
1044//      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
1045// TODO: not sure if I want to set dQ in x or y direction...
1046        // the short dimension is the 8mm tubes, use this direction as dQ?
1047        // 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]
1048        // WRT the beam center. use qx or qy directly. Still not happy with this way...
1049
1050
1051        qBin_qxqy[] =  p*delQ   
1052        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
1053
1054        iBin_qxqy = 0
1055        iBin2_qxqy = 0
1056        eBin_qxqy = 0
1057        eBin2D_qxqy = 0
1058        nBin_qxqy = 0   //number of intensities added to each bin
1059
1060// now there are situations of:
1061// 1 panel
1062// 2 panels
1063// 4 panels
1064//
1065// this needs to be a double loop now...
1066// TODO:
1067// -- the iErr wave and accumulation of error is NOT CALCULATED CORRECTLY YET
1068// -- the solid angle per pixel is not yet implemented.
1069//    it will be present for WORK data other than RAW, but not for RAW
1070
1071// if any of the masks don't exist, display the error, and proceed with the averaging, using all data
1072        if(maskMissing == 1)
1073                Print "Mask file not found for at least one detector - so all data is used"
1074        endif
1075
1076        Variable mask_val
1077// use set 1 (no number) only
1078        if(nSets >= 1)
1079                xDim=DimSize(inten,0)
1080                yDim=DimSize(inten,1)
1081       
1082                for(ii=0;ii<xDim;ii+=1)
1083                        for(jj=0;jj<yDim;jj+=1)
1084                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1085                                qVal = qTotal[ii][jj]
1086                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1087                                val = inten[ii][jj]
1088                               
1089                                if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
1090                                        mask_val = 0
1091                                else
1092                                        mask_val = mask[ii][jj]
1093                                endif
1094                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1095                                        iBin_qxqy[binIndex] += val
1096                                        iBin2_qxqy[binIndex] += val*val
1097                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1098                                        nBin_qxqy[binIndex] += 1
1099                                endif
1100                        endfor
1101                endfor
1102               
1103        endif
1104
1105// add in set 2 (set 1 already done)
1106        if(nSets >= 2)
1107                xDim=DimSize(inten2,0)
1108                yDim=DimSize(inten2,1)
1109       
1110                for(ii=0;ii<xDim;ii+=1)
1111                        for(jj=0;jj<yDim;jj+=1)
1112                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1113                                qVal = qTotal2[ii][jj]
1114                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1115                                val = inten2[ii][jj]
1116                               
1117                                if(isVCALC || maskMissing)
1118                                        mask_val = 0
1119                                else
1120                                        mask_val = mask2[ii][jj]
1121                                endif
1122                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1123                                        iBin_qxqy[binIndex] += val
1124                                        iBin2_qxqy[binIndex] += val*val
1125                                        eBin2D_qxqy[binIndex] += iErr2[ii][jj]*iErr2[ii][jj]
1126                                        nBin_qxqy[binIndex] += 1
1127                                endif
1128                        endfor
1129                endfor
1130               
1131        endif
1132
1133// add in set 3 and 4 (set 1 and 2 already done)
1134        if(nSets == 4)
1135                xDim=DimSize(inten3,0)
1136                yDim=DimSize(inten3,1)
1137       
1138                for(ii=0;ii<xDim;ii+=1)
1139                        for(jj=0;jj<yDim;jj+=1)
1140                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1141                                qVal = qTotal3[ii][jj]
1142                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1143                                val = inten3[ii][jj]
1144                               
1145                                if(isVCALC || maskMissing)
1146                                        mask_val = 0
1147                                else
1148                                        mask_val = mask3[ii][jj]
1149                                endif
1150                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1151                                        iBin_qxqy[binIndex] += val
1152                                        iBin2_qxqy[binIndex] += val*val
1153                                        eBin2D_qxqy[binIndex] += iErr3[ii][jj]*iErr3[ii][jj]
1154                                        nBin_qxqy[binIndex] += 1
1155                                endif
1156                        endfor
1157                endfor
1158               
1159               
1160                xDim=DimSize(inten4,0)
1161                yDim=DimSize(inten4,1)
1162       
1163                for(ii=0;ii<xDim;ii+=1)
1164                        for(jj=0;jj<yDim;jj+=1)
1165                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1166                                qVal = qTotal4[ii][jj]
1167                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1168                                val = inten4[ii][jj]
1169                               
1170                                if(isVCALC || maskMissing)
1171                                        mask_val = 0
1172                                else
1173                                        mask_val = mask4[ii][jj]
1174                                endif
1175                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1176                                        iBin_qxqy[binIndex] += val
1177                                        iBin2_qxqy[binIndex] += val*val
1178                                        eBin2D_qxqy[binIndex] += iErr4[ii][jj]*iErr4[ii][jj]
1179                                        nBin_qxqy[binIndex] += 1
1180                                endif
1181                        endfor
1182                endfor
1183               
1184        endif
1185
1186
1187// after looping through all of the data on the panels, calculate errors on I(q),
1188// just like in CircSectAve.ipf
1189// TODO:
1190// -- 2D Errors were NOT properly acculumated above, so this loop of calculations is NOT MEANINGFUL (yet)
1191// x- the error on the 1D intensity, is correctly calculated as the standard error of the mean.
1192        for(ii=0;ii<nq;ii+=1)
1193                if(nBin_qxqy[ii] == 0)
1194                        //no pixels in annuli, data unknown
1195                        iBin_qxqy[ii] = 0
1196                        eBin_qxqy[ii] = 1
1197                        eBin2D_qxqy[ii] = NaN
1198                else
1199                        if(nBin_qxqy[ii] <= 1)
1200                                //need more than one pixel to determine error
1201                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1202                                eBin_qxqy[ii] = 1
1203                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1204                        else
1205                                //assume that the intensity in each pixel in annuli is normally distributed about mean...
1206                                //  -- this is correctly calculating the error as the standard error of the mean, as
1207                                //    was always done for SANS as well.
1208                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1209                                avesq = iBin_qxqy[ii]^2
1210                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1211                                var = aveisq-avesq
1212                                if(var<=0)
1213                                        eBin_qxqy[ii] = 1e-6
1214                                else
1215                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1216                                endif
1217                                // and calculate as it is propagated pixel-by-pixel
1218                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1219                        endif
1220                endif
1221        endfor
1222       
1223        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo
1224       
1225        // find the last non-zero point, working backwards
1226        val=nq
1227        do
1228                val -= 1
1229        while((nBin_qxqy[val] == 0) && val > 0)
1230       
1231//      print val, nBin_qxqy[val]
1232        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1233
1234        if(val == 0)
1235                // all the points were deleted
1236                return(0)
1237        endif
1238       
1239       
1240        // since the beam center is not always on the detector, many of the low Q bins will have zero pixels
1241        // find the first non-zero point, working forwards
1242        val = -1
1243        do
1244                val += 1
1245        while(nBin_qxqy[val] == 0)     
1246        DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1247
1248        // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
1249        val = numpnts(nBin_qxqy)-1
1250        do
1251                if(nBin_qxqy[val] == 0)
1252                        DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1253                endif
1254                val -= 1
1255        while(val>0)
1256       
1257        SetDataFolder root:
1258       
1259        return(0)
1260End
1261
1262
Note: See TracBrowser for help on using the repository browser.