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

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

minor changes, fast ways to patch mis-information in the headers has been added to the VSANS menu for easier access.

added comments here and there

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