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

Last change on this file since 1051 was 1051, checked in by srkline, 5 years ago

a lot of little changes:

changed the name of the Raw Data display procedure file (removed test)

lots of bug fixes, moving items from the macros menu to proper locations, getting the file status to display properly, some error checking, and cleaning up a few TODO items.

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