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

Last change on this file since 1134 was 1134, checked in by srkline, 4 years ago

more changes to VCALC functionality and bug fixes for VCALC

File size: 47.3 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        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        tmpSig = sqrt(tmpInten)         // corrected based on John's memo, from 8/9/99
162
163        tmpInten += gnoise(tmpSig)
164        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        qTot = VC_CalcQval(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
208        qx = VC_CalcQX(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
209        qy = VC_CalcQY(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
210        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//
413Function 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//
440Function 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//
470Function 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//
499Function 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        Variable isVCALC
710        if(cmpstr(folderStr,"VCALC") == 0)
711                isVCALC = 1
712        endif
713       
714        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
715        String instPath = ":entry:instrument:detector_"
716               
717        if(isVCALC)
718                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)               // 2D detector data
719        else
720                Wave inten = V_getDetectorDataW(folderStr,detStr)
721        endif
722
723        Wave qx = $(folderPath+instPath+detStr+":qx_"+detStr)
724        Wave qy = $(folderPath+instPath+detStr+":qy_"+detStr)
725       
726        Variable xDim,yDim,delQ
727       
728        xDim=DimSize(inten,0)
729        yDim=DimSize(inten,1)
730       
731        if(xDim<yDim)
732                delQ = abs(qx[0][0] - qx[1][0])/2
733        else
734                delQ = abs(qy[0][1] - qy[0][0])/2
735        endif
736
737// TODOHIGHRES
738// -- is this how I want to handle the too-fine resolution of 1x1 binning?
739        NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
740
741        if(cmpstr(detStr,"B") == 0 && gHighResBinning == 1)
742                delQ = 4*delQ
743                Print "Reset delta Q for binning the back detector to 4x pix = ",delQ
744        endif
745       
746        // set the global
747        Variable/G $(folderPath+instPath+detStr+":gDelQ_"+detStr) = delQ
748//      Print "SET delQ = ",delQ," for ",type
749       
750        return(delQ)
751end
752
753
754//TODO -- need a switch here to dispatch to the averaging type
755Proc VC_BinQxQy_to_1D(folderStr,type)
756        String folderStr
757        String type
758//      Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4)
759//      Prompt type,"detector identifier"
760
761
762        VC_fDoBinning_QxQy2D(folderStr, type)
763
764
765/// this is for a tall, narrow slit mode       
766//      VC_fBinDetector_byRows(folderStr,type)
767       
768End
769
770
771// folderStr is RAW, VCALC, SAM, etc.
772// type is "B", "FL" for single binning, "FLR", or "MLRTB" or similar if multiple panels are combined
773//
774Proc VC_Graph_1D_detType(folderStr,type)
775        String folderStr,type
776       
777        SetDataFolder $("root:Packages:NIST:VSANS:"+folderStr)
778       
779        Display $("iBin_qxqy"+"_"+type) vs $("qBin_qxqy"+"_"+type)
780        ModifyGraph mirror=2,grid=1,log=1
781        ModifyGraph mode=4,marker=19,msize=2
782//      ErrorBars/T=0 iBin_qxqy Y,wave=(eBin2D_qxqy,eBin2D_qxqy)                // for simulations, I don't have 2D uncertainty
783        ErrorBars/T=0 $("iBin_qxqy"+"_"+type) Y,wave=($("eBin_qxqy"+"_"+type),$("eBin_qxqy"+"_"+type))
784        legend
785       
786        SetDataFolder root:
787
788End
789
790
791
792//////////
793//
794//              Function that bins a 2D detctor panel into I(q) based on the q-value of the pixel
795//              - each pixel QxQyQz has been calculated beforehand
796//              - if multiple panels are selected to be combined, it is done here during the binning
797//              - the setting of deltaQ step is still a little suspect (TODO)
798//
799//
800// see the equivalent function in PlotUtils2D_v40.ipf
801//
802//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
803//
804// this has been modified to accept different detector panels and to take arrays
805// -- type = FL or FR or...other panel identifiers
806//
807// TODO "iErr" is not always defined correctly since it doesn't really apply here for data that is not 2D simulation
808//
809//
810// updated Feb2016 to take new folder structure
811// TODO
812// -- VERIFY
813// -- figure out what the best location is to put the averaged data? currently @ top level of WORK folder
814//    but this is a lousy choice.
815// x- binning is now Mask-aware. If mask is not present, all data is used. If data is from VCALC, all data is used
816// x- Where do I put the solid angle correction? In here as a weight for each point, or later on as
817//    a blanket correction (matrix multiply) for an entire panel? (Solid Angle correction is done in the
818//    step where data is added to a WORK file (see Raw_to_Work())
819//
820//
821// TODO:
822// -- some of the input parameters for the resolution calcuation are either assumed (apOff) or are currently
823//    hard-wired. these need to be corrected before even the pinhole resolution is correct
824// x- resolution calculation is in the correct place. The calculation is done per-panel (specified by TYPE),
825//    and then the unwanted points can be discarded (all 6 columns) as the data is trimmed and concatenated
826//    is separate functions that are resolution-aware.
827//
828//
829// folderStr = WORK folder, type = the binning type (may include multiple detectors)
830Function VC_fDoBinning_QxQy2D(folderStr,type,collimationStr)
831        String folderStr,type,collimationStr
832       
833        Variable nSets = 0
834        Variable xDim,yDim
835        Variable ii,jj
836        Variable qVal,nq,var,avesq,aveisq
837        Variable binIndex,val,isVCALC=0,maskMissing
838
839        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
840        String instPath = ":entry:instrument:detector_"
841        String detStr
842               
843        if(cmpstr(folderStr,"VCALC") == 0)
844                isVCALC = 1
845        endif
846       
847        detStr = type
848
849// now switch on the type to determine which waves to declare and create
850// since there may be more than one panel to step through. There may be two, there may be four
851//
852
853// TODO:
854// -- Solid_Angle -- waves will be present for WORK data other than RAW, but not for RAW
855//
856// assume that the mask files are missing unless we can find them. If VCALC data,
857//  then the Mask is missing by definition
858        maskMissing = 1
859
860        strswitch(type) // string switch
861//              case "FL":              // execute if case matches expression
862//              case "FR":
863//                      detStr = type
864//                      if(isVCALC)
865//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
866//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation
867//                      else
868//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
869//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
870//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
871//                              if(WaveExists(mask) == 1)
872//                                      maskMissing = 0
873//                              endif
874//                             
875//                      endif
876//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
877//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
878//                      nSets = 1
879//                      break   
880                                                               
881//              case "FT":             
882//              case "FB":
883//                      detStr = type
884//                      if(isVCALC)
885//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
886//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
887//                      else
888//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
889//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
890//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
891//                              if(WaveExists(mask) == 1)
892//                                      maskMissing = 0
893//                              endif
894//                      endif
895//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
896//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
897//                      nSets = 1
898//                      break
899                       
900//              case "ML":             
901//              case "MR":
902//                      detStr = type
903//                      if(isVCALC)
904//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
905//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
906//                      else
907//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
908//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
909//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
910//                              if(WaveExists(mask) == 1)
911//                                      maskMissing = 0
912//                              endif
913//                      endif   
914//                      //TODO:
915//                      // -- decide on the proper deltaQ for binning. either nominal value for LR, or one
916//                      //    determined specifically for that panel (currently using one tube width as deltaQ)
917//                      // -- this is repeated multiple times in this switch
918//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
919//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
920//                      nSets = 1
921//                      break   
922                                       
923//              case "MT":             
924//              case "MB":
925//                      detStr = type
926//                      if(isVCALC)
927//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
928//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
929//                      else
930//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
931//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
932//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
933//                              if(WaveExists(mask) == 1)
934//                                      maskMissing = 0
935//                              endif
936//                      endif   
937//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
938//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
939//                      nSets = 1
940//                      break   
941
942// only one panel, simply pick that panel and move on out of the switch
943                case "FL":
944                case "FR":
945                case "FT":
946                case "FB":
947                case "ML":
948                case "MR":
949                case "MT":
950                case "MB":                     
951                case "B":       
952                        if(isVCALC)
953                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
954                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
955                        else
956                                Wave inten = V_getDetectorDataW(folderStr,detStr)
957                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
958                        endif   
959                        Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
960                        if(WaveExists(mask) == 1)
961                                maskMissing = 0
962                        endif
963                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
964                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values 
965                        nSets = 1
966                        break   
967                       
968                case "FLR":
969                // detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
970                // TODO
971                // -- see if I can un-hard-wire some of this below when more than one panel is combined
972                        if(isVCALC)
973                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
974                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
975                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
976                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
977                        else
978                                Wave inten = V_getDetectorDataW(folderStr,"FL")
979                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
980                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
981                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
982                        endif   
983                        Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
984                        Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
985                        if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
986                                maskMissing = 0
987                        endif
988                       
989                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
990                       
991                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
992                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
993               
994                        nSets = 2
995                        break                   
996               
997                case "FTB":
998                        if(isVCALC)
999                                WAVE inten = $(folderPath+instPath+"FT"+":det_"+"FT")
1000                                WAVE/Z iErr = $("iErr_"+"FT")                   // 2D errors -- may not exist, especially for simulation               
1001                                WAVE inten2 = $(folderPath+instPath+"FB"+":det_"+"FB")
1002                                WAVE/Z iErr2 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
1003                        else
1004                                Wave inten = V_getDetectorDataW(folderStr,"FT")
1005                                Wave iErr = V_getDetectorDataErrW(folderStr,"FT")
1006                                Wave inten2 = V_getDetectorDataW(folderStr,"FB")
1007                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FB")
1008                        endif   
1009                        Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
1010                        Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
1011                        if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1012                                maskMissing = 0
1013                        endif
1014                       
1015                        NVAR delQ = $(folderPath+instPath+"FT"+":gDelQ_FT")
1016                       
1017                        Wave qTotal = $(folderPath+instPath+"FT"+":qTot_"+"FT")                 // 2D q-values 
1018                        Wave qTotal2 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
1019       
1020                        nSets = 2
1021                        break           
1022               
1023                case "FLRTB":
1024                        if(isVCALC)
1025                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
1026                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
1027                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
1028                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
1029                                WAVE inten3 = $(folderPath+instPath+"FT"+":det_"+"FT")
1030                                WAVE/Z iErr3 = $("iErr_"+"FT")                  // 2D errors -- may not exist, especially for simulation               
1031                                WAVE inten4 = $(folderPath+instPath+"FB"+":det_"+"FB")
1032                                WAVE/Z iErr4 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
1033                        else
1034                                Wave inten = V_getDetectorDataW(folderStr,"FL")
1035                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
1036                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
1037                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
1038                                Wave inten3 = V_getDetectorDataW(folderStr,"FT")
1039                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"FT")
1040                                Wave inten4 = V_getDetectorDataW(folderStr,"FB")
1041                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"FB")
1042                               
1043                        endif   
1044                        Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
1045                        Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
1046                        Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
1047                        Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
1048                        if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
1049                                maskMissing = 0
1050                        endif
1051                       
1052                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
1053                       
1054                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
1055                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
1056                        Wave qTotal3 = $(folderPath+instPath+"FT"+":qTot_"+"FT")                        // 2D q-values 
1057                        Wave qTotal4 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
1058               
1059                        nSets = 4
1060                        break           
1061                       
1062                case "MLR":
1063                        if(isVCALC)
1064                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
1065                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
1066                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
1067                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
1068                        else
1069                                Wave inten = V_getDetectorDataW(folderStr,"ML")
1070                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
1071                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
1072                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
1073                               
1074                        endif   
1075                        Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
1076                        Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
1077                        if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1078                                maskMissing = 0
1079                        endif
1080                       
1081                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
1082                       
1083                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
1084                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
1085               
1086                        nSets = 2
1087                        break                   
1088               
1089                case "MTB":
1090                        if(isVCALC)
1091                                WAVE inten = $(folderPath+instPath+"MT"+":det_"+"MT")
1092                                WAVE/Z iErr = $("iErr_"+"MT")                   // 2D errors -- may not exist, especially for simulation               
1093                                WAVE inten2 = $(folderPath+instPath+"MB"+":det_"+"MB")
1094                                WAVE/Z iErr2 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
1095                        else
1096                                Wave inten = V_getDetectorDataW(folderStr,"MT")
1097                                Wave iErr = V_getDetectorDataErrW(folderStr,"MT")
1098                                Wave inten2 = V_getDetectorDataW(folderStr,"MB")
1099                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MB")
1100                               
1101                        endif   
1102                        Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
1103                        Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
1104                        if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1105                                maskMissing = 0
1106                        endif
1107                       
1108                        NVAR delQ = $(folderPath+instPath+"MT"+":gDelQ_MT")
1109                       
1110                        Wave qTotal = $(folderPath+instPath+"MT"+":qTot_"+"MT")                 // 2D q-values 
1111                        Wave qTotal2 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
1112               
1113                        nSets = 2
1114                        break                           
1115               
1116                case "MLRTB":
1117                        if(isVCALC)
1118                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
1119                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
1120                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
1121                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
1122                                WAVE inten3 = $(folderPath+instPath+"MT"+":det_"+"MT")
1123                                WAVE/Z iErr3 = $("iErr_"+"MT")                  // 2D errors -- may not exist, especially for simulation               
1124                                WAVE inten4 = $(folderPath+instPath+"MB"+":det_"+"MB")
1125                                WAVE/Z iErr4 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
1126                        else
1127                                Wave inten = V_getDetectorDataW(folderStr,"ML")
1128                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
1129                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
1130                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
1131                                Wave inten3 = V_getDetectorDataW(folderStr,"MT")
1132                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"MT")
1133                                Wave inten4 = V_getDetectorDataW(folderStr,"MB")
1134                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"MB")
1135                               
1136                        endif
1137                        Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
1138                        Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
1139                        Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
1140                        Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
1141                        if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
1142                                maskMissing = 0
1143                        endif
1144                               
1145                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
1146                       
1147                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
1148                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
1149                        Wave qTotal3 = $(folderPath+instPath+"MT"+":qTot_"+"MT")                        // 2D q-values 
1150                        Wave qTotal4 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
1151               
1152                        nSets = 4
1153                        break                                                                   
1154                                       
1155                default:
1156                        nSets = 0                                                       
1157                        Print "ERROR   ---- type is not recognized "
1158        endswitch
1159
1160//      Print "delQ = ",delQ," for ",type
1161
1162        if(nSets == 0)
1163                SetDataFolder root:
1164                return(0)
1165        endif
1166
1167
1168// RAW data is currently read in and the 2D error wave is correctly generated
1169// 2D error is propagated through all reduction steps, but I have not
1170// verified that it is an exact duplication of the 1D error
1171//
1172//
1173//
1174// IF ther is no 2D error wave present for some reason, make a fake one
1175        if(WaveExists(iErr)==0  && WaveExists(inten) != 0)
1176                Duplicate/O inten,iErr
1177                Wave iErr=iErr
1178//              iErr = 1+sqrt(inten+0.75)                       // can't use this -- it applies to counts, not intensity (already a count rate...)
1179                iErr = sqrt(inten+0.75)                 // TODO -- here I'm just using some fictional value
1180        endif
1181        if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0)
1182                Duplicate/O inten2,iErr2
1183                Wave iErr2=iErr2
1184//              iErr2 = 1+sqrt(inten2+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1185                iErr2 = sqrt(inten2+0.75)                       // TODO -- here I'm just using some fictional value
1186        endif
1187        if(WaveExists(iErr3)==0  && WaveExists(inten3) != 0)
1188                Duplicate/O inten3,iErr3
1189                Wave iErr3=iErr3
1190//              iErr3 = 1+sqrt(inten3+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1191                iErr3 = sqrt(inten3+0.75)                       // TODO -- here I'm just using some fictional value
1192        endif
1193        if(WaveExists(iErr4)==0  && WaveExists(inten4) != 0)
1194                Duplicate/O inten4,iErr4
1195                Wave iErr4=iErr4
1196//              iErr4 = 1+sqrt(inten4+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1197                iErr4 = sqrt(inten4+0.75)                       // TODO -- here I'm just using some fictional value
1198        endif
1199
1200        // TODO -- nq will need to be larger, once the back detector is installed
1201        //
1202        // note that the back panel of 320x320 (1mm res) results in 447 data points!
1203        // - so I upped nq to 600
1204
1205        if(cmpstr(type,"B") == 0)
1206                nq = 8000
1207        else
1208                nq=600
1209        endif
1210
1211//******TODO****** -- where to put the averaged data -- right now, folderStr is forced to ""   
1212//      SetDataFolder $("root:"+folderStr)              //should already be here, but make sure...     
1213        Make/O/D/N=(nq)  $(folderPath+":"+"iBin_qxqy"+"_"+type)
1214        Make/O/D/N=(nq)  $(folderPath+":"+"qBin_qxqy"+"_"+type)
1215        Make/O/D/N=(nq)  $(folderPath+":"+"nBin_qxqy"+"_"+type)
1216        Make/O/D/N=(nq)  $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1217        Make/O/D/N=(nq)  $(folderPath+":"+"eBin_qxqy"+"_"+type)
1218        Make/O/D/N=(nq)  $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1219       
1220        Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+type)
1221        Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type)
1222        Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy"+"_"+type)
1223        Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1224        Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy"+"_"+type)
1225        Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1226       
1227       
1228//      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
1229// TODO: not sure if I want to set dQ in x or y direction...
1230        // the short dimension is the 8mm tubes, use this direction as dQ?
1231        // 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]
1232        // WRT the beam center. use qx or qy directly. Still not happy with this way...
1233
1234
1235        qBin_qxqy[] =  p*delQ   
1236        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
1237
1238        iBin_qxqy = 0
1239        iBin2_qxqy = 0
1240        eBin_qxqy = 0
1241        eBin2D_qxqy = 0
1242        nBin_qxqy = 0   //number of intensities added to each bin
1243
1244// now there are situations of:
1245// 1 panel
1246// 2 panels
1247// 4 panels
1248//
1249// this needs to be a double loop now...
1250// TODO:
1251// -- the iErr (=2D) wave and accumulation of error is NOT CALCULATED CORRECTLY YET
1252// -- verify the 2D error propagation by reducing it to 1D error
1253//
1254//
1255// The 1D error does not use iErr, and IS CALCULATED CORRECTLY
1256//
1257// x- the solid angle per pixel will be present for WORK data other than RAW, but not for RAW
1258
1259//
1260// if any of the masks don't exist, display the error, and proceed with the averaging, using all data
1261        if(maskMissing == 1)
1262                Print "Mask file not found for at least one detector - so all data is used"
1263        endif
1264       
1265        NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
1266        if(gIgnoreDetB && cmpstr(type,"B") == 0)
1267                maskMissing = 1
1268                Print "Mask skipped for B due to possible mismatch (Panel B ignored in preferences)"
1269        endif
1270
1271        Variable mask_val
1272// use set 1 (no number) only
1273        if(nSets >= 1)
1274                xDim=DimSize(inten,0)
1275                yDim=DimSize(inten,1)
1276       
1277                for(ii=0;ii<xDim;ii+=1)
1278                        for(jj=0;jj<yDim;jj+=1)
1279                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1280                                qVal = qTotal[ii][jj]
1281                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1282                                val = inten[ii][jj]
1283                               
1284//                              if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
1285                                if(maskMissing)         // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
1286                                        mask_val = 0
1287                                else
1288                                        mask_val = mask[ii][jj]
1289                                endif
1290                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1291                                        iBin_qxqy[binIndex] += val
1292                                        iBin2_qxqy[binIndex] += val*val
1293                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1294                                        nBin_qxqy[binIndex] += 1
1295                                endif
1296                        endfor
1297                endfor
1298               
1299        endif
1300
1301// add in set 2 (set 1 already done)
1302        if(nSets >= 2)
1303                xDim=DimSize(inten2,0)
1304                yDim=DimSize(inten2,1)
1305       
1306                for(ii=0;ii<xDim;ii+=1)
1307                        for(jj=0;jj<yDim;jj+=1)
1308                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1309                                qVal = qTotal2[ii][jj]
1310                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1311                                val = inten2[ii][jj]
1312                               
1313//                              if(isVCALC || maskMissing)
1314                                if(maskMissing)
1315                                        mask_val = 0
1316                                else
1317                                        mask_val = mask2[ii][jj]
1318                                endif
1319                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1320                                        iBin_qxqy[binIndex] += val
1321                                        iBin2_qxqy[binIndex] += val*val
1322                                        eBin2D_qxqy[binIndex] += iErr2[ii][jj]*iErr2[ii][jj]
1323                                        nBin_qxqy[binIndex] += 1
1324                                endif
1325                        endfor
1326                endfor
1327               
1328        endif
1329
1330// add in set 3 and 4 (set 1 and 2 already done)
1331        if(nSets == 4)
1332                xDim=DimSize(inten3,0)
1333                yDim=DimSize(inten3,1)
1334       
1335                for(ii=0;ii<xDim;ii+=1)
1336                        for(jj=0;jj<yDim;jj+=1)
1337                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1338                                qVal = qTotal3[ii][jj]
1339                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1340                                val = inten3[ii][jj]
1341                               
1342//                              if(isVCALC || maskMissing)
1343                                if(maskMissing)
1344                                        mask_val = 0
1345                                else
1346                                        mask_val = mask3[ii][jj]
1347                                endif
1348                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1349                                        iBin_qxqy[binIndex] += val
1350                                        iBin2_qxqy[binIndex] += val*val
1351                                        eBin2D_qxqy[binIndex] += iErr3[ii][jj]*iErr3[ii][jj]
1352                                        nBin_qxqy[binIndex] += 1
1353                                endif
1354                        endfor
1355                endfor
1356               
1357               
1358                xDim=DimSize(inten4,0)
1359                yDim=DimSize(inten4,1)
1360       
1361                for(ii=0;ii<xDim;ii+=1)
1362                        for(jj=0;jj<yDim;jj+=1)
1363                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1364                                qVal = qTotal4[ii][jj]
1365                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1366                                val = inten4[ii][jj]
1367                               
1368//                              if(isVCALC || maskMissing)
1369                                if(maskMissing)
1370                                        mask_val = 0
1371                                else
1372                                        mask_val = mask4[ii][jj]
1373                                endif
1374                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1375                                        iBin_qxqy[binIndex] += val
1376                                        iBin2_qxqy[binIndex] += val*val
1377                                        eBin2D_qxqy[binIndex] += iErr4[ii][jj]*iErr4[ii][jj]
1378                                        nBin_qxqy[binIndex] += 1
1379                                endif
1380                        endfor
1381                endfor
1382               
1383        endif
1384
1385
1386// after looping through all of the data on the panels, calculate errors on I(q),
1387// just like in CircSectAve.ipf
1388// TODO:
1389// -- 2D Errors were (maybe) properly acculumated through reduction, so this loop of calculations is NOT VERIFIED (yet)
1390// x- the error on the 1D intensity, is correctly calculated as the standard error of the mean.
1391        for(ii=0;ii<nq;ii+=1)
1392                if(nBin_qxqy[ii] == 0)
1393                        //no pixels in annuli, data unknown
1394                        iBin_qxqy[ii] = 0
1395                        eBin_qxqy[ii] = 1
1396                        eBin2D_qxqy[ii] = NaN
1397                else
1398                        if(nBin_qxqy[ii] <= 1)
1399                                //need more than one pixel to determine error
1400                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1401                                eBin_qxqy[ii] = 1
1402                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1403                        else
1404                                //assume that the intensity in each pixel in annuli is normally distributed about mean...
1405                                //  -- this is correctly calculating the error as the standard error of the mean, as
1406                                //    was always done for SANS as well.
1407                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1408                                avesq = iBin_qxqy[ii]^2
1409                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1410                                var = aveisq-avesq
1411                                if(var<=0)
1412                                        eBin_qxqy[ii] = 1e-6
1413                                else
1414                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1415                                endif
1416                                // and calculate as it is propagated pixel-by-pixel
1417                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1418                        endif
1419                endif
1420        endfor
1421       
1422        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo
1423       
1424        // find the last non-zero point, working backwards
1425        val=nq
1426        do
1427                val -= 1
1428        while((nBin_qxqy[val] == 0) && val > 0)
1429       
1430//      print val, nBin_qxqy[val]
1431        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1432
1433        if(val == 0)
1434                // all the points were deleted
1435                return(0)
1436        endif
1437       
1438       
1439        // since the beam center is not always on the detector, many of the low Q bins will have zero pixels
1440        // find the first non-zero point, working forwards
1441        val = -1
1442        do
1443                val += 1
1444        while(nBin_qxqy[val] == 0)     
1445        DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1446
1447        // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
1448        val = numpnts(nBin_qxqy)-1
1449        do
1450                if(nBin_qxqy[val] == 0)
1451                        DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1452                endif
1453                val -= 1
1454        while(val>0)
1455
1456// utility function to remove NaN values from the waves
1457        V_RemoveNaNsQIS(qBin_qxqy, iBin_qxqy, eBin_qxqy)
1458
1459       
1460        // TODO:
1461        // -- This is where I calculate the resolution in SANS (see CircSectAve)
1462        // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
1463        // -- from the top of the function, folderStr = work folder, type = "FLRTB" or other type of averaging
1464        //
1465        nq = numpnts(qBin_qxqy)
1466        Make/O/D/N=(nq)  $(folderPath+":"+"sigmaQ_qxqy"+"_"+type)
1467        Make/O/D/N=(nq)  $(folderPath+":"+"qBar_qxqy"+"_"+type)
1468        Make/O/D/N=(nq)  $(folderPath+":"+"fSubS_qxqy"+"_"+type)
1469        Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+type)
1470        Wave qbar = $(folderPath+":"+"qBar_qxqy_"+type)
1471        Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+type)
1472        Variable ret1,ret2,ret3                 
1473
1474// all of the different collimation conditions are handled within the V_getResolution function
1475// which is responsible for switching based on the different collimation types (white beam, slit, Xtal, etc)
1476// to calculate the correct resolution, or fill the waves with the correct "flags"
1477//
1478
1479// For white beam data, the wavelength distribution can't be represented as a gaussian, but all of the other
1480//  geometric corrections still apply. Passing zero for the lambdaWidth will return the geometry contribution,
1481//  as long as the wavelength can be handled separately. It appears to be correct to do as a double integral,
1482//  with the inner(lambda) calculated first, then the outer(geometry).
1483//
1484
1485// possible values are:
1486//
1487// pinhole
1488// pinhole_whiteBeam
1489// convergingPinholes
1490//
1491// *slit data should be reduced using the slit routine, not here, proceed but warn
1492// narrowSlit
1493// narrowSlit_whiteBeam
1494//
1495        ii=0
1496        do
1497                V_getResolution(qBin_qxqy[ii],folderStr,type,collimationStr,ret1,ret2,ret3)
1498                sigmaq[ii] = ret1       
1499                qbar[ii] = ret2
1500                fsubs[ii] = ret3       
1501                ii+=1
1502        while(ii<nq)
1503
1504
1505        SetDataFolder root:
1506       
1507        return(0)
1508End
1509
1510
Note: See TracBrowser for help on using the repository browser.