source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_DetectorCorrections.ipf @ 1024

Last change on this file since 1024 was 1024, checked in by srkline, 6 years ago

minor changes to prefix functions with "V_" to avoid conflicts with non-VSANS functions.

File size: 36.0 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2#pragma version=1.0
3#pragma IgorVersion=6.1
4
5
6
7//
8// functions to apply corrections to the detector panels
9//
10// these are meant to be called by the procedures that convert "raw" data to
11// "adjusted" or corrected data sets
12//
13// may be relocated in the future
14//
15
16
17
18//
19// detector dead time
20//
21// input is the data array (N tubes x M pixels)
22// input of N x 1 array of dead time values
23//
24// output is the corrected counts in data, overwriting the input data
25//
26// Note that the equation in Roe (eqn 2.15, p. 63) looks different, but it is really the
27// same old equation, just written in a more complex form.
28//
29// TODO
30// -- verify the direction of the tubes and indexing
31// x- decide on the appropriate functional form for the tubes
32// x- need count time as input
33// -- be sure I'm working in the right data folder
34// -- clean up when done
35// -- calculate + return the error contribution?
36// -- verify the error propagation
37Function V_DeadTimeCorrectionTubes(dataW,data_errW,dtW,ctTime)
38        Wave dataW,data_errW,dtW
39        Variable ctTime
40       
41        // do I count on the orientation as an input, or do I just figure it out on my own?
42        String orientation
43        Variable dimX,dimY
44        dimX = DimSize(dataW,0)
45        dimY = DimSize(dataw,1)
46        if(dimX > dimY)
47                orientation = "horizontal"
48        else
49                orientation = "vertical"
50        endif
51       
52        // sum the counts in each tube and divide by time for total cr per tube
53        Variable npt
54       
55        if(cmpstr(orientation,"vertical")==0)
56                //      this is data dimensioned as (Ntubes,Npix)
57               
58                MatrixOp/O sumTubes = sumRows(dataW)            // n x 1 result
59                sumTubes /= ctTime              //now count rate per tube
60               
61                dataW[][] = dataW[p][q]/(1-sumTubes[p]*dtW[p])          //correct the data
62                data_errW[][] = data_errW[p][q]/(1-sumTubes[p]*dtW[p])          // propagate the error wave
63
64        elseif(cmpstr(orientation,"horizontal")==0)
65        //      this is data (horizontal) dimensioned as (Npix,Ntubes)
66
67                MatrixOp/O sumTubes = sumCols(dataW)            // 1 x m result
68                sumTubes /= ctTime
69               
70                dataW[][] = dataW[p][q]/(1-sumTubes[q]*dtW[q])
71                data_errW[][] = data_errW[p][q]/(1-sumTubes[q]*dtW[q])
72       
73        else           
74                DoAlert 0,"Orientation not correctly passed in DeadTimeCorrectionTubes(). No correction done."
75        endif
76       
77        return(0)
78end
79
80// test function
81Function V_testDTCor()
82
83        String detStr = ""
84        String fname = "RAW"
85        Variable ctTime
86       
87        detStr = "FR"
88        Wave w = V_getDetectorDataW(fname,detStr)
89        Wave w_err = V_getDetectorDataErrW(fname,detStr)
90        Wave w_dt = V_getDetector_deadtime(fname,detStr)
91
92        ctTime = V_getCount_time(fname)
93       
94//      ctTime = 10
95        V_DeadTimeCorrectionTubes(w,w_err,w_dt,ctTime)
96
97End
98
99
100//
101// Non-linear data correction
102//
103// DOES NOT modify the data, only calculates the spatial relationship
104//
105// input is the data array (N tubes x M pixels)
106// input of N x M array of quadratic coefficients
107//
108// output is wave of corrected real space distance corresponding to each pixel of the data
109//
110//
111// TODO
112// -- UNITS!!!! currently this is mm, which certainly doesn't match anything else!!!
113// -- verify the direction of the tubes and indexing
114// -- be sure I'm working in the right data folder
115// -- clean up when done
116// -- calculate + return the error contribution?
117// -- do I want this to return a wave?
118// -- do I need to write a separate function that returns the distance wave for later calculations?
119// -- do I want to make the distance array 3D to keep the x and y dims together? Calculate them all right now?
120// -- what else do I need to pass to the function? (fname=folder? detStr?)
121// -- need a separate block or function to handle "B" detector which will be ? different
122//
123//
124Function V_NonLinearCorrection(dataW,coefW,tube_width,detStr,destPath)
125        Wave dataW,coefW
126        Variable tube_width
127        String detStr,destPath
128       
129         
130        // do I count on the orientation as an input, or do I just figure it out on my own?
131        String orientation
132        Variable dimX,dimY
133        dimX = DimSize(dataW,0)
134        dimY = DimSize(dataW,1)
135        if(dimX > dimY)
136                orientation = "horizontal"
137        else
138                orientation = "vertical"
139        endif
140
141        // make a wave of the same dimensions, in the same data folder for the distance
142        // ?? or a 3D wave?
143        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
144        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
145        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
146        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
147       
148        // then per tube, do the quadratic calculation to get the real space distance along the tube
149        // the distance perpendicular to the tube is n*(8.4mm) per tube index
150       
151        if(cmpstr(orientation,"vertical")==0)
152                //      this is data dimensioned as (Ntubes,Npix)
153                data_realDistX[][] = tube_width*p
154                data_realDistY[][] = coefW[0][p] + coefW[1][p]*q + coefW[2][p]*q*q
155       
156        elseif(cmpstr(orientation,"horizontal")==0)
157                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
158                data_realDistX[][] = coefW[0][q] + coefW[1][q]*p + coefW[2][q]*p*p
159                data_realDistY[][] = tube_width*q
160
161        else           
162                DoAlert 0,"Orientation not correctly passed in NonLinearCorrection(). No correction done."
163        endif
164       
165        return(0)
166end
167
168// TODO:
169// -- the cal_x and y coefficients are totally fake
170// -- the wave assignment may not be correct.. so beware
171//
172//
173Function V_NonLinearCorrection_B(folder,detStr,destPath)
174        String folder,detStr,destPath
175
176        if(cmpstr(detStr,"B") != 0)
177                return(0)
178        endif
179       
180        // do I count on the orientation as an input, or do I just figure it out on my own?
181        Variable dimX,dimY
182       
183        Wave dataW = V_getDetectorDataW(folder,detStr)
184       
185        dimX = DimSize(dataW,0)
186        dimY = DimSize(dataW,1)
187
188        // make a wave of the same dimensions, in the same data folder for the distance
189        // ?? or a 3D wave?
190        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
191        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
192        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
193        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
194       
195       
196        Wave cal_x = V_getDet_cal_x(folder,detStr)
197        Wave cal_y = V_getDet_cal_y(folder,detStr)
198       
199        data_realDistX[][] = cal_x[0]*p
200        data_realDistY[][] = cal_y[0]*q
201       
202        return(0)
203end
204
205
206//
207//
208// TODO
209// -- VERIFY the calculations
210// -- verify where this needs to be done (if the beam center is changed)
211// -- then the q-calculation needs to be re-done
212// -- the position along the tube length is referenced to tube[0], for no particular reason
213//    It may be better to take an average? but [0] is an ASSUMPTION
214// -- distance along tube is simple interpolation, or do I use the coefficients to
215//    calculate the actual value
216//
217// -- distance in the lateral direction is based on tube width, which is a fixed parameter
218//
219//
220Function V_ConvertBeamCtr_to_mm(folder,detStr,destPath)
221        String folder,detStr,destPath
222       
223        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
224        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")     
225
226        String orientation
227        Variable dimX,dimY,xCtr,yCtr
228        dimX = DimSize(data_realDistX,0)
229        dimY = DimSize(data_realDistX,1)
230        if(dimX > dimY)
231                orientation = "horizontal"
232        else
233                orientation = "vertical"
234        endif
235       
236        xCtr = V_getDet_beam_center_x(folder,detStr)
237        yCtr = V_getDet_beam_center_y(folder,detStr)   
238       
239        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
240        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
241        WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
242        WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
243
244        Variable tube_width = V_getDet_tubeWidth(folder,detStr)
245
246//
247        if(cmpstr(orientation,"vertical")==0)
248                //      this is data dimensioned as (Ntubes,Npix)
249//              data_realDistX[][] = tube_width*p
250//              data_realDistY[][] = coefW[0][p] + coefW[1][p]*q + coefW[2][p]*q*q
251                x_mm[0] = tube_width*xCtr
252                y_mm[0] = data_realDistY[0][yCtr]
253        else
254                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
255//              data_realDistX[][] = coefW[0][q] + coefW[1][q]*p + coefW[2][q]*p*p
256//              data_realDistY[][] = tube_width*q
257                x_mm[0] = data_realDistX[xCtr][0]
258                y_mm[0] = tube_width*yCtr
259        endif
260               
261        return(0)
262end
263
264
265//
266//
267// TODO
268// -- VERIFY the calculations
269// -- verify where this needs to be done (if the beam center is changed)
270// -- then the q-calculation needs to be re-done
271// -- the position along the tube length is referenced to tube[0], for no particular reason
272//    It may be better to take an average? but [0] is an ASSUMPTION
273// -- distance along tube is simple interpolation
274//
275// -- distance in the lateral direction is based on tube width, which is well known
276//
277//
278Function V_ConvertBeamCtr_to_mmB(folder,detStr,destPath)
279        String folder,detStr,destPath
280       
281        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
282        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")     
283       
284        Variable xCtr,yCtr
285        xCtr = V_getDet_beam_center_x(folder,detStr)
286        yCtr = V_getDet_beam_center_y(folder,detStr)   
287       
288        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
289        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
290        WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
291        WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
292
293        x_mm[0] = data_realDistX[xCtr][0]
294        y_mm[0] = data_realDistY[0][yCtr]
295               
296        return(0)
297end
298
299
300
301
302
303
304/////
305//
306// non-linear corrections to the tube pixels
307// - returns the distance in mm (although this may change)
308//
309// c0,c1,c2,pix
310// c0-c2 are the fit coefficients
311// pix is the test pixel
312//
313// returns the distance in mm (relative to ctr pixel)
314// ctr is the center pixel, as defined when fitting to quadratic was done
315//
316Function V_TubePixel_to_mm(c0,c1,c2,pix)
317        Variable c0,c1,c2,pix
318       
319        Variable dist
320        dist = c0 + c1*pix + c2*pix*pix
321       
322        return(dist)
323End
324//
325////
326
327// TODO
328// get rid of this in the real data
329//
330// TESTING ONLY
331Proc V_MakeFakeCalibrationWaves()
332        // make these in the RAW data folder, before converting to a work folder
333        // - then they will be "found" by get()
334        // -- only for the tube, not the Back det
335       
336//      DoAlert 0, "re-do this and do a better job of filling the fake calibration data"
337
338        DoAlert 0, "Calibration waves are read in from the data file"
339       
340//      V_fMakeFakeCalibrationWaves()
341End
342
343
344
345// TODO
346// get rid of this in the real data
347//
348// TESTING ONLY
349//
350// orientation does not matter, there are 48 tubes in each bank
351// so dimension (3,48) for everything.
352//
353// -- but the orientation does indicate TB vs LR, which has implications for
354//  the (fictional) dimension of the pixel along the tube axis, at least as far
355// as for making the fake coefficients.
356//
357Function V_fMakeFakeCalibrationWaves()
358
359        Variable ii,pixSize
360        String detStr,fname="RAW",orientation
361       
362        for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
363                detStr = StringFromList(ii, ksDetectorListNoB, ";")
364//              Wave w = V_getDetectorDataW(fname,detStr)
365                Make/O/D/N=(3,48) $("root:Packages:NIST:VSANS:RAW:entry:instrument:detector_"+detStr+":spatial_calibration")
366                Wave calib = $("root:Packages:NIST:VSANS:RAW:entry:instrument:detector_"+detStr+":spatial_calibration")
367                // !!!! this overwrites what is there
368
369                orientation = V_getDet_tubeOrientation(fname,detStr)
370                if(cmpstr(orientation,"vertical")==0)
371                //      this is vertical tube data dimensioned as (Ntubes,Npix)
372                        pixSize = 8.4           //V_getDet_y_pixel_size(fname,detStr)
373                       
374                elseif(cmpstr(orientation,"horizontal")==0)
375                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
376                        pixSize = 4                     //V_getDet_x_pixel_size(fname,detStr)
377                       
378                else           
379                        DoAlert 0,"Orientation not correctly passed in NonLinearCorrection(). No correction done."
380                endif
381               
382                calib[0][] = -(128/2)*pixSize                   //approx (n/2)*pixSixe
383                calib[1][] = pixSize
384                calib[2][] = 2e-4
385               
386        endfor
387       
388        // now fake calibration for "B"
389        Wave cal_x = V_getDet_cal_x("RAW","B")
390        Wave cal_y = V_getDet_cal_y("RAW","B")
391       
392        cal_x = 1               // mm, ignore the other 2 values
393        cal_y = 1               // mm
394        return(0)
395End
396
397//
398// TODO:
399// -- MUST VERIFY the definition of SDD and how (if) setback is written to the data files
400// -- currently I'm assuming that the SDD is the "nominal" value which is correct for the
401//    L/R panels, but is not correct for the T/B panels (must add in the setback)
402//
403//
404//
405// data_realDistX, Y must be previously generated from running NonLinearCorrection()
406//
407// call with:
408// fname as the WORK folder, "RAW"
409// detStr = detector String, "FL"
410// destPath = path to destination WORK folder ("root:Packages:NIST:VSANS:"+folder)
411//
412Function V_Detector_CalcQVals(fname,detStr,destPath)
413        String fname,detStr,destPath
414
415        String orientation
416        Variable xCtr,yCtr,lambda,sdd
417       
418// get all of the geometry information 
419        orientation = V_getDet_tubeOrientation(fname,detStr)
420//      sdd = V_getDet_distance(fname,detStr)           //[cm]
421//      sdd += V_getDet_TBSetback(fname,detStr)/10              // written in [mm], convert to [cm], returns 0 for L/R/B panels
422
423        sdd = V_getDet_ActualDistance(fname,detStr)             //sdd derived, including setback [cm]
424        sdd/=100                // sdd reported in cm, pass in m
425        // this is the ctr in pixels
426//      xCtr = V_getDet_beam_center_x(fname,detStr)
427//      yCtr = V_getDet_beam_center_y(fname,detStr)
428        // this is ctr in mm
429        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
430        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
431        lambda = V_getWavelength(fname)
432        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
433        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
434
435// make the new waves
436        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
437        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
438        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
439        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
440        Wave qTot = $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
441        Wave qx = $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
442        Wave qy = $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
443        Wave qz = $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
444
445// calculate all of the q-values
446        qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
447        qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
448        qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
449        qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
450       
451       
452        return(0)
453End
454
455
456//function to calculate the overall q-value, given all of the necesary trig inputs
457//
458// TODO:
459// -- verify the calculation (accuracy - in all input conditions)
460// -- verify the units of everything here, it's currently all jumbled and wrong... and repeated
461// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
462//    to each pixel
463//
464//sdd is in meters
465//wavelength is in Angstroms
466//
467//returned magnitude of Q is in 1/Angstroms
468//
469Function V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
470        Variable xaxval,yaxval,xctr,yctr,sdd,lam
471        Wave distX,distY
472       
473        Variable dx,dy,qval,two_theta,dist
474               
475        sdd *=100               //convert to cm
476        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
477        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
478        dist = sqrt(dx^2 + dy^2)
479       
480        dist /= 10  // convert mm to cm
481       
482        two_theta = atan(dist/sdd)
483
484        qval = 4*Pi/lam*sin(two_theta/2)
485       
486        return qval
487End
488
489//calculates just the q-value in the x-direction on the detector
490// TODO:
491// -- verify the calculation (accuracy - in all input conditions)
492// -- verify the units of everything here, it's currently all jumbled and wrong... and repeated
493// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
494//    to each pixel
495//
496//
497// this properly accounts for qz
498//
499Function V_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
500        Variable xaxval,yaxval,xctr,yctr,sdd,lam
501        Wave distX,distY
502
503        Variable qx,qval,phi,dx,dy,dist,two_theta
504       
505        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
506       
507        sdd *=100               //convert to cm
508        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
509        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
510        phi = V_FindPhi(dx,dy)
511       
512        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
513        dist = sqrt(dx^2 + dy^2)
514        dist /= 10  // convert mm to cm
515
516        two_theta = atan(dist/sdd)
517
518        qx = qval*cos(two_theta/2)*cos(phi)
519       
520        return qx
521End
522
523//calculates just the q-value in the y-direction on the detector
524// TODO:
525// -- verify the calculation (accuracy - in all input conditions)
526// -- verify the units of everything here, it's currently all jumbled and wrong... and repeated
527// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
528//    to each pixel
529//
530//
531// this properly accounts for qz
532//
533Function V_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
534        Variable xaxval,yaxval,xctr,yctr,sdd,lam
535        Wave distX,distY
536
537        Variable qy,qval,phi,dx,dy,dist,two_theta
538       
539        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
540       
541        sdd *=100               //convert to cm
542        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
543        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
544        phi = V_FindPhi(dx,dy)
545       
546        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
547        dist = sqrt(dx^2 + dy^2)
548        dist /= 10  // convert mm to cm
549
550        two_theta = atan(dist/sdd)
551
552        qy = qval*cos(two_theta/2)*sin(phi)
553       
554        return qy
555End
556
557//calculates just the q-value in the z-direction on the detector
558// TODO:
559// -- verify the calculation (accuracy - in all input conditions)
560// -- verify the units of everything here, it's currently all jumbled and wrong... and repeated
561// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
562//    to each pixel
563//
564// not actually used for any calculations, but here for completeness if anyone asks, or for 2D data export
565//
566// this properly accounts for qz
567//
568Function V_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
569        Variable xaxval,yaxval,xctr,yctr,sdd,lam
570        Wave distX,distY
571
572        Variable qz,qval,phi,dx,dy,dist,two_theta
573       
574        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
575       
576        sdd *=100               //convert to cm
577        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
578        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
579       
580        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
581        dist = sqrt(dx^2 + dy^2)
582        dist /= 10  // convert mm to cm
583
584        two_theta = atan(dist/sdd)
585
586        qz = qval*sin(two_theta/2)
587       
588        return qz
589End
590
591
592//
593// TODO -- VERIFY calculations
594// -- This is the actual solid angle per pixel, not a ratio vs. some "unit SA"
595//    Do I just correct for the different area vs. the "nominal" central area?
596// -- decide how to implement - either directly change the data values (as was done in the past)
597//    or use this as a weighting for when the data is binned to I(q). In the second method, 2D data
598//    would need this to be applied before exporting
599// -- do I keep a wave note indicating that this correction has been applied to the data
600//    so that it can be "un-applied"?
601// -- do I calculate theta from geometry directly, or get it from Q (Assuming it's present?)
602//    (probably just from geometry, since I need SDD and dx and dy values...)
603//
604//
605Function V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
606        Wave w,w_err
607        String fname,detStr,destPath
608
609        Variable sdd,xCtr,yCtr,lambda
610
611// get all of the geometry information 
612//      orientation = V_getDet_tubeOrientation(fname,detStr)
613        sdd = V_getDet_ActualDistance(fname,detStr)
614        sdd/=100                // sdd in cm, pass in m
615
616        // this is ctr in mm
617        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
618        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
619        lambda = V_getWavelength(fname)
620       
621        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
622       
623        Wave data_realDistX = data_realDistX
624        Wave data_realDistY = data_realDistY
625
626        Duplicate/O w solid_angle,tmp_theta,tmp_dist            //in the current df
627
628//// calculate the scattering angle
629//      dx = (distX - xctr)             //delta x in mm
630//      dy = (distY - yctr)             //delta y in mm
631        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
632       
633        tmp_dist /= 10  // convert mm to cm
634        sdd *=100               //convert to cm
635
636        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle
637
638        Variable ii,jj,numx,numy,dx,dy
639        numx = DimSize(tmp_theta,0)
640        numy = DimSize(tmp_theta,1)
641       
642        for(ii=0        ;ii<numx;ii+=1)
643                for(jj=0;jj<numy;jj+=1)
644                       
645                        if(ii==0)               //do a forward difference if ii==0
646                                dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj])        //delta x for the pixel
647                        else
648                                dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj])        //delta x for the pixel
649                        endif
650                       
651                       
652                        if(jj==0)
653                                dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj])        //delta y for the pixel
654                        else
655                                dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1])        //delta y for the pixel
656                        endif
657       
658                        dx /= 10
659                        dy /= 10                // convert mm to cm (since sdd is in cm)
660                        solid_angle[ii][jj] = dx*dy             //this is in cm^2
661                endfor
662        endfor
663       
664        // to cover up any issues w/negative dx or dy
665        solid_angle = abs(solid_angle)
666       
667        // solid_angle correction
668        // == dx*dy*cos^3/sdd^2
669        solid_angle *= (cos(tmp_theta))^3
670        solid_angle /= sdd^2
671       
672        // Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!)
673        w /= solid_angle
674       
675       
676        // TODO:
677        // correctly apply the correction to the error wave (assume a perfect value?)
678        // w_err /= solid_angle         //is this correct??
679
680// TODO -- clean up after I'm satisfied computations are correct               
681//      KillWaves/Z tmp_theta,tmp_dist
682       
683        return(0)
684end
685
686
687////////////
688// TODO: all of below is untested code
689//   copied from SANS
690//
691//
692// TODO :
693//   -- DoAlert 0,"This has not yet been updated for VSANS"
694//
695//performs solid angle and non-linear detector corrections to raw data as it is "added" to a work folder
696//function is called by Raw_to_work() and Add_raw_to_work() functions
697//works on the actual data array, assumes that is is already on LINEAR scale
698//
699Function V_DetCorr(data,data_err,realsread,doEfficiency,doTrans)
700        Wave data,data_err,realsread
701        Variable doEfficiency,doTrans
702
703        DoAlert 0,"This has not yet been updated for VSANS"
704       
705        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,xx0,yy0
706        Variable ii,jj,dtdist,dtdis2
707        Variable xi,xd,yd,rad,ratio,domega,xy
708        Variable lambda,trans,trans_err,lat_err,tmp_err,lat_corr
709       
710//      Print "...doing jacobian and non-linear corrections"
711
712        NVAR pixelsX = root:myGlobals:gNPixelsX
713        NVAR pixelsY = root:myGlobals:gNPixelsY
714       
715        //set up values to send to auxiliary trig functions
716        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
717        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
718
719        x0 = realsread[16]
720        y0 = realsread[17]
721        sx = realsread[10]
722        sx3 = realsread[11]
723        sy = realsread[13]
724        sy3 = realsread[14]
725       
726        dtdist = 1000*realsread[18]     //sdd in mm
727        dtdis2 = dtdist^2
728       
729        lambda = realsRead[26]
730        trans = RealsRead[4]
731        trans_err = RealsRead[41]               //new, March 2011
732       
733
734        //waves to contain repeated function calls
735        Make/O/N=(pixelsX) fyy,xx,yy            //Assumes square detector !!!
736        ii=0
737        do
738                xi = ii
739//              fyy[ii] = dc_fy(ii+1,sy,sy3,ycenter)
740//              xx[ii] = dc_fxn(ii+1,sx,sx3,xcenter)
741//              yy[ii] = dc_fym(ii+1,sy,sy3,ycenter)
742                ii+=1
743        while(ii<pixelsX)
744       
745        Make/O/N=(pixelsX,pixelsY) SolidAngle           // testing only
746       
747        ii=0
748        do
749                xi = ii
750//              xd = dc_fx(ii+1,sx,sx3,xcenter)-xx0
751                jj=0
752                do
753                        yd = fyy[jj]-yy0
754                        //rad is the distance of pixel ij from the sample
755                        //domega is the ratio of the solid angle of pixel ij versus center pixel
756                        // product xy = 1 for a detector with a linear spatial response (modern Ordela)
757                        // solid angle calculated, dW^3 >=1, so multiply data to raise measured values to correct values.
758                        rad = sqrt(dtdis2 + xd^2 + yd^2)
759                        domega = rad/dtdist
760                        ratio = domega^3
761                        xy = xx[ii]*yy[jj]
762                       
763                        data[ii][jj] *= xy*ratio
764                       
765                        solidAngle[ii][jj] = xy*ratio           //testing only 
766                        data_err[ii][jj] *= xy*ratio                    //error propagation assumes that SA and Jacobian are exact, so simply scale error
767                       
768                       
769                        // correction factor for detector efficiency JBG memo det_eff_cor2.doc 3/20/07
770                        // correction inserted 11/2007 SRK
771                        // large angle detector efficiency is >= 1 and will "bump up" the measured value at the highest angles
772                        // so divide here to get the correct answer (5/22/08 SRK)
773                        if(doEfficiency)
774//                              data[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
775//                              data_err[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
776//                              solidAngle[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)           //testing only
777                        endif
778                       
779                        // large angle transmission calculation is <= 1 and will "bump down" the measured value at the highest angles
780                        // so divide here to get the correct answer
781                        if(doTrans)
782                       
783                                if(trans<0.1 && ii==0 && jj==0)
784                                        Print "***transmission is less than 0.1*** and is a significant correction"
785                                endif
786                               
787                                if(trans==0)
788                                        if(ii==0 && jj==0)
789                                                Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
790                                        endif
791                                        trans = 1
792                                endif
793                               
794                                // pass in the transmission error, and the error in the correction is returned as the last parameter
795                                lat_corr = V_LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,lat_err)           //moved from 1D avg SRK 11/2007
796                                data[ii][jj] /= lat_corr                        //divide by the correction factor
797                                //
798                                //
799                                //
800                                // relative errors add in quadrature
801                                tmp_err = (data_err[ii][jj]/lat_corr)^2 + (lat_err/lat_corr)^2*data[ii][jj]*data[ii][jj]/lat_corr^2
802                                tmp_err = sqrt(tmp_err)
803                               
804                                data_err[ii][jj] = tmp_err
805                               
806//                              solidAngle[ii][jj] = lat_err
807
808                               
809                                //solidAngle[ii][jj] = LargeAngleTransmissionCorr(trans,dtdist,xd,yd)           //testing only
810                        endif
811                       
812                        jj+=1
813                while(jj<pixelsX)
814                ii+=1
815        while(ii<pixelsX)
816       
817        //clean up waves
818       
819        Return(0)
820End
821
822
823
824
825// DIVIDE the intensity by this correction to get the right answer
826// TODO:
827//   --         DoAlert 0,"This has not yet been updated for VSANS"
828//
829//
830Function V_LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,err)
831        Variable trans,dtdist,xd,yd,trans_err,&err
832
833        DoAlert 0,"This has not yet been updated for VSANS"
834       
835        //angle dependent transmission correction
836        Variable uval,arg,cos_th,correction,theta
837       
838        ////this section is the trans_correct() VAX routine
839//      if(trans<0.1)
840//              Print "***transmission is less than 0.1*** and is a significant correction"
841//      endif
842//      if(trans==0)
843//              Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
844//              trans = 1
845//      endif
846       
847        theta = atan( (sqrt(xd^2 + yd^2))/dtdist )              //theta at the input pixel
848       
849        //optical thickness
850        uval = -ln(trans)               //use natural logarithm
851        cos_th = cos(theta)
852        arg = (1-cos_th)/cos_th
853       
854        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
855        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120
856        // OR
857        if((uval<0.01) || (cos_th>0.99))       
858                //small arg, approx correction
859                correction= 1-0.5*uval*arg
860        else
861                //large arg, exact correction
862                correction = (1-exp(-uval*arg))/(uval*arg)
863        endif
864
865        Variable tmp
866       
867        if(trans == 1)
868                err = 0         //no correction, no error
869        else
870                //sigT, calculated from the Taylor expansion
871                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3)
872                tmp *= tmp
873                tmp *= trans_err^2
874                tmp = sqrt(tmp)         //sigT
875               
876                err = tmp
877        endif
878       
879//      Printf "trans error = %g\r",trans_err
880//      Printf "correction = %g +/- %g\r", correction, err
881       
882        //end of transmission/pathlength correction
883
884        return(correction)
885end
886
887
888//
889// TODO:
890//   --         DoAlert 0,"This has not yet been updated for VSANS"
891//
892//test procedure, not called anymore
893Proc V_AbsoluteScaling(type,c0,c1,c2,c3,c4,c5)
894        String type
895        Variable c0=1,c1=0.1,c2=0.95,c3=0.1,c4=1,c5=32.0
896        Prompt type,"WORK data type",popup,"CAL;COR;SAM"
897        Prompt c0, "Sample Transmission"
898        Prompt c1, "Sample Thickness (cm)"
899        Prompt c2, "Standard Transmission"
900        Prompt c3, "Standard Thickness (cm)"
901        Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)"
902        Prompt c5, "Standard Cross-Section (cm-1)"
903
904        Variable err
905        //call the function to do the math
906        //data from "type" will be scaled and deposited in ABS
907        err = V_Absolute_Scale(type,c0,c1,c2,c3,c4,c5)
908       
909        if(err)
910                Abort "Error in Absolute_Scale()"
911        endif
912       
913        //contents are always dumped to ABS
914        type = "ABS"
915       
916        String newTitle = "WORK_"+type
917        DoWindow/F SANS_Data
918        DoWindow/T SANS_Data, newTitle
919        KillStrings/Z newTitle
920       
921        //need to update the display with "data" from the correct dataFolder
922        //reset the current display type to "type"
923        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
924        gCurDispType = Type     
925       
926        V_fRawWindowHook()
927       
928End
929
930//
931// TODO:
932//   --         DoAlert 0,"This has not yet been updated for VSANS"
933//
934//s_ is the standard
935//w_ is the "work" file
936//both are work files and should already be normalized to 10^8 monitor counts
937Function V_Absolute_Scale(type,w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err)
938        String type
939        Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
940
941        DoAlert 0,"This has not yet been updated for VSANS"
942               
943        //convert the "type" data to absolute scale using the given standard information
944        //copying the "type" waves to ABS
945       
946        //check for existence of data, rescale to linear if needed
947        String destPath
948        //check for "type"
949        if(WaveExists($("root:Packages:NIST:"+Type + ":data")) == 0)
950                Print "There is no work file in "+type+"--Aborting"
951                Return(1)               //error condition
952        Endif
953        //check for log-scaling of the "type" data and adjust if necessary
954        destPath = "root:Packages:NIST:"+Type
955        NVAR gIsLogScale = $(destPath + ":gIsLogScale")
956        if(gIsLogScale)
957                Duplicate/O $(destPath + ":linear_data") $(destPath + ":data")//back to linear scale
958                Variable/G $(destPath + ":gIsLogScale")=0       //the "type" data is not logscale anymore
959        endif
960       
961        //copy "oldtype" information to ABS
962        //overwriting out the old contents of the ABS folder (/O option in Duplicate)
963        //copy over the waves data,vlegend,text,integers,reals(read)
964
965        String oldType= "root:Packages:NIST:"+type              //this is where the data to be absoluted is
966        //copy from current dir (type) to ABS, defined by destPath
967        Duplicate/O $(oldType + ":data"),$"root:Packages:NIST:ABS:data"
968        Duplicate/O $(oldType + ":linear_data"),$"root:Packages:NIST:ABS:linear_data"
969        Duplicate/O $(oldType + ":linear_data_error"),$"root:Packages:NIST:ABS:linear_data_error"
970//      Duplicate/O $(oldType + ":vlegend"),$"root:Packages:NIST:ABS:vlegend"
971        Duplicate/O $(oldType + ":textread"),$"root:Packages:NIST:ABS:textread"
972        Duplicate/O $(oldType + ":integersread"),$"root:Packages:NIST:ABS:integersread"
973        Duplicate/O $(oldType + ":realsread"),$"root:Packages:NIST:ABS:realsread"
974        //need to save a copy of filelist string too (from the current type folder)
975        SVAR oldFileList = $(oldType + ":fileList")
976        //need to copy filelist string too
977        String/G $"root:Packages:NIST:ABS:fileList" = oldFileList
978       
979        //now switch to ABS folder
980        //make appropriate wave references
981        WAVE data=$"root:Packages:NIST:ABS:linear_data"                                 // these wave references point to the "type" data in ABS
982        WAVE data_err=$"root:Packages:NIST:ABS:linear_data_error"                                       // these wave references point to the "type" data in ABS
983        WAVE data_copy=$"root:Packages:NIST:ABS:data"                                   // just for display
984        WAVE/T textread=$"root:Packages:NIST:ABS:textread"                      //that are to be directly operated on
985        WAVE integersread=$"root:Packages:NIST:ABS:integersread"
986        WAVE realsread=$"root:Packages:NIST:ABS:realsread"
987        Variable/G $"root:Packages:NIST:ABS:gIsLogscale"=0                      //make new flag in ABS folder, data is linear scale
988       
989        //do the actual absolute scaling here, modifying the data in ABS
990        Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
991       
992        w_moncount = realsread[0]               //monitor count in "type"
993        if(w_moncount == 0)
994                //zero monitor counts will give divide by zero ---
995                DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
996                Return(1)               //report error
997        Endif
998       
999        //calculate scale factor
1000        Variable scale,trans_err
1001        s1 = defmon/realsread[0]                //[0] is monitor count (s1 should be 1)
1002        s2 = s_thick/w_thick
1003        s3 = s_trans/w_trans
1004        s4 = s_cross/s_izero
1005       
1006        // kappa comes in as s_izero, so be sure to use 1/kappa_err
1007       
1008        data *= s1*s2*s3*s4
1009       
1010        scale = s1*s2*s3*s4
1011        trans_err = realsRead[41]
1012       
1013//      print scale
1014//      print data[0][0]
1015       
1016        data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1017
1018//      print data_err[0][0]
1019       
1020// keep "data" in sync with linear_data
1021        data_copy = data
1022       
1023        //********* 15APR02
1024        // DO NOt correct for atenuators here - the COR step already does this, putting all of the data one equal
1025        // footing (zero atten) before doing the subtraction.
1026        //
1027        //Print "ABS data multiplied by  ",s1*s2*s3*s4/attenFactor
1028       
1029        //update the ABS header information
1030        textread[1] = date() + " " + time()             //date + time stamp
1031       
1032        Return (0) //no error
1033End
1034
1035
1036//
1037// TODO:
1038//   --         DoAlert 0,"This has not yet been updated for VSANS"
1039//
1040//
1041// match the attenuation of the RAW data to the "type" data
1042// so that they can be properly added
1043//
1044// are the attenuator numbers the same? if so exit
1045//
1046// if not, find the attenuator number for type
1047// - find both attenuation factors
1048//
1049// rescale the raw data to match the ratio of the two attenuation factors
1050// -- adjust the detector count (rw)
1051// -- the linear data
1052//
1053//
1054Function V_Adjust_RAW_Attenuation(type)
1055        String type
1056
1057        DoAlert 0,"This has not yet been updated for VSANS"
1058       
1059        WAVE rw=$("root:Packages:NIST:RAW:realsread")
1060        WAVE linear_data=$("root:Packages:NIST:RAW:linear_data")
1061        WAVE data=$("root:Packages:NIST:RAW:data")
1062        WAVE data_err=$("root:Packages:NIST:RAW:linear_data_error")
1063        WAVE/T tw = $("root:Packages:NIST:RAW:textRead")
1064       
1065        WAVE dest_reals=$("root:Packages:NIST:"+type+":realsread")
1066
1067        Variable dest_atten,raw_atten,tol
1068        Variable lambda,raw_atten_err,raw_AttenFactor,dest_attenFactor,dest_atten_err
1069        String fileStr
1070
1071        dest_atten = dest_reals[3]
1072        raw_atten = rw[3]
1073       
1074        tol = 0.1               // within 0.1 atten units is OK
1075        if(abs(dest_atten - raw_atten) < tol )
1076                return(0)
1077        endif
1078
1079        fileStr = tw[3]
1080        lambda = rw[26]
1081        // TODO access correct values
1082        raw_AttenFactor = 1//AttenuationFactor(fileStr,lambda,raw_atten,raw_atten_err)
1083        dest_AttenFactor = 1//AttenuationFactor(fileStr,lambda,dest_atten,dest_atten_err)
1084               
1085        rw[2] *= dest_AttenFactor/raw_AttenFactor
1086        linear_data *= dest_AttenFactor/raw_AttenFactor
1087       
1088        // to keep "data" and linear_data in sync
1089        data = linear_data
1090       
1091        return(0)
1092End
1093
1094//
1095// TODO:
1096//   --         DoAlert 0,"This has not yet been updated for VSANS"
1097//
1098//************************
1099//unused testing procedure, may not be up-to-date with other procedures
1100//check before re-implementing
1101//
1102Proc V_DIV_a_Workfile(type)
1103        String type
1104        Prompt type,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
1105       
1106        //macro will take whatever is in SELECTED folder and DIVide it by the current
1107        //contents of the DIV folder - the function will check for existence
1108        //before proceeding
1109
1110        Abort "This has not yet been updated for VSANS"
1111       
1112        Variable err
1113        err = V_DIVCorrection(type)             //returns err = 1 if data doesn't exist in specified folders
1114       
1115        if(err)
1116                Abort "error in V_DIVCorrection()"
1117        endif
1118       
1119        //contents are NOT always dumped to CAL, but are in the new type folder
1120       
1121        String newTitle = "WORK_"+type
1122        DoWindow/F VSANS_Data
1123        DoWindow/T VSANS_Data, newTitle
1124        KillStrings/Z newTitle
1125       
1126        //need to update the display with "data" from the correct dataFolder
1127        //reset the current displaytype to "type"
1128        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=Type
1129       
1130        V_UpdateDisplayInformation(type)
1131       
1132End
1133
1134
1135//
1136// TODO:
1137//   --         DoAlert 0,"This has not yet been updated for VSANS"
1138//   -- how is the error propagation handled?
1139//
1140//function will divide the contents of "workType" folder with the contents of
1141//the DIV folder + detStr
1142// all data is linear scale for the calculation
1143//
1144Function V_DIVCorrection(data,data_err,detStr,workType)
1145        Wave data,data_err
1146        String detStr,workType
1147       
1148        //check for existence of data in type and DIV
1149        // if the desired data doesn't exist, let the user know, and abort
1150        String destPath=""
1151
1152        if(WaveExists(data) == 0)
1153                Print "The data wave does not exist in V_DIVCorrection()"
1154                Return(1)               //error condition
1155        Endif
1156       
1157        //check for DIV
1158        // if the DIV workfile doesn't exist, let the user know,and abort
1159
1160        WAVE/Z div_data = $("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")
1161        if(WaveExists(div_data) == 0)
1162                Print "The DIV wave does not exist in V_DIVCorrection()"
1163                Return(1)               //error condition
1164        Endif
1165        //files exist, proceed
1166
1167        data /= div_data
1168       
1169        data_err /= div_data
1170       
1171        Return(0)
1172End
1173
1174
1175//////////////////////////
Note: See TracBrowser for help on using the repository browser.