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

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

Significant changes to the base READ of individual data fields from data files. Now, if the field requested is from a WORK file, and it does not exist, an error condition is returned (or a null wave). Calling procedures are responsible for handling errors. This prevents a string of open file dialogs if fields are missing from a file if they were never in the file to begin with (like sensor logs, polarization hardware, etc.)

New get/write calls were added for the updated temperature sensor fields.

group_ID is now only in the sample block, not the duplicated in the reduction block, and is correctly a number not a string.

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