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

Last change on this file since 1192 was 1192, checked in by srkline, 3 years ago

added comments in NCNR_Utils identifying the proper instrument prefix values for each instrument

Saving individual data sets from VSANS will now properly ignore the back detector (if preference set) and will handle all of the current binning types.

In the preferences, the bin width (as a multiple of pixels) is now active. This setting is a multiple of the natural deltaQ of one tube width.

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