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

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

LOTS of changes to accommodate the beam center being reported in cm rather than pixels. Required a lot of changes to VCALC (to fill in simulated data), and in the reading and interpreting of data for display, and most importantly, the calculation of q.

There may still be a few residual bugs with this. I am still re-testing with new sample data sets.

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