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

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

many minor changes after real VSANS data collected.

additional procedures added to allow easy correction of the incorrect header information from NICE.

Most notable addition is the pinhole resolution added to the calculation and the I(q) output. White beam is also treated (incorrectly) as a gaussian distrivution, but the results of smeared fitting look to be quite good.

Trimming and sorting routines are now (pinhole) resolution aware.

File identification routines have been updated to use the proper definitions of "purpose" and "intent". Both fields are now in the catalog, to allow for better sorting.

File size: 44.9 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,900,1e-6,0.01,qTot[p][q])   
64                        break
65                case "Debye":
66                        tmpInten = VC_Debye(10,300,0.0001,qTot[p][q])
67                        break
68                case "Sphere":
69                        tmpInten = VC_SphereForm(1,60,1e-6,0.001,qTot[p][q])   
70                        break
71                case "AgBeh":
72                        tmpInten = VC_BroadPeak(1e-11,3,20,100.0,0.1,3,0.1,qTot[p][q])
73                        break
74                case "Vycor":
75                        tmpInten = VC_BroadPeak(1e-9,3,20,500.0,0.015,3,0.1,qTot[p][q])
76                        break   
77                case "Empty Cell":
78                        tmpInten = VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
79                        break
80                case "Blocked Beam":
81                        tmpInten = VC_BlockedBeam(0.01,qTot[p][q])
82                        break
83                case "Debye +":
84                        tmpInten = VC_Debye(10,300,0.0001,qTot[p][q])
85                        addEmpBgd = 1
86                        break
87                case "AgBeh +":
88                        tmpInten = VC_BroadPeak(1e-11,3,20,100.0,0.1,3,0.1,qTot[p][q])
89                        addEmpBgd = 1
90                        break
91                case "Empty Cell +":
92                        tmpInten = VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
93                        tmpInten += VC_BlockedBeam(0.01,qTot[p][q])
94                        break
95                default:
96                        tmpInten = VC_Debye(10,300,0.1,qTot[p][q])
97        endswitch
98
99
100        if(addEmpBgd == 1)
101                tmpInten += VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
102                tmpInten += VC_BlockedBeam(0.01,qTot[p][q])
103        endif
104
105       
106// 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 tmpCalib
267                tmpCalib[0] = 1
268                tmpCalib[1] = 1
269                tmpcalib[2] = 10000
270        endif
271       
272//      Wave w_calib = V_getDetTube_spatialCalib("VCALC",detStr)
273        Variable tube_width = 8.4                       // TODO: UNITS!!! Hard-wired value in [mm]
274        if(cmpstr(detStr,"B") == 0)
275                V_NonLinearCorrection_B("VCALC",data,tmpCalib,tmpCalib,detStr,destPath)
276        else
277                V_NonLinearCorrection("VCALC",data,tmpCalib,tube_width,detStr,destPath)
278        endif
279                               
280        Wave/Z data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
281        Wave/Z data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
282        NVAR gUseNonLinearDet = root:Packages:NIST:VSANS:VCALC:gUseNonLinearDet
283
284        if(kBCTR_CM)
285                if(gUseNonLinearDet && WaveExists(data_realDistX) && WaveExists(data_realDistY))
286                        // beam ctr is in cm already
287
288                        // calculate all of the q-values
289                        qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
290                        qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
291                        qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
292                        qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
293               
294        //              Print "det, x_mm, y_mm ",detStr,num2str(newX),num2str(newY)
295        //              Print "det, x_pix, y_pix ",detStr,num2str(xCtr),num2str(yCtr)
296                else
297                        // do the q-calculation using linear detector
298                        //VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
299                        qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
300                        qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
301                        qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
302                        qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
303                endif   
304       
305       
306        else
307        // using the old calculation with beam center in pixels
308                if(gUseNonLinearDet && WaveExists(data_realDistX) && WaveExists(data_realDistY))
309                        // convert the beam centers to mm
310//                      String orientation
311                        Variable newX,newY
312                        dimX = DimSize(data_realDistX,0)
313                        dimY = DimSize(data_realDistX,1)
314                        if(dimX > dimY)
315                                orientation = "horizontal"
316                        else
317                                orientation = "vertical"
318                        endif
319                       
320               
321                //
322                        if(cmpstr(orientation,"vertical")==0)
323                                //      this is data dimensioned as (Ntubes,Npix)
324                                newX = tube_width*xCtr
325                                newY = data_realDistY[0][yCtr]
326                        else
327                                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
328                                newX = data_realDistX[xCtr][0]
329                                newY = tube_width*yCtr
330                        endif   
331       
332                        //if detector "B", different calculation for the centers (not tubes)
333                        if(cmpstr(detStr,"B")==0)
334                                newX = data_realDistX[xCtr][0]
335                                newY = data_realDistY[0][yCtr]
336                                //newX = xCtr
337                                //newY = yCtr
338                        endif           
339                                       
340                        // calculate all of the q-values
341                        qTot = V_CalcQval(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
342                        qx = V_CalcQX(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
343                        qy = V_CalcQY(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
344                        qz = V_CalcQZ(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
345               
346        //              Print "det, x_mm, y_mm ",detStr,num2str(newX),num2str(newY)
347        //              Print "det, x_pix, y_pix ",detStr,num2str(xCtr),num2str(yCtr)
348                else
349                        // do the q-calculation using linear detector
350                        VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
351                endif   
352       
353        endif
354       
355        KillWaves/Z tmpCalib
356       
357        return(0)
358End
359
360
361//////////////////////
362// NOTE: The Q calculations are different than what is in GaussUtils in that they take into
363// accout the different x/y pixel sizes and the beam center not being on the detector -
364// off a different edge for each LRTB type
365/////////////////////
366
367//function to calculate the overall q-value, given all of the necesary trig inputs
368//and are in detector coordinates (1,128) rather than axis values
369//the pixel locations need not be integers, reals are ok inputs
370//sdd is in [cm]
371//wavelength is in Angstroms
372//
373//returned magnitude of Q is in 1/Angstroms
374//
375//
376Function VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
377        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
378       
379        Variable dx,dy,qval,two_theta,dist
380               
381
382        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
383        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
384        dist = sqrt(dx^2 + dy^2)
385       
386        two_theta = atan(dist/sdd)
387
388        qval = 4*Pi/lam*sin(two_theta/2)
389       
390        return qval
391End
392
393//calculates just the q-value in the x-direction on the detector
394//input/output is the same as CalcQval()
395//ALL inputs are in detector coordinates
396//
397//sdd is in [cm]
398//wavelength is in Angstroms
399//
400// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
401// now properly accounts for qz
402//
403Function VC_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
404        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
405
406        Variable qx,qval,phi,dx,dy,dist,two_theta
407       
408        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
409       
410//      sdd *=100               //convert to cm
411        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
412        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
413        phi = V_FindPhi(dx,dy)
414       
415        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
416        dist = sqrt(dx^2 + dy^2)
417        two_theta = atan(dist/sdd)
418
419        qx = qval*cos(two_theta/2)*cos(phi)
420       
421        return qx
422End
423
424//calculates just the q-value in the y-direction on the detector
425//input/output is the same as CalcQval()
426//ALL inputs are in detector coordinates
427//sdd is in [cm]
428//wavelength is in Angstroms
429//
430// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
431// now properly accounts for qz
432//
433Function VC_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
434        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
435       
436        Variable dy,qval,dx,phi,qy,dist,two_theta
437       
438        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
439       
440//      sdd *=100               //convert to cm
441        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
442        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
443        phi = V_FindPhi(dx,dy)
444       
445        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
446        dist = sqrt(dx^2 + dy^2)
447        two_theta = atan(dist/sdd)
448       
449        qy = qval*cos(two_theta/2)*sin(phi)
450       
451        return qy
452End
453
454//calculates just the z-component of the q-vector, not measured on the detector
455//input/output is the same as CalcQval()
456//ALL inputs are in detector coordinates
457//sdd is in [cm]
458//wavelength is in Angstroms
459//
460// not actually used, but here for completeness if anyone asks
461//
462Function VC_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
463        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
464       
465        Variable dy,qval,dx,phi,qz,dist,two_theta
466       
467        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
468       
469//      sdd *=100               //convert to cm
470       
471        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
472        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
473        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
474        dist = sqrt(dx^2 + dy^2)
475        two_theta = atan(dist/sdd)
476       
477        qz = qval*sin(two_theta/2)
478       
479        return qz
480End
481
482//phi is defined from +x axis, proceeding CCW around [0,2Pi]
483Threadsafe Function V_FindPhi(vx,vy)
484        variable vx,vy
485       
486        variable phi
487       
488        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
489       
490        // special cases
491        if(vx==0 && vy > 0)
492                return(pi/2)
493        endif
494        if(vx==0 && vy < 0)
495                return(3*pi/2)
496        endif
497        if(vx >= 0 && vy == 0)
498                return(0)
499        endif
500        if(vx < 0 && vy == 0)
501                return(pi)
502        endif
503       
504       
505        if(vx > 0 && vy > 0)
506                return(phi)
507        endif
508        if(vx < 0 && vy > 0)
509                return(phi + pi)
510        endif
511        if(vx < 0 && vy < 0)
512                return(phi + pi)
513        endif
514        if( vx > 0 && vy < 0)
515                return(phi + 2*pi)
516        endif
517       
518        return(phi)
519end
520
521Function VC_SphereForm(scale,radius,delrho,bkg,x)                               
522        Variable scale,radius,delrho,bkg
523        Variable x
524       
525        // variables are:                                                       
526        //[0] scale
527        //[1] radius (A)
528        //[2] delrho (A-2)
529        //[3] background (cm-1)
530       
531//      Variable scale,radius,delrho,bkg                               
532//      scale = w[0]
533//      radius = w[1]
534//      delrho = w[2]
535//      bkg = w[3]
536       
537       
538        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
539        // and is rescaled to give [=] cm^-1
540       
541        Variable bes,f,vol,f2
542        //
543        //handle q==0 separately
544        If(x==0)
545                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
546                return(f)
547        Endif
548       
549//      bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
550       
551        bes = 3*sqrt(pi/(2*x*radius))*BesselJ(1.5,x*radius)/(x*radius)
552       
553        vol = 4*pi/3*radius^3
554        f = vol*bes*delrho              // [=] A
555        // normalize to single particle volume, convert to 1/cm
556        f2 = f * f / vol * 1.0e8                // [=] 1/cm
557       
558        return (scale*f2+bkg)   // Scale, then add in the background
559       
560End
561
562Function VC_Debye(scale,rg,bkg,x)
563        Variable scale,rg,bkg
564        Variable x
565       
566        // variables are:
567        //[0] scale factor
568        //[1] radius of gyration [A]
569        //[2] background        [cm-1]
570       
571        // calculates (scale*debye)+bkg
572        Variable Pq,qr2
573       
574        qr2=(x*rg)^2
575        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
576       
577        //scale
578        Pq *= scale
579        // then add in the background
580        return (Pq+bkg)
581End
582
583// a sum of a power law and debye to approximate the scattering from a real empty cell
584//
585//      make/O/D coef_ECEmp = {2.2e-8,3.346,0.0065,9.0,0.016}
586//
587Function VC_EC_Empirical(aa,mm,scale,rg,bkg,x)
588        Variable aa,mm,scale,rg,bkg
589        Variable x
590       
591        // variables are:
592        //[0] = A
593        //[1] = power m
594        //[2] scale factor
595        //[3] radius of gyration [A]
596        //[4] background        [cm-1]
597       
598        Variable Iq
599       
600        // calculates (scale*debye)+bkg
601        Variable Pq,qr2
602       
603//      if(x*Rg < 1e-3)         //added Oct 2008 to avoid numerical errors at low arg values
604//              return(scale+bkg)
605//      endif
606       
607        Iq = aa*x^-mm
608       
609        qr2=(x*rg)^2
610        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
611       
612        //scale
613        Pq *= scale
614        // then add the terms up
615        return (Iq + Pq + bkg)
616End
617
618// blocked beam
619//
620Function VC_BlockedBeam(bkg,x)
621        Variable bkg
622        Variable x
623       
624        return (bkg)
625End
626
627
628//
629// a broad peak to simulate silver behenate or vycor
630//
631// peak @ 0.1 ~ AgBeh
632//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}               
633//
634//
635// peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
636//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}             
637//
638//
639Function VC_BroadPeak(aa,nn,cc,LL,Qzero,mm,bgd,x)
640        Variable aa,nn,cc,LL,Qzero,mm,bgd
641        Variable x
642       
643        // variables are:                                                       
644        //[0] Porod term scaling
645        //[1] Porod exponent
646        //[2] Lorentzian term scaling
647        //[3] Lorentzian screening length [A]
648        //[4] peak location [1/A]
649        //[5] Lorentzian exponent
650        //[6] background
651       
652//      local variables
653        Variable inten, qval
654//      x is the q-value for the calculation
655        qval = x
656//      do the calculation and return the function value
657       
658        inten = aa/(qval)^nn + cc/(1 + (abs(qval-Qzero)*LL)^mm) + bgd
659
660        Return (inten)
661       
662End
663
664//
665// updated to new folder structure Feb 2016
666// folderStr = RAW,SAM, VCALC or other
667// detStr is the panel identifer "ML", etc.
668//
669Function SetDeltaQ(folderStr,detStr)
670        String folderStr,detStr
671
672        Variable isVCALC
673        if(cmpstr(folderStr,"VCALC") == 0)
674                isVCALC = 1
675        endif
676       
677        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
678        String instPath = ":entry:instrument:detector_"
679               
680        if(isVCALC)
681                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)               // 2D detector data
682        else
683                Wave inten = V_getDetectorDataW(folderStr,detStr)
684        endif
685
686        Wave qx = $(folderPath+instPath+detStr+":qx_"+detStr)
687        Wave qy = $(folderPath+instPath+detStr+":qy_"+detStr)
688       
689        Variable xDim,yDim,delQ
690       
691        xDim=DimSize(inten,0)
692        yDim=DimSize(inten,1)
693       
694        if(xDim<yDim)
695                delQ = abs(qx[0][0] - qx[1][0])/2
696        else
697                delQ = abs(qy[0][1] - qy[0][0])/2
698        endif
699       
700        // set the global
701        Variable/G $(folderPath+instPath+detStr+":gDelQ_"+detStr) = delQ
702//      Print "SET delQ = ",delQ," for ",type
703       
704        return(delQ)
705end
706
707
708//TODO -- need a switch here to dispatch to the averaging type
709Proc VC_BinQxQy_to_1D(folderStr,type)
710        String folderStr
711        String type
712//      Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4)
713//      Prompt type,"detector identifier"
714
715
716        VC_fDoBinning_QxQy2D(folderStr, type)
717
718
719/// this is for a tall, narrow slit mode       
720//      VC_fBinDetector_byRows(folderStr,type)
721       
722End
723
724
725// folderStr is RAW, VCALC, SAM, etc.
726// type is "B", "FL" for single binning, "FLR", or "MLRTB" or similar if multiple panels are combined
727//
728Proc VC_Graph_1D_detType(folderStr,type)
729        String folderStr,type
730       
731        SetDataFolder $("root:Packages:NIST:VSANS:"+folderStr)
732       
733        Display $("iBin_qxqy"+"_"+type) vs $("qBin_qxqy"+"_"+type)
734        ModifyGraph mirror=2,grid=1,log=1
735        ModifyGraph mode=4,marker=19,msize=2
736//      ErrorBars/T=0 iBin_qxqy Y,wave=(eBin2D_qxqy,eBin2D_qxqy)                // for simulations, I don't have 2D uncertainty
737        ErrorBars/T=0 $("iBin_qxqy"+"_"+type) Y,wave=($("eBin_qxqy"+"_"+type),$("eBin_qxqy"+"_"+type))
738        legend
739       
740        SetDataFolder root:
741
742End
743
744
745
746//////////
747//
748//              Function that bins a 2D detctor panel into I(q) based on the q-value of the pixel
749//              - each pixel QxQyQz has been calculated beforehand
750//              - if multiple panels are selected to be combined, it is done here during the binning
751//              - the setting of deltaQ step is still a little suspect (TODO)
752//
753//
754// see the equivalent function in PlotUtils2D_v40.ipf
755//
756//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
757//
758// this has been modified to accept different detector panels and to take arrays
759// -- type = FL or FR or...other panel identifiers
760//
761// TODO "iErr" is not always defined correctly since it doesn't really apply here for data that is not 2D simulation
762//
763//
764// updated Feb2016 to take new folder structure
765// TODO
766// -- VERIFY
767// -- figure out what the best location is to put the averaged data? currently @ top level of WORK folder
768//    but this is a lousy choice.
769// x- binning is now Mask-aware. If mask is not present, all data is used. If data is from VCALC, all data is used
770// x- Where do I put the solid angle correction? In here as a weight for each point, or later on as
771//    a blanket correction (matrix multiply) for an entire panel? (Solid Angle correction is done in the
772//    step where data is added to a WORK file (see Raw_to_Work())
773//
774//
775//
776// folderStr = WORK folder, type = the binning type (may include multiple detectors)
777Function VC_fDoBinning_QxQy2D(folderStr,type)
778        String folderStr,type
779       
780        Variable nSets = 0
781        Variable xDim,yDim
782        Variable ii,jj
783        Variable qVal,nq,var,avesq,aveisq
784        Variable binIndex,val,isVCALC=0,maskMissing
785
786        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
787        String instPath = ":entry:instrument:detector_"
788        String detStr
789               
790        if(cmpstr(folderStr,"VCALC") == 0)
791                isVCALC = 1
792        endif
793       
794// now switch on the type to determine which waves to declare and create
795// since there may be more than one panel to step through. There may be two, there may be four
796//
797
798// TODO:
799// -- Solid_Angle -- waves will be present for WORK data other than RAW, but not for RAW
800//
801// assume that the mask files are missing unless we can find them. If VCALC data,
802//  then the Mask is missing by definition
803        maskMissing = 1
804
805        strswitch(type) // string switch
806                case "FL":              // execute if case matches expression
807                case "FR":
808                        detStr = type
809                        if(isVCALC)
810                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
811                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation
812                        else
813                                Wave inten = V_getDetectorDataW(folderStr,detStr)
814                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
815                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
816                                if(WaveExists(mask) == 1)
817                                        maskMissing = 0
818                                endif
819                               
820                        endif
821                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
822                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
823                        nSets = 1
824                        break   
825                                                               
826                case "FT":             
827                case "FB":
828                        detStr = type
829                        if(isVCALC)
830                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
831                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
832                        else
833                                Wave inten = V_getDetectorDataW(folderStr,detStr)
834                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
835                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
836                                if(WaveExists(mask) == 1)
837                                        maskMissing = 0
838                                endif
839                        endif
840                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
841                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
842                        nSets = 1
843                        break
844                       
845                case "ML":             
846                case "MR":
847                        detStr = type
848                        if(isVCALC)
849                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
850                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
851                        else
852                                Wave inten = V_getDetectorDataW(folderStr,detStr)
853                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
854                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
855                                if(WaveExists(mask) == 1)
856                                        maskMissing = 0
857                                endif
858                        endif   
859                        //TODO:
860                        // -- decide on the proper deltaQ for binning. either nominal value for LR, or one
861                        //    determined specifically for that panel (currently using one tube width as deltaQ)
862                        // -- this is repeated multiple times in this switch
863                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
864//                      NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
865                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
866                        nSets = 1
867                        break   
868                                       
869                case "MT":             
870                case "MB":
871                        detStr = type
872                        if(isVCALC)
873                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
874                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
875                        else
876                                Wave inten = V_getDetectorDataW(folderStr,detStr)
877                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
878                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
879                                if(WaveExists(mask) == 1)
880                                        maskMissing = 0
881                                endif
882                        endif   
883                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
884                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
885                        nSets = 1
886                        break   
887                                       
888                case "B":       
889                        detStr = type
890                        if(isVCALC)
891                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
892                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
893                        else
894                                Wave inten = V_getDetectorDataW(folderStr,detStr)
895                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
896                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
897                                if(WaveExists(mask) == 1)
898                                        maskMissing = 0
899                                endif
900                        endif   
901                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_B")
902                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values 
903                        nSets = 1
904                        break   
905                       
906                case "FLR":
907                // detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
908                // TODO
909                // -- see if I can un-hard-wire some of this below when more than one panel is combined
910                        if(isVCALC)
911                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
912                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
913                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
914                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
915                        else
916                                Wave inten = V_getDetectorDataW(folderStr,"FL")
917                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
918                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
919                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
920                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
921                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
922                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
923                                        maskMissing = 0
924                                endif
925                        endif   
926                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
927                       
928                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
929                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
930               
931                        nSets = 2
932                        break                   
933               
934                case "FTB":
935                        if(isVCALC)
936                                WAVE inten = $(folderPath+instPath+"FT"+":det_"+"FT")
937                                WAVE/Z iErr = $("iErr_"+"FT")                   // 2D errors -- may not exist, especially for simulation               
938                                WAVE inten2 = $(folderPath+instPath+"FB"+":det_"+"FB")
939                                WAVE/Z iErr2 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
940                        else
941                                Wave inten = V_getDetectorDataW(folderStr,"FT")
942                                Wave iErr = V_getDetectorDataErrW(folderStr,"FT")
943                                Wave inten2 = V_getDetectorDataW(folderStr,"FB")
944                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FB")
945                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
946                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
947                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
948                                        maskMissing = 0
949                                endif
950                        endif   
951                        NVAR delQ = $(folderPath+instPath+"FT"+":gDelQ_FT")
952                       
953                        Wave qTotal = $(folderPath+instPath+"FT"+":qTot_"+"FT")                 // 2D q-values 
954                        Wave qTotal2 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
955       
956                        nSets = 2
957                        break           
958               
959                case "FLRTB":
960                        if(isVCALC)
961                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
962                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
963                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
964                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
965                                WAVE inten3 = $(folderPath+instPath+"FT"+":det_"+"FT")
966                                WAVE/Z iErr3 = $("iErr_"+"FT")                  // 2D errors -- may not exist, especially for simulation               
967                                WAVE inten4 = $(folderPath+instPath+"FB"+":det_"+"FB")
968                                WAVE/Z iErr4 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
969                        else
970                                Wave inten = V_getDetectorDataW(folderStr,"FL")
971                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
972                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
973                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
974                                Wave inten3 = V_getDetectorDataW(folderStr,"FT")
975                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"FT")
976                                Wave inten4 = V_getDetectorDataW(folderStr,"FB")
977                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"FB")
978                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
979                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
980                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
981                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
982                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
983                                        maskMissing = 0
984                                endif
985                        endif   
986                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
987                       
988                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
989                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
990                        Wave qTotal3 = $(folderPath+instPath+"FT"+":qTot_"+"FT")                        // 2D q-values 
991                        Wave qTotal4 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
992               
993                        nSets = 4
994                        break           
995                       
996                case "MLR":
997                        if(isVCALC)
998                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
999                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
1000                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
1001                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
1002                        else
1003                                Wave inten = V_getDetectorDataW(folderStr,"ML")
1004                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
1005                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
1006                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
1007                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
1008                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
1009                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1010                                        maskMissing = 0
1011                                endif
1012                        endif   
1013                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
1014                       
1015                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
1016                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
1017               
1018                        nSets = 2
1019                        break                   
1020               
1021                case "MTB":
1022                        if(isVCALC)
1023                                WAVE inten = $(folderPath+instPath+"MT"+":det_"+"MT")
1024                                WAVE/Z iErr = $("iErr_"+"MT")                   // 2D errors -- may not exist, especially for simulation               
1025                                WAVE inten2 = $(folderPath+instPath+"MB"+":det_"+"MB")
1026                                WAVE/Z iErr2 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
1027                        else
1028                                Wave inten = V_getDetectorDataW(folderStr,"MT")
1029                                Wave iErr = V_getDetectorDataErrW(folderStr,"MT")
1030                                Wave inten2 = V_getDetectorDataW(folderStr,"MB")
1031                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MB")
1032                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
1033                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
1034                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1035                                        maskMissing = 0
1036                                endif
1037                        endif   
1038                        NVAR delQ = $(folderPath+instPath+"MT"+":gDelQ_MT")
1039                       
1040                        Wave qTotal = $(folderPath+instPath+"MT"+":qTot_"+"MT")                 // 2D q-values 
1041                        Wave qTotal2 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
1042               
1043                        nSets = 2
1044                        break                           
1045               
1046                case "MLRTB":
1047                        if(isVCALC)
1048                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
1049                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
1050                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
1051                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
1052                                WAVE inten3 = $(folderPath+instPath+"MT"+":det_"+"MT")
1053                                WAVE/Z iErr3 = $("iErr_"+"MT")                  // 2D errors -- may not exist, especially for simulation               
1054                                WAVE inten4 = $(folderPath+instPath+"MB"+":det_"+"MB")
1055                                WAVE/Z iErr4 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
1056                        else
1057                                Wave inten = V_getDetectorDataW(folderStr,"ML")
1058                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
1059                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
1060                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
1061                                Wave inten3 = V_getDetectorDataW(folderStr,"MT")
1062                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"MT")
1063                                Wave inten4 = V_getDetectorDataW(folderStr,"MB")
1064                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"MB")
1065                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
1066                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
1067                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
1068                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
1069                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
1070                                        maskMissing = 0
1071                                endif
1072                        endif   
1073                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
1074                       
1075                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
1076                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
1077                        Wave qTotal3 = $(folderPath+instPath+"MT"+":qTot_"+"MT")                        // 2D q-values 
1078                        Wave qTotal4 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
1079               
1080                        nSets = 4
1081                        break                                                                   
1082                                       
1083                default:
1084                        nSets = 0                                                       // optional default expression executed
1085                        Print "ERROR   ---- type is not recognized "
1086        endswitch
1087
1088//      Print "delQ = ",delQ," for ",type
1089
1090        if(nSets == 0)
1091                SetDataFolder root:
1092                return(0)
1093        endif
1094
1095
1096//TODO: properly define the 2D errors here - I'll have this if I do the simulation
1097// -- need to propagate the 2D errors up to this point
1098//
1099        if(WaveExists(iErr)==0  && WaveExists(inten) != 0)
1100                Duplicate/O inten,iErr
1101                Wave iErr=iErr
1102//              iErr = 1+sqrt(inten+0.75)                       // can't use this -- it applies to counts, not intensity (already a count rate...)
1103                iErr = sqrt(inten+0.75)                 // TODO -- here I'm just using some fictional value
1104        endif
1105        if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0)
1106                Duplicate/O inten2,iErr2
1107                Wave iErr2=iErr2
1108//              iErr2 = 1+sqrt(inten2+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1109                iErr2 = sqrt(inten2+0.75)                       // TODO -- here I'm just using some fictional value
1110        endif
1111        if(WaveExists(iErr3)==0  && WaveExists(inten3) != 0)
1112                Duplicate/O inten3,iErr3
1113                Wave iErr3=iErr3
1114//              iErr3 = 1+sqrt(inten3+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1115                iErr3 = sqrt(inten3+0.75)                       // TODO -- here I'm just using some fictional value
1116        endif
1117        if(WaveExists(iErr4)==0  && WaveExists(inten4) != 0)
1118                Duplicate/O inten4,iErr4
1119                Wave iErr4=iErr4
1120//              iErr4 = 1+sqrt(inten4+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1121                iErr4 = sqrt(inten4+0.75)                       // TODO -- here I'm just using some fictional value
1122        endif
1123
1124        // TODO -- nq will need to be larger, once the back detector is installed
1125        //
1126        // note that the back panel of 320x320 (1mm res) results in 447 data points!
1127        // - so I upped nq to 600
1128
1129        nq = 600
1130
1131//******TODO****** -- where to put the averaged data -- right now, folderStr is forced to ""   
1132//      SetDataFolder $("root:"+folderStr)              //should already be here, but make sure...     
1133        Make/O/D/N=(nq)  $(folderPath+":"+"iBin_qxqy"+"_"+type)
1134        Make/O/D/N=(nq)  $(folderPath+":"+"qBin_qxqy"+"_"+type)
1135        Make/O/D/N=(nq)  $(folderPath+":"+"nBin_qxqy"+"_"+type)
1136        Make/O/D/N=(nq)  $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1137        Make/O/D/N=(nq)  $(folderPath+":"+"eBin_qxqy"+"_"+type)
1138        Make/O/D/N=(nq)  $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1139       
1140        Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+type)
1141        Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type)
1142        Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy"+"_"+type)
1143        Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1144        Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy"+"_"+type)
1145        Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1146       
1147       
1148//      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
1149// TODO: not sure if I want to set dQ in x or y direction...
1150        // the short dimension is the 8mm tubes, use this direction as dQ?
1151        // 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]
1152        // WRT the beam center. use qx or qy directly. Still not happy with this way...
1153
1154
1155        qBin_qxqy[] =  p*delQ   
1156        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
1157
1158        iBin_qxqy = 0
1159        iBin2_qxqy = 0
1160        eBin_qxqy = 0
1161        eBin2D_qxqy = 0
1162        nBin_qxqy = 0   //number of intensities added to each bin
1163
1164// now there are situations of:
1165// 1 panel
1166// 2 panels
1167// 4 panels
1168//
1169// this needs to be a double loop now...
1170// TODO:
1171// -- the iErr (=2D) wave and accumulation of error is NOT CALCULATED CORRECTLY YET
1172// -- the solid angle per pixel is not completely implemented.
1173//    it will be present for WORK data other than RAW, but not for RAW
1174
1175// if any of the masks don't exist, display the error, and proceed with the averaging, using all data
1176        if(maskMissing == 1)
1177                Print "Mask file not found for at least one detector - so all data is used"
1178        endif
1179
1180        Variable mask_val
1181// use set 1 (no number) only
1182        if(nSets >= 1)
1183                xDim=DimSize(inten,0)
1184                yDim=DimSize(inten,1)
1185       
1186                for(ii=0;ii<xDim;ii+=1)
1187                        for(jj=0;jj<yDim;jj+=1)
1188                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1189                                qVal = qTotal[ii][jj]
1190                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1191                                val = inten[ii][jj]
1192                               
1193                                if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
1194                                        mask_val = 0
1195                                else
1196                                        mask_val = mask[ii][jj]
1197                                endif
1198                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1199                                        iBin_qxqy[binIndex] += val
1200                                        iBin2_qxqy[binIndex] += val*val
1201                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1202                                        nBin_qxqy[binIndex] += 1
1203                                endif
1204                        endfor
1205                endfor
1206               
1207        endif
1208
1209// add in set 2 (set 1 already done)
1210        if(nSets >= 2)
1211                xDim=DimSize(inten2,0)
1212                yDim=DimSize(inten2,1)
1213       
1214                for(ii=0;ii<xDim;ii+=1)
1215                        for(jj=0;jj<yDim;jj+=1)
1216                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1217                                qVal = qTotal2[ii][jj]
1218                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1219                                val = inten2[ii][jj]
1220                               
1221                                if(isVCALC || maskMissing)
1222                                        mask_val = 0
1223                                else
1224                                        mask_val = mask2[ii][jj]
1225                                endif
1226                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1227                                        iBin_qxqy[binIndex] += val
1228                                        iBin2_qxqy[binIndex] += val*val
1229                                        eBin2D_qxqy[binIndex] += iErr2[ii][jj]*iErr2[ii][jj]
1230                                        nBin_qxqy[binIndex] += 1
1231                                endif
1232                        endfor
1233                endfor
1234               
1235        endif
1236
1237// add in set 3 and 4 (set 1 and 2 already done)
1238        if(nSets == 4)
1239                xDim=DimSize(inten3,0)
1240                yDim=DimSize(inten3,1)
1241       
1242                for(ii=0;ii<xDim;ii+=1)
1243                        for(jj=0;jj<yDim;jj+=1)
1244                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1245                                qVal = qTotal3[ii][jj]
1246                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1247                                val = inten3[ii][jj]
1248                               
1249                                if(isVCALC || maskMissing)
1250                                        mask_val = 0
1251                                else
1252                                        mask_val = mask3[ii][jj]
1253                                endif
1254                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1255                                        iBin_qxqy[binIndex] += val
1256                                        iBin2_qxqy[binIndex] += val*val
1257                                        eBin2D_qxqy[binIndex] += iErr3[ii][jj]*iErr3[ii][jj]
1258                                        nBin_qxqy[binIndex] += 1
1259                                endif
1260                        endfor
1261                endfor
1262               
1263               
1264                xDim=DimSize(inten4,0)
1265                yDim=DimSize(inten4,1)
1266       
1267                for(ii=0;ii<xDim;ii+=1)
1268                        for(jj=0;jj<yDim;jj+=1)
1269                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1270                                qVal = qTotal4[ii][jj]
1271                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1272                                val = inten4[ii][jj]
1273                               
1274                                if(isVCALC || maskMissing)
1275                                        mask_val = 0
1276                                else
1277                                        mask_val = mask4[ii][jj]
1278                                endif
1279                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1280                                        iBin_qxqy[binIndex] += val
1281                                        iBin2_qxqy[binIndex] += val*val
1282                                        eBin2D_qxqy[binIndex] += iErr4[ii][jj]*iErr4[ii][jj]
1283                                        nBin_qxqy[binIndex] += 1
1284                                endif
1285                        endfor
1286                endfor
1287               
1288        endif
1289
1290
1291// after looping through all of the data on the panels, calculate errors on I(q),
1292// just like in CircSectAve.ipf
1293// TODO:
1294// -- 2D Errors were NOT properly acculumated through reduction, so this loop of calculations is NOT MEANINGFUL (yet)
1295// x- the error on the 1D intensity, is correctly calculated as the standard error of the mean.
1296        for(ii=0;ii<nq;ii+=1)
1297                if(nBin_qxqy[ii] == 0)
1298                        //no pixels in annuli, data unknown
1299                        iBin_qxqy[ii] = 0
1300                        eBin_qxqy[ii] = 1
1301                        eBin2D_qxqy[ii] = NaN
1302                else
1303                        if(nBin_qxqy[ii] <= 1)
1304                                //need more than one pixel to determine error
1305                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1306                                eBin_qxqy[ii] = 1
1307                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1308                        else
1309                                //assume that the intensity in each pixel in annuli is normally distributed about mean...
1310                                //  -- this is correctly calculating the error as the standard error of the mean, as
1311                                //    was always done for SANS as well.
1312                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1313                                avesq = iBin_qxqy[ii]^2
1314                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1315                                var = aveisq-avesq
1316                                if(var<=0)
1317                                        eBin_qxqy[ii] = 1e-6
1318                                else
1319                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1320                                endif
1321                                // and calculate as it is propagated pixel-by-pixel
1322                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1323                        endif
1324                endif
1325        endfor
1326       
1327        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo
1328       
1329        // find the last non-zero point, working backwards
1330        val=nq
1331        do
1332                val -= 1
1333        while((nBin_qxqy[val] == 0) && val > 0)
1334       
1335//      print val, nBin_qxqy[val]
1336        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1337
1338        if(val == 0)
1339                // all the points were deleted
1340                return(0)
1341        endif
1342       
1343       
1344        // since the beam center is not always on the detector, many of the low Q bins will have zero pixels
1345        // find the first non-zero point, working forwards
1346        val = -1
1347        do
1348                val += 1
1349        while(nBin_qxqy[val] == 0)     
1350        DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1351
1352        // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
1353        val = numpnts(nBin_qxqy)-1
1354        do
1355                if(nBin_qxqy[val] == 0)
1356                        DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1357                endif
1358                val -= 1
1359        while(val>0)
1360       
1361        // TODO:
1362        // -- This is where I calculate the resolution in SANS (see CircSectAve)
1363        // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
1364        // -- from the top of the function, folderStr = work folder, type = "FLRTB" or other type of averaging
1365        //
1366        nq = numpnts(qBin_qxqy)
1367        Make/O/D/N=(nq)  $(folderPath+":"+"sigmaQ_qxqy"+"_"+type)
1368        Make/O/D/N=(nq)  $(folderPath+":"+"qBar_qxqy"+"_"+type)
1369        Make/O/D/N=(nq)  $(folderPath+":"+"fSubS_qxqy"+"_"+type)
1370        Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+type)
1371        Wave qbar = $(folderPath+":"+"qBar_qxqy_"+type)
1372        Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+type)
1373                               
1374
1375        Variable ret1,ret2,ret3
1376        Variable lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses
1377
1378// TODO: check the units of all of the inputs
1379
1380// lambda = wavelength [A]
1381        lambda = V_getWavelength(folderStr)
1382       
1383// lambdaWidth = [dimensionless]
1384        lambdaWidth = V_getWavelength_spread(folderStr)
1385       
1386// DDet = detector pixel resolution [cm]        **assumes square pixel
1387        DDet = 0.8              // TODO -- this is hard-wired
1388       
1389// apOff = sample aperture to sample distance [cm]
1390        apOff = 10              // TODO -- this is hard-wired
1391       
1392// S1 = source aperture diameter [mm]
1393        S1 = str2num(V_getSourceAp_size(folderStr))
1394       
1395// S2 = sample aperture diameter [mm]
1396        S2 = V_getSampleAp2_size(folderStr)*10          // sample ap 1 or 2? the "external" may not exist?
1397       
1398// L1 = source to sample distance [m]
1399        L1 = V_getSourceAp_distance(folderStr)/100
1400
1401// L2 = sample to detector distance [m]
1402        L2 = V_getDet_ActualDistance(folderStr,type[0]+"L")/100         //TODO "L" panel is hard wired and is WRONG
1403       
1404// BS = beam stop diameter [mm]
1405        //BS = V_getBeamStopC2_size(folderStr)          // TODO: what are the units? which BS is in? carr2, carr3, back, none?
1406        BS = 25.4                       //TODO hard-wired value
1407       
1408// del_r = step size [mm] = binWidth*(mm/pixel)
1409        del_r = 1*8                     // TODO: 8mm/pixel hard-wired
1410
1411// usingLenses = flag for lenses = 0 if no lenses, non-zero if lenses are in-beam
1412        usingLenses = 0
1413
1414
1415Print "Resolution lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses"
1416Print lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses
1417
1418        ii=0
1419        do
1420                V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3)
1421                sigmaq[ii] = ret1       
1422                qbar[ii] = ret2
1423                fsubs[ii] = ret3       
1424                ii+=1
1425        while(ii<nq)
1426       
1427       
1428       
1429        SetDataFolder root:
1430       
1431        return(0)
1432End
1433
1434
Note: See TracBrowser for help on using the repository browser.