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

Last change on this file since 1249 was 1249, checked in by srkline, 2 years ago

minor changes - can't find the diff command in TortoiseSVN...

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