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

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

added procedures to output QxQy_ASCII data. Each panel is output into its own file. Output format is the same as for 2D SANS data, including the 2D resolution function. However, reading the data back in with the SANS analysis macros currently does not redimension the data back to the matrix correctly as it assumes a square detector.

I will need to add the X and Y dimensions of each panel into the header, and then make use of these values when they are read in. Also, writing the QxQy? data is quick for the M and F panels ( 0.8 s) but is extremely slow for the back, High-res panel (120 s) since there are 1.1.E6 points there vs. 6144 pts. I'll need to find a way to speed this operation up.

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