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

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

many additions.

Moved unused igor/nexus testing files to Vx_ prefix since they're garbage. Pulled out the useful bits for mask and div R/W and moved those to theire appropriate procedures.

Testing the simple correction of data, only tested basic subtraction. All of it still needs to be verified since I don't have any real header numbers and units yet.

Adjusted the columns on the file catalog to be more appropriate, and added a hook to allow loading of raw data files directly from the table and a popup contextual menu. May add more functionality to it later.

Corrected how the 1D data is plotted so that it correctly uses the binning type. I(q) save now also uses the binning as specified.

File size: 40.1 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        nq = 600
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//******TODO****** -- where to put the averaged data -- right now, folderStr is forced to ""   
1052//      SetDataFolder $("root:"+folderStr)              //should already be here, but make sure...     
1053        Make/O/D/N=(nq)  $(folderPath+":"+"iBin_qxqy"+"_"+type)
1054        Make/O/D/N=(nq)  $(folderPath+":"+"qBin_qxqy"+"_"+type)
1055        Make/O/D/N=(nq)  $(folderPath+":"+"nBin_qxqy"+"_"+type)
1056        Make/O/D/N=(nq)  $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1057        Make/O/D/N=(nq)  $(folderPath+":"+"eBin_qxqy"+"_"+type)
1058        Make/O/D/N=(nq)  $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1059       
1060        Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+type)
1061        Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type)
1062        Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy"+"_"+type)
1063        Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1064        Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy"+"_"+type)
1065        Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1066       
1067       
1068//      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
1069// TODO: not sure if I want to set dQ in x or y direction...
1070        // the short dimension is the 8mm tubes, use this direction as dQ?
1071        // 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]
1072        // WRT the beam center. use qx or qy directly. Still not happy with this way...
1073
1074
1075        qBin_qxqy[] =  p*delQ   
1076        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
1077
1078        iBin_qxqy = 0
1079        iBin2_qxqy = 0
1080        eBin_qxqy = 0
1081        eBin2D_qxqy = 0
1082        nBin_qxqy = 0   //number of intensities added to each bin
1083
1084// now there are situations of:
1085// 1 panel
1086// 2 panels
1087// 4 panels
1088//
1089// this needs to be a double loop now...
1090// TODO:
1091// -- the iErr wave and accumulation of error is NOT CALCULATED CORRECTLY YET
1092// -- the solid angle per pixel is not yet implemented.
1093//    it will be present for WORK data other than RAW, but not for RAW
1094
1095// if any of the masks don't exist, display the error, and proceed with the averaging, using all data
1096        if(maskMissing == 1)
1097                Print "Mask file not found for at least one detector - so all data is used"
1098        endif
1099
1100        Variable mask_val
1101// use set 1 (no number) only
1102        if(nSets >= 1)
1103                xDim=DimSize(inten,0)
1104                yDim=DimSize(inten,1)
1105       
1106                for(ii=0;ii<xDim;ii+=1)
1107                        for(jj=0;jj<yDim;jj+=1)
1108                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1109                                qVal = qTotal[ii][jj]
1110                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1111                                val = inten[ii][jj]
1112                               
1113                                if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
1114                                        mask_val = 0
1115                                else
1116                                        mask_val = mask[ii][jj]
1117                                endif
1118                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1119                                        iBin_qxqy[binIndex] += val
1120                                        iBin2_qxqy[binIndex] += val*val
1121                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1122                                        nBin_qxqy[binIndex] += 1
1123                                endif
1124                        endfor
1125                endfor
1126               
1127        endif
1128
1129// add in set 2 (set 1 already done)
1130        if(nSets >= 2)
1131                xDim=DimSize(inten2,0)
1132                yDim=DimSize(inten2,1)
1133       
1134                for(ii=0;ii<xDim;ii+=1)
1135                        for(jj=0;jj<yDim;jj+=1)
1136                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1137                                qVal = qTotal2[ii][jj]
1138                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1139                                val = inten2[ii][jj]
1140                               
1141                                if(isVCALC || maskMissing)
1142                                        mask_val = 0
1143                                else
1144                                        mask_val = mask2[ii][jj]
1145                                endif
1146                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1147                                        iBin_qxqy[binIndex] += val
1148                                        iBin2_qxqy[binIndex] += val*val
1149                                        eBin2D_qxqy[binIndex] += iErr2[ii][jj]*iErr2[ii][jj]
1150                                        nBin_qxqy[binIndex] += 1
1151                                endif
1152                        endfor
1153                endfor
1154               
1155        endif
1156
1157// add in set 3 and 4 (set 1 and 2 already done)
1158        if(nSets == 4)
1159                xDim=DimSize(inten3,0)
1160                yDim=DimSize(inten3,1)
1161       
1162                for(ii=0;ii<xDim;ii+=1)
1163                        for(jj=0;jj<yDim;jj+=1)
1164                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1165                                qVal = qTotal3[ii][jj]
1166                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1167                                val = inten3[ii][jj]
1168                               
1169                                if(isVCALC || maskMissing)
1170                                        mask_val = 0
1171                                else
1172                                        mask_val = mask3[ii][jj]
1173                                endif
1174                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1175                                        iBin_qxqy[binIndex] += val
1176                                        iBin2_qxqy[binIndex] += val*val
1177                                        eBin2D_qxqy[binIndex] += iErr3[ii][jj]*iErr3[ii][jj]
1178                                        nBin_qxqy[binIndex] += 1
1179                                endif
1180                        endfor
1181                endfor
1182               
1183               
1184                xDim=DimSize(inten4,0)
1185                yDim=DimSize(inten4,1)
1186       
1187                for(ii=0;ii<xDim;ii+=1)
1188                        for(jj=0;jj<yDim;jj+=1)
1189                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1190                                qVal = qTotal4[ii][jj]
1191                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1192                                val = inten4[ii][jj]
1193                               
1194                                if(isVCALC || maskMissing)
1195                                        mask_val = 0
1196                                else
1197                                        mask_val = mask4[ii][jj]
1198                                endif
1199                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1200                                        iBin_qxqy[binIndex] += val
1201                                        iBin2_qxqy[binIndex] += val*val
1202                                        eBin2D_qxqy[binIndex] += iErr4[ii][jj]*iErr4[ii][jj]
1203                                        nBin_qxqy[binIndex] += 1
1204                                endif
1205                        endfor
1206                endfor
1207               
1208        endif
1209
1210
1211// after looping through all of the data on the panels, calculate errors on I(q),
1212// just like in CircSectAve.ipf
1213// TODO:
1214// -- 2D Errors were NOT properly acculumated above, so this loop of calculations is NOT MEANINGFUL (yet)
1215// x- the error on the 1D intensity, is correctly calculated as the standard error of the mean.
1216        for(ii=0;ii<nq;ii+=1)
1217                if(nBin_qxqy[ii] == 0)
1218                        //no pixels in annuli, data unknown
1219                        iBin_qxqy[ii] = 0
1220                        eBin_qxqy[ii] = 1
1221                        eBin2D_qxqy[ii] = NaN
1222                else
1223                        if(nBin_qxqy[ii] <= 1)
1224                                //need more than one pixel to determine error
1225                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1226                                eBin_qxqy[ii] = 1
1227                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1228                        else
1229                                //assume that the intensity in each pixel in annuli is normally distributed about mean...
1230                                //  -- this is correctly calculating the error as the standard error of the mean, as
1231                                //    was always done for SANS as well.
1232                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1233                                avesq = iBin_qxqy[ii]^2
1234                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1235                                var = aveisq-avesq
1236                                if(var<=0)
1237                                        eBin_qxqy[ii] = 1e-6
1238                                else
1239                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1240                                endif
1241                                // and calculate as it is propagated pixel-by-pixel
1242                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1243                        endif
1244                endif
1245        endfor
1246       
1247        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo
1248       
1249        // find the last non-zero point, working backwards
1250        val=nq
1251        do
1252                val -= 1
1253        while((nBin_qxqy[val] == 0) && val > 0)
1254       
1255//      print val, nBin_qxqy[val]
1256        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1257
1258        if(val == 0)
1259                // all the points were deleted
1260                return(0)
1261        endif
1262       
1263       
1264        // since the beam center is not always on the detector, many of the low Q bins will have zero pixels
1265        // find the first non-zero point, working forwards
1266        val = -1
1267        do
1268                val += 1
1269        while(nBin_qxqy[val] == 0)     
1270        DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1271
1272        // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
1273        val = numpnts(nBin_qxqy)-1
1274        do
1275                if(nBin_qxqy[val] == 0)
1276                        DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1277                endif
1278                val -= 1
1279        while(val>0)
1280       
1281        SetDataFolder root:
1282       
1283        return(0)
1284End
1285
1286
Note: See TracBrowser for help on using the repository browser.