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

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

lots of changes to plotting of q-values, generating fake data with non-linear corrections, masking of data, etc.

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