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

Last change on this file since 999 was 999, checked in by srkline, 7 years ago

changes to a few analysis models to make these Igor 7-ready

adding mask editing utilities

many changes to event mode for easier processing of split lists

updated event mode help file

+ lots more!

File size: 39.4 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// TODO: hard wired for a sphere - change this to allow minimal selections and altering of coefficients
16// TODO: 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        String funcStr = VCALC_getModelFunctionStr()
53        strswitch(funcStr)
54                case "Big Debye":
55                        tmpInten = VC_Debye(100,3000,0.0001,qTot[p][q])
56                        break
57                case "Big Sphere":
58                        tmpInten = VC_SphereForm(1,900,1e-6,0.01,qTot[p][q])   
59                        break
60                case "Debye":
61                        tmpInten = VC_Debye(10,300,0.0001,qTot[p][q])
62                        break
63                case "Sphere":
64                        tmpInten = VC_SphereForm(1,60,1e-6,0.001,qTot[p][q])   
65                        break
66                case "AgBeh":
67                        tmpInten = VC_BroadPeak(1e-9,3,20,100.0,0.1,3,0.1,qTot[p][q])
68                        break
69                case "Vycor":
70                        tmpInten = VC_BroadPeak(1e-9,3,20,500.0,0.015,3,0.1,qTot[p][q])
71                        break   
72                case "Empty Cell":
73                        tmpInten = VC_EC_Empirical(2.2e-8,3.346,0.0065,9.0,0.016,qTot[p][q])
74                        break
75                case "Blocked Beam":
76                        tmpInten = VC_BlockedBeam(1,qTot[p][q])
77                        break
78                default:
79                        tmpInten = VC_Debye(10,300,0.1,qTot[p][q])
80        endswitch
81
82// TODO: this is faked to get around the singularity at the center of the back detector
83//
84//
85        if(cmpstr(type,"B") == 0)
86                Variable nx,ny,px,py
87                nx = VCALC_get_nPix_X(type)
88                ny = VCALC_get_nPix_Y(type)
89                px = trunc(nx/2)
90                py = trunc(ny/2)
91               
92                tmpInten[px][py] = (tmpInten[px][py+1] + tmpInten[px][py-1])/2
93        endif
94
95
96
97///////////////
98//      // calculate the scattering cross section simply to be able to estimate the transmission
99//      Variable sig_sas=0
100//     
101//      // remember that the random deviate is the coherent portion ONLY - the incoherent background is
102//      // subtracted before the calculation.
103//      CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
104//
105//      if(sig_sas > 100)
106//              DoAlert 0,"SAS cross section > 100. Estimates of multiple scattering are unreliable. Choosing a model with a well-defined Rg may help"
107//      endif           
108//
109//      // calculate the multiple scattering fraction for display (10/2009)
110//      Variable ii,nMax=10,tau
111//      mScat=0
112//      tau = thick*sig_sas
113//      // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
114//      // neutrons out of those that were scattered
115//      for(ii=2;ii<nMax;ii+=1)
116//              mScat += tau^(ii)/factorial(ii)
117////            print tau^(ii)/factorial(ii)
118//      endfor
119//      estTrans = exp(-1*thick*sig_sas)                //thickness and sigma both in units of cm
120//      mscat *= (estTrans)/(1-estTrans)
121//
122////    if(mScat > 0.1)         //  Display warning
123//
124//      Print "Sig_sas = ",sig_sas
125////////////////////
126       
127        prob_i = trans*thick*pixSizeX*pixSizeY/(sdd)^2*tmpInten                 //probability of a neutron in q-bin(i)
128               
129        tmpInten = (imon)*prob_i                //tmpInten is not the model calculation anymore!!
130
131
132/// **** can I safely assume a Gaussian error in the count rate??
133        tmpSig = sqrt(tmpInten)         // corrected based on John's memo, from 8/9/99
134
135        tmpInten += gnoise(tmpSig)
136        tmpInten = (tmpInten[p][q] < 0) ? 0 : tmpInten[p][q]                    // MAR 2013 -- is this the right thing to do
137        tmpInten = trunc(tmpInten)
138               
139       
140        det = tmpInten
141
142// if I want "absolute" scale -- then I lose the integer nature of the detector (but keep the random)
143//      det /= trans*thick*pixSizeX*pixSizeY/(sdd)^2*imon
144
145       
146        KillWaves/Z tmpInten,tmpSig,prob_i     
147        SetDataFolder root:
148
149        return(0)
150End
151
152
153// For a given detector panel, calculate the q-values
154// -work with everything as arrays
155// Input needed:
156// detector data
157// detector type (LRTB?)
158// beam center (may be off the detector)
159// SDD
160// lambda
161//
162// pixel dimensions for detector type (global constants)
163// - data dimensions read directly from array
164//
165// --What is calculated:
166// array of Q
167// array of qx,qy,qz
168// array of error already exists
169//
170//
171// -- sdd in meters
172// -- lambda in Angstroms
173Function VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
174        Wave data,qTot,qx,qy,qz
175        Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
176               
177        // loop over the array and calculate the values - this is done as a wave assignment
178// TODO -- be sure that it's p,q -- or maybe p+1,q+1 as used in WriteQIS.ipf   
179        qTot = VC_CalcQval(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
180        qx = VC_CalcQX(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
181        qy = VC_CalcQY(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
182        qz = VC_CalcQZ(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
183       
184        return(0)
185End
186
187
188// for testing, a version that will calculate the q-arrays for VCALC based on whatever nonlinear coefficients
189// exist in the RAW data folder
190//
191// reverts to the "regular" linear detector if waves not found or a flag is set
192//
193// need to convert the beam center to mm
194// need to call the VSANS V_CalcQval routines (these use the real-space distance, not pixel dims)
195//
196// TODO:
197// -- tube width is hard-wired in
198//
199//
200Function VC_Detector_2Q_NonLin(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY,detStr)
201        Wave data,qTot,qx,qy,qz
202        Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
203        String detStr
204       
205        String destPath = "root:Packages:NIST:VSANS:RAW"
206                       
207        Wave/Z data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
208        Wave/Z data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
209        NVAR gUseNonLinearDet = root:Packages:NIST:VSANS:VCALC:gUseNonLinearDet
210       
211        if(gUseNonLinearDet && WaveExists(data_realDistX) && WaveExists(data_realDistY))
212                // convert the beam centers to mm
213                String orientation
214                Variable dimX,dimY,newX,newY
215                dimX = DimSize(data_realDistX,0)
216                dimY = DimSize(data_realDistX,1)
217                if(dimX > dimY)
218                        orientation = "horizontal"
219                else
220                        orientation = "vertical"
221                endif
222               
223                Variable tube_width = 8.4               //mm
224       
225        //
226                if(cmpstr(orientation,"vertical")==0)
227                        //      this is data dimensioned as (Ntubes,Npix)
228                        newX = tube_width*xCtr
229                        newY = data_realDistY[0][yCtr]
230                else
231                        //      this is data (horizontal) dimensioned as (Npix,Ntubes)
232                        newX = data_realDistX[xCtr][0]
233                        newY = tube_width*yCtr
234                endif   
235
236                //if detector "B", different calculation for the centers (not tubes)
237                if(cmpstr(detStr,"B")==0)
238                        newX = data_realDistX[xCtr][0]
239                        newY = data_realDistY[0][yCtr]
240                        //newX = xCtr
241                        //newY = yCtr
242                endif           
243                               
244                // calculate all of the q-values
245                qTot = V_CalcQval(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
246                qx = V_CalcQX(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
247                qy = V_CalcQY(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
248                qz = V_CalcQZ(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
249       
250//              Print "det, x_mm, y_mm ",detStr,num2str(newX),num2str(newY)
251//              Print "det, x_pix, y_pix ",detStr,num2str(xCtr),num2str(yCtr)
252        else
253                // do the q-calculation using linear detector
254                VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
255        endif
256       
257        return(0)
258End
259
260
261//////////////////////
262// NOTE: The Q calculations are different than what is in GaussUtils in that they take into
263// accout the different x/y pixel sizes and the beam center not being on the detector -
264// off a different edge for each LRTB type
265/////////////////////
266
267//function to calculate the overall q-value, given all of the necesary trig inputs
268//and are in detector coordinates (1,128) rather than axis values
269//the pixel locations need not be integers, reals are ok inputs
270//sdd is in meters
271//wavelength is in Angstroms
272//
273//returned magnitude of Q is in 1/Angstroms
274//
275Function VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
276        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
277       
278        Variable dx,dy,qval,two_theta,dist
279               
280        sdd *=100               //convert to cm
281        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
282        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
283        dist = sqrt(dx^2 + dy^2)
284       
285        two_theta = atan(dist/sdd)
286
287        qval = 4*Pi/lam*sin(two_theta/2)
288       
289        return qval
290End
291
292//calculates just the q-value in the x-direction on the detector
293//input/output is the same as CalcQval()
294//ALL inputs are in detector coordinates
295//
296//sdd is in meters
297//wavelength is in Angstroms
298//
299// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
300// now properly accounts for qz
301//
302Function VC_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
303        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
304
305        Variable qx,qval,phi,dx,dy,dist,two_theta
306       
307        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
308       
309        sdd *=100               //convert to cm
310        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
311        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
312        phi = V_FindPhi(dx,dy)
313       
314        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
315        dist = sqrt(dx^2 + dy^2)
316        two_theta = atan(dist/sdd)
317
318        qx = qval*cos(two_theta/2)*cos(phi)
319       
320        return qx
321End
322
323//calculates just the q-value in the y-direction on the detector
324//input/output is the same as CalcQval()
325//ALL inputs are in detector coordinates
326//sdd is in meters
327//wavelength is in Angstroms
328//
329// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
330// now properly accounts for qz
331//
332Function VC_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
333        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
334       
335        Variable dy,qval,dx,phi,qy,dist,two_theta
336       
337        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
338       
339        sdd *=100               //convert to cm
340        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
341        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
342        phi = V_FindPhi(dx,dy)
343       
344        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
345        dist = sqrt(dx^2 + dy^2)
346        two_theta = atan(dist/sdd)
347       
348        qy = qval*cos(two_theta/2)*sin(phi)
349       
350        return qy
351End
352
353//calculates just the z-component of the q-vector, not measured on the detector
354//input/output is the same as CalcQval()
355//ALL inputs are in detector coordinates
356//sdd is in meters
357//wavelength is in Angstroms
358//
359// not actually used, but here for completeness if anyone asks
360//
361Function VC_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
362        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
363       
364        Variable dy,qval,dx,phi,qz,dist,two_theta
365       
366        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
367       
368        sdd *=100               //convert to cm
369       
370        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
371        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
372        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
373        dist = sqrt(dx^2 + dy^2)
374        two_theta = atan(dist/sdd)
375       
376        qz = qval*sin(two_theta/2)
377       
378        return qz
379End
380
381//phi is defined from +x axis, proceeding CCW around [0,2Pi]
382Threadsafe Function V_FindPhi(vx,vy)
383        variable vx,vy
384       
385        variable phi
386       
387        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
388       
389        // special cases
390        if(vx==0 && vy > 0)
391                return(pi/2)
392        endif
393        if(vx==0 && vy < 0)
394                return(3*pi/2)
395        endif
396        if(vx >= 0 && vy == 0)
397                return(0)
398        endif
399        if(vx < 0 && vy == 0)
400                return(pi)
401        endif
402       
403       
404        if(vx > 0 && vy > 0)
405                return(phi)
406        endif
407        if(vx < 0 && vy > 0)
408                return(phi + pi)
409        endif
410        if(vx < 0 && vy < 0)
411                return(phi + pi)
412        endif
413        if( vx > 0 && vy < 0)
414                return(phi + 2*pi)
415        endif
416       
417        return(phi)
418end
419
420Function VC_SphereForm(scale,radius,delrho,bkg,x)                               
421        Variable scale,radius,delrho,bkg
422        Variable x
423       
424        // variables are:                                                       
425        //[0] scale
426        //[1] radius (A)
427        //[2] delrho (A-2)
428        //[3] background (cm-1)
429       
430//      Variable scale,radius,delrho,bkg                               
431//      scale = w[0]
432//      radius = w[1]
433//      delrho = w[2]
434//      bkg = w[3]
435       
436       
437        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
438        // and is rescaled to give [=] cm^-1
439       
440        Variable bes,f,vol,f2
441        //
442        //handle q==0 separately
443        If(x==0)
444                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
445                return(f)
446        Endif
447       
448//      bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
449       
450        bes = 3*sqrt(pi/(2*x*radius))*BesselJ(1.5,x*radius)/(x*radius)
451       
452        vol = 4*pi/3*radius^3
453        f = vol*bes*delrho              // [=] A
454        // normalize to single particle volume, convert to 1/cm
455        f2 = f * f / vol * 1.0e8                // [=] 1/cm
456       
457        return (scale*f2+bkg)   // Scale, then add in the background
458       
459End
460
461Function VC_Debye(scale,rg,bkg,x)
462        Variable scale,rg,bkg
463        Variable x
464       
465        // variables are:
466        //[0] scale factor
467        //[1] radius of gyration [A]
468        //[2] background        [cm-1]
469       
470        // calculates (scale*debye)+bkg
471        Variable Pq,qr2
472       
473        qr2=(x*rg)^2
474        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
475       
476        //scale
477        Pq *= scale
478        // then add in the background
479        return (Pq+bkg)
480End
481
482// a sum of a power law and debye to approximate the scattering from a real empty cell
483//
484//      make/O/D coef_ECEmp = {2.2e-8,3.346,0.0065,9.0,0.016}
485//
486Function VC_EC_Empirical(aa,mm,scale,rg,bkg,x)
487        Variable aa,mm,scale,rg,bkg
488        Variable x
489       
490        // variables are:
491        //[0] = A
492        //[1] = power m
493        //[2] scale factor
494        //[3] radius of gyration [A]
495        //[4] background        [cm-1]
496       
497        Variable Iq
498       
499        // calculates (scale*debye)+bkg
500        Variable Pq,qr2
501       
502//      if(x*Rg < 1e-3)         //added Oct 2008 to avoid numerical errors at low arg values
503//              return(scale+bkg)
504//      endif
505       
506        Iq = aa*x^-mm
507       
508        qr2=(x*rg)^2
509        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
510       
511        //scale
512        Pq *= scale
513        // then add the terms up
514        return (Iq + Pq + bkg)
515End
516
517// blocked beam
518//
519Function VC_BlockedBeam(bkg,x)
520        Variable bkg
521        Variable x
522       
523        return (bkg)
524End
525
526
527//
528// a broad peak to simulate silver behenate or vycor
529//
530// peak @ 0.1 ~ AgBeh
531//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}               
532//
533//
534// peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
535//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}             
536//
537//
538Function VC_BroadPeak(aa,nn,cc,LL,Qzero,mm,bgd,x)
539        Variable aa,nn,cc,LL,Qzero,mm,bgd
540        Variable x
541       
542        // variables are:                                                       
543        //[0] Porod term scaling
544        //[1] Porod exponent
545        //[2] Lorentzian term scaling
546        //[3] Lorentzian screening length [A]
547        //[4] peak location [1/A]
548        //[5] Lorentzian exponent
549        //[6] background
550       
551//      local variables
552        Variable inten, qval
553//      x is the q-value for the calculation
554        qval = x
555//      do the calculation and return the function value
556       
557        inten = aa/(qval)^nn + cc/(1 + (abs(qval-Qzero)*LL)^mm) + bgd
558
559        Return (inten)
560       
561End
562
563//
564// updated to new folder structure Feb 2016
565// folderStr = RAW,SAM, VCALC or other
566// detStr is the panel identifer "ML", etc.
567//
568Function SetDeltaQ(folderStr,detStr)
569        String folderStr,detStr
570
571        Variable isVCALC
572        if(cmpstr(folderStr,"VCALC") == 0)
573                isVCALC = 1
574        endif
575       
576        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
577        String instPath = ":entry:instrument:detector_"
578               
579        if(isVCALC)
580                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)               // 2D detector data
581        else
582                Wave inten = V_getDetectorDataW(folderStr,detStr)
583        endif
584
585        Wave qx = $(folderPath+instPath+detStr+":qx_"+detStr)
586        Wave qy = $(folderPath+instPath+detStr+":qy_"+detStr)
587       
588        Variable xDim,yDim,delQ
589       
590        xDim=DimSize(inten,0)
591        yDim=DimSize(inten,1)
592       
593        if(xDim<yDim)
594                delQ = abs(qx[0][0] - qx[1][0])/2
595        else
596                delQ = abs(qy[0][1] - qy[0][0])/2
597        endif
598       
599        // set the global
600        Variable/G $(folderPath+instPath+detStr+":gDelQ_"+detStr) = delQ
601//      Print "SET delQ = ",delQ," for ",type
602       
603        return(delQ)
604end
605
606
607//TODO -- need a switch here to dispatch to the averaging type
608Proc VC_BinQxQy_to_1D(folderStr,type)
609        String folderStr
610        String type
611//      Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4)
612//      Prompt type,"detector identifier"
613
614
615        VC_fDoBinning_QxQy2D(folderStr, type)
616
617
618/// this is for a tall, narrow slit mode       
619//      VC_fBinDetector_byRows(folderStr,type)
620       
621End
622
623
624// folderStr is RAW, VCALC, SAM, etc.
625// type is "B", "FL" for single binning, "FLR", or "MLRTB" or similar if multiple panels are combined
626//
627Proc VC_Graph_1D_detType(folderStr,type)
628        String folderStr,type
629       
630        SetDataFolder $("root:Packages:NIST:VSANS:"+folderStr)
631       
632        Display $("iBin_qxqy"+"_"+type) vs $("qBin_qxqy"+"_"+type)
633        ModifyGraph mirror=2,grid=1,log=1
634        ModifyGraph mode=4,marker=19,msize=2
635//      ErrorBars/T=0 iBin_qxqy Y,wave=(eBin2D_qxqy,eBin2D_qxqy)                // for simulations, I don't have 2D uncertainty
636        ErrorBars/T=0 $("iBin_qxqy"+"_"+type) Y,wave=($("eBin_qxqy"+"_"+type),$("eBin_qxqy"+"_"+type))
637        legend
638       
639        SetDataFolder root:
640
641End
642
643
644
645//////////
646//
647//              Function that bins a 2D detctor panel into I(q) based on the q-value of the pixel
648//              - each pixel QxQyQz has been calculated beforehand
649//              - if multiple panels are selected to be combined, it is done here during the binning
650//              - the setting of deltaQ step is still a little suspect (TODO)
651//
652//
653// see the equivalent function in PlotUtils2D_v40.ipf
654//
655//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
656//
657// this has been modified to accept different detector panels and to take arrays
658// -- type = FL or FR or...other panel identifiers
659//
660// TODO "iErr" is not always defined correctly since it doesn't really apply here for data that is not 2D simulation
661//
662//
663// updated Feb2016 to take new folder structure
664// TODO
665// -- VERIFY
666// -- figure out what the best location is to put the averaged data? currently @ top level of WORK folder
667//    but this is a lousy choice.
668// x- binning is now Mask-aware. If mask is not present, all data is used. If data is from VCALC, all data is used
669// -- Where do I put the solid angle correction? In here as a weight for each point, or later on as
670//    a blanket correction (matrix multiply) for an entire panel?
671//
672//
673//
674// folderStr = WORK folder, type = the binning type (may include multiple detectors)
675Function VC_fDoBinning_QxQy2D(folderStr,type)
676        String folderStr,type
677       
678        Variable nSets = 0
679        Variable xDim,yDim
680        Variable ii,jj
681        Variable qVal,nq,var,avesq,aveisq
682        Variable binIndex,val,isVCALC=0,maskMissing
683
684        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
685        String instPath = ":entry:instrument:detector_"
686        String detStr
687               
688        if(cmpstr(folderStr,"VCALC") == 0)
689                isVCALC = 1
690        endif
691       
692// now switch on the type to determine which waves to declare and create
693// since there may be more than one panel to step through. There may be two, there may be four
694//
695
696// TODO:
697// -- Solid_Angle -- waves will be present for WORK data other than RAW, but not for RAW
698//
699// assume that the mask files are missing unless we can find them. If VCALC data,
700//  then the Mask is missing by definition
701        maskMissing = 1
702
703        strswitch(type) // string switch
704                case "FL":              // execute if case matches expression
705                case "FR":
706                        detStr = type
707                        if(isVCALC)
708                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
709                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation
710                        else
711                                Wave inten = V_getDetectorDataW(folderStr,detStr)
712                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
713                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
714                                if(WaveExists(mask) == 1)
715                                        maskMissing = 0
716                                endif
717                               
718                        endif
719                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
720                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
721                        nSets = 1
722                        break   
723                                                               
724                case "FT":             
725                case "FB":
726                        detStr = type
727                        if(isVCALC)
728                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
729                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
730                        else
731                                Wave inten = V_getDetectorDataW(folderStr,detStr)
732                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
733                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
734                                if(WaveExists(mask) == 1)
735                                        maskMissing = 0
736                                endif
737                        endif
738                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
739                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
740                        nSets = 1
741                        break
742                       
743                case "ML":             
744                case "MR":
745                        detStr = type
746                        if(isVCALC)
747                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
748                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
749                        else
750                                Wave inten = V_getDetectorDataW(folderStr,detStr)
751                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
752                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
753                                if(WaveExists(mask) == 1)
754                                        maskMissing = 0
755                                endif
756                        endif   
757                        //TODO:
758                        // -- decide on the proper deltaQ for binning. either nominal value for LR, or one
759                        //    determined specifically for that panel (currently using one tube width as deltaQ)
760                        // -- this is repeated multiple times in this switch
761                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
762//                      NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
763                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
764                        nSets = 1
765                        break   
766                                       
767                case "MT":             
768                case "MB":
769                        detStr = type
770                        if(isVCALC)
771                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
772                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
773                        else
774                                Wave inten = V_getDetectorDataW(folderStr,detStr)
775                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
776                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
777                                if(WaveExists(mask) == 1)
778                                        maskMissing = 0
779                                endif
780                        endif   
781                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
782                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
783                        nSets = 1
784                        break   
785                                       
786                case "B":       
787                        detStr = type
788                        if(isVCALC)
789                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
790                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
791                        else
792                                Wave inten = V_getDetectorDataW(folderStr,detStr)
793                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
794                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
795                                if(WaveExists(mask) == 1)
796                                        maskMissing = 0
797                                endif
798                        endif   
799                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_B")
800                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values 
801                        nSets = 1
802                        break   
803                       
804                case "FLR":
805                // detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
806                // TODO
807                // -- see if I can un-hard-wire some of this below when more than one panel is combined
808                        if(isVCALC)
809                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
810                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
811                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
812                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
813                        else
814                                Wave inten = V_getDetectorDataW(folderStr,"FL")
815                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
816                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
817                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
818                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
819                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
820                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
821                                        maskMissing = 0
822                                endif
823                        endif   
824                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
825                       
826                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
827                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
828               
829                        nSets = 2
830                        break                   
831               
832                case "FTB":
833                        if(isVCALC)
834                                WAVE inten = $(folderPath+instPath+"FT"+":det_"+"FT")
835                                WAVE/Z iErr = $("iErr_"+"FT")                   // 2D errors -- may not exist, especially for simulation               
836                                WAVE inten2 = $(folderPath+instPath+"FB"+":det_"+"FB")
837                                WAVE/Z iErr2 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
838                        else
839                                Wave inten = V_getDetectorDataW(folderStr,"FT")
840                                Wave iErr = V_getDetectorDataErrW(folderStr,"FT")
841                                Wave inten2 = V_getDetectorDataW(folderStr,"FB")
842                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FB")
843                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
844                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
845                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
846                                        maskMissing = 0
847                                endif
848                        endif   
849                        NVAR delQ = $(folderPath+instPath+"FT"+":gDelQ_FT")
850                       
851                        Wave qTotal = $(folderPath+instPath+"FT"+":qTot_"+"FT")                 // 2D q-values 
852                        Wave qTotal2 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
853       
854                        nSets = 2
855                        break           
856               
857                case "FLRTB":
858                        if(isVCALC)
859                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
860                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
861                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
862                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
863                                WAVE inten3 = $(folderPath+instPath+"FT"+":det_"+"FT")
864                                WAVE/Z iErr3 = $("iErr_"+"FT")                  // 2D errors -- may not exist, especially for simulation               
865                                WAVE inten4 = $(folderPath+instPath+"FB"+":det_"+"FB")
866                                WAVE/Z iErr4 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
867                        else
868                                Wave inten = V_getDetectorDataW(folderStr,"FL")
869                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
870                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
871                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
872                                Wave inten3 = V_getDetectorDataW(folderStr,"FT")
873                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"FT")
874                                Wave inten4 = V_getDetectorDataW(folderStr,"FB")
875                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"FB")
876                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
877                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
878                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
879                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
880                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
881                                        maskMissing = 0
882                                endif
883                        endif   
884                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
885                       
886                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
887                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
888                        Wave qTotal3 = $(folderPath+instPath+"FT"+":qTot_"+"FT")                        // 2D q-values 
889                        Wave qTotal4 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
890               
891                        nSets = 4
892                        break           
893                       
894                case "MLR":
895                        if(isVCALC)
896                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
897                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
898                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
899                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
900                        else
901                                Wave inten = V_getDetectorDataW(folderStr,"ML")
902                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
903                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
904                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
905                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
906                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
907                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
908                                        maskMissing = 0
909                                endif
910                        endif   
911                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
912                       
913                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
914                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
915               
916                        nSets = 2
917                        break                   
918               
919                case "MTB":
920                        if(isVCALC)
921                                WAVE inten = $(folderPath+instPath+"MT"+":det_"+"MT")
922                                WAVE/Z iErr = $("iErr_"+"MT")                   // 2D errors -- may not exist, especially for simulation               
923                                WAVE inten2 = $(folderPath+instPath+"MB"+":det_"+"MB")
924                                WAVE/Z iErr2 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
925                        else
926                                Wave inten = V_getDetectorDataW(folderStr,"MT")
927                                Wave iErr = V_getDetectorDataErrW(folderStr,"MT")
928                                Wave inten2 = V_getDetectorDataW(folderStr,"MB")
929                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MB")
930                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
931                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
932                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
933                                        maskMissing = 0
934                                endif
935                        endif   
936                        NVAR delQ = $(folderPath+instPath+"MT"+":gDelQ_MT")
937                       
938                        Wave qTotal = $(folderPath+instPath+"MT"+":qTot_"+"MT")                 // 2D q-values 
939                        Wave qTotal2 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
940               
941                        nSets = 2
942                        break                           
943               
944                case "MLRTB":
945                        if(isVCALC)
946                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
947                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
948                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
949                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
950                                WAVE inten3 = $(folderPath+instPath+"MT"+":det_"+"MT")
951                                WAVE/Z iErr3 = $("iErr_"+"MT")                  // 2D errors -- may not exist, especially for simulation               
952                                WAVE inten4 = $(folderPath+instPath+"MB"+":det_"+"MB")
953                                WAVE/Z iErr4 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
954                        else
955                                Wave inten = V_getDetectorDataW(folderStr,"ML")
956                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
957                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
958                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
959                                Wave inten3 = V_getDetectorDataW(folderStr,"MT")
960                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"MT")
961                                Wave inten4 = V_getDetectorDataW(folderStr,"MB")
962                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"MB")
963                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
964                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
965                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
966                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
967                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
968                                        maskMissing = 0
969                                endif
970                        endif   
971                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
972                       
973                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
974                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
975                        Wave qTotal3 = $(folderPath+instPath+"MT"+":qTot_"+"MT")                        // 2D q-values 
976                        Wave qTotal4 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
977               
978                        nSets = 4
979                        break                                                                   
980                                       
981                default:
982                        nSets = 0                                                       // optional default expression executed
983                        Print "ERROR   ---- type is not recognized "
984        endswitch
985
986//      Print "delQ = ",delQ," for ",type
987
988        if(nSets == 0)
989                SetDataFolder root:
990                return(0)
991        endif
992
993
994//TODO: properly define the errors here - I'll have this if I do the simulation
995// -- need to propagate the errors up to this point
996//
997        if(WaveExists(iErr)==0  && WaveExists(inten) != 0)
998                Duplicate/O inten,iErr
999                Wave iErr=iErr
1000//              iErr = 1+sqrt(inten+0.75)                       // can't use this -- it applies to counts, not intensity (already a count rate...)
1001                iErr = sqrt(inten+0.75)                 // TODO -- here I'm just using some fictional value
1002        endif
1003        if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0)
1004                Duplicate/O inten2,iErr2
1005                Wave iErr2=iErr2
1006//              iErr2 = 1+sqrt(inten2+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1007                iErr2 = sqrt(inten2+0.75)                       // TODO -- here I'm just using some fictional value
1008        endif
1009        if(WaveExists(iErr3)==0  && WaveExists(inten3) != 0)
1010                Duplicate/O inten3,iErr3
1011                Wave iErr3=iErr3
1012//              iErr3 = 1+sqrt(inten3+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1013                iErr3 = sqrt(inten3+0.75)                       // TODO -- here I'm just using some fictional value
1014        endif
1015        if(WaveExists(iErr4)==0  && WaveExists(inten4) != 0)
1016                Duplicate/O inten4,iErr4
1017                Wave iErr4=iErr4
1018//              iErr4 = 1+sqrt(inten4+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1019                iErr4 = sqrt(inten4+0.75)                       // TODO -- here I'm just using some fictional value
1020        endif
1021
1022        nq = 600
1023
1024        // note that the back panel of 320x320 (1mm res) results in 447 data points!
1025        // - so I upped nq to 600
1026
1027//******TODO****** -- where to put the averaged data -- right now, folderStr is forced to ""   
1028//      SetDataFolder $("root:"+folderStr)              //should already be here, but make sure...     
1029        Make/O/D/N=(nq)  $(folderPath+":"+"iBin_qxqy"+"_"+type)
1030        Make/O/D/N=(nq)  $(folderPath+":"+"qBin_qxqy"+"_"+type)
1031        Make/O/D/N=(nq)  $(folderPath+":"+"nBin_qxqy"+"_"+type)
1032        Make/O/D/N=(nq)  $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1033        Make/O/D/N=(nq)  $(folderPath+":"+"eBin_qxqy"+"_"+type)
1034        Make/O/D/N=(nq)  $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1035       
1036        Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+type)
1037        Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type)
1038        Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy"+"_"+type)
1039        Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1040        Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy"+"_"+type)
1041        Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1042       
1043       
1044//      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
1045// TODO: not sure if I want to set dQ in x or y direction...
1046        // the short dimension is the 8mm tubes, use this direction as dQ?
1047        // 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]
1048        // WRT the beam center. use qx or qy directly. Still not happy with this way...
1049
1050
1051        qBin_qxqy[] =  p*delQ   
1052        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
1053
1054        iBin_qxqy = 0
1055        iBin2_qxqy = 0
1056        eBin_qxqy = 0
1057        eBin2D_qxqy = 0
1058        nBin_qxqy = 0   //number of intensities added to each bin
1059
1060// now there are situations of:
1061// 1 panel
1062// 2 panels
1063// 4 panels
1064//
1065// this needs to be a double loop now...
1066// TODO:
1067// -- the iErr wave and accumulation of error is NOT CALCULATED CORRECTLY YET
1068// -- the solid angle per pixel is not yet implemented.
1069//    it will be present for WORK data other than RAW, but not for RAW
1070
1071// if any of the masks don't exist, display the error, and proceed with the averaging, using all data
1072        if(maskMissing == 1)
1073                Print "Mask file not found for at least one detector - so all data is used"
1074        endif
1075
1076        Variable mask_val
1077// use set 1 (no number) only
1078        if(nSets >= 1)
1079                xDim=DimSize(inten,0)
1080                yDim=DimSize(inten,1)
1081       
1082                for(ii=0;ii<xDim;ii+=1)
1083                        for(jj=0;jj<yDim;jj+=1)
1084                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1085                                qVal = qTotal[ii][jj]
1086                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1087                                val = inten[ii][jj]
1088                               
1089                                if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
1090                                        mask_val = 0
1091                                else
1092                                        mask_val = mask[ii][jj]
1093                                endif
1094                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1095                                        iBin_qxqy[binIndex] += val
1096                                        iBin2_qxqy[binIndex] += val*val
1097                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1098                                        nBin_qxqy[binIndex] += 1
1099                                endif
1100                        endfor
1101                endfor
1102               
1103        endif
1104
1105// add in set 2 (set 1 already done)
1106        if(nSets >= 2)
1107                xDim=DimSize(inten2,0)
1108                yDim=DimSize(inten2,1)
1109       
1110                for(ii=0;ii<xDim;ii+=1)
1111                        for(jj=0;jj<yDim;jj+=1)
1112                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1113                                qVal = qTotal2[ii][jj]
1114                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1115                                val = inten2[ii][jj]
1116                               
1117                                if(isVCALC || maskMissing)
1118                                        mask_val = 0
1119                                else
1120                                        mask_val = mask2[ii][jj]
1121                                endif
1122                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1123                                        iBin_qxqy[binIndex] += val
1124                                        iBin2_qxqy[binIndex] += val*val
1125                                        eBin2D_qxqy[binIndex] += iErr2[ii][jj]*iErr2[ii][jj]
1126                                        nBin_qxqy[binIndex] += 1
1127                                endif
1128                        endfor
1129                endfor
1130               
1131        endif
1132
1133// add in set 3 and 4 (set 1 and 2 already done)
1134        if(nSets == 4)
1135                xDim=DimSize(inten3,0)
1136                yDim=DimSize(inten3,1)
1137       
1138                for(ii=0;ii<xDim;ii+=1)
1139                        for(jj=0;jj<yDim;jj+=1)
1140                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1141                                qVal = qTotal3[ii][jj]
1142                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1143                                val = inten3[ii][jj]
1144                               
1145                                if(isVCALC || maskMissing)
1146                                        mask_val = 0
1147                                else
1148                                        mask_val = mask3[ii][jj]
1149                                endif
1150                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1151                                        iBin_qxqy[binIndex] += val
1152                                        iBin2_qxqy[binIndex] += val*val
1153                                        eBin2D_qxqy[binIndex] += iErr3[ii][jj]*iErr3[ii][jj]
1154                                        nBin_qxqy[binIndex] += 1
1155                                endif
1156                        endfor
1157                endfor
1158               
1159               
1160                xDim=DimSize(inten4,0)
1161                yDim=DimSize(inten4,1)
1162       
1163                for(ii=0;ii<xDim;ii+=1)
1164                        for(jj=0;jj<yDim;jj+=1)
1165                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1166                                qVal = qTotal4[ii][jj]
1167                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1168                                val = inten4[ii][jj]
1169                               
1170                                if(isVCALC || maskMissing)
1171                                        mask_val = 0
1172                                else
1173                                        mask_val = mask4[ii][jj]
1174                                endif
1175                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1176                                        iBin_qxqy[binIndex] += val
1177                                        iBin2_qxqy[binIndex] += val*val
1178                                        eBin2D_qxqy[binIndex] += iErr4[ii][jj]*iErr4[ii][jj]
1179                                        nBin_qxqy[binIndex] += 1
1180                                endif
1181                        endfor
1182                endfor
1183               
1184        endif
1185
1186
1187// after looping through all of the data on the panels, calculate errors on I(q),
1188// just like in CircSectAve.ipf
1189// TODO:
1190// -- Errors were NOT properly acculumated above, so this loop of calculations is NOT MEANINGFUL (yet)
1191        for(ii=0;ii<nq;ii+=1)
1192                if(nBin_qxqy[ii] == 0)
1193                        //no pixels in annuli, data unknown
1194                        iBin_qxqy[ii] = 0
1195                        eBin_qxqy[ii] = 1
1196                        eBin2D_qxqy[ii] = NaN
1197                else
1198                        if(nBin_qxqy[ii] <= 1)
1199                                //need more than one pixel to determine error
1200                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1201                                eBin_qxqy[ii] = 1
1202                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1203                        else
1204                                //assume that the intensity in each pixel in annuli is normally distributed about mean...
1205                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1206                                avesq = iBin_qxqy[ii]^2
1207                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1208                                var = aveisq-avesq
1209                                if(var<=0)
1210                                        eBin_qxqy[ii] = 1e-6
1211                                else
1212                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1213                                endif
1214                                // and calculate as it is propagated pixel-by-pixel
1215                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1216                        endif
1217                endif
1218        endfor
1219       
1220        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo
1221       
1222        // find the last non-zero point, working backwards
1223        val=nq
1224        do
1225                val -= 1
1226        while((nBin_qxqy[val] == 0) && val > 0)
1227       
1228//      print val, nBin_qxqy[val]
1229        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1230
1231        if(val == 0)
1232                // all the points were deleted
1233                return(0)
1234        endif
1235       
1236       
1237        // since the beam center is not always on the detector, many of the low Q bins will have zero pixels
1238        // find the first non-zero point, working forwards
1239        val = -1
1240        do
1241                val += 1
1242        while(nBin_qxqy[val] == 0)     
1243        DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1244
1245        // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
1246        val = numpnts(nBin_qxqy)-1
1247        do
1248                if(nBin_qxqy[val] == 0)
1249                        DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1250                endif
1251                val -= 1
1252        while(val>0)
1253       
1254        SetDataFolder root:
1255       
1256        return(0)
1257End
1258
1259
Note: See TracBrowser for help on using the repository browser.