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

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

changes to streamline the data plotting of 1D data, in preparation for different modes of combining detector panels. Also will allow better integration with protocols to combine 1D data, which can now be part of the protocol.

Other changes, but I can't remember whtat they were...

File size: 40.2 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        //
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// TODO: 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// -- 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?
695//
696//
697//
698// folderStr = WORK folder, type = the binning type (may include multiple detectors)
699Function VC_fDoBinning_QxQy2D(folderStr,type)
700        String folderStr,type
701       
702        Variable nSets = 0
703        Variable xDim,yDim
704        Variable ii,jj
705        Variable qVal,nq,var,avesq,aveisq
706        Variable binIndex,val,isVCALC=0,maskMissing
707
708        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
709        String instPath = ":entry:instrument:detector_"
710        String detStr
711               
712        if(cmpstr(folderStr,"VCALC") == 0)
713                isVCALC = 1
714        endif
715       
716// now switch on the type to determine which waves to declare and create
717// since there may be more than one panel to step through. There may be two, there may be four
718//
719
720// TODO:
721// -- Solid_Angle -- waves will be present for WORK data other than RAW, but not for RAW
722//
723// assume that the mask files are missing unless we can find them. If VCALC data,
724//  then the Mask is missing by definition
725        maskMissing = 1
726
727        strswitch(type) // string switch
728                case "FL":              // execute if case matches expression
729                case "FR":
730                        detStr = type
731                        if(isVCALC)
732                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
733                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation
734                        else
735                                Wave inten = V_getDetectorDataW(folderStr,detStr)
736                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
737                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
738                                if(WaveExists(mask) == 1)
739                                        maskMissing = 0
740                                endif
741                               
742                        endif
743                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
744                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
745                        nSets = 1
746                        break   
747                                                               
748                case "FT":             
749                case "FB":
750                        detStr = type
751                        if(isVCALC)
752                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
753                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
754                        else
755                                Wave inten = V_getDetectorDataW(folderStr,detStr)
756                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
757                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
758                                if(WaveExists(mask) == 1)
759                                        maskMissing = 0
760                                endif
761                        endif
762                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
763                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
764                        nSets = 1
765                        break
766                       
767                case "ML":             
768                case "MR":
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                        //TODO:
782                        // -- decide on the proper deltaQ for binning. either nominal value for LR, or one
783                        //    determined specifically for that panel (currently using one tube width as deltaQ)
784                        // -- this is repeated multiple times in this switch
785                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
786//                      NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
787                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
788                        nSets = 1
789                        break   
790                                       
791                case "MT":             
792                case "MB":
793                        detStr = type
794                        if(isVCALC)
795                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
796                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
797                        else
798                                Wave inten = V_getDetectorDataW(folderStr,detStr)
799                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
800                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
801                                if(WaveExists(mask) == 1)
802                                        maskMissing = 0
803                                endif
804                        endif   
805                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
806                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
807                        nSets = 1
808                        break   
809                                       
810                case "B":       
811                        detStr = type
812                        if(isVCALC)
813                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
814                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
815                        else
816                                Wave inten = V_getDetectorDataW(folderStr,detStr)
817                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
818                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
819                                if(WaveExists(mask) == 1)
820                                        maskMissing = 0
821                                endif
822                        endif   
823                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_B")
824                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values 
825                        nSets = 1
826                        break   
827                       
828                case "FLR":
829                // detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
830                // TODO
831                // -- see if I can un-hard-wire some of this below when more than one panel is combined
832                        if(isVCALC)
833                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
834                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
835                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
836                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
837                        else
838                                Wave inten = V_getDetectorDataW(folderStr,"FL")
839                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
840                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
841                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
842                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
843                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
844                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
845                                        maskMissing = 0
846                                endif
847                        endif   
848                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
849                       
850                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
851                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
852               
853                        nSets = 2
854                        break                   
855               
856                case "FTB":
857                        if(isVCALC)
858                                WAVE inten = $(folderPath+instPath+"FT"+":det_"+"FT")
859                                WAVE/Z iErr = $("iErr_"+"FT")                   // 2D errors -- may not exist, especially for simulation               
860                                WAVE inten2 = $(folderPath+instPath+"FB"+":det_"+"FB")
861                                WAVE/Z iErr2 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
862                        else
863                                Wave inten = V_getDetectorDataW(folderStr,"FT")
864                                Wave iErr = V_getDetectorDataErrW(folderStr,"FT")
865                                Wave inten2 = V_getDetectorDataW(folderStr,"FB")
866                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FB")
867                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
868                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
869                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
870                                        maskMissing = 0
871                                endif
872                        endif   
873                        NVAR delQ = $(folderPath+instPath+"FT"+":gDelQ_FT")
874                       
875                        Wave qTotal = $(folderPath+instPath+"FT"+":qTot_"+"FT")                 // 2D q-values 
876                        Wave qTotal2 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
877       
878                        nSets = 2
879                        break           
880               
881                case "FLRTB":
882                        if(isVCALC)
883                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
884                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
885                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
886                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
887                                WAVE inten3 = $(folderPath+instPath+"FT"+":det_"+"FT")
888                                WAVE/Z iErr3 = $("iErr_"+"FT")                  // 2D errors -- may not exist, especially for simulation               
889                                WAVE inten4 = $(folderPath+instPath+"FB"+":det_"+"FB")
890                                WAVE/Z iErr4 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
891                        else
892                                Wave inten = V_getDetectorDataW(folderStr,"FL")
893                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
894                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
895                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
896                                Wave inten3 = V_getDetectorDataW(folderStr,"FT")
897                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"FT")
898                                Wave inten4 = V_getDetectorDataW(folderStr,"FB")
899                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"FB")
900                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
901                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
902                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
903                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
904                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
905                                        maskMissing = 0
906                                endif
907                        endif   
908                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
909                       
910                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
911                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
912                        Wave qTotal3 = $(folderPath+instPath+"FT"+":qTot_"+"FT")                        // 2D q-values 
913                        Wave qTotal4 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
914               
915                        nSets = 4
916                        break           
917                       
918                case "MLR":
919                        if(isVCALC)
920                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
921                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
922                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
923                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
924                        else
925                                Wave inten = V_getDetectorDataW(folderStr,"ML")
926                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
927                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
928                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
929                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
930                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
931                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
932                                        maskMissing = 0
933                                endif
934                        endif   
935                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
936                       
937                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
938                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
939               
940                        nSets = 2
941                        break                   
942               
943                case "MTB":
944                        if(isVCALC)
945                                WAVE inten = $(folderPath+instPath+"MT"+":det_"+"MT")
946                                WAVE/Z iErr = $("iErr_"+"MT")                   // 2D errors -- may not exist, especially for simulation               
947                                WAVE inten2 = $(folderPath+instPath+"MB"+":det_"+"MB")
948                                WAVE/Z iErr2 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
949                        else
950                                Wave inten = V_getDetectorDataW(folderStr,"MT")
951                                Wave iErr = V_getDetectorDataErrW(folderStr,"MT")
952                                Wave inten2 = V_getDetectorDataW(folderStr,"MB")
953                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MB")
954                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
955                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
956                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
957                                        maskMissing = 0
958                                endif
959                        endif   
960                        NVAR delQ = $(folderPath+instPath+"MT"+":gDelQ_MT")
961                       
962                        Wave qTotal = $(folderPath+instPath+"MT"+":qTot_"+"MT")                 // 2D q-values 
963                        Wave qTotal2 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
964               
965                        nSets = 2
966                        break                           
967               
968                case "MLRTB":
969                        if(isVCALC)
970                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
971                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
972                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
973                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
974                                WAVE inten3 = $(folderPath+instPath+"MT"+":det_"+"MT")
975                                WAVE/Z iErr3 = $("iErr_"+"MT")                  // 2D errors -- may not exist, especially for simulation               
976                                WAVE inten4 = $(folderPath+instPath+"MB"+":det_"+"MB")
977                                WAVE/Z iErr4 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
978                        else
979                                Wave inten = V_getDetectorDataW(folderStr,"ML")
980                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
981                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
982                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
983                                Wave inten3 = V_getDetectorDataW(folderStr,"MT")
984                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"MT")
985                                Wave inten4 = V_getDetectorDataW(folderStr,"MB")
986                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"MB")
987                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
988                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
989                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
990                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
991                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
992                                        maskMissing = 0
993                                endif
994                        endif   
995                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
996                       
997                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
998                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
999                        Wave qTotal3 = $(folderPath+instPath+"MT"+":qTot_"+"MT")                        // 2D q-values 
1000                        Wave qTotal4 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
1001               
1002                        nSets = 4
1003                        break                                                                   
1004                                       
1005                default:
1006                        nSets = 0                                                       // optional default expression executed
1007                        Print "ERROR   ---- type is not recognized "
1008        endswitch
1009
1010//      Print "delQ = ",delQ," for ",type
1011
1012        if(nSets == 0)
1013                SetDataFolder root:
1014                return(0)
1015        endif
1016
1017
1018//TODO: properly define the errors here - I'll have this if I do the simulation
1019// -- need to propagate the errors up to this point
1020//
1021        if(WaveExists(iErr)==0  && WaveExists(inten) != 0)
1022                Duplicate/O inten,iErr
1023                Wave iErr=iErr
1024//              iErr = 1+sqrt(inten+0.75)                       // can't use this -- it applies to counts, not intensity (already a count rate...)
1025                iErr = sqrt(inten+0.75)                 // TODO -- here I'm just using some fictional value
1026        endif
1027        if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0)
1028                Duplicate/O inten2,iErr2
1029                Wave iErr2=iErr2
1030//              iErr2 = 1+sqrt(inten2+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1031                iErr2 = sqrt(inten2+0.75)                       // TODO -- here I'm just using some fictional value
1032        endif
1033        if(WaveExists(iErr3)==0  && WaveExists(inten3) != 0)
1034                Duplicate/O inten3,iErr3
1035                Wave iErr3=iErr3
1036//              iErr3 = 1+sqrt(inten3+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1037                iErr3 = sqrt(inten3+0.75)                       // TODO -- here I'm just using some fictional value
1038        endif
1039        if(WaveExists(iErr4)==0  && WaveExists(inten4) != 0)
1040                Duplicate/O inten4,iErr4
1041                Wave iErr4=iErr4
1042//              iErr4 = 1+sqrt(inten4+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1043                iErr4 = sqrt(inten4+0.75)                       // TODO -- here I'm just using some fictional value
1044        endif
1045
1046        // TODO -- nq will need to be larger, once the back detector is installed
1047        //
1048        // note that the back panel of 320x320 (1mm res) results in 447 data points!
1049        // - so I upped nq to 600
1050
1051        nq = 600
1052
1053//******TODO****** -- where to put the averaged data -- right now, folderStr is forced to ""   
1054//      SetDataFolder $("root:"+folderStr)              //should already be here, but make sure...     
1055        Make/O/D/N=(nq)  $(folderPath+":"+"iBin_qxqy"+"_"+type)
1056        Make/O/D/N=(nq)  $(folderPath+":"+"qBin_qxqy"+"_"+type)
1057        Make/O/D/N=(nq)  $(folderPath+":"+"nBin_qxqy"+"_"+type)
1058        Make/O/D/N=(nq)  $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1059        Make/O/D/N=(nq)  $(folderPath+":"+"eBin_qxqy"+"_"+type)
1060        Make/O/D/N=(nq)  $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1061       
1062        Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+type)
1063        Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type)
1064        Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy"+"_"+type)
1065        Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1066        Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy"+"_"+type)
1067        Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1068       
1069       
1070//      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
1071// TODO: not sure if I want to set dQ in x or y direction...
1072        // the short dimension is the 8mm tubes, use this direction as dQ?
1073        // 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]
1074        // WRT the beam center. use qx or qy directly. Still not happy with this way...
1075
1076
1077        qBin_qxqy[] =  p*delQ   
1078        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
1079
1080        iBin_qxqy = 0
1081        iBin2_qxqy = 0
1082        eBin_qxqy = 0
1083        eBin2D_qxqy = 0
1084        nBin_qxqy = 0   //number of intensities added to each bin
1085
1086// now there are situations of:
1087// 1 panel
1088// 2 panels
1089// 4 panels
1090//
1091// this needs to be a double loop now...
1092// TODO:
1093// -- the iErr wave and accumulation of error is NOT CALCULATED CORRECTLY YET
1094// -- the solid angle per pixel is not yet implemented.
1095//    it will be present for WORK data other than RAW, but not for RAW
1096
1097// if any of the masks don't exist, display the error, and proceed with the averaging, using all data
1098        if(maskMissing == 1)
1099                Print "Mask file not found for at least one detector - so all data is used"
1100        endif
1101
1102        Variable mask_val
1103// use set 1 (no number) only
1104        if(nSets >= 1)
1105                xDim=DimSize(inten,0)
1106                yDim=DimSize(inten,1)
1107       
1108                for(ii=0;ii<xDim;ii+=1)
1109                        for(jj=0;jj<yDim;jj+=1)
1110                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1111                                qVal = qTotal[ii][jj]
1112                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1113                                val = inten[ii][jj]
1114                               
1115                                if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
1116                                        mask_val = 0
1117                                else
1118                                        mask_val = mask[ii][jj]
1119                                endif
1120                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1121                                        iBin_qxqy[binIndex] += val
1122                                        iBin2_qxqy[binIndex] += val*val
1123                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1124                                        nBin_qxqy[binIndex] += 1
1125                                endif
1126                        endfor
1127                endfor
1128               
1129        endif
1130
1131// add in set 2 (set 1 already done)
1132        if(nSets >= 2)
1133                xDim=DimSize(inten2,0)
1134                yDim=DimSize(inten2,1)
1135       
1136                for(ii=0;ii<xDim;ii+=1)
1137                        for(jj=0;jj<yDim;jj+=1)
1138                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1139                                qVal = qTotal2[ii][jj]
1140                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1141                                val = inten2[ii][jj]
1142                               
1143                                if(isVCALC || maskMissing)
1144                                        mask_val = 0
1145                                else
1146                                        mask_val = mask2[ii][jj]
1147                                endif
1148                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1149                                        iBin_qxqy[binIndex] += val
1150                                        iBin2_qxqy[binIndex] += val*val
1151                                        eBin2D_qxqy[binIndex] += iErr2[ii][jj]*iErr2[ii][jj]
1152                                        nBin_qxqy[binIndex] += 1
1153                                endif
1154                        endfor
1155                endfor
1156               
1157        endif
1158
1159// add in set 3 and 4 (set 1 and 2 already done)
1160        if(nSets == 4)
1161                xDim=DimSize(inten3,0)
1162                yDim=DimSize(inten3,1)
1163       
1164                for(ii=0;ii<xDim;ii+=1)
1165                        for(jj=0;jj<yDim;jj+=1)
1166                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1167                                qVal = qTotal3[ii][jj]
1168                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1169                                val = inten3[ii][jj]
1170                               
1171                                if(isVCALC || maskMissing)
1172                                        mask_val = 0
1173                                else
1174                                        mask_val = mask3[ii][jj]
1175                                endif
1176                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1177                                        iBin_qxqy[binIndex] += val
1178                                        iBin2_qxqy[binIndex] += val*val
1179                                        eBin2D_qxqy[binIndex] += iErr3[ii][jj]*iErr3[ii][jj]
1180                                        nBin_qxqy[binIndex] += 1
1181                                endif
1182                        endfor
1183                endfor
1184               
1185               
1186                xDim=DimSize(inten4,0)
1187                yDim=DimSize(inten4,1)
1188       
1189                for(ii=0;ii<xDim;ii+=1)
1190                        for(jj=0;jj<yDim;jj+=1)
1191                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1192                                qVal = qTotal4[ii][jj]
1193                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1194                                val = inten4[ii][jj]
1195                               
1196                                if(isVCALC || maskMissing)
1197                                        mask_val = 0
1198                                else
1199                                        mask_val = mask4[ii][jj]
1200                                endif
1201                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1202                                        iBin_qxqy[binIndex] += val
1203                                        iBin2_qxqy[binIndex] += val*val
1204                                        eBin2D_qxqy[binIndex] += iErr4[ii][jj]*iErr4[ii][jj]
1205                                        nBin_qxqy[binIndex] += 1
1206                                endif
1207                        endfor
1208                endfor
1209               
1210        endif
1211
1212
1213// after looping through all of the data on the panels, calculate errors on I(q),
1214// just like in CircSectAve.ipf
1215// TODO:
1216// -- 2D Errors were NOT properly acculumated above, so this loop of calculations is NOT MEANINGFUL (yet)
1217// x- the error on the 1D intensity, is correctly calculated as the standard error of the mean.
1218        for(ii=0;ii<nq;ii+=1)
1219                if(nBin_qxqy[ii] == 0)
1220                        //no pixels in annuli, data unknown
1221                        iBin_qxqy[ii] = 0
1222                        eBin_qxqy[ii] = 1
1223                        eBin2D_qxqy[ii] = NaN
1224                else
1225                        if(nBin_qxqy[ii] <= 1)
1226                                //need more than one pixel to determine error
1227                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1228                                eBin_qxqy[ii] = 1
1229                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1230                        else
1231                                //assume that the intensity in each pixel in annuli is normally distributed about mean...
1232                                //  -- this is correctly calculating the error as the standard error of the mean, as
1233                                //    was always done for SANS as well.
1234                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1235                                avesq = iBin_qxqy[ii]^2
1236                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1237                                var = aveisq-avesq
1238                                if(var<=0)
1239                                        eBin_qxqy[ii] = 1e-6
1240                                else
1241                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1242                                endif
1243                                // and calculate as it is propagated pixel-by-pixel
1244                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1245                        endif
1246                endif
1247        endfor
1248       
1249        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo
1250       
1251        // find the last non-zero point, working backwards
1252        val=nq
1253        do
1254                val -= 1
1255        while((nBin_qxqy[val] == 0) && val > 0)
1256       
1257//      print val, nBin_qxqy[val]
1258        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1259
1260        if(val == 0)
1261                // all the points were deleted
1262                return(0)
1263        endif
1264       
1265       
1266        // since the beam center is not always on the detector, many of the low Q bins will have zero pixels
1267        // find the first non-zero point, working forwards
1268        val = -1
1269        do
1270                val += 1
1271        while(nBin_qxqy[val] == 0)     
1272        DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1273
1274        // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
1275        val = numpnts(nBin_qxqy)-1
1276        do
1277                if(nBin_qxqy[val] == 0)
1278                        DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1279                endif
1280                val -= 1
1281        while(val>0)
1282       
1283        SetDataFolder root:
1284       
1285        return(0)
1286End
1287
1288
Note: See TracBrowser for help on using the repository browser.