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

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

New dimensions added for the back detector. many functions neede to be updated to accomodate these changes. Beam center is handled in the same way (in cm, not pixels) as other panels even though this panel is like the 2D detectors on SANS.

Still missing is the real values for caibration, pixel size, dead time, etc. that are yet to be measured.

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