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

Last change on this file since 957 was 955, checked in by srkline, 8 years ago

some reorganization of the r/w routines to generate HDF test files for SANS and VSANS (all are housed together for testing). also some reorganization of the detector binning routines to get functions grouped in more logical locations.
File size: 29.6 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2
3/////////////////////////
4//
5// Utility functions to:
6//              calculate Q, Qx, Qy, Qz
7//              fill the detector panels with simulated data (the model functions are here)
8//              bin the 2D detector to 1D I(Q) based on Q and deltaQ (bin width)
9//
10/////////////////////////
11
12
13
14
15// 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)
17Function FillPanel_wModelData(det,qTot,type)
18        Wave det,qTot
19        String type
20
21        SetDataFolder root:Packages:NIST:VSANS:VCALC:Front
22
23        // q-values and detector arrays already allocated and calculated
24        Duplicate/O det tmpInten,tmpSig,prob_i
25               
26        Variable imon,trans,thick,sdd,pixSizeX,pixSizeY,sdd_offset
27
28        //imon = V_BeamIntensity()*CountTime
29        imon = VCALC_getImon()          //TODO: currently from the panel, not calculated
30        trans = 0.8
31        thick = 0.1
32       
33        // need SDD
34        // need pixel dimensions
35        // nominal sdd in meters, offset in mm, want result in cm !
36        sdd = VCALC_getSDD(type)*100    +  VCALC_getTopBottomSDDOffset(type) / 10               // result is sdd in [cm]
37
38        pixSizeX = VCALC_getPixSizeX(type)              // cm
39        pixSizeY = VCALC_getPixSizeY(type)
40       
41       
42        //?? pick the function from a popup on the panel? (bypass the analysis panel, or maybe it's better to
43        //  keep the panel to keep people used to using it.)
44        // peak @ 0.1 ~ AgBeh
45        //      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}               
46        //
47        // peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
48        //      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}             
49        String funcStr = VCALC_getModelFunctionStr()
50        strswitch(funcStr)
51                case "Big Debye":
52                        tmpInten = V_Debye(10,3000,0.0001,qTot[p][q])
53                        break
54                case "Big Sphere":
55                        tmpInten = V_SphereForm(1,900,1e-6,0.01,qTot[p][q])     
56                        break
57                case "Debye":
58                        tmpInten = V_Debye(10,300,0.0001,qTot[p][q])
59                        break
60                case "Sphere":
61                        tmpInten = V_SphereForm(1,60,1e-6,0.001,qTot[p][q])     
62                        break
63                case "AgBeh":
64                        tmpInten = V_BroadPeak(1e-9,3,20,100.0,0.1,3,0.1,qTot[p][q])
65                        break
66                case "Vycor":
67                        tmpInten = V_BroadPeak(1e-9,3,20,500.0,0.015,3,0.1,qTot[p][q])
68                        break   
69                case "Empty Cell":
70                        tmpInten = V_EC_Empirical(2.2e-8,3.346,0.0065,9.0,0.016,qTot[p][q])
71                        break
72                case "Blocked Beam":
73                        tmpInten = V_BlockedBeam(1,qTot[p][q])
74                        break
75                default:
76                        tmpInten = V_Debye(10,300,0.1,qTot[p][q])
77        endswitch
78
79
80///////////////
81//      // calculate the scattering cross section simply to be able to estimate the transmission
82//      Variable sig_sas=0
83//     
84//      // remember that the random deviate is the coherent portion ONLY - the incoherent background is
85//      // subtracted before the calculation.
86//      CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
87//
88//      if(sig_sas > 100)
89//              DoAlert 0,"SAS cross section > 100. Estimates of multiple scattering are unreliable. Choosing a model with a well-defined Rg may help"
90//      endif           
91//
92//      // calculate the multiple scattering fraction for display (10/2009)
93//      Variable ii,nMax=10,tau
94//      mScat=0
95//      tau = thick*sig_sas
96//      // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
97//      // neutrons out of those that were scattered
98//      for(ii=2;ii<nMax;ii+=1)
99//              mScat += tau^(ii)/factorial(ii)
100////            print tau^(ii)/factorial(ii)
101//      endfor
102//      estTrans = exp(-1*thick*sig_sas)                //thickness and sigma both in units of cm
103//      mscat *= (estTrans)/(1-estTrans)
104//
105////    if(mScat > 0.1)         //  Display warning
106//
107//      Print "Sig_sas = ",sig_sas
108////////////////////
109       
110        prob_i = trans*thick*pixSizeX*pixSizeY/(sdd)^2*tmpInten                 //probability of a neutron in q-bin(i)
111               
112        tmpInten = (imon)*prob_i                //tmpInten is not the model calculation anymore!!
113
114
115/// **** can I safely assume a Gaussian error in the count rate??
116        tmpSig = sqrt(tmpInten)         // corrected based on John's memo, from 8/9/99
117
118        tmpInten += gnoise(tmpSig)
119        tmpInten = (tmpInten[p][q] < 0) ? 0 : tmpInten[p][q]                    // MAR 2013 -- is this the right thing to do
120        tmpInten = trunc(tmpInten)
121               
122       
123        det = tmpInten
124
125// if I want "absolute" scale -- then I lose the integer nature of the detector (but keep the random)
126//      det /= trans*thick*pixSizeX*pixSizeY/(sdd)^2*imon
127
128       
129        KillWaves/Z tmpInten,tmpSig,prob_i     
130        SetDataFolder root:
131
132        return(0)
133End
134
135
136// For a given detector panel, calculate the q-values
137// -work with everything as arrays
138// Input needed:
139// detector data
140// detector type (LRTB?)
141// beam center (may be off the detector)
142// SDD
143// lambda
144//
145// pixel dimensions for detector type (global constants)
146// - data dimensions read directly from array
147//
148// --What is calculated:
149// array of Q
150// array of qx,qy,qz
151// array of error already exists
152//
153//
154// -- sdd in meters
155// -- lambda in Angstroms
156Function V_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
157        Wave data,qTot,qx,qy,qz
158        Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
159               
160        // loop over the array and calculate the values - this is done as a wave assignment
161// TODO -- be sure that it's p,q -- or maybe p+1,q+1 as used in WriteQIS.ipf   
162        qTot = V_CalcQval(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
163        qx = V_CalcQX(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
164        qy = V_CalcQY(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
165        qz = V_CalcQZ(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
166       
167        return(0)
168End
169
170
171//////////////////////
172// NOTE: The Q calculations are different than what is in GaussUtils in that they take into
173// accout the different x/y pixel sizes and the beam center not being on the detector -
174// off a different edge for each LRTB type
175/////////////////////
176
177//function to calculate the overall q-value, given all of the necesary trig inputs
178//and are in detector coordinates (1,128) rather than axis values
179//the pixel locations need not be integers, reals are ok inputs
180//sdd is in meters
181//wavelength is in Angstroms
182//
183//returned magnitude of Q is in 1/Angstroms
184//
185Function V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
186        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
187       
188        Variable dx,dy,qval,two_theta,dist
189               
190        sdd *=100               //convert to cm
191        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
192        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
193        dist = sqrt(dx^2 + dy^2)
194       
195        two_theta = atan(dist/sdd)
196
197        qval = 4*Pi/lam*sin(two_theta/2)
198       
199        return qval
200End
201
202//calculates just the q-value in the x-direction on the detector
203//input/output is the same as CalcQval()
204//ALL inputs are in detector coordinates
205//
206//sdd is in meters
207//wavelength is in Angstroms
208//
209// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
210// now properly accounts for qz
211//
212Function V_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
213        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
214
215        Variable qx,qval,phi,dx,dy,dist,two_theta
216       
217        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
218       
219        sdd *=100               //convert to cm
220        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
221        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
222        phi = V_FindPhi(dx,dy)
223       
224        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
225        dist = sqrt(dx^2 + dy^2)
226        two_theta = atan(dist/sdd)
227
228        qx = qval*cos(two_theta/2)*cos(phi)
229       
230        return qx
231End
232
233//calculates just the q-value in the y-direction on the detector
234//input/output is the same as CalcQval()
235//ALL inputs are in detector coordinates
236//sdd is in meters
237//wavelength is in Angstroms
238//
239// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
240// now properly accounts for qz
241//
242Function V_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
243        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
244       
245        Variable dy,qval,dx,phi,qy,dist,two_theta
246       
247        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
248       
249        sdd *=100               //convert to cm
250        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
251        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
252        phi = V_FindPhi(dx,dy)
253       
254        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
255        dist = sqrt(dx^2 + dy^2)
256        two_theta = atan(dist/sdd)
257       
258        qy = qval*cos(two_theta/2)*sin(phi)
259       
260        return qy
261End
262
263//calculates just the z-component of the q-vector, not measured on the detector
264//input/output is the same as CalcQval()
265//ALL inputs are in detector coordinates
266//sdd is in meters
267//wavelength is in Angstroms
268//
269// not actually used, but here for completeness if anyone asks
270//
271Function V_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
272        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
273       
274        Variable dy,qval,dx,phi,qz,dist,two_theta
275       
276        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
277       
278        sdd *=100               //convert to cm
279       
280        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
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        two_theta = atan(dist/sdd)
285       
286        qz = qval*sin(two_theta/2)
287       
288        return qz
289End
290
291//phi is defined from +x axis, proceeding CCW around [0,2Pi]
292Threadsafe Function V_FindPhi(vx,vy)
293        variable vx,vy
294       
295        variable phi
296       
297        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
298       
299        // special cases
300        if(vx==0 && vy > 0)
301                return(pi/2)
302        endif
303        if(vx==0 && vy < 0)
304                return(3*pi/2)
305        endif
306        if(vx >= 0 && vy == 0)
307                return(0)
308        endif
309        if(vx < 0 && vy == 0)
310                return(pi)
311        endif
312       
313       
314        if(vx > 0 && vy > 0)
315                return(phi)
316        endif
317        if(vx < 0 && vy > 0)
318                return(phi + pi)
319        endif
320        if(vx < 0 && vy < 0)
321                return(phi + pi)
322        endif
323        if( vx > 0 && vy < 0)
324                return(phi + 2*pi)
325        endif
326       
327        return(phi)
328end
329
330Function V_SphereForm(scale,radius,delrho,bkg,x)                               
331        Variable scale,radius,delrho,bkg
332        Variable x
333       
334        // variables are:                                                       
335        //[0] scale
336        //[1] radius (A)
337        //[2] delrho (A-2)
338        //[3] background (cm-1)
339       
340//      Variable scale,radius,delrho,bkg                               
341//      scale = w[0]
342//      radius = w[1]
343//      delrho = w[2]
344//      bkg = w[3]
345       
346       
347        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
348        // and is rescaled to give [=] cm^-1
349       
350        Variable bes,f,vol,f2
351        //
352        //handle q==0 separately
353        If(x==0)
354                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
355                return(f)
356        Endif
357       
358//      bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
359       
360        bes = 3*sqrt(pi/(2*x*radius))*BesselJ(1.5,x*radius)/(x*radius)
361       
362        vol = 4*pi/3*radius^3
363        f = vol*bes*delrho              // [=] A
364        // normalize to single particle volume, convert to 1/cm
365        f2 = f * f / vol * 1.0e8                // [=] 1/cm
366       
367        return (scale*f2+bkg)   // Scale, then add in the background
368       
369End
370
371Function V_Debye(scale,rg,bkg,x)
372        Variable scale,rg,bkg
373        Variable x
374       
375        // variables are:
376        //[0] scale factor
377        //[1] radius of gyration [A]
378        //[2] background        [cm-1]
379       
380        // calculates (scale*debye)+bkg
381        Variable Pq,qr2
382       
383        qr2=(x*rg)^2
384        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
385       
386        //scale
387        Pq *= scale
388        // then add in the background
389        return (Pq+bkg)
390End
391
392// a sum of a power law and debye to approximate the scattering from a real empty cell
393//
394//      make/O/D coef_ECEmp = {2.2e-8,3.346,0.0065,9.0,0.016}
395//
396Function V_EC_Empirical(aa,mm,scale,rg,bkg,x)
397        Variable aa,mm,scale,rg,bkg
398        Variable x
399       
400        // variables are:
401        //[0] = A
402        //[1] = power m
403        //[2] scale factor
404        //[3] radius of gyration [A]
405        //[4] background        [cm-1]
406       
407        Variable Iq
408       
409        // calculates (scale*debye)+bkg
410        Variable Pq,qr2
411       
412//      if(x*Rg < 1e-3)         //added Oct 2008 to avoid numerical errors at low arg values
413//              return(scale+bkg)
414//      endif
415       
416        Iq = aa*x^-mm
417       
418        qr2=(x*rg)^2
419        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
420       
421        //scale
422        Pq *= scale
423        // then add the terms up
424        return (Iq + Pq + bkg)
425End
426
427// blocked beam
428//
429Function V_BlockedBeam(bkg,x)
430        Variable bkg
431        Variable x
432       
433        return (bkg)
434End
435
436
437//
438// a broad peak to simulate silver behenate or vycor
439//
440// peak @ 0.1 ~ AgBeh
441//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}               
442//
443//
444// peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
445//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}             
446//
447//
448Function V_BroadPeak(aa,nn,cc,LL,Qzero,mm,bgd,x)
449        Variable aa,nn,cc,LL,Qzero,mm,bgd
450        Variable x
451       
452        // variables are:                                                       
453        //[0] Porod term scaling
454        //[1] Porod exponent
455        //[2] Lorentzian term scaling
456        //[3] Lorentzian screening length [A]
457        //[4] peak location [1/A]
458        //[5] Lorentzian exponent
459        //[6] background
460       
461//      local variables
462        Variable inten, qval
463//      x is the q-value for the calculation
464        qval = x
465//      do the calculation and return the function value
466       
467        inten = aa/(qval)^nn + cc/(1 + (abs(qval-Qzero)*LL)^mm) + bgd
468
469        Return (inten)
470       
471End
472
473
474
475Function SetDeltaQ(folderStr,type)
476        String folderStr,type
477
478        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + folderStr + ":det_"+type)            // 2D detector data
479       
480        Variable xDim,yDim,delQ
481       
482        xDim=DimSize(inten,0)
483        yDim=DimSize(inten,1)
484       
485        if(xDim<yDim)
486                WAVE qx = $("root:Packages:NIST:VSANS:VCALC:" + folderStr + ":qx_"+type)
487                delQ = abs(qx[0][0] - qx[1][0])/2
488        else
489                WAVE qy = $("root:Packages:NIST:VSANS:VCALC:" + folderStr + ":qy_"+type)
490                delQ = abs(qy[0][1] - qy[0][0])/2
491        endif
492       
493        // set the global
494        Variable/G $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_"+type) = delQ
495//      Print "SET delQ = ",delQ," for ",type
496       
497        return(0)
498end
499
500
501//TODO -- need a switch here to dispatch to the averaging type
502Proc V_BinQxQy_to_1D(folderStr,type)
503        String folderStr
504        String type
505//      Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4)
506//      Prompt type,"detector identifier"
507
508
509        V_fDoBinning_QxQy2D(folderStr, type)
510
511
512/// this is for a tall, narrow slit mode       
513//      V_fBinDetector_byRows(folderStr,type)
514       
515End
516
517
518Proc V_Graph_1D_detType(folderStr,type)
519        String folderStr,type
520       
521        SetDataFolder root:Packages:NIST:VSANS:VCALC
522       
523        Display $("iBin_qxqy"+"_"+type) vs $("qBin_qxqy"+"_"+type)
524        ModifyGraph mirror=2,grid=1,log=1
525        ModifyGraph mode=4,marker=19,msize=2
526//      ErrorBars/T=0 iBin_qxqy Y,wave=(eBin2D_qxqy,eBin2D_qxqy)                // for simulations, I don't have 2D uncertainty
527        ErrorBars/T=0 $("iBin_qxqy"+"_"+type) Y,wave=($("eBin_qxqy"+"_"+type),$("eBin_qxqy"+"_"+type))
528        legend
529       
530        SetDataFolder root:
531
532End
533
534
535
536//////////
537//
538//              Function that bins a 2D detctor panel into I(q) based on the q-value of the pixel
539//              - each pixel QxQyQz has been calculated beforehand
540//              - if multiple panels are selected to be combined, it is done here during the binning
541//              - the setting of deltaQ step is still a little suspect (TODO)
542//
543//
544// see the equivalent function in PlotUtils2D_v40.ipf
545//
546//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
547//
548// this has been modified to accept different detector panels and to take arrays
549// -- type = FL or FR or...other panel identifiers
550//
551// TODO "iErr" is not always defined correctly since it doesn't really apply here for data that is not 2D simulation
552//
553Function V_fDoBinning_QxQy2D(folderStr,type)
554        String folderStr,type
555       
556        Variable nSets = 0
557        Variable xDim,yDim
558        Variable ii,jj
559        Variable qVal,nq,var,avesq,aveisq
560        Variable binIndex,val
561
562       
563        SetDataFolder root:Packages:NIST:VSANS:VCALC
564       
565// now switch on the type to determine which waves to declare and create
566// since there may be more than one panel to step through. There may be two, there may be four
567//
568
569        strswitch(type) // string switch
570                case "FL":              // execute if case matches expression
571                case "FR":
572                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_FL")
573                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+type)              // 2D detector data
574                        WAVE/Z iErr = $("iErr_"+type)                   // 2D errors -- may not exist, especially for simulation
575                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+type)                     // 2D q-values
576                        nSets = 1
577                        break   
578                                                               
579                case "FT":             
580                case "FB":
581                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_FT")
582                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+type)              // 2D detector data
583                        WAVE/Z iErr = $("iErr_"+type)                   // 2D errors -- may not exist, especially for simulation
584                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+type)                     // 2D q-values
585                        nSets = 1
586                        break
587                       
588                case "ML":             
589                case "MR":             
590                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_ML")
591                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+type)             // 2D detector data
592                        WAVE/Z iErr = $("iErr_"+type)                   // 2D errors -- may not exist, especially for simulation
593                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+type)                    // 2D q-values
594                        nSets = 1
595                        break   
596                                       
597                case "MT":             
598                case "MB":             
599                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_MT")
600                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+type)             // 2D detector data
601                        WAVE/Z iErr = $("iErr_"+type)                   // 2D errors -- may not exist, especially for simulation
602                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+type)                    // 2D q-values
603                        nSets = 1
604                        break   
605                                       
606                case "B":               
607                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_B")
608                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Back" + ":det_"+type)               // 2D detector data
609                        WAVE/Z iErr = $("iErr_"+type)                   // 2D errors -- may not exist, especially for simulation
610                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Back" +":qTot_"+type)                      // 2D q-values
611                        nSets = 1
612                        break   
613                       
614                case "FLR":
615                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_FL")
616                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FL")              // 2D detector data
617                        WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation
618                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FL")                     // 2D q-values
619                        WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FR")             // 2D detector data
620                        WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation
621                        Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FR")                    // 2D q-values
622                        nSets = 2
623                        break                   
624               
625                case "FTB":
626                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_FT")
627                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FT")              // 2D detector data
628                        WAVE/Z iErr = $("iErr_"+"FT")                   // 2D errors -- may not exist, especially for simulation
629                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FT")                     // 2D q-values
630                        WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FB")             // 2D detector data
631                        WAVE/Z iErr2 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation
632                        Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FB")                    // 2D q-values
633                        nSets = 2
634                        break           
635               
636                case "FLRTB":
637                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_FL")
638                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FL")              // 2D detector data
639                        WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation
640                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FL")                     // 2D q-values
641                        WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FR")             // 2D detector data
642                        WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation
643                        Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FR")                    // 2D q-values
644                        WAVE inten3 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FT")             // 2D detector data
645                        WAVE/Z iErr3 = $("iErr_"+"FT")                  // 2D errors -- may not exist, especially for simulation
646                        Wave qTotal3 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FT")                    // 2D q-values
647                        WAVE inten4 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FB")             // 2D detector data
648                        WAVE/Z iErr4 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation
649                        Wave qTotal4 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FB")                    // 2D q-values
650                        nSets = 4
651                        break           
652                       
653
654                case "MLR":
655                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_ML")
656                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"ML")             // 2D detector data
657                        WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation
658                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"ML")                    // 2D q-values
659                        WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MR")            // 2D detector data
660                        WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation
661                        Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MR")                   // 2D q-values
662                        nSets = 2
663                        break                   
664               
665                case "MTB":
666                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_MT")
667                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MT")             // 2D detector data
668                        WAVE/Z iErr = $("iErr_"+"MT")                   // 2D errors -- may not exist, especially for simulation
669                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MT")                    // 2D q-values
670                        WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MB")            // 2D detector data
671                        WAVE/Z iErr2 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation
672                        Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MB")                   // 2D q-values
673                        nSets = 2
674                        break                           
675               
676                case "MLRTB":
677                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_ML")
678                        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"ML")             // 2D detector data
679                        WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation
680                        Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"ML")                    // 2D q-values
681                        WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MR")            // 2D detector data
682                        WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation
683                        Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MR")                   // 2D q-values
684                        WAVE inten3 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MT")            // 2D detector data
685                        WAVE/Z iErr3 = $("iErr_"+"MT")                  // 2D errors -- may not exist, especially for simulation
686                        Wave qTotal3 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MT")                   // 2D q-values
687                        WAVE inten4 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MB")            // 2D detector data
688                        WAVE/Z iErr4 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation
689                        Wave qTotal4 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MB")                   // 2D q-values
690                        nSets = 4
691                        break                                                                   
692                                       
693                default:
694                        nSets = 0                                                       // optional default expression executed
695                        Print "ERROR   ---- type is not recognized "
696        endswitch
697
698//      Print "delQ = ",delQ," for ",type
699
700        if(nSets == 0)
701                return(0)
702        endif
703
704
705//TODO: properly define the errors here - I'll have this if I do the simulation
706        if(WaveExists(iErr)==0  && WaveExists(inten) != 0)
707                Duplicate/O inten,iErr
708                Wave iErr=iErr
709//              iErr = 1+sqrt(inten+0.75)                       // can't use this -- it applies to counts, not intensity (already a count rate...)
710                iErr = sqrt(inten+0.75)                 // TODO -- here I'm just using some fictional value
711        endif
712        if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0)
713                Duplicate/O inten2,iErr2
714                Wave iErr2=iErr2
715//              iErr2 = 1+sqrt(inten2+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
716                iErr2 = sqrt(inten2+0.75)                       // TODO -- here I'm just using some fictional value
717        endif
718        if(WaveExists(iErr3)==0  && WaveExists(inten3) != 0)
719                Duplicate/O inten3,iErr3
720                Wave iErr3=iErr3
721//              iErr3 = 1+sqrt(inten3+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
722                iErr3 = sqrt(inten3+0.75)                       // TODO -- here I'm just using some fictional value
723        endif
724        if(WaveExists(iErr4)==0  && WaveExists(inten4) != 0)
725                Duplicate/O inten4,iErr4
726                Wave iErr4=iErr4
727//              iErr4 = 1+sqrt(inten4+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
728                iErr4 = sqrt(inten4+0.75)                       // TODO -- here I'm just using some fictional value
729        endif
730
731        nq = 600
732
733        // note that the back panel of 320x320 (1mm res) results in 447 data points!
734        // - so I upped nq to 600
735
736//******TODO****** -- where to put the averaged data -- right now, folderStr is forced to ""   
737//      SetDataFolder $("root:"+folderStr)              //should already be here, but make sure...     
738        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"iBin_qxqy"+"_"+type)
739        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"qBin_qxqy"+"_"+type)
740        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"nBin_qxqy"+"_"+type)
741        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"iBin2_qxqy"+"_"+type)
742        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"eBin_qxqy"+"_"+type)
743        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"eBin2D_qxqy"+"_"+type)
744       
745        Wave iBin_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"iBin_qxqy_"+type)
746        Wave qBin_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"qBin_qxqy"+"_"+type)
747        Wave nBin_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"nBin_qxqy"+"_"+type)
748        Wave iBin2_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"iBin2_qxqy"+"_"+type)
749        Wave eBin_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"eBin_qxqy"+"_"+type)
750        Wave eBin2D_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"eBin2D_qxqy"+"_"+type)
751       
752       
753//      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
754// TODO: not sure if I want to set dQ in x or y direction...
755        // the short dimension is the 8mm tubes, use this direction as dQ?
756        // 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]
757        // WRT the beam center. use qx or qy directly. Still not happy with this way...
758
759
760        qBin_qxqy[] =  p*delQ   
761        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
762
763        iBin_qxqy = 0
764        iBin2_qxqy = 0
765        eBin_qxqy = 0
766        eBin2D_qxqy = 0
767        nBin_qxqy = 0   //number of intensities added to each bin
768
769// now there are situations of:
770// 1 panel
771// 2 panels
772// 4 panels
773//
774// this needs to be a double loop now...
775
776// use set 1 (no number) only
777        if(nSets >= 1)
778                xDim=DimSize(inten,0)
779                yDim=DimSize(inten,1)
780       
781                for(ii=0;ii<xDim;ii+=1)
782                        for(jj=0;jj<yDim;jj+=1)
783                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
784                                qVal = qTotal[ii][jj]
785                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
786                                val = inten[ii][jj]
787                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
788                                        iBin_qxqy[binIndex] += val
789                                        iBin2_qxqy[binIndex] += val*val
790                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
791                                        nBin_qxqy[binIndex] += 1
792                                endif
793                        endfor
794                endfor
795               
796        endif
797
798// add in set 2 (set 1 already done)
799        if(nSets >= 2)
800                xDim=DimSize(inten2,0)
801                yDim=DimSize(inten2,1)
802       
803                for(ii=0;ii<xDim;ii+=1)
804                        for(jj=0;jj<yDim;jj+=1)
805                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
806                                qVal = qTotal2[ii][jj]
807                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
808                                val = inten2[ii][jj]
809                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
810                                        iBin_qxqy[binIndex] += val
811                                        iBin2_qxqy[binIndex] += val*val
812                                        eBin2D_qxqy[binIndex] += iErr2[ii][jj]*iErr2[ii][jj]
813                                        nBin_qxqy[binIndex] += 1
814                                endif
815                        endfor
816                endfor
817               
818        endif
819
820// add in set 3 and 4 (set 1 and 2already done)
821        if(nSets == 4)
822                xDim=DimSize(inten3,0)
823                yDim=DimSize(inten3,1)
824       
825                for(ii=0;ii<xDim;ii+=1)
826                        for(jj=0;jj<yDim;jj+=1)
827                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
828                                qVal = qTotal3[ii][jj]
829                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
830                                val = inten3[ii][jj]
831                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
832                                        iBin_qxqy[binIndex] += val
833                                        iBin2_qxqy[binIndex] += val*val
834                                        eBin2D_qxqy[binIndex] += iErr3[ii][jj]*iErr3[ii][jj]
835                                        nBin_qxqy[binIndex] += 1
836                                endif
837                        endfor
838                endfor
839               
840               
841                xDim=DimSize(inten4,0)
842                yDim=DimSize(inten4,1)
843       
844                for(ii=0;ii<xDim;ii+=1)
845                        for(jj=0;jj<yDim;jj+=1)
846                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
847                                qVal = qTotal4[ii][jj]
848                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
849                                val = inten4[ii][jj]
850                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
851                                        iBin_qxqy[binIndex] += val
852                                        iBin2_qxqy[binIndex] += val*val
853                                        eBin2D_qxqy[binIndex] += iErr4[ii][jj]*iErr4[ii][jj]
854                                        nBin_qxqy[binIndex] += 1
855                                endif
856                        endfor
857                endfor
858               
859        endif
860
861
862// after looping through all of the data on the panels, calculate errors on I(q),
863// just like in CircSectAve.ipf
864        for(ii=0;ii<nq;ii+=1)
865                if(nBin_qxqy[ii] == 0)
866                        //no pixels in annuli, data unknown
867                        iBin_qxqy[ii] = 0
868                        eBin_qxqy[ii] = 1
869                        eBin2D_qxqy[ii] = NaN
870                else
871                        if(nBin_qxqy[ii] <= 1)
872                                //need more than one pixel to determine error
873                                iBin_qxqy[ii] /= nBin_qxqy[ii]
874                                eBin_qxqy[ii] = 1
875                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
876                        else
877                                //assume that the intensity in each pixel in annuli is normally distributed about mean...
878                                iBin_qxqy[ii] /= nBin_qxqy[ii]
879                                avesq = iBin_qxqy[ii]^2
880                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
881                                var = aveisq-avesq
882                                if(var<=0)
883                                        eBin_qxqy[ii] = 1e-6
884                                else
885                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
886                                endif
887                                // and calculate as it is propagated pixel-by-pixel
888                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
889                        endif
890                endif
891        endfor
892       
893        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo
894       
895        // find the last non-zero point, working backwards
896        val=nq
897        do
898                val -= 1
899        while((nBin_qxqy[val] == 0) && val > 0)
900       
901//      print val, nBin_qxqy[val]
902        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
903
904        if(val == 0)
905                // all the points were deleted
906                return(0)
907        endif
908       
909       
910        // since the beam center is not always on the detector, many of the low Q bins will have zero pixels
911        // find the first non-zero point, working forwards
912        val = -1
913        do
914                val += 1
915        while(nBin_qxqy[val] == 0)     
916        DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
917
918        // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
919        val = numpnts(nBin_qxqy)-1
920        do
921                if(nBin_qxqy[val] == 0)
922                        DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
923                endif
924                val -= 1
925        while(val>0)
926       
927        SetDataFolder root:
928       
929        return(0)
930End
931
932
Note: See TracBrowser for help on using the repository browser.