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

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

Added the angle dependent transmission correction to the data correction in the raw_to_work step, in 2D

added a testing file that can generate fake event data, read, write, and decode it. Read is based on GBLoadWave. Hoepfully I'll not need to write an XOP. manipulation of the 64 bit words are done with simple bit shifts and logic.

also added are a number of error checking routines to improve behavior when wave, folders, etc. are missing.

File size: 35.3 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
796//                              lat_corr = V_LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,lat_err)           //moved from 1D avg SRK 11/2007
797
798                                data[ii][jj] /= lat_corr                        //divide by the correction factor
799                                //
800                                //
801                                //
802                                // relative errors add in quadrature
803                                tmp_err = (data_err[ii][jj]/lat_corr)^2 + (lat_err/lat_corr)^2*data[ii][jj]*data[ii][jj]/lat_corr^2
804                                tmp_err = sqrt(tmp_err)
805                               
806                                data_err[ii][jj] = tmp_err
807                               
808//                              solidAngle[ii][jj] = lat_err
809
810                               
811                                //solidAngle[ii][jj] = LargeAngleTransmissionCorr(trans,dtdist,xd,yd)           //testing only
812                        endif
813                       
814                        jj+=1
815                while(jj<pixelsX)
816                ii+=1
817        while(ii<pixelsX)
818       
819        //clean up waves
820       
821        Return(0)
822End
823
824
825
826//
827// DIVIDE the intensity by this correction to get the right answer
828// TODO:
829//   --         DoAlert 0,"This has not yet been updated for VSANS"
830//
831//
832
833// Apply the large angle transmssion correction as the data is converted to WORK
834// so that whether the data is saved as 2D or 1D, the correction has properly been done.
835//
836// This is, however, a SAMPLE dependent calculation, not purely instrument geometry.
837//
838Function V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
839        Wave w,w_err
840        String fname,detStr,destPath
841
842        Variable sdd,xCtr,yCtr,trans,trans_err,uval
843
844// get all of the geometry information 
845//      orientation = V_getDet_tubeOrientation(fname,detStr)
846        sdd = V_getDet_ActualDistance(fname,detStr)
847        sdd/=100                // sdd in cm, pass in m
848
849        // this is ctr in mm
850        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
851        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
852        trans = V_getSampleTransmission(fname)
853        trans_err = V_getSampleTransError(fname)
854       
855        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
856       
857        Wave data_realDistX = data_realDistX
858        Wave data_realDistY = data_realDistY
859
860        Duplicate/O w lat_corr,tmp_theta,tmp_dist,lat_err,tmp_err               //in the current df
861
862//// calculate the scattering angle
863//      dx = (distX - xctr)             //delta x in mm
864//      dy = (distY - yctr)             //delta y in mm
865        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
866       
867        tmp_dist /= 10  // convert mm to cm
868        sdd *=100               //convert to cm
869
870        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle
871
872        Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp
873        numx = DimSize(tmp_theta,0)
874        numy = DimSize(tmp_theta,1)
875       
876       
877        //optical thickness
878        uval = -ln(trans)               //use natural logarithm
879       
880        for(ii=0        ;ii<numx;ii+=1)
881                for(jj=0;jj<numy;jj+=1)
882                       
883                        cos_th = cos(tmp_theta[ii][jj])
884                        arg = (1-cos_th)/cos_th
885                       
886                        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
887                        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120
888                        // OR
889                        if((uval<0.01) || (cos_th>0.99))       
890                                //small arg, approx correction
891                                lat_corr[ii][jj] = 1-0.5*uval*arg
892                        else
893                                //large arg, exact correction
894                                lat_corr[ii][jj] = (1-exp(-uval*arg))/(uval*arg)
895                        endif
896                         
897                        // TODO
898                        // -- properly calculate and apply the 2D error propagation
899                        if(trans == 1)
900                                lat_err[ii][jj] = 0             //no correction, no error
901                        else
902                                //sigT, calculated from the Taylor expansion
903                                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3)
904                                tmp *= tmp
905                                tmp *= trans_err^2
906                                tmp = sqrt(tmp)         //sigT
907                               
908                                lat_err[ii][jj] = tmp
909                        endif
910                         
911 
912                endfor
913        endfor
914       
915
916       
917        // Here it is! Apply the correction to the intensity (divide -- to get the proper correction)
918        w /= lat_corr
919
920        // relative errors add in quadrature to the current 2D error
921        tmp_err = (w_err/lat_corr)^2 + (lat_err/lat_corr)^2*w*w/lat_corr^2
922        tmp_err = sqrt(tmp_err)
923       
924        w_err = tmp_err
925       
926        // TODO:
927        // correctly apply the correction to the error wave (assume a perfect value?)
928        // w_err /= tmp         //is this correct??
929
930        // TODO -- clean up after I'm satisfied computations are correct               
931        KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err
932       
933        return(0)
934end
935
936
937//
938// TODO:
939//   --         DoAlert 0,"This has not yet been updated for VSANS"
940//
941//test procedure, not called anymore
942Proc V_AbsoluteScaling(type,c0,c1,c2,c3,c4,c5,I_err)
943        String type
944        Variable c0=1,c1=0.1,c2=0.95,c3=0.1,c4=1,c5=32.0,I_err=0.32
945        Prompt type,"WORK data type",popup,"CAL;COR;SAM"
946        Prompt c0, "Sample Transmission"
947        Prompt c1, "Sample Thickness (cm)"
948        Prompt c2, "Standard Transmission"
949        Prompt c3, "Standard Thickness (cm)"
950        Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)"
951        Prompt c5, "Standard Cross-Section (cm-1)"
952        Prompt I_err, "error in I(q=0) (one std dev)"
953
954        Variable err
955        //call the function to do the math
956        //data from "type" will be scaled and deposited in ABS
957        err = V_Absolute_Scale(type,c0,c1,c2,c3,c4,c5,I_err)
958       
959        if(err)
960                Abort "Error in V_Absolute_Scale()"
961        endif
962       
963        //contents are always dumped to ABS
964        type = "ABS"
965       
966        //need to update the display with "data" from the correct dataFolder
967        //reset the current display type to "type"
968        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
969        gCurDispType = Type     
970       
971        V_UpdateDisplayInformation(Type)
972       
973End
974
975//
976// TODO:
977//
978// kappa comes in as s_izero, so be sure to use 1/kappa_err
979//
980//convert the "type" data to absolute scale using the given standard information
981//s_ is the standard
982//w_ is the "work" file
983//both are work files and should already be normalized to 10^8 monitor counts
984Function V_Absolute_Scale(type,w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err)
985        String type
986        Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
987
988
989        Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
990        Variable scale,trans_err
991        Variable err,ii
992        String detStr
993       
994        // be sure that the starting data exists
995        err = V_WorkDataExists(type)
996        if(err==1)
997                return(err)
998        endif
999               
1000        //copy from current dir (type) to ABS
1001        V_CopyHDFToWorkFolder(type,"ABS")       
1002
1003       
1004        w_moncount = V_getMonitorCount(type)            //monitor count in "type"
1005        if(w_moncount == 0)
1006                //zero monitor counts will give divide by zero ---
1007                DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
1008                Return(1)               //report error
1009        Endif
1010       
1011        //calculate scale factor
1012        s1 = defmon/w_moncount          // monitor count (s1 should be 1)
1013        s2 = s_thick/w_thick
1014        s3 = s_trans/w_trans
1015        s4 = s_cross/s_izero
1016        scale = s1*s2*s3*s4
1017
1018        trans_err = V_getSampleTransError(type)
1019       
1020        // kappa comes in as s_izero, so be sure to use 1/kappa_err
1021
1022        // and now loop through all of the detectors
1023        //do the actual absolute scaling here, modifying the data in ABS
1024        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
1025                detStr = StringFromList(ii, ksDetectorListAll, ";")
1026                Wave data = V_getDetectorDataW("ABS",detStr)
1027                Wave data_err = V_getDetectorDataErrW("ABS",detStr)
1028               
1029                data *= s1*s2*s3*s4
1030                data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1031        endfor
1032       
1033        //********* 15APR02
1034        // DO NOT correct for atenuators here - the COR step already does this, putting all of the data one equal
1035        // footing (zero atten) before doing the subtraction.
1036       
1037        Return (0) //no error
1038End
1039
1040
1041//
1042// TODO:
1043//   --         DoAlert 0,"This has not yet been updated for VSANS"
1044//
1045//
1046// match the attenuation of the RAW data to the "type" data
1047// so that they can be properly added
1048//
1049// are the attenuator numbers the same? if so exit
1050//
1051// if not, find the attenuator number for type
1052// - find both attenuation factors
1053//
1054// rescale the raw data to match the ratio of the two attenuation factors
1055// -- adjust the detector count (rw)
1056// -- the linear data
1057//
1058//
1059Function V_Adjust_RAW_Attenuation(type)
1060        String type
1061
1062        DoAlert 0,"This has not yet been updated for VSANS"
1063       
1064        WAVE rw=$("root:Packages:NIST:RAW:realsread")
1065        WAVE linear_data=$("root:Packages:NIST:RAW:linear_data")
1066        WAVE data=$("root:Packages:NIST:RAW:data")
1067        WAVE data_err=$("root:Packages:NIST:RAW:linear_data_error")
1068        WAVE/T tw = $("root:Packages:NIST:RAW:textRead")
1069       
1070        WAVE dest_reals=$("root:Packages:NIST:"+type+":realsread")
1071
1072        Variable dest_atten,raw_atten,tol
1073        Variable lambda,raw_atten_err,raw_AttenFactor,dest_attenFactor,dest_atten_err
1074        String fileStr
1075
1076        dest_atten = dest_reals[3]
1077        raw_atten = rw[3]
1078       
1079        tol = 0.1               // within 0.1 atten units is OK
1080        if(abs(dest_atten - raw_atten) < tol )
1081                return(0)
1082        endif
1083
1084        fileStr = tw[3]
1085        lambda = rw[26]
1086        // TODO access correct values
1087        raw_AttenFactor = 1//AttenuationFactor(fileStr,lambda,raw_atten,raw_atten_err)
1088        dest_AttenFactor = 1//AttenuationFactor(fileStr,lambda,dest_atten,dest_atten_err)
1089               
1090        rw[2] *= dest_AttenFactor/raw_AttenFactor
1091        linear_data *= dest_AttenFactor/raw_AttenFactor
1092       
1093        // to keep "data" and linear_data in sync
1094        data = linear_data
1095       
1096        return(0)
1097End
1098
1099//
1100// TODO:
1101//   --         DoAlert 0,"This has not yet been updated for VSANS"
1102//
1103//************************
1104//unused testing procedure, may not be up-to-date with other procedures
1105//check before re-implementing
1106//
1107Proc V_DIV_a_Workfile(type)
1108        String type
1109        Prompt type,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
1110       
1111        //macro will take whatever is in SELECTED folder and DIVide it by the current
1112        //contents of the DIV folder - the function will check for existence
1113        //before proceeding
1114
1115        Abort "This has not yet been updated for VSANS"
1116       
1117        Variable err
1118        err = V_DIVCorrection(type)             //returns err = 1 if data doesn't exist in specified folders
1119       
1120        if(err)
1121                Abort "error in V_DIVCorrection()"
1122        endif
1123       
1124        //contents are NOT always dumped to CAL, but are in the new type folder
1125       
1126        String newTitle = "WORK_"+type
1127        DoWindow/F VSANS_Data
1128        DoWindow/T VSANS_Data, newTitle
1129        KillStrings/Z newTitle
1130       
1131        //need to update the display with "data" from the correct dataFolder
1132        //reset the current displaytype to "type"
1133        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=Type
1134       
1135        V_UpdateDisplayInformation(type)
1136       
1137End
1138
1139
1140//
1141// TODO:
1142//   --         DoAlert 0,"This has not yet been updated for VSANS"
1143//   -- how is the error propagation handled?
1144//
1145//function will divide the contents of "workType" folder with the contents of
1146//the DIV folder + detStr
1147// all data is linear scale for the calculation
1148//
1149Function V_DIVCorrection(data,data_err,detStr,workType)
1150        Wave data,data_err
1151        String detStr,workType
1152       
1153        //check for existence of data in type and DIV
1154        // if the desired data doesn't exist, let the user know, and abort
1155        String destPath=""
1156
1157        if(WaveExists(data) == 0)
1158                Print "The data wave does not exist in V_DIVCorrection()"
1159                Return(1)               //error condition
1160        Endif
1161       
1162        //check for DIV
1163        // if the DIV workfile doesn't exist, let the user know,and abort
1164
1165        WAVE/Z div_data = $("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")
1166        if(WaveExists(div_data) == 0)
1167                Print "The DIV wave does not exist in V_DIVCorrection()"
1168                Return(1)               //error condition
1169        Endif
1170        //files exist, proceed
1171
1172        data /= div_data
1173       
1174        data_err /= div_data
1175       
1176        Return(0)
1177End
1178
1179
1180//////////////////////////
Note: See TracBrowser for help on using the repository browser.