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

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

added a few corrections to the reduction:

Added sample apertures to the patch panel so that they can be corrected

A flag is now written to the data files if the "flip" has been done, and it will refuse to flip again. This flag can be reset if something goes wrong.

Multiple reduce now allows run numbers to be entered as is for SANS,

Filtering of files for the protocol panel should be better behaved now.

File size: 47.0 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// TODO:
776// -- some of the input parameters for the resolution calcuation are either assumed (apOff) or are currently
777//    hard-wired. these need to be corrected before even the pinhole resolution is correct
778// x- resolution calculation is in the correct place. The calculation is done per-panel (specified by TYPE),
779//    and then the unwanted points can be discarded (all 6 columns) as the data is trimmed and concatenated
780//    is separate functions that are resolution-aware.
781//
782//
783// folderStr = WORK folder, type = the binning type (may include multiple detectors)
784Function VC_fDoBinning_QxQy2D(folderStr,type)
785        String folderStr,type
786       
787        Variable nSets = 0
788        Variable xDim,yDim
789        Variable ii,jj
790        Variable qVal,nq,var,avesq,aveisq
791        Variable binIndex,val,isVCALC=0,maskMissing
792
793        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
794        String instPath = ":entry:instrument:detector_"
795        String detStr
796               
797        if(cmpstr(folderStr,"VCALC") == 0)
798                isVCALC = 1
799        endif
800       
801        detStr = type
802
803// now switch on the type to determine which waves to declare and create
804// since there may be more than one panel to step through. There may be two, there may be four
805//
806
807// TODO:
808// -- Solid_Angle -- waves will be present for WORK data other than RAW, but not for RAW
809//
810// assume that the mask files are missing unless we can find them. If VCALC data,
811//  then the Mask is missing by definition
812        maskMissing = 1
813
814        strswitch(type) // string switch
815//              case "FL":              // execute if case matches expression
816//              case "FR":
817//                      detStr = type
818//                      if(isVCALC)
819//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
820//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation
821//                      else
822//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
823//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
824//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
825//                              if(WaveExists(mask) == 1)
826//                                      maskMissing = 0
827//                              endif
828//                             
829//                      endif
830//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
831//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
832//                      nSets = 1
833//                      break   
834                                                               
835//              case "FT":             
836//              case "FB":
837//                      detStr = type
838//                      if(isVCALC)
839//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
840//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
841//                      else
842//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
843//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
844//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
845//                              if(WaveExists(mask) == 1)
846//                                      maskMissing = 0
847//                              endif
848//                      endif
849//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
850//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
851//                      nSets = 1
852//                      break
853                       
854//              case "ML":             
855//              case "MR":
856//                      detStr = type
857//                      if(isVCALC)
858//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
859//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
860//                      else
861//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
862//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
863//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
864//                              if(WaveExists(mask) == 1)
865//                                      maskMissing = 0
866//                              endif
867//                      endif   
868//                      //TODO:
869//                      // -- decide on the proper deltaQ for binning. either nominal value for LR, or one
870//                      //    determined specifically for that panel (currently using one tube width as deltaQ)
871//                      // -- this is repeated multiple times in this switch
872//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
873//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
874//                      nSets = 1
875//                      break   
876                                       
877//              case "MT":             
878//              case "MB":
879//                      detStr = type
880//                      if(isVCALC)
881//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
882//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
883//                      else
884//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
885//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
886//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
887//                              if(WaveExists(mask) == 1)
888//                                      maskMissing = 0
889//                              endif
890//                      endif   
891//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
892//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
893//                      nSets = 1
894//                      break   
895
896// only one panel, simply pick that panel and move on out of the switch
897                case "FL":
898                case "FR":
899                case "FT":
900                case "FB":
901                case "ML":
902                case "MR":
903                case "MT":
904                case "MB":                     
905                case "B":       
906                        if(isVCALC)
907                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
908                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
909                        else
910                                Wave inten = V_getDetectorDataW(folderStr,detStr)
911                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
912                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
913                                if(WaveExists(mask) == 1)
914                                        maskMissing = 0
915                                endif
916                        endif   
917                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
918                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values 
919                        nSets = 1
920                        break   
921                       
922                case "FLR":
923                // detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
924                // TODO
925                // -- see if I can un-hard-wire some of this below when more than one panel is combined
926                        if(isVCALC)
927                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
928                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
929                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
930                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
931                        else
932                                Wave inten = V_getDetectorDataW(folderStr,"FL")
933                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
934                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
935                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
936                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
937                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
938                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
939                                        maskMissing = 0
940                                endif
941                        endif   
942                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
943                       
944                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
945                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
946               
947                        nSets = 2
948                        break                   
949               
950                case "FTB":
951                        if(isVCALC)
952                                WAVE inten = $(folderPath+instPath+"FT"+":det_"+"FT")
953                                WAVE/Z iErr = $("iErr_"+"FT")                   // 2D errors -- may not exist, especially for simulation               
954                                WAVE inten2 = $(folderPath+instPath+"FB"+":det_"+"FB")
955                                WAVE/Z iErr2 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
956                        else
957                                Wave inten = V_getDetectorDataW(folderStr,"FT")
958                                Wave iErr = V_getDetectorDataErrW(folderStr,"FT")
959                                Wave inten2 = V_getDetectorDataW(folderStr,"FB")
960                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FB")
961                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
962                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
963                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
964                                        maskMissing = 0
965                                endif
966                        endif   
967                        NVAR delQ = $(folderPath+instPath+"FT"+":gDelQ_FT")
968                       
969                        Wave qTotal = $(folderPath+instPath+"FT"+":qTot_"+"FT")                 // 2D q-values 
970                        Wave qTotal2 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
971       
972                        nSets = 2
973                        break           
974               
975                case "FLRTB":
976                        if(isVCALC)
977                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
978                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
979                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
980                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
981                                WAVE inten3 = $(folderPath+instPath+"FT"+":det_"+"FT")
982                                WAVE/Z iErr3 = $("iErr_"+"FT")                  // 2D errors -- may not exist, especially for simulation               
983                                WAVE inten4 = $(folderPath+instPath+"FB"+":det_"+"FB")
984                                WAVE/Z iErr4 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
985                        else
986                                Wave inten = V_getDetectorDataW(folderStr,"FL")
987                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
988                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
989                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
990                                Wave inten3 = V_getDetectorDataW(folderStr,"FT")
991                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"FT")
992                                Wave inten4 = V_getDetectorDataW(folderStr,"FB")
993                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"FB")
994                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
995                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
996                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
997                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
998                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
999                                        maskMissing = 0
1000                                endif
1001                        endif   
1002                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
1003                       
1004                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
1005                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
1006                        Wave qTotal3 = $(folderPath+instPath+"FT"+":qTot_"+"FT")                        // 2D q-values 
1007                        Wave qTotal4 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
1008               
1009                        nSets = 4
1010                        break           
1011                       
1012                case "MLR":
1013                        if(isVCALC)
1014                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
1015                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
1016                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
1017                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
1018                        else
1019                                Wave inten = V_getDetectorDataW(folderStr,"ML")
1020                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
1021                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
1022                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
1023                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
1024                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
1025                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1026                                        maskMissing = 0
1027                                endif
1028                        endif   
1029                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
1030                       
1031                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
1032                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
1033               
1034                        nSets = 2
1035                        break                   
1036               
1037                case "MTB":
1038                        if(isVCALC)
1039                                WAVE inten = $(folderPath+instPath+"MT"+":det_"+"MT")
1040                                WAVE/Z iErr = $("iErr_"+"MT")                   // 2D errors -- may not exist, especially for simulation               
1041                                WAVE inten2 = $(folderPath+instPath+"MB"+":det_"+"MB")
1042                                WAVE/Z iErr2 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
1043                        else
1044                                Wave inten = V_getDetectorDataW(folderStr,"MT")
1045                                Wave iErr = V_getDetectorDataErrW(folderStr,"MT")
1046                                Wave inten2 = V_getDetectorDataW(folderStr,"MB")
1047                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MB")
1048                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
1049                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
1050                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1051                                        maskMissing = 0
1052                                endif
1053                        endif   
1054                        NVAR delQ = $(folderPath+instPath+"MT"+":gDelQ_MT")
1055                       
1056                        Wave qTotal = $(folderPath+instPath+"MT"+":qTot_"+"MT")                 // 2D q-values 
1057                        Wave qTotal2 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
1058               
1059                        nSets = 2
1060                        break                           
1061               
1062                case "MLRTB":
1063                        if(isVCALC)
1064                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
1065                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
1066                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
1067                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
1068                                WAVE inten3 = $(folderPath+instPath+"MT"+":det_"+"MT")
1069                                WAVE/Z iErr3 = $("iErr_"+"MT")                  // 2D errors -- may not exist, especially for simulation               
1070                                WAVE inten4 = $(folderPath+instPath+"MB"+":det_"+"MB")
1071                                WAVE/Z iErr4 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
1072                        else
1073                                Wave inten = V_getDetectorDataW(folderStr,"ML")
1074                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
1075                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
1076                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
1077                                Wave inten3 = V_getDetectorDataW(folderStr,"MT")
1078                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"MT")
1079                                Wave inten4 = V_getDetectorDataW(folderStr,"MB")
1080                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"MB")
1081                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
1082                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
1083                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
1084                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
1085                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
1086                                        maskMissing = 0
1087                                endif
1088                        endif   
1089                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
1090                       
1091                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
1092                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
1093                        Wave qTotal3 = $(folderPath+instPath+"MT"+":qTot_"+"MT")                        // 2D q-values 
1094                        Wave qTotal4 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
1095               
1096                        nSets = 4
1097                        break                                                                   
1098                                       
1099                default:
1100                        nSets = 0                                                       
1101                        Print "ERROR   ---- type is not recognized "
1102        endswitch
1103
1104//      Print "delQ = ",delQ," for ",type
1105
1106        if(nSets == 0)
1107                SetDataFolder root:
1108                return(0)
1109        endif
1110
1111
1112// RAW data is currently read in and the 2D error wave is correctly generated
1113// 2D error is propagated through all reduction steps, but I have not
1114// verified that it is an exact duplication of the 1D error
1115//
1116//
1117//
1118// IF ther is no 2D error wave present for some reason, make a fake one
1119        if(WaveExists(iErr)==0  && WaveExists(inten) != 0)
1120                Duplicate/O inten,iErr
1121                Wave iErr=iErr
1122//              iErr = 1+sqrt(inten+0.75)                       // can't use this -- it applies to counts, not intensity (already a count rate...)
1123                iErr = sqrt(inten+0.75)                 // TODO -- here I'm just using some fictional value
1124        endif
1125        if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0)
1126                Duplicate/O inten2,iErr2
1127                Wave iErr2=iErr2
1128//              iErr2 = 1+sqrt(inten2+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1129                iErr2 = sqrt(inten2+0.75)                       // TODO -- here I'm just using some fictional value
1130        endif
1131        if(WaveExists(iErr3)==0  && WaveExists(inten3) != 0)
1132                Duplicate/O inten3,iErr3
1133                Wave iErr3=iErr3
1134//              iErr3 = 1+sqrt(inten3+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1135                iErr3 = sqrt(inten3+0.75)                       // TODO -- here I'm just using some fictional value
1136        endif
1137        if(WaveExists(iErr4)==0  && WaveExists(inten4) != 0)
1138                Duplicate/O inten4,iErr4
1139                Wave iErr4=iErr4
1140//              iErr4 = 1+sqrt(inten4+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1141                iErr4 = sqrt(inten4+0.75)                       // TODO -- here I'm just using some fictional value
1142        endif
1143
1144        // TODO -- nq will need to be larger, once the back detector is installed
1145        //
1146        // note that the back panel of 320x320 (1mm res) results in 447 data points!
1147        // - so I upped nq to 600
1148
1149        nq = 600
1150
1151//******TODO****** -- where to put the averaged data -- right now, folderStr is forced to ""   
1152//      SetDataFolder $("root:"+folderStr)              //should already be here, but make sure...     
1153        Make/O/D/N=(nq)  $(folderPath+":"+"iBin_qxqy"+"_"+type)
1154        Make/O/D/N=(nq)  $(folderPath+":"+"qBin_qxqy"+"_"+type)
1155        Make/O/D/N=(nq)  $(folderPath+":"+"nBin_qxqy"+"_"+type)
1156        Make/O/D/N=(nq)  $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1157        Make/O/D/N=(nq)  $(folderPath+":"+"eBin_qxqy"+"_"+type)
1158        Make/O/D/N=(nq)  $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1159       
1160        Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+type)
1161        Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type)
1162        Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy"+"_"+type)
1163        Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1164        Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy"+"_"+type)
1165        Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1166       
1167       
1168//      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
1169// TODO: not sure if I want to set dQ in x or y direction...
1170        // the short dimension is the 8mm tubes, use this direction as dQ?
1171        // 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]
1172        // WRT the beam center. use qx or qy directly. Still not happy with this way...
1173
1174
1175        qBin_qxqy[] =  p*delQ   
1176        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
1177
1178        iBin_qxqy = 0
1179        iBin2_qxqy = 0
1180        eBin_qxqy = 0
1181        eBin2D_qxqy = 0
1182        nBin_qxqy = 0   //number of intensities added to each bin
1183
1184// now there are situations of:
1185// 1 panel
1186// 2 panels
1187// 4 panels
1188//
1189// this needs to be a double loop now...
1190// TODO:
1191// -- the iErr (=2D) wave and accumulation of error is NOT CALCULATED CORRECTLY YET
1192// -- verify the 2D error propagation by reducing it to 1D error
1193//
1194//
1195// The 1D error does not use iErr, and IS CALCULATED CORRECTLY
1196//
1197// x- the solid angle per pixel will be present for WORK data other than RAW, but not for RAW
1198
1199//
1200// if any of the masks don't exist, display the error, and proceed with the averaging, using all data
1201        if(maskMissing == 1)
1202                Print "Mask file not found for at least one detector - so all data is used"
1203        endif
1204
1205        Variable mask_val
1206// use set 1 (no number) only
1207        if(nSets >= 1)
1208                xDim=DimSize(inten,0)
1209                yDim=DimSize(inten,1)
1210       
1211                for(ii=0;ii<xDim;ii+=1)
1212                        for(jj=0;jj<yDim;jj+=1)
1213                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1214                                qVal = qTotal[ii][jj]
1215                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1216                                val = inten[ii][jj]
1217                               
1218                                if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
1219                                        mask_val = 0
1220                                else
1221                                        mask_val = mask[ii][jj]
1222                                endif
1223                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1224                                        iBin_qxqy[binIndex] += val
1225                                        iBin2_qxqy[binIndex] += val*val
1226                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1227                                        nBin_qxqy[binIndex] += 1
1228                                endif
1229                        endfor
1230                endfor
1231               
1232        endif
1233
1234// add in set 2 (set 1 already done)
1235        if(nSets >= 2)
1236                xDim=DimSize(inten2,0)
1237                yDim=DimSize(inten2,1)
1238       
1239                for(ii=0;ii<xDim;ii+=1)
1240                        for(jj=0;jj<yDim;jj+=1)
1241                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1242                                qVal = qTotal2[ii][jj]
1243                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1244                                val = inten2[ii][jj]
1245                               
1246                                if(isVCALC || maskMissing)
1247                                        mask_val = 0
1248                                else
1249                                        mask_val = mask2[ii][jj]
1250                                endif
1251                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1252                                        iBin_qxqy[binIndex] += val
1253                                        iBin2_qxqy[binIndex] += val*val
1254                                        eBin2D_qxqy[binIndex] += iErr2[ii][jj]*iErr2[ii][jj]
1255                                        nBin_qxqy[binIndex] += 1
1256                                endif
1257                        endfor
1258                endfor
1259               
1260        endif
1261
1262// add in set 3 and 4 (set 1 and 2 already done)
1263        if(nSets == 4)
1264                xDim=DimSize(inten3,0)
1265                yDim=DimSize(inten3,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 = qTotal3[ii][jj]
1271                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1272                                val = inten3[ii][jj]
1273                               
1274                                if(isVCALC || maskMissing)
1275                                        mask_val = 0
1276                                else
1277                                        mask_val = mask3[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] += iErr3[ii][jj]*iErr3[ii][jj]
1283                                        nBin_qxqy[binIndex] += 1
1284                                endif
1285                        endfor
1286                endfor
1287               
1288               
1289                xDim=DimSize(inten4,0)
1290                yDim=DimSize(inten4,1)
1291       
1292                for(ii=0;ii<xDim;ii+=1)
1293                        for(jj=0;jj<yDim;jj+=1)
1294                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1295                                qVal = qTotal4[ii][jj]
1296                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1297                                val = inten4[ii][jj]
1298                               
1299                                if(isVCALC || maskMissing)
1300                                        mask_val = 0
1301                                else
1302                                        mask_val = mask4[ii][jj]
1303                                endif
1304                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1305                                        iBin_qxqy[binIndex] += val
1306                                        iBin2_qxqy[binIndex] += val*val
1307                                        eBin2D_qxqy[binIndex] += iErr4[ii][jj]*iErr4[ii][jj]
1308                                        nBin_qxqy[binIndex] += 1
1309                                endif
1310                        endfor
1311                endfor
1312               
1313        endif
1314
1315
1316// after looping through all of the data on the panels, calculate errors on I(q),
1317// just like in CircSectAve.ipf
1318// TODO:
1319// -- 2D Errors were (maybe) properly acculumated through reduction, so this loop of calculations is NOT VERIFIED (yet)
1320// x- the error on the 1D intensity, is correctly calculated as the standard error of the mean.
1321        for(ii=0;ii<nq;ii+=1)
1322                if(nBin_qxqy[ii] == 0)
1323                        //no pixels in annuli, data unknown
1324                        iBin_qxqy[ii] = 0
1325                        eBin_qxqy[ii] = 1
1326                        eBin2D_qxqy[ii] = NaN
1327                else
1328                        if(nBin_qxqy[ii] <= 1)
1329                                //need more than one pixel to determine error
1330                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1331                                eBin_qxqy[ii] = 1
1332                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1333                        else
1334                                //assume that the intensity in each pixel in annuli is normally distributed about mean...
1335                                //  -- this is correctly calculating the error as the standard error of the mean, as
1336                                //    was always done for SANS as well.
1337                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1338                                avesq = iBin_qxqy[ii]^2
1339                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1340                                var = aveisq-avesq
1341                                if(var<=0)
1342                                        eBin_qxqy[ii] = 1e-6
1343                                else
1344                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1345                                endif
1346                                // and calculate as it is propagated pixel-by-pixel
1347                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1348                        endif
1349                endif
1350        endfor
1351       
1352        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo
1353       
1354        // find the last non-zero point, working backwards
1355        val=nq
1356        do
1357                val -= 1
1358        while((nBin_qxqy[val] == 0) && val > 0)
1359       
1360//      print val, nBin_qxqy[val]
1361        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1362
1363        if(val == 0)
1364                // all the points were deleted
1365                return(0)
1366        endif
1367       
1368       
1369        // since the beam center is not always on the detector, many of the low Q bins will have zero pixels
1370        // find the first non-zero point, working forwards
1371        val = -1
1372        do
1373                val += 1
1374        while(nBin_qxqy[val] == 0)     
1375        DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1376
1377        // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
1378        val = numpnts(nBin_qxqy)-1
1379        do
1380                if(nBin_qxqy[val] == 0)
1381                        DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1382                endif
1383                val -= 1
1384        while(val>0)
1385       
1386        // TODO:
1387        // -- This is where I calculate the resolution in SANS (see CircSectAve)
1388        // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
1389        // -- from the top of the function, folderStr = work folder, type = "FLRTB" or other type of averaging
1390        //
1391        nq = numpnts(qBin_qxqy)
1392        Make/O/D/N=(nq)  $(folderPath+":"+"sigmaQ_qxqy"+"_"+type)
1393        Make/O/D/N=(nq)  $(folderPath+":"+"qBar_qxqy"+"_"+type)
1394        Make/O/D/N=(nq)  $(folderPath+":"+"fSubS_qxqy"+"_"+type)
1395        Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+type)
1396        Wave qbar = $(folderPath+":"+"qBar_qxqy_"+type)
1397        Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+type)
1398                               
1399
1400        Variable ret1,ret2,ret3
1401        Variable lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses
1402
1403// TODO: check the units of all of the inputs
1404
1405// lambda = wavelength [A]
1406        lambda = V_getWavelength(folderStr)
1407       
1408// lambdaWidth = [dimensionless]
1409        lambdaWidth = V_getWavelength_spread(folderStr)
1410       
1411// DDet = detector pixel resolution [cm]        **assumes square pixel
1412        DDet = 0.8              // TODO -- this is hard-wired
1413       
1414// apOff = sample aperture to sample distance [cm]
1415        apOff = 10              // TODO -- this is hard-wired
1416       
1417// S1 = source aperture diameter [mm]
1418        S1 = str2num(V_getSourceAp_size(folderStr))
1419       
1420// S2 = sample aperture diameter [cm]
1421// as of 3/2018, the "internal" sample aperture is not in use, only the external
1422// TODO : verify the units on the Ap2 (external)
1423// sample aperture 1(internal) is set to report "12.7 mm" as a STRING
1424// sample aperture 2(external) reports the number typed in...
1425//
1426// so I'm trusting [cm] is in the file
1427        S2 = V_getSampleAp2_size(folderStr)*10          // sample ap 1 or 2? 2 = the "external", convert to [cm]
1428       
1429// L1 = source to sample distance [m]
1430        L1 = V_getSourceAp_distance(folderStr)/100
1431
1432// L2 = sample to detector distance [m]
1433// take the first two characters of the "type" to get the correct distance.
1434// if the type is say, MLRTB, then the implicit assumption in combining all four panels is that the resolution
1435// is not an issue for the slightly different distances.
1436        L2 = V_getDet_ActualDistance(folderStr,type[0,1])/100           //convert cm to m
1437       
1438// BS = beam stop diameter [mm]
1439        //BS = V_getBeamStopC2_size(folderStr)          // TODO: what are the units? which BS is in? carr2, carr3, back, none?
1440        BS = 25.4                       //TODO hard-wired value
1441       
1442// del_r = step size [mm] = binWidth*(mm/pixel)
1443        del_r = 1*8                     // TODO: 8mm/pixel hard-wired
1444
1445// usingLenses = flag for lenses = 0 if no lenses, non-zero if lenses are in-beam
1446        usingLenses = 0
1447
1448if(cmpstr(detStr,"FL")==0)
1449        Print "(FL) Resolution lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses"
1450        Print lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses
1451endif
1452
1453
1454// TODO:
1455// this is the point where I need to switch on the different collimation types (white beam, slit, Xtal, etc)
1456// to calculate the correct resolution, or fill the waves with the correct "flags"
1457//
1458
1459// For white beam data, the wavelength distribution can't be represented as a gaussian, but all of the other
1460//  geometric corrections still apply. Passing zero for the lambdaWidth will return the grometry contribution,
1461//  as long as the wavelength can be handled separately. It appears to be correct to do as a double integral,
1462// with the inner(lambda) calculated first, then the outer(geometry).
1463//
1464
1465
1466        ii=0
1467        do
1468                V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3)
1469                sigmaq[ii] = ret1       
1470                qbar[ii] = ret2
1471                fsubs[ii] = ret3       
1472                ii+=1
1473        while(ii<nq)
1474       
1475       
1476       
1477        SetDataFolder root:
1478       
1479        return(0)
1480End
1481
1482
Note: See TracBrowser for help on using the repository browser.