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

Last change on this file since 1108 was 1108, checked in by srkline, 4 years ago

changes mostly to VCALC to add in what bits of information I have about the instrument dimensions. Added in stubs (based on NG3 SANS) for the beam intensity. Added in a preset condition for Front+Middle. Still need a more uniform way to do this.

File size: 51.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// x- hard wired for a sphere - change this to allow minimal selections and altering of coefficients
16// x- add the "fake" 2D simulation to fill the panels which are then later averaged as I(Q)
17//
18// NOTE - this is a VCALC only routine, so it's not been made completely generic
19//
20Function FillPanel_wModelData(det,qTot,type)
21        Wave det,qTot
22        String type
23
24//      SetDataFolder root:Packages:NIST:VSANS:VCALC:Front
25
26        // q-values and detector arrays already allocated and calculated
27        Duplicate/O det tmpInten,tmpSig,prob_i
28               
29        Variable imon,trans,thick,sdd,pixSizeX,pixSizeY,sdd_offset
30
31        //imon = VC_BeamIntensity()*CountTime
32        imon = VCALC_getImon()          //TODO: currently from the panel, not calculated
33        trans = 0.8
34        thick = 0.1
35       
36        // need SDD
37        // need pixel dimensions
38        // nominal sdd in cm, offset in cm, want result in cm !
39        sdd = VCALC_getSDD(type)        +  VCALC_getTopBottomSDDSetback(type)           // 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,2000,1e-6,0.01,qTot[p][q])   
64                        break
65                case "Debye":
66                        tmpInten = VC_Debye(10,300,0.0001,qTot[p][q])
67                        break
68                case "Sphere":
69                        tmpInten = VC_SphereForm(1,60,1e-6,0.001,qTot[p][q])   
70                        break
71                case "AgBeh":
72                        tmpInten = VC_BroadPeak(1e-11,3,20,100.0,0.1,3,0.1,qTot[p][q])
73                        break
74                case "Vycor":
75                        tmpInten = VC_BroadPeak(1e-9,3,20,500.0,0.015,3,0.1,qTot[p][q])
76                        break   
77                case "Empty Cell":
78                        tmpInten = VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
79                        break
80                case "Blocked Beam":
81                        tmpInten = VC_BlockedBeam(0.01,qTot[p][q])
82                        break
83                case "Debye +":
84                        tmpInten = VC_Debye(10,300,0.0001,qTot[p][q])
85                        addEmpBgd = 1
86                        break
87                case "AgBeh +":
88                        tmpInten = VC_BroadPeak(1e-11,3,20,100.0,0.1,3,0.1,qTot[p][q])
89                        addEmpBgd = 1
90                        break
91                case "Empty Cell +":
92                        tmpInten = VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
93                        tmpInten += VC_BlockedBeam(0.01,qTot[p][q])
94                        break
95                default:
96                        tmpInten = VC_Debye(10,300,0.1,qTot[p][q])
97        endswitch
98
99
100        if(addEmpBgd == 1)
101                tmpInten += VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
102                tmpInten += VC_BlockedBeam(0.01,qTot[p][q])
103        endif
104
105       
106// x- this is faked to get around the singularity at the center of the back detector
107//
108//
109        if(cmpstr(type,"B") == 0)
110                Variable nx,ny,px,py
111                nx = VCALC_get_nPix_X(type)
112                ny = VCALC_get_nPix_Y(type)
113                px = trunc(nx/2)
114                py = trunc(ny/2)
115               
116                tmpInten[px][py] = (tmpInten[px][py+1] + tmpInten[px][py-1])/2
117        endif
118
119
120
121///////////////
122//      // calculate the scattering cross section simply to be able to estimate the transmission
123//      Variable sig_sas=0
124//     
125//      // remember that the random deviate is the coherent portion ONLY - the incoherent background is
126//      // subtracted before the calculation.
127//      CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
128//
129//      if(sig_sas > 100)
130//              DoAlert 0,"SAS cross section > 100. Estimates of multiple scattering are unreliable. Choosing a model with a well-defined Rg may help"
131//      endif           
132//
133//      // calculate the multiple scattering fraction for display (10/2009)
134//      Variable ii,nMax=10,tau
135//      mScat=0
136//      tau = thick*sig_sas
137//      // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
138//      // neutrons out of those that were scattered
139//      for(ii=2;ii<nMax;ii+=1)
140//              mScat += tau^(ii)/factorial(ii)
141////            print tau^(ii)/factorial(ii)
142//      endfor
143//      estTrans = exp(-1*thick*sig_sas)                //thickness and sigma both in units of cm
144//      mscat *= (estTrans)/(1-estTrans)
145//
146////    if(mScat > 0.1)         //  Display warning
147//
148//      Print "Sig_sas = ",sig_sas
149////////////////////
150       
151        prob_i = trans*thick*pixSizeX*pixSizeY/(sdd)^2*tmpInten                 //probability of a neutron in q-bin(i)
152               
153        tmpInten = (imon)*prob_i                //tmpInten is not the model calculation anymore!!
154
155
156/// **** can I safely assume a Gaussian error in the count rate??
157        tmpSig = sqrt(tmpInten)         // corrected based on John's memo, from 8/9/99
158
159        tmpInten += gnoise(tmpSig)
160        tmpInten = (tmpInten[p][q] < 0) ? 0 : tmpInten[p][q]                    // MAR 2013 -- is this the right thing to do
161        tmpInten = trunc(tmpInten)
162               
163       
164        det = tmpInten
165
166// if I want "absolute" scale -- then I lose the integer nature of the detector (but keep the random)
167//      det /= trans*thick*pixSizeX*pixSizeY/(sdd)^2*imon
168
169       
170        KillWaves/Z tmpInten,tmpSig,prob_i     
171        SetDataFolder root:
172
173        return(0)
174End
175
176
177// For a given detector panel, calculate the q-values
178// -work with everything as arrays
179// Input needed:
180// detector data
181// detector type (LRTB?)
182// beam center (may be off the detector)
183// SDD
184// lambda
185//
186// pixel dimensions for detector type (global constants)
187// - data dimensions read directly from array
188//
189// --What is calculated:
190// array of Q
191// array of qx,qy,qz
192// array of error already exists
193//
194//
195// -- sdd in cm
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 call the VSANS V_CalcQval routines (these use the real-space distance, not pixel dims)
218//
219// ***** everything passed in is [cm], except for wavelength [A]
220//
221// ****  TODO :: calibration constants are still in [mm]
222//
223//
224// TODO:
225// -- tube width is hard-wired in
226//
227//
228Function VC_Detector_2Q_NonLin(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY,detStr)
229        Wave data,qTot,qx,qy,qz
230        Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
231        String detStr
232       
233        String destPath = "root:Packages:NIST:VSANS:VCALC"
234       
235        // be sure that the real distance waves exist
236        // TODO -- this may not be the best location?
237
238// calibration waves do not exist yet, so make some fake ones   '
239        // do I count on the orientation as an input, or do I just figure it out on my own?
240        String orientation
241        Variable dimX,dimY
242        dimX = DimSize(data,0)
243        dimY = DimSize(data,1)
244        if(dimX > dimY)
245                orientation = "horizontal"
246        else
247                orientation = "vertical"
248        endif
249       
250        if(cmpstr(orientation,"vertical")==0)
251                Make/O/D/N=(3,48) tmpCalib
252                // for the "tall" L/R banks
253                tmpCalib[0][] = -512
254                tmpCalib[1][] = 8
255                tmpCalib[2][] = 0
256        else
257                Make/O/D/N=(3,48) tmpCalib
258                // for the "short" T/B banks
259                tmpCalib[0][] = -256
260                tmpCalib[1][] = 4
261                tmpCalib[2][] = 0
262        endif
263        // override if back panel
264        if(cmpstr(detStr,"B") == 0)
265                // and for the back detector "B"
266                Make/O/D/N=3 tmpCalibX,tmpCalibY
267                tmpCalibX[0] = VCALC_getPixSizeX(detStr)                        // pixel size in [cm]  VCALC_getPixSizeX(detStr) is [cm]
268                tmpCalibX[1] = 1
269                tmpcalibX[2] = 10000
270                tmpCalibY[0] = VCALC_getPixSizeY(detStr)                        // pixel size in [cm]  VCALC_getPixSizeX(detStr) is [cm]
271                tmpCalibY[1] = 1
272                tmpcalibY[2] = 10000
273        endif
274       
275//      Wave w_calib = V_getDetTube_spatialCalib("VCALC",detStr)
276        Variable tube_width = 8.4                       // TODO: UNITS!!! Hard-wired value in [mm]
277        if(cmpstr(detStr,"B") == 0)
278                V_NonLinearCorrection_B("VCALC",data,tmpCalibX,tmpCalibY,detStr,destPath)
279                // beam center is in pixels, so use the old routine
280                V_ConvertBeamCtr_to_mmB("VCALC","B",destPath)
281        else
282                V_NonLinearCorrection("VCALC",data,tmpCalib,tube_width,detStr,destPath)
283        endif
284                               
285        Wave/Z data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
286        Wave/Z data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
287        NVAR gUseNonLinearDet = root:Packages:NIST:VSANS:VCALC:gUseNonLinearDet
288
289        if(kBCTR_CM && cmpstr(detStr,"B") != 0)
290                if(gUseNonLinearDet && WaveExists(data_realDistX) && WaveExists(data_realDistY))
291                        // beam ctr is in cm already
292
293                        // calculate all of the q-values
294                        qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
295                        qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
296                        qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
297                        qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
298               
299        //              Print "det, x_mm, y_mm ",detStr,num2str(newX),num2str(newY)
300        //              Print "det, x_pix, y_pix ",detStr,num2str(xCtr),num2str(yCtr)
301                else
302                        // do the q-calculation using linear detector
303                        //VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
304                        qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
305                        qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
306                        qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
307                        qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
308                endif   
309       
310       
311        else
312        // using the old calculation with beam center in pixels
313                if(gUseNonLinearDet && WaveExists(data_realDistX) && WaveExists(data_realDistY))
314                        // convert the beam centers to mm
315//                      String orientation
316                        Variable newX,newY
317                        dimX = DimSize(data_realDistX,0)
318                        dimY = DimSize(data_realDistX,1)
319                        if(dimX > dimY)
320                                orientation = "horizontal"
321                        else
322                                orientation = "vertical"
323                        endif
324                       
325               
326                //
327                        if(cmpstr(orientation,"vertical")==0)
328                                //      this is data dimensioned as (Ntubes,Npix)
329                                newX = tube_width*xCtr
330                                newY = data_realDistY[0][yCtr]
331                        else
332                                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
333                                newX = data_realDistX[xCtr][0]
334                                newY = tube_width*yCtr
335                        endif   
336       
337                        //if detector "B", different calculation for the centers (not tubes)
338                        if(cmpstr(detStr,"B")==0)
339                                newX = data_realDistX[xCtr][0]
340                                newY = data_realDistY[0][yCtr]
341                                //newX = xCtr
342                                //newY = yCtr
343                        endif           
344                                       
345                        // calculate all of the q-values
346                        qTot = V_CalcQval(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
347                        qx = V_CalcQX(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
348                        qy = V_CalcQY(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
349                        qz = V_CalcQZ(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
350               
351        //              Print "det, x_mm, y_mm ",detStr,num2str(newX),num2str(newY)
352        //              Print "det, x_pix, y_pix ",detStr,num2str(xCtr),num2str(yCtr)
353                else
354                        // do the q-calculation using linear detector
355                        VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
356                endif   
357       
358        endif
359       
360        KillWaves/Z tmpCalib,tmpCalibX,tmpCalibY
361       
362        return(0)
363End
364
365
366//////////////////////
367// NOTE: The Q calculations are different than what is in GaussUtils in that they take into
368// accout the different x/y pixel sizes and the beam center not being on the detector -
369// off a different edge for each LRTB type
370/////////////////////
371
372//function to calculate the overall q-value, given all of the necesary trig inputs
373//and are in detector coordinates (1,128) rather than axis values
374//the pixel locations need not be integers, reals are ok inputs
375//sdd is in [cm]
376//wavelength is in Angstroms
377//
378//returned magnitude of Q is in 1/Angstroms
379//
380//
381Function VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
382        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
383       
384        Variable dx,dy,qval,two_theta,dist
385               
386
387        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
388        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
389        dist = sqrt(dx^2 + dy^2)
390       
391        two_theta = atan(dist/sdd)
392
393        qval = 4*Pi/lam*sin(two_theta/2)
394       
395        return qval
396End
397
398//calculates just the q-value in the x-direction on the detector
399//input/output is the same as CalcQval()
400//ALL inputs are in detector coordinates
401//
402//sdd is in [cm]
403//wavelength is in Angstroms
404//
405// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
406// now properly accounts for qz
407//
408Function VC_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
409        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
410
411        Variable qx,qval,phi,dx,dy,dist,two_theta
412       
413        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
414       
415//      sdd *=100               //convert to cm
416        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
417        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
418        phi = V_FindPhi(dx,dy)
419       
420        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
421        dist = sqrt(dx^2 + dy^2)
422        two_theta = atan(dist/sdd)
423
424        qx = qval*cos(two_theta/2)*cos(phi)
425       
426        return qx
427End
428
429//calculates just the q-value in the y-direction on the detector
430//input/output is the same as CalcQval()
431//ALL inputs are in detector coordinates
432//sdd is in [cm]
433//wavelength is in Angstroms
434//
435// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
436// now properly accounts for qz
437//
438Function VC_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
439        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
440       
441        Variable dy,qval,dx,phi,qy,dist,two_theta
442       
443        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
444       
445//      sdd *=100               //convert to cm
446        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
447        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
448        phi = V_FindPhi(dx,dy)
449       
450        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
451        dist = sqrt(dx^2 + dy^2)
452        two_theta = atan(dist/sdd)
453       
454        qy = qval*cos(two_theta/2)*sin(phi)
455       
456        return qy
457End
458
459//calculates just the z-component of the q-vector, not measured on the detector
460//input/output is the same as CalcQval()
461//ALL inputs are in detector coordinates
462//sdd is in [cm]
463//wavelength is in Angstroms
464//
465// not actually used, but here for completeness if anyone asks
466//
467Function VC_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
468        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
469       
470        Variable dy,qval,dx,phi,qz,dist,two_theta
471       
472        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
473       
474//      sdd *=100               //convert to cm
475       
476        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
477        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
478        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
479        dist = sqrt(dx^2 + dy^2)
480        two_theta = atan(dist/sdd)
481       
482        qz = qval*sin(two_theta/2)
483       
484        return qz
485End
486
487//phi is defined from +x axis, proceeding CCW around [0,2Pi]
488Threadsafe Function V_FindPhi(vx,vy)
489        variable vx,vy
490       
491        variable phi
492       
493        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
494       
495        // special cases
496        if(vx==0 && vy > 0)
497                return(pi/2)
498        endif
499        if(vx==0 && vy < 0)
500                return(3*pi/2)
501        endif
502        if(vx >= 0 && vy == 0)
503                return(0)
504        endif
505        if(vx < 0 && vy == 0)
506                return(pi)
507        endif
508       
509       
510        if(vx > 0 && vy > 0)
511                return(phi)
512        endif
513        if(vx < 0 && vy > 0)
514                return(phi + pi)
515        endif
516        if(vx < 0 && vy < 0)
517                return(phi + pi)
518        endif
519        if( vx > 0 && vy < 0)
520                return(phi + 2*pi)
521        endif
522       
523        return(phi)
524end
525
526Function VC_SphereForm(scale,radius,delrho,bkg,x)                               
527        Variable scale,radius,delrho,bkg
528        Variable x
529       
530        // variables are:                                                       
531        //[0] scale
532        //[1] radius (A)
533        //[2] delrho (A-2)
534        //[3] background (cm-1)
535       
536//      Variable scale,radius,delrho,bkg                               
537//      scale = w[0]
538//      radius = w[1]
539//      delrho = w[2]
540//      bkg = w[3]
541       
542       
543        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
544        // and is rescaled to give [=] cm^-1
545       
546        Variable bes,f,vol,f2
547        //
548        //handle q==0 separately
549        If(x==0)
550                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
551                return(f)
552        Endif
553       
554//      bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
555       
556        bes = 3*sqrt(pi/(2*x*radius))*BesselJ(1.5,x*radius)/(x*radius)
557       
558        vol = 4*pi/3*radius^3
559        f = vol*bes*delrho              // [=] A
560        // normalize to single particle volume, convert to 1/cm
561        f2 = f * f / vol * 1.0e8                // [=] 1/cm
562       
563        return (scale*f2+bkg)   // Scale, then add in the background
564       
565End
566
567Function VC_Debye(scale,rg,bkg,x)
568        Variable scale,rg,bkg
569        Variable x
570       
571        // variables are:
572        //[0] scale factor
573        //[1] radius of gyration [A]
574        //[2] background        [cm-1]
575       
576        // calculates (scale*debye)+bkg
577        Variable Pq,qr2
578       
579        qr2=(x*rg)^2
580        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
581       
582        //scale
583        Pq *= scale
584        // then add in the background
585        return (Pq+bkg)
586End
587
588// a sum of a power law and debye to approximate the scattering from a real empty cell
589//
590//      make/O/D coef_ECEmp = {2.2e-8,3.346,0.0065,9.0,0.016}
591//
592Function VC_EC_Empirical(aa,mm,scale,rg,bkg,x)
593        Variable aa,mm,scale,rg,bkg
594        Variable x
595       
596        // variables are:
597        //[0] = A
598        //[1] = power m
599        //[2] scale factor
600        //[3] radius of gyration [A]
601        //[4] background        [cm-1]
602       
603        Variable Iq
604       
605        // calculates (scale*debye)+bkg
606        Variable Pq,qr2
607       
608//      if(x*Rg < 1e-3)         //added Oct 2008 to avoid numerical errors at low arg values
609//              return(scale+bkg)
610//      endif
611       
612        Iq = aa*x^-mm
613       
614        qr2=(x*rg)^2
615        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
616       
617        //scale
618        Pq *= scale
619        // then add the terms up
620        return (Iq + Pq + bkg)
621End
622
623// blocked beam
624//
625Function VC_BlockedBeam(bkg,x)
626        Variable bkg
627        Variable x
628       
629        return (bkg)
630End
631
632
633//
634// a broad peak to simulate silver behenate or vycor
635//
636// peak @ 0.1 ~ AgBeh
637//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}               
638//
639//
640// peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
641//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}             
642//
643//
644Function VC_BroadPeak(aa,nn,cc,LL,Qzero,mm,bgd,x)
645        Variable aa,nn,cc,LL,Qzero,mm,bgd
646        Variable x
647       
648        // variables are:                                                       
649        //[0] Porod term scaling
650        //[1] Porod exponent
651        //[2] Lorentzian term scaling
652        //[3] Lorentzian screening length [A]
653        //[4] peak location [1/A]
654        //[5] Lorentzian exponent
655        //[6] background
656       
657//      local variables
658        Variable inten, qval
659//      x is the q-value for the calculation
660        qval = x
661//      do the calculation and return the function value
662       
663        inten = aa/(qval)^nn + cc/(1 + (abs(qval-Qzero)*LL)^mm) + bgd
664
665        Return (inten)
666       
667End
668
669//
670// updated to new folder structure Feb 2016
671// folderStr = RAW,SAM, VCALC or other
672// detStr is the panel identifer "ML", etc.
673//
674Function V_SetDeltaQ(folderStr,detStr)
675        String folderStr,detStr
676
677        Variable isVCALC
678        if(cmpstr(folderStr,"VCALC") == 0)
679                isVCALC = 1
680        endif
681       
682        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
683        String instPath = ":entry:instrument:detector_"
684               
685        if(isVCALC)
686                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)               // 2D detector data
687        else
688                Wave inten = V_getDetectorDataW(folderStr,detStr)
689        endif
690
691        Wave qx = $(folderPath+instPath+detStr+":qx_"+detStr)
692        Wave qy = $(folderPath+instPath+detStr+":qy_"+detStr)
693       
694        Variable xDim,yDim,delQ
695       
696        xDim=DimSize(inten,0)
697        yDim=DimSize(inten,1)
698       
699        if(xDim<yDim)
700                delQ = abs(qx[0][0] - qx[1][0])/2
701        else
702                delQ = abs(qy[0][1] - qy[0][0])/2
703        endif
704       
705        // set the global
706        Variable/G $(folderPath+instPath+detStr+":gDelQ_"+detStr) = delQ
707//      Print "SET delQ = ",delQ," for ",type
708       
709        return(delQ)
710end
711
712
713//TODO -- need a switch here to dispatch to the averaging type
714Proc VC_BinQxQy_to_1D(folderStr,type)
715        String folderStr
716        String type
717//      Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4)
718//      Prompt type,"detector identifier"
719
720
721        VC_fDoBinning_QxQy2D(folderStr, type)
722
723
724/// this is for a tall, narrow slit mode       
725//      VC_fBinDetector_byRows(folderStr,type)
726       
727End
728
729
730// folderStr is RAW, VCALC, SAM, etc.
731// type is "B", "FL" for single binning, "FLR", or "MLRTB" or similar if multiple panels are combined
732//
733Proc VC_Graph_1D_detType(folderStr,type)
734        String folderStr,type
735       
736        SetDataFolder $("root:Packages:NIST:VSANS:"+folderStr)
737       
738        Display $("iBin_qxqy"+"_"+type) vs $("qBin_qxqy"+"_"+type)
739        ModifyGraph mirror=2,grid=1,log=1
740        ModifyGraph mode=4,marker=19,msize=2
741//      ErrorBars/T=0 iBin_qxqy Y,wave=(eBin2D_qxqy,eBin2D_qxqy)                // for simulations, I don't have 2D uncertainty
742        ErrorBars/T=0 $("iBin_qxqy"+"_"+type) Y,wave=($("eBin_qxqy"+"_"+type),$("eBin_qxqy"+"_"+type))
743        legend
744       
745        SetDataFolder root:
746
747End
748
749
750
751//////////
752//
753//              Function that bins a 2D detctor panel into I(q) based on the q-value of the pixel
754//              - each pixel QxQyQz has been calculated beforehand
755//              - if multiple panels are selected to be combined, it is done here during the binning
756//              - the setting of deltaQ step is still a little suspect (TODO)
757//
758//
759// see the equivalent function in PlotUtils2D_v40.ipf
760//
761//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
762//
763// this has been modified to accept different detector panels and to take arrays
764// -- type = FL or FR or...other panel identifiers
765//
766// TODO "iErr" is not always defined correctly since it doesn't really apply here for data that is not 2D simulation
767//
768//
769// updated Feb2016 to take new folder structure
770// TODO
771// -- VERIFY
772// -- figure out what the best location is to put the averaged data? currently @ top level of WORK folder
773//    but this is a lousy choice.
774// x- binning is now Mask-aware. If mask is not present, all data is used. If data is from VCALC, all data is used
775// x- Where do I put the solid angle correction? In here as a weight for each point, or later on as
776//    a blanket correction (matrix multiply) for an entire panel? (Solid Angle correction is done in the
777//    step where data is added to a WORK file (see Raw_to_Work())
778//
779//
780// TODO:
781// -- some of the input parameters for the resolution calcuation are either assumed (apOff) or are currently
782//    hard-wired. these need to be corrected before even the pinhole resolution is correct
783// x- resolution calculation is in the correct place. The calculation is done per-panel (specified by TYPE),
784//    and then the unwanted points can be discarded (all 6 columns) as the data is trimmed and concatenated
785//    is separate functions that are resolution-aware.
786//
787//
788// folderStr = WORK folder, type = the binning type (may include multiple detectors)
789Function VC_fDoBinning_QxQy2D(folderStr,type,collimationStr)
790        String folderStr,type,collimationStr
791       
792        Variable nSets = 0
793        Variable xDim,yDim
794        Variable ii,jj
795        Variable qVal,nq,var,avesq,aveisq
796        Variable binIndex,val,isVCALC=0,maskMissing
797
798        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
799        String instPath = ":entry:instrument:detector_"
800        String detStr
801               
802        if(cmpstr(folderStr,"VCALC") == 0)
803                isVCALC = 1
804        endif
805       
806        detStr = type
807
808// now switch on the type to determine which waves to declare and create
809// since there may be more than one panel to step through. There may be two, there may be four
810//
811
812// TODO:
813// -- Solid_Angle -- waves will be present for WORK data other than RAW, but not for RAW
814//
815// assume that the mask files are missing unless we can find them. If VCALC data,
816//  then the Mask is missing by definition
817        maskMissing = 1
818
819        strswitch(type) // string switch
820//              case "FL":              // execute if case matches expression
821//              case "FR":
822//                      detStr = type
823//                      if(isVCALC)
824//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
825//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation
826//                      else
827//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
828//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
829//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
830//                              if(WaveExists(mask) == 1)
831//                                      maskMissing = 0
832//                              endif
833//                             
834//                      endif
835//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
836//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
837//                      nSets = 1
838//                      break   
839                                                               
840//              case "FT":             
841//              case "FB":
842//                      detStr = type
843//                      if(isVCALC)
844//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
845//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
846//                      else
847//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
848//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
849//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
850//                              if(WaveExists(mask) == 1)
851//                                      maskMissing = 0
852//                              endif
853//                      endif
854//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
855//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
856//                      nSets = 1
857//                      break
858                       
859//              case "ML":             
860//              case "MR":
861//                      detStr = type
862//                      if(isVCALC)
863//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
864//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
865//                      else
866//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
867//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
868//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
869//                              if(WaveExists(mask) == 1)
870//                                      maskMissing = 0
871//                              endif
872//                      endif   
873//                      //TODO:
874//                      // -- decide on the proper deltaQ for binning. either nominal value for LR, or one
875//                      //    determined specifically for that panel (currently using one tube width as deltaQ)
876//                      // -- this is repeated multiple times in this switch
877//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
878//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
879//                      nSets = 1
880//                      break   
881                                       
882//              case "MT":             
883//              case "MB":
884//                      detStr = type
885//                      if(isVCALC)
886//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
887//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
888//                      else
889//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
890//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
891//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
892//                              if(WaveExists(mask) == 1)
893//                                      maskMissing = 0
894//                              endif
895//                      endif   
896//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
897//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
898//                      nSets = 1
899//                      break   
900
901// only one panel, simply pick that panel and move on out of the switch
902                case "FL":
903                case "FR":
904                case "FT":
905                case "FB":
906                case "ML":
907                case "MR":
908                case "MT":
909                case "MB":                     
910                case "B":       
911                        if(isVCALC)
912                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
913                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
914                        else
915                                Wave inten = V_getDetectorDataW(folderStr,detStr)
916                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
917                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
918                                if(WaveExists(mask) == 1)
919                                        maskMissing = 0
920                                endif
921                        endif   
922                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
923                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values 
924                        nSets = 1
925                        break   
926                       
927                case "FLR":
928                // detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
929                // TODO
930                // -- see if I can un-hard-wire some of this below when more than one panel is combined
931                        if(isVCALC)
932                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
933                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
934                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
935                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
936                        else
937                                Wave inten = V_getDetectorDataW(folderStr,"FL")
938                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
939                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
940                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
941                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
942                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
943                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
944                                        maskMissing = 0
945                                endif
946                        endif   
947                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
948                       
949                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
950                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
951               
952                        nSets = 2
953                        break                   
954               
955                case "FTB":
956                        if(isVCALC)
957                                WAVE inten = $(folderPath+instPath+"FT"+":det_"+"FT")
958                                WAVE/Z iErr = $("iErr_"+"FT")                   // 2D errors -- may not exist, especially for simulation               
959                                WAVE inten2 = $(folderPath+instPath+"FB"+":det_"+"FB")
960                                WAVE/Z iErr2 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
961                        else
962                                Wave inten = V_getDetectorDataW(folderStr,"FT")
963                                Wave iErr = V_getDetectorDataErrW(folderStr,"FT")
964                                Wave inten2 = V_getDetectorDataW(folderStr,"FB")
965                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FB")
966                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
967                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
968                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
969                                        maskMissing = 0
970                                endif
971                        endif   
972                        NVAR delQ = $(folderPath+instPath+"FT"+":gDelQ_FT")
973                       
974                        Wave qTotal = $(folderPath+instPath+"FT"+":qTot_"+"FT")                 // 2D q-values 
975                        Wave qTotal2 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
976       
977                        nSets = 2
978                        break           
979               
980                case "FLRTB":
981                        if(isVCALC)
982                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
983                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
984                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
985                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
986                                WAVE inten3 = $(folderPath+instPath+"FT"+":det_"+"FT")
987                                WAVE/Z iErr3 = $("iErr_"+"FT")                  // 2D errors -- may not exist, especially for simulation               
988                                WAVE inten4 = $(folderPath+instPath+"FB"+":det_"+"FB")
989                                WAVE/Z iErr4 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
990                        else
991                                Wave inten = V_getDetectorDataW(folderStr,"FL")
992                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
993                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
994                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
995                                Wave inten3 = V_getDetectorDataW(folderStr,"FT")
996                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"FT")
997                                Wave inten4 = V_getDetectorDataW(folderStr,"FB")
998                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"FB")
999                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
1000                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
1001                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
1002                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
1003                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
1004                                        maskMissing = 0
1005                                endif
1006                        endif   
1007                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
1008                       
1009                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
1010                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
1011                        Wave qTotal3 = $(folderPath+instPath+"FT"+":qTot_"+"FT")                        // 2D q-values 
1012                        Wave qTotal4 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
1013               
1014                        nSets = 4
1015                        break           
1016                       
1017                case "MLR":
1018                        if(isVCALC)
1019                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
1020                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
1021                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
1022                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
1023                        else
1024                                Wave inten = V_getDetectorDataW(folderStr,"ML")
1025                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
1026                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
1027                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
1028                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
1029                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
1030                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1031                                        maskMissing = 0
1032                                endif
1033                        endif   
1034                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
1035                       
1036                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
1037                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
1038               
1039                        nSets = 2
1040                        break                   
1041               
1042                case "MTB":
1043                        if(isVCALC)
1044                                WAVE inten = $(folderPath+instPath+"MT"+":det_"+"MT")
1045                                WAVE/Z iErr = $("iErr_"+"MT")                   // 2D errors -- may not exist, especially for simulation               
1046                                WAVE inten2 = $(folderPath+instPath+"MB"+":det_"+"MB")
1047                                WAVE/Z iErr2 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
1048                        else
1049                                Wave inten = V_getDetectorDataW(folderStr,"MT")
1050                                Wave iErr = V_getDetectorDataErrW(folderStr,"MT")
1051                                Wave inten2 = V_getDetectorDataW(folderStr,"MB")
1052                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MB")
1053                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
1054                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
1055                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1056                                        maskMissing = 0
1057                                endif
1058                        endif   
1059                        NVAR delQ = $(folderPath+instPath+"MT"+":gDelQ_MT")
1060                       
1061                        Wave qTotal = $(folderPath+instPath+"MT"+":qTot_"+"MT")                 // 2D q-values 
1062                        Wave qTotal2 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
1063               
1064                        nSets = 2
1065                        break                           
1066               
1067                case "MLRTB":
1068                        if(isVCALC)
1069                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
1070                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
1071                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
1072                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
1073                                WAVE inten3 = $(folderPath+instPath+"MT"+":det_"+"MT")
1074                                WAVE/Z iErr3 = $("iErr_"+"MT")                  // 2D errors -- may not exist, especially for simulation               
1075                                WAVE inten4 = $(folderPath+instPath+"MB"+":det_"+"MB")
1076                                WAVE/Z iErr4 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
1077                        else
1078                                Wave inten = V_getDetectorDataW(folderStr,"ML")
1079                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
1080                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
1081                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
1082                                Wave inten3 = V_getDetectorDataW(folderStr,"MT")
1083                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"MT")
1084                                Wave inten4 = V_getDetectorDataW(folderStr,"MB")
1085                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"MB")
1086                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
1087                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
1088                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
1089                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
1090                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
1091                                        maskMissing = 0
1092                                endif
1093                        endif   
1094                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
1095                       
1096                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
1097                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
1098                        Wave qTotal3 = $(folderPath+instPath+"MT"+":qTot_"+"MT")                        // 2D q-values 
1099                        Wave qTotal4 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
1100               
1101                        nSets = 4
1102                        break                                                                   
1103                                       
1104                default:
1105                        nSets = 0                                                       
1106                        Print "ERROR   ---- type is not recognized "
1107        endswitch
1108
1109//      Print "delQ = ",delQ," for ",type
1110
1111        if(nSets == 0)
1112                SetDataFolder root:
1113                return(0)
1114        endif
1115
1116
1117// RAW data is currently read in and the 2D error wave is correctly generated
1118// 2D error is propagated through all reduction steps, but I have not
1119// verified that it is an exact duplication of the 1D error
1120//
1121//
1122//
1123// IF ther is no 2D error wave present for some reason, make a fake one
1124        if(WaveExists(iErr)==0  && WaveExists(inten) != 0)
1125                Duplicate/O inten,iErr
1126                Wave iErr=iErr
1127//              iErr = 1+sqrt(inten+0.75)                       // can't use this -- it applies to counts, not intensity (already a count rate...)
1128                iErr = sqrt(inten+0.75)                 // TODO -- here I'm just using some fictional value
1129        endif
1130        if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0)
1131                Duplicate/O inten2,iErr2
1132                Wave iErr2=iErr2
1133//              iErr2 = 1+sqrt(inten2+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1134                iErr2 = sqrt(inten2+0.75)                       // TODO -- here I'm just using some fictional value
1135        endif
1136        if(WaveExists(iErr3)==0  && WaveExists(inten3) != 0)
1137                Duplicate/O inten3,iErr3
1138                Wave iErr3=iErr3
1139//              iErr3 = 1+sqrt(inten3+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1140                iErr3 = sqrt(inten3+0.75)                       // TODO -- here I'm just using some fictional value
1141        endif
1142        if(WaveExists(iErr4)==0  && WaveExists(inten4) != 0)
1143                Duplicate/O inten4,iErr4
1144                Wave iErr4=iErr4
1145//              iErr4 = 1+sqrt(inten4+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1146                iErr4 = sqrt(inten4+0.75)                       // TODO -- here I'm just using some fictional value
1147        endif
1148
1149        // TODO -- nq will need to be larger, once the back detector is installed
1150        //
1151        // note that the back panel of 320x320 (1mm res) results in 447 data points!
1152        // - so I upped nq to 600
1153
1154        if(cmpstr(type,"B") == 0)
1155                nq = 8000
1156        else
1157                nq=600
1158        endif
1159
1160//******TODO****** -- where to put the averaged data -- right now, folderStr is forced to ""   
1161//      SetDataFolder $("root:"+folderStr)              //should already be here, but make sure...     
1162        Make/O/D/N=(nq)  $(folderPath+":"+"iBin_qxqy"+"_"+type)
1163        Make/O/D/N=(nq)  $(folderPath+":"+"qBin_qxqy"+"_"+type)
1164        Make/O/D/N=(nq)  $(folderPath+":"+"nBin_qxqy"+"_"+type)
1165        Make/O/D/N=(nq)  $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1166        Make/O/D/N=(nq)  $(folderPath+":"+"eBin_qxqy"+"_"+type)
1167        Make/O/D/N=(nq)  $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1168       
1169        Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+type)
1170        Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type)
1171        Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy"+"_"+type)
1172        Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1173        Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy"+"_"+type)
1174        Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1175       
1176       
1177//      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
1178// TODO: not sure if I want to set dQ in x or y direction...
1179        // the short dimension is the 8mm tubes, use this direction as dQ?
1180        // 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]
1181        // WRT the beam center. use qx or qy directly. Still not happy with this way...
1182
1183
1184        qBin_qxqy[] =  p*delQ   
1185        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
1186
1187        iBin_qxqy = 0
1188        iBin2_qxqy = 0
1189        eBin_qxqy = 0
1190        eBin2D_qxqy = 0
1191        nBin_qxqy = 0   //number of intensities added to each bin
1192
1193// now there are situations of:
1194// 1 panel
1195// 2 panels
1196// 4 panels
1197//
1198// this needs to be a double loop now...
1199// TODO:
1200// -- the iErr (=2D) wave and accumulation of error is NOT CALCULATED CORRECTLY YET
1201// -- verify the 2D error propagation by reducing it to 1D error
1202//
1203//
1204// The 1D error does not use iErr, and IS CALCULATED CORRECTLY
1205//
1206// x- the solid angle per pixel will be present for WORK data other than RAW, but not for RAW
1207
1208//
1209// if any of the masks don't exist, display the error, and proceed with the averaging, using all data
1210        if(maskMissing == 1)
1211                Print "Mask file not found for at least one detector - so all data is used"
1212        endif
1213       
1214        NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
1215        if(gIgnoreDetB && cmpstr(type,"B") == 0)
1216                maskMissing = 1
1217                Print "Mask skipped for B due to possible mismatch (Panel B ignored in preferences)"
1218        endif
1219
1220        Variable mask_val
1221// use set 1 (no number) only
1222        if(nSets >= 1)
1223                xDim=DimSize(inten,0)
1224                yDim=DimSize(inten,1)
1225       
1226                for(ii=0;ii<xDim;ii+=1)
1227                        for(jj=0;jj<yDim;jj+=1)
1228                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1229                                qVal = qTotal[ii][jj]
1230                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1231                                val = inten[ii][jj]
1232                               
1233                                if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
1234                                        mask_val = 0
1235                                else
1236                                        mask_val = mask[ii][jj]
1237                                endif
1238                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1239                                        iBin_qxqy[binIndex] += val
1240                                        iBin2_qxqy[binIndex] += val*val
1241                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1242                                        nBin_qxqy[binIndex] += 1
1243                                endif
1244                        endfor
1245                endfor
1246               
1247        endif
1248
1249// add in set 2 (set 1 already done)
1250        if(nSets >= 2)
1251                xDim=DimSize(inten2,0)
1252                yDim=DimSize(inten2,1)
1253       
1254                for(ii=0;ii<xDim;ii+=1)
1255                        for(jj=0;jj<yDim;jj+=1)
1256                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1257                                qVal = qTotal2[ii][jj]
1258                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1259                                val = inten2[ii][jj]
1260                               
1261                                if(isVCALC || maskMissing)
1262                                        mask_val = 0
1263                                else
1264                                        mask_val = mask2[ii][jj]
1265                                endif
1266                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1267                                        iBin_qxqy[binIndex] += val
1268                                        iBin2_qxqy[binIndex] += val*val
1269                                        eBin2D_qxqy[binIndex] += iErr2[ii][jj]*iErr2[ii][jj]
1270                                        nBin_qxqy[binIndex] += 1
1271                                endif
1272                        endfor
1273                endfor
1274               
1275        endif
1276
1277// add in set 3 and 4 (set 1 and 2 already done)
1278        if(nSets == 4)
1279                xDim=DimSize(inten3,0)
1280                yDim=DimSize(inten3,1)
1281       
1282                for(ii=0;ii<xDim;ii+=1)
1283                        for(jj=0;jj<yDim;jj+=1)
1284                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1285                                qVal = qTotal3[ii][jj]
1286                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1287                                val = inten3[ii][jj]
1288                               
1289                                if(isVCALC || maskMissing)
1290                                        mask_val = 0
1291                                else
1292                                        mask_val = mask3[ii][jj]
1293                                endif
1294                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1295                                        iBin_qxqy[binIndex] += val
1296                                        iBin2_qxqy[binIndex] += val*val
1297                                        eBin2D_qxqy[binIndex] += iErr3[ii][jj]*iErr3[ii][jj]
1298                                        nBin_qxqy[binIndex] += 1
1299                                endif
1300                        endfor
1301                endfor
1302               
1303               
1304                xDim=DimSize(inten4,0)
1305                yDim=DimSize(inten4,1)
1306       
1307                for(ii=0;ii<xDim;ii+=1)
1308                        for(jj=0;jj<yDim;jj+=1)
1309                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1310                                qVal = qTotal4[ii][jj]
1311                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1312                                val = inten4[ii][jj]
1313                               
1314                                if(isVCALC || maskMissing)
1315                                        mask_val = 0
1316                                else
1317                                        mask_val = mask4[ii][jj]
1318                                endif
1319                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1320                                        iBin_qxqy[binIndex] += val
1321                                        iBin2_qxqy[binIndex] += val*val
1322                                        eBin2D_qxqy[binIndex] += iErr4[ii][jj]*iErr4[ii][jj]
1323                                        nBin_qxqy[binIndex] += 1
1324                                endif
1325                        endfor
1326                endfor
1327               
1328        endif
1329
1330
1331// after looping through all of the data on the panels, calculate errors on I(q),
1332// just like in CircSectAve.ipf
1333// TODO:
1334// -- 2D Errors were (maybe) properly acculumated through reduction, so this loop of calculations is NOT VERIFIED (yet)
1335// x- the error on the 1D intensity, is correctly calculated as the standard error of the mean.
1336        for(ii=0;ii<nq;ii+=1)
1337                if(nBin_qxqy[ii] == 0)
1338                        //no pixels in annuli, data unknown
1339                        iBin_qxqy[ii] = 0
1340                        eBin_qxqy[ii] = 1
1341                        eBin2D_qxqy[ii] = NaN
1342                else
1343                        if(nBin_qxqy[ii] <= 1)
1344                                //need more than one pixel to determine error
1345                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1346                                eBin_qxqy[ii] = 1
1347                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1348                        else
1349                                //assume that the intensity in each pixel in annuli is normally distributed about mean...
1350                                //  -- this is correctly calculating the error as the standard error of the mean, as
1351                                //    was always done for SANS as well.
1352                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1353                                avesq = iBin_qxqy[ii]^2
1354                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1355                                var = aveisq-avesq
1356                                if(var<=0)
1357                                        eBin_qxqy[ii] = 1e-6
1358                                else
1359                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1360                                endif
1361                                // and calculate as it is propagated pixel-by-pixel
1362                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1363                        endif
1364                endif
1365        endfor
1366       
1367        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo
1368       
1369        // find the last non-zero point, working backwards
1370        val=nq
1371        do
1372                val -= 1
1373        while((nBin_qxqy[val] == 0) && val > 0)
1374       
1375//      print val, nBin_qxqy[val]
1376        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1377
1378        if(val == 0)
1379                // all the points were deleted
1380                return(0)
1381        endif
1382       
1383       
1384        // since the beam center is not always on the detector, many of the low Q bins will have zero pixels
1385        // find the first non-zero point, working forwards
1386        val = -1
1387        do
1388                val += 1
1389        while(nBin_qxqy[val] == 0)     
1390        DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1391
1392        // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
1393        val = numpnts(nBin_qxqy)-1
1394        do
1395                if(nBin_qxqy[val] == 0)
1396                        DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1397                endif
1398                val -= 1
1399        while(val>0)
1400
1401// utility function to remove NaN values from the waves
1402        V_RemoveNaNsQIS(qBin_qxqy, iBin_qxqy, eBin_qxqy)
1403
1404       
1405        // TODO:
1406        // -- This is where I calculate the resolution in SANS (see CircSectAve)
1407        // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
1408        // -- from the top of the function, folderStr = work folder, type = "FLRTB" or other type of averaging
1409        //
1410        nq = numpnts(qBin_qxqy)
1411        Make/O/D/N=(nq)  $(folderPath+":"+"sigmaQ_qxqy"+"_"+type)
1412        Make/O/D/N=(nq)  $(folderPath+":"+"qBar_qxqy"+"_"+type)
1413        Make/O/D/N=(nq)  $(folderPath+":"+"fSubS_qxqy"+"_"+type)
1414        Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+type)
1415        Wave qbar = $(folderPath+":"+"qBar_qxqy_"+type)
1416        Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+type)
1417                               
1418
1419        Variable ret1,ret2,ret3
1420        Variable lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses
1421
1422// TODO: check the units of all of the inputs
1423
1424// lambda = wavelength [A]
1425        lambda = V_getWavelength(folderStr)
1426       
1427// lambdaWidth = [dimensionless]
1428        lambdaWidth = V_getWavelength_spread(folderStr)
1429       
1430// DDet = detector pixel resolution [cm]        **assumes square pixel
1431        // V_getDet_pixel_fwhm_x(folderStr,detStr)
1432        // V_getDet_pixel_fwhm_y(folderStr,detStr)
1433//      DDet = 0.8              // TODO -- this is hard-wired
1434
1435        if(strlen(type) == 1)
1436                // it's "B"
1437                DDet = V_getDet_pixel_fwhm_x(folderStr,type)            // value is already in cm
1438        else
1439                DDet = V_getDet_pixel_fwhm_x(folderStr,type[0,1])               // value is already in cm
1440        endif
1441               
1442// apOff = sample aperture to sample distance [cm]
1443        apOff = 10              // TODO -- this is hard-wired
1444       
1445// S1 = source aperture diameter [mm]
1446// may be either circle or rectangle
1447        String s1_shape="",bs_shape=""
1448        Variable width,height,equiv_S1,equiv_bs
1449       
1450       
1451        s1_shape = V_getSourceAp_shape(folderStr)
1452        if(cmpstr(s1_shape,"CIRCLE") == 0)
1453                S1 = str2num(V_getSourceAp_size(folderStr))
1454        else
1455                S1 = V_getSourceAp_height(folderStr)            // TODO: need the width or at least an equivalent diameter
1456        endif
1457       
1458       
1459// S2 = sample aperture diameter [cm]
1460// as of 3/2018, the "internal" sample aperture is not in use, only the external
1461// TODO : verify the units on the Ap2 (external)
1462// sample aperture 1(internal) is set to report "12.7 mm" as a STRING
1463// sample aperture 2(external) reports the number typed in...
1464//
1465// so I'm trusting [cm] is in the file
1466        S2 = V_getSampleAp2_size(folderStr)*10          // sample ap 1 or 2? 2 = the "external", convert to [mm]
1467       
1468// L1 = source to sample distance [m]
1469        L1 = V_getSourceAp_distance(folderStr)/100
1470
1471// L2 = sample to detector distance [m]
1472// take the first two characters of the "type" to get the correct distance.
1473// if the type is say, MLRTB, then the implicit assumption in combining all four panels is that the resolution
1474// is not an issue for the slightly different distances.
1475        if(strlen(type) == 1)
1476                // it's "B"
1477                L2 = V_getDet_ActualDistance(folderStr,type)/100                //convert cm to m
1478        else
1479                L2 = V_getDet_ActualDistance(folderStr,type[0,1])/100           //convert cm to m
1480        endif
1481       
1482// BS = beam stop diameter [mm]
1483//TODO:? which BS is in? carr2, carr3, none?
1484// -- need to check the detector, num_beamstops field, then description, then shape/size or shape/height and shape/width
1485//
1486// TODO: the values in the file are incorrect!!! BS = 1000 mm diameter!!!
1487        BS = V_DeduceBeamstopDiameter(folderStr,type)           //returns diameter in [mm]
1488//      BS = V_getBeamStopC2_size(folderStr)            // Units are [mm]
1489//      BS = 25.4                       //TODO hard-wired value
1490       
1491//      bs_shape = V_getBeamStopC2_shape(folderStr)
1492//      if(cmpstr(s1_shape,"CIRCLE") == 0)
1493//              bs = V_getBeamStopC2_size(folderStr)
1494//      else
1495//              bs = V_getBeamStopC2_height(folderStr) 
1496//      endif
1497
1498
1499       
1500// del_r = step size [mm] = binWidth*(mm/pixel)
1501        del_r = 1*DDet*10               // TODO: this is probably not the correct value
1502
1503// usingLenses = flag for lenses = 0 if no lenses, non-zero if lenses are in-beam
1504        usingLenses = 0
1505
1506if(cmpstr(detStr,"FL")==0)
1507        Print "(FL) Resolution lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses"
1508        Print lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses
1509endif
1510
1511
1512// TODO:
1513// this is the point where I need to switch on the different collimation types (white beam, slit, Xtal, etc)
1514// to calculate the correct resolution, or fill the waves with the correct "flags"
1515//
1516
1517// For white beam data, the wavelength distribution can't be represented as a gaussian, but all of the other
1518//  geometric corrections still apply. Passing zero for the lambdaWidth will return the geometry contribution,
1519//  as long as the wavelength can be handled separately. It appears to be correct to do as a double integral,
1520//  with the inner(lambda) calculated first, then the outer(geometry).
1521//
1522
1523// possible values are:
1524//
1525// pinhole
1526// pinhole_whiteBeam
1527// convergingPinholes
1528//
1529// *slit data should be reduced using the slit routine, not here, proceed but warn
1530// narrowSlit
1531// narrowSlit_whiteBeam
1532
1533        if(cmpstr(collimationStr,"pinhole") == 0)
1534
1535                ii=0
1536                do
1537                        V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3)
1538                        sigmaq[ii] = ret1       
1539                        qbar[ii] = ret2
1540                        fsubs[ii] = ret3       
1541                        ii+=1
1542                while(ii<nq)
1543       
1544        endif
1545       
1546
1547        if(cmpstr(collimationStr,"pinhole_whiteBeam") == 0)
1548
1549//              set lambdaWidth == 0 so that the gaussian resolution calculates only the geometry contribution.
1550// the white beam distribution will need to be flagged some other way
1551//
1552                lambdaWidth = 0
1553               
1554                ii=0
1555                do
1556                        V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3)
1557                        sigmaq[ii] = ret1       
1558                        qbar[ii] = ret2
1559                        fsubs[ii] = ret3       
1560                        ii+=1
1561                while(ii<nq)
1562       
1563        endif
1564
1565        if(cmpstr(collimationStr,"convergingPinholes") == 0)
1566
1567//              set usingLenses == 1 so that the Gaussian resolution calculation will be for a focus condition
1568//
1569                usingLenses = 1
1570               
1571                ii=0
1572                do
1573                        V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3)
1574                        sigmaq[ii] = ret1       
1575                        qbar[ii] = ret2
1576                        fsubs[ii] = ret3       
1577                        ii+=1
1578                while(ii<nq)
1579       
1580        endif
1581
1582
1583// should not end up here, except for odd testing cases
1584        if(cmpstr(collimationStr,"narrowSlit") == 0)
1585
1586                Print "??? Slit data is being averaged as pinhole - reset the AVERAGE parameters in the protocol ???"
1587                ii=0
1588                do
1589                        V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3)
1590                        sigmaq[ii] = ret1       
1591                        qbar[ii] = ret2
1592                        fsubs[ii] = ret3       
1593                        ii+=1
1594                while(ii<nq)
1595       
1596        endif
1597       
1598// should not end up here, except for odd testing cases
1599        if(cmpstr(collimationStr,"narrowSlit_whiteBeam") == 0)
1600
1601//              set lambdaWidth == 0 so that the gaussian resolution calculates only the geometry contribution.
1602// the white beam distribution will need to be flagged some other way
1603//
1604                Print "??? Slit data is being averaged as pinhole - reset the AVERAGE parameters in the protocol ???"
1605
1606                lambdaWidth = 0
1607               
1608                ii=0
1609                do
1610                        V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3)
1611                        sigmaq[ii] = ret1       
1612                        qbar[ii] = ret2
1613                        fsubs[ii] = ret3       
1614                        ii+=1
1615                while(ii<nq)
1616       
1617        endif
1618
1619
1620
1621               
1622        SetDataFolder root:
1623       
1624        return(0)
1625End
1626
1627
Note: See TracBrowser for help on using the repository browser.