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

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

Significant changes to the base READ of individual data fields from data files. Now, if the field requested is from a WORK file, and it does not exist, an error condition is returned (or a null wave). Calling procedures are responsible for handling errors. This prevents a string of open file dialogs if fields are missing from a file if they were never in the file to begin with (like sensor logs, polarization hardware, etc.)

New get/write calls were added for the updated temperature sensor fields.

group_ID is now only in the sample block, not the duplicated in the reduction block, and is correctly a number not a string.

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// x- 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// x- be sure I'm working in the right data folder (all waves are passed in)
34// x- clean up when done
35// x- calculate + return the error contribution?
36// x- 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// x- verify the direction of the tubes and indexing
114// x- be sure I'm working in the right data folder (it is passed in, and the full path is used)
115// x- clean up when done
116// x- calculate + return the error contribution? (there is none for this operation)
117// x- do I want this to return a wave? (no, default names are generated)
118// x- do I need to write a separate function that returns the distance wave for later calculations?
119// x- do I want to make the distance array 3D to keep the x and y dims together? Calculate them all right now?
120// x- what else do I need to pass to the function? (fname=folder? detStr?)
121// y- (yes,see below) 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//
272// -- not much is known about the "B" detector, so this
273//    all hinges on the non-linear corrections being done correctly for that detector
274//
275//
276Function V_ConvertBeamCtr_to_mmB(folder,detStr,destPath)
277        String folder,detStr,destPath
278       
279        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
280        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")     
281       
282        Variable xCtr,yCtr
283        xCtr = V_getDet_beam_center_x(folder,detStr)
284        yCtr = V_getDet_beam_center_y(folder,detStr)   
285       
286        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
287        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
288        WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
289        WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
290
291        x_mm[0] = data_realDistX[xCtr][0]
292        y_mm[0] = data_realDistY[0][yCtr]
293               
294        return(0)
295end
296
297
298
299
300
301
302/////
303//
304// non-linear corrections to the tube pixels
305// - returns the distance in mm (although this may change)
306//
307// c0,c1,c2,pix
308// c0-c2 are the fit coefficients
309// pix is the test pixel
310//
311// returns the distance in mm (relative to ctr pixel)
312// ctr is the center pixel, as defined when fitting to quadratic was done
313//
314Function V_TubePixel_to_mm(c0,c1,c2,pix)
315        Variable c0,c1,c2,pix
316       
317        Variable dist
318        dist = c0 + c1*pix + c2*pix*pix
319       
320        return(dist)
321End
322//
323////
324
325// TODO
326// get rid of this in the real data
327//
328// TESTING ONLY
329Proc V_MakeFakeCalibrationWaves()
330        // make these in the RAW data folder, before converting to a work folder
331        // - then they will be "found" by get()
332        // -- only for the tube, not the Back det
333       
334//      DoAlert 0, "re-do this and do a better job of filling the fake calibration data"
335
336        DoAlert 0, "Calibration waves are read in from the data file"
337       
338//      V_fMakeFakeCalibrationWaves()
339End
340
341
342
343// TODO
344// get rid of this in the real data
345//
346// TESTING ONLY
347//
348// orientation does not matter, there are 48 tubes in each bank
349// so dimension (3,48) for everything.
350//
351// -- but the orientation does indicate TB vs LR, which has implications for
352//  the (fictional) dimension of the pixel along the tube axis, at least as far
353// as for making the fake coefficients.
354//
355Function V_fMakeFakeCalibrationWaves()
356
357        Variable ii,pixSize
358        String detStr,fname="RAW",orientation
359       
360        for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
361                detStr = StringFromList(ii, ksDetectorListNoB, ";")
362//              Wave w = V_getDetectorDataW(fname,detStr)
363                Make/O/D/N=(3,48) $("root:Packages:NIST:VSANS:RAW:entry:instrument:detector_"+detStr+":spatial_calibration")
364                Wave calib = $("root:Packages:NIST:VSANS:RAW:entry:instrument:detector_"+detStr+":spatial_calibration")
365                // !!!! this overwrites what is there
366
367                orientation = V_getDet_tubeOrientation(fname,detStr)
368                if(cmpstr(orientation,"vertical")==0)
369                //      this is vertical tube data dimensioned as (Ntubes,Npix)
370                        pixSize = 8.4           //V_getDet_y_pixel_size(fname,detStr)
371                       
372                elseif(cmpstr(orientation,"horizontal")==0)
373                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
374                        pixSize = 4                     //V_getDet_x_pixel_size(fname,detStr)
375                       
376                else           
377                        DoAlert 0,"Orientation not correctly passed in NonLinearCorrection(). No correction done."
378                endif
379               
380                calib[0][] = -(128/2)*pixSize                   //approx (n/2)*pixSixe
381                calib[1][] = pixSize
382                calib[2][] = 2e-4
383               
384        endfor
385       
386        // now fake calibration for "B"
387        Wave cal_x = V_getDet_cal_x("RAW","B")
388        Wave cal_y = V_getDet_cal_y("RAW","B")
389       
390        cal_x = 1               // mm, ignore the other 2 values
391        cal_y = 1               // mm
392        return(0)
393End
394
395//
396// TODO:
397// -- MUST VERIFY the definition of SDD and how (if) setback is written to the data files
398// -- currently I'm assuming that the SDD is the "nominal" value which is correct for the
399//    L/R panels, but is not correct for the T/B panels (must add in the setback)
400//
401//
402//
403// data_realDistX, Y must be previously generated from running NonLinearCorrection()
404//
405// call with:
406// fname as the WORK folder, "RAW"
407// detStr = detector String, "FL"
408// destPath = path to destination WORK folder ("root:Packages:NIST:VSANS:"+folder)
409//
410Function V_Detector_CalcQVals(fname,detStr,destPath)
411        String fname,detStr,destPath
412
413        String orientation
414        Variable xCtr,yCtr,lambda,sdd
415       
416// get all of the geometry information 
417        orientation = V_getDet_tubeOrientation(fname,detStr)
418//      sdd = V_getDet_distance(fname,detStr)           //[cm]
419//      sdd += V_getDet_TBSetback(fname,detStr)/10              // written in [mm], convert to [cm], returns 0 for L/R/B panels
420
421        sdd = V_getDet_ActualDistance(fname,detStr)             //sdd derived, including setback [cm]
422        sdd/=100                // sdd reported in cm, pass in m
423        // this is the ctr in pixels
424//      xCtr = V_getDet_beam_center_x(fname,detStr)
425//      yCtr = V_getDet_beam_center_y(fname,detStr)
426        // this is ctr in mm
427        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
428        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
429        lambda = V_getWavelength(fname)
430        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
431        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
432
433// make the new waves
434        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
435        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
436        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
437        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
438        Wave qTot = $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
439        Wave qx = $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
440        Wave qy = $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
441        Wave qz = $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
442
443// calculate all of the q-values
444        qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
445        qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
446        qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
447        qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
448       
449       
450        return(0)
451End
452
453
454//function to calculate the overall q-value, given all of the necesary trig inputs
455//
456// TODO:
457// -- verify the calculation (accuracy - in all input conditions)
458// -- verify the units of everything here, it's currently all jumbled and wrong... and repeated
459// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
460//    to each pixel
461//
462//sdd is in meters
463//wavelength is in Angstroms
464//
465//returned magnitude of Q is in 1/Angstroms
466//
467Function V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
468        Variable xaxval,yaxval,xctr,yctr,sdd,lam
469        Wave distX,distY
470       
471        Variable dx,dy,qval,two_theta,dist
472               
473        sdd *=100               //convert to cm
474        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
475        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
476        dist = sqrt(dx^2 + dy^2)
477       
478        dist /= 10  // convert mm to cm
479       
480        two_theta = atan(dist/sdd)
481
482        qval = 4*Pi/lam*sin(two_theta/2)
483       
484        return qval
485End
486
487//calculates just the q-value in the x-direction on the detector
488// TODO:
489// -- verify the calculation (accuracy - in all input conditions)
490// -- verify the units of everything here, it's currently all jumbled and wrong... and repeated
491// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
492//    to each pixel
493//
494//
495// this properly accounts for qz
496//
497Function V_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
498        Variable xaxval,yaxval,xctr,yctr,sdd,lam
499        Wave distX,distY
500
501        Variable qx,qval,phi,dx,dy,dist,two_theta
502       
503        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
504       
505        sdd *=100               //convert to cm
506        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
507        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
508        phi = V_FindPhi(dx,dy)
509       
510        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
511        dist = sqrt(dx^2 + dy^2)
512        dist /= 10  // convert mm to cm
513
514        two_theta = atan(dist/sdd)
515
516        qx = qval*cos(two_theta/2)*cos(phi)
517       
518        return qx
519End
520
521//calculates just the q-value in the y-direction on the detector
522// TODO:
523// -- verify the calculation (accuracy - in all input conditions)
524// -- verify the units of everything here, it's currently all jumbled and wrong... and repeated
525// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
526//    to each pixel
527//
528//
529// this properly accounts for qz
530//
531Function V_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
532        Variable xaxval,yaxval,xctr,yctr,sdd,lam
533        Wave distX,distY
534
535        Variable qy,qval,phi,dx,dy,dist,two_theta
536       
537        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
538       
539        sdd *=100               //convert to cm
540        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
541        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
542        phi = V_FindPhi(dx,dy)
543       
544        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
545        dist = sqrt(dx^2 + dy^2)
546        dist /= 10  // convert mm to cm
547
548        two_theta = atan(dist/sdd)
549
550        qy = qval*cos(two_theta/2)*sin(phi)
551       
552        return qy
553End
554
555//calculates just the q-value in the z-direction on the detector
556// TODO:
557// -- verify the calculation (accuracy - in all input conditions)
558// -- verify the units of everything here, it's currently all jumbled and wrong... and repeated
559// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
560//    to each pixel
561//
562// not actually used for any calculations, but here for completeness if anyone asks, or for 2D data export
563//
564// this properly accounts for qz
565//
566Function V_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
567        Variable xaxval,yaxval,xctr,yctr,sdd,lam
568        Wave distX,distY
569
570        Variable qz,qval,phi,dx,dy,dist,two_theta
571       
572        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
573       
574        sdd *=100               //convert to cm
575        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
576        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
577       
578        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
579        dist = sqrt(dx^2 + dy^2)
580        dist /= 10  // convert mm to cm
581
582        two_theta = atan(dist/sdd)
583
584        qz = qval*sin(two_theta/2)
585       
586        return qz
587End
588
589
590//
591// TODO -- VERIFY calculations
592// -- This is the actual solid angle per pixel, not a ratio vs. some "unit SA"
593//    Do I just correct for the different area vs. the "nominal" central area?
594// -- decide how to implement - either directly change the data values (as was done in the past)
595//    or use this as a weighting for when the data is binned to I(q). In the second method, 2D data
596//    would need this to be applied before exporting
597// -- do I keep a wave note indicating that this correction has been applied to the data
598//    so that it can be "un-applied"?
599// -- do I calculate theta from geometry directly, or get it from Q (Assuming it's present?)
600//    (probably just from geometry, since I need SDD and dx and dy values...)
601//
602//
603Function V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
604        Wave w,w_err
605        String fname,detStr,destPath
606
607        Variable sdd,xCtr,yCtr,lambda
608
609// get all of the geometry information 
610//      orientation = V_getDet_tubeOrientation(fname,detStr)
611        sdd = V_getDet_ActualDistance(fname,detStr)
612        sdd/=100                // sdd in cm, pass in m
613
614        // this is ctr in mm
615        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
616        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
617        lambda = V_getWavelength(fname)
618       
619        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
620       
621        Wave data_realDistX = data_realDistX
622        Wave data_realDistY = data_realDistY
623
624        Duplicate/O w solid_angle,tmp_theta,tmp_dist            //in the current df
625
626//// calculate the scattering angle
627//      dx = (distX - xctr)             //delta x in mm
628//      dy = (distY - yctr)             //delta y in mm
629        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
630       
631        tmp_dist /= 10  // convert mm to cm
632        sdd *=100               //convert to cm
633
634        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle
635
636        Variable ii,jj,numx,numy,dx,dy
637        numx = DimSize(tmp_theta,0)
638        numy = DimSize(tmp_theta,1)
639       
640        for(ii=0        ;ii<numx;ii+=1)
641                for(jj=0;jj<numy;jj+=1)
642                       
643                        if(ii==0)               //do a forward difference if ii==0
644                                dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj])        //delta x for the pixel
645                        else
646                                dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj])        //delta x for the pixel
647                        endif
648                       
649                       
650                        if(jj==0)
651                                dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj])        //delta y for the pixel
652                        else
653                                dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1])        //delta y for the pixel
654                        endif
655       
656                        dx /= 10
657                        dy /= 10                // convert mm to cm (since sdd is in cm)
658                        solid_angle[ii][jj] = dx*dy             //this is in cm^2
659                endfor
660        endfor
661       
662        // to cover up any issues w/negative dx or dy
663        solid_angle = abs(solid_angle)
664       
665        // solid_angle correction
666        // == dx*dy*cos^3/sdd^2
667        solid_angle *= (cos(tmp_theta))^3
668        solid_angle /= sdd^2
669       
670        // Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!)
671        w /= solid_angle
672       
673       
674        // TODO:
675        // correctly apply the correction to the error wave (assume a perfect value?)
676        w_err /= solid_angle            //is this correct??
677
678// TODO -- clean up after I'm satisfied computations are correct               
679//      KillWaves/Z tmp_theta,tmp_dist
680       
681        return(0)
682end
683
684
685////////////
686// TODO: all of below is untested code
687//   copied from SANS
688//
689//
690// TODO :
691//   -- DoAlert 0,"This has not yet been updated for VSANS"
692//
693//performs solid angle and non-linear detector corrections to raw data as it is "added" to a work folder
694//function is called by Raw_to_work() and Add_raw_to_work() functions
695//works on the actual data array, assumes that is is already on LINEAR scale
696//
697Function V_DetCorr(data,data_err,realsread,doEfficiency,doTrans)
698        Wave data,data_err,realsread
699        Variable doEfficiency,doTrans
700
701        DoAlert 0,"This has not yet been updated for VSANS"
702       
703        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,xx0,yy0
704        Variable ii,jj,dtdist,dtdis2
705        Variable xi,xd,yd,rad,ratio,domega,xy
706        Variable lambda,trans,trans_err,lat_err,tmp_err,lat_corr
707       
708//      Print "...doing jacobian and non-linear corrections"
709
710        NVAR pixelsX = root:myGlobals:gNPixelsX
711        NVAR pixelsY = root:myGlobals:gNPixelsY
712       
713        //set up values to send to auxiliary trig functions
714        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
715        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
716
717        x0 = realsread[16]
718        y0 = realsread[17]
719        sx = realsread[10]
720        sx3 = realsread[11]
721        sy = realsread[13]
722        sy3 = realsread[14]
723       
724        dtdist = 1000*realsread[18]     //sdd in mm
725        dtdis2 = dtdist^2
726       
727        lambda = realsRead[26]
728        trans = RealsRead[4]
729        trans_err = RealsRead[41]               //new, March 2011
730       
731
732        //waves to contain repeated function calls
733        Make/O/N=(pixelsX) fyy,xx,yy            //Assumes square detector !!!
734        ii=0
735        do
736                xi = ii
737//              fyy[ii] = dc_fy(ii+1,sy,sy3,ycenter)
738//              xx[ii] = dc_fxn(ii+1,sx,sx3,xcenter)
739//              yy[ii] = dc_fym(ii+1,sy,sy3,ycenter)
740                ii+=1
741        while(ii<pixelsX)
742       
743        Make/O/N=(pixelsX,pixelsY) SolidAngle           // testing only
744       
745        ii=0
746        do
747                xi = ii
748//              xd = dc_fx(ii+1,sx,sx3,xcenter)-xx0
749                jj=0
750                do
751                        yd = fyy[jj]-yy0
752                        //rad is the distance of pixel ij from the sample
753                        //domega is the ratio of the solid angle of pixel ij versus center pixel
754                        // product xy = 1 for a detector with a linear spatial response (modern Ordela)
755                        // solid angle calculated, dW^3 >=1, so multiply data to raise measured values to correct values.
756                        rad = sqrt(dtdis2 + xd^2 + yd^2)
757                        domega = rad/dtdist
758                        ratio = domega^3
759                        xy = xx[ii]*yy[jj]
760                       
761                        data[ii][jj] *= xy*ratio
762                       
763                        solidAngle[ii][jj] = xy*ratio           //testing only 
764                        data_err[ii][jj] *= xy*ratio                    //error propagation assumes that SA and Jacobian are exact, so simply scale error
765                       
766                       
767                        // correction factor for detector efficiency JBG memo det_eff_cor2.doc 3/20/07
768                        // correction inserted 11/2007 SRK
769                        // large angle detector efficiency is >= 1 and will "bump up" the measured value at the highest angles
770                        // so divide here to get the correct answer (5/22/08 SRK)
771                        if(doEfficiency)
772//                              data[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
773//                              data_err[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
774//                              solidAngle[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)           //testing only
775                        endif
776                       
777                        // large angle transmission calculation is <= 1 and will "bump down" the measured value at the highest angles
778                        // so divide here to get the correct answer
779                        if(doTrans)
780                       
781                                if(trans<0.1 && ii==0 && jj==0)
782                                        Print "***transmission is less than 0.1*** and is a significant correction"
783                                endif
784                               
785                                if(trans==0)
786                                        if(ii==0 && jj==0)
787                                                Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
788                                        endif
789                                        trans = 1
790                                endif
791                               
792                                // pass in the transmission error, and the error in the correction is returned as the last parameter
793
794//                              lat_corr = V_LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,lat_err)           //moved from 1D avg SRK 11/2007
795
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//
830
831// Apply the large angle transmssion correction as the data is converted to WORK
832// so that whether the data is saved as 2D or 1D, the correction has properly been done.
833//
834// This is, however, a SAMPLE dependent calculation, not purely instrument geometry.
835//
836Function V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
837        Wave w,w_err
838        String fname,detStr,destPath
839
840        Variable sdd,xCtr,yCtr,trans,trans_err,uval
841
842// get all of the geometry information 
843//      orientation = V_getDet_tubeOrientation(fname,detStr)
844        sdd = V_getDet_ActualDistance(fname,detStr)
845        sdd/=100                // sdd in cm, pass in m
846
847        // this is ctr in mm
848        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
849        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
850        trans = V_getSampleTransmission(fname)
851        trans_err = V_getSampleTransError(fname)
852       
853        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
854       
855        Wave data_realDistX = data_realDistX
856        Wave data_realDistY = data_realDistY
857
858        Duplicate/O w lat_corr,tmp_theta,tmp_dist,lat_err,tmp_err               //in the current df
859
860//// calculate the scattering angle
861//      dx = (distX - xctr)             //delta x in mm
862//      dy = (distY - yctr)             //delta y in mm
863        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
864       
865        tmp_dist /= 10  // convert mm to cm
866        sdd *=100               //convert to cm
867
868        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle
869
870        Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp
871        numx = DimSize(tmp_theta,0)
872        numy = DimSize(tmp_theta,1)
873       
874       
875        //optical thickness
876        uval = -ln(trans)               //use natural logarithm
877       
878        for(ii=0        ;ii<numx;ii+=1)
879                for(jj=0;jj<numy;jj+=1)
880                       
881                        cos_th = cos(tmp_theta[ii][jj])
882                        arg = (1-cos_th)/cos_th
883                       
884                        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
885                        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120
886                        // OR
887                        if((uval<0.01) || (cos_th>0.99))       
888                                //small arg, approx correction
889                                lat_corr[ii][jj] = 1-0.5*uval*arg
890                        else
891                                //large arg, exact correction
892                                lat_corr[ii][jj] = (1-exp(-uval*arg))/(uval*arg)
893                        endif
894                         
895                        // TODO
896                        // -- properly calculate and apply the 2D error propagation
897                        if(trans == 1)
898                                lat_err[ii][jj] = 0             //no correction, no error
899                        else
900                                //sigT, calculated from the Taylor expansion
901                                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3)
902                                tmp *= tmp
903                                tmp *= trans_err^2
904                                tmp = sqrt(tmp)         //sigT
905                               
906                                lat_err[ii][jj] = tmp
907                        endif
908                         
909 
910                endfor
911        endfor
912       
913
914       
915        // Here it is! Apply the correction to the intensity (divide -- to get the proper correction)
916        w /= lat_corr
917
918        // relative errors add in quadrature to the current 2D error
919        tmp_err = (w_err/lat_corr)^2 + (lat_err/lat_corr)^2*w*w/lat_corr^2
920        tmp_err = sqrt(tmp_err)
921       
922        w_err = tmp_err
923       
924        // TODO:
925        // correctly apply the correction to the error wave (assume a perfect value?)
926        // w_err /= tmp         //is this correct??
927
928        // TODO -- clean up after I'm satisfied computations are correct               
929        KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err
930       
931        return(0)
932end
933
934
935//
936// TODO:
937//   --         DoAlert 0,"This has not yet been updated for VSANS"
938//
939//test procedure, not called anymore
940Proc V_AbsoluteScaling(type,c0,c1,c2,c3,c4,c5,I_err)
941        String type
942        Variable c0=1,c1=0.1,c2=0.95,c3=0.1,c4=1,c5=32.0,I_err=0.32
943        Prompt type,"WORK data type",popup,"CAL;COR;SAM"
944        Prompt c0, "Sample Transmission"
945        Prompt c1, "Sample Thickness (cm)"
946        Prompt c2, "Standard Transmission"
947        Prompt c3, "Standard Thickness (cm)"
948        Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)"
949        Prompt c5, "Standard Cross-Section (cm-1)"
950        Prompt I_err, "error in I(q=0) (one std dev)"
951
952        Variable err
953        //call the function to do the math
954        //data from "type" will be scaled and deposited in ABS
955        err = V_Absolute_Scale(type,c0,c1,c2,c3,c4,c5,I_err)
956       
957        if(err)
958                Abort "Error in V_Absolute_Scale()"
959        endif
960       
961        //contents are always dumped to ABS
962        type = "ABS"
963       
964        //need to update the display with "data" from the correct dataFolder
965        //reset the current display type to "type"
966        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
967        gCurDispType = Type     
968       
969        V_UpdateDisplayInformation(Type)
970       
971End
972
973//
974// TODO:
975//
976// kappa comes in as s_izero, so be sure to use 1/kappa_err
977//
978//convert the "type" data to absolute scale using the given standard information
979//s_ is the standard
980//w_ is the "work" file
981//both are work files and should already be normalized to 10^8 monitor counts
982Function V_Absolute_Scale(type,w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err)
983        String type
984        Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
985
986
987        Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
988        Variable scale,trans_err
989        Variable err,ii
990        String detStr
991       
992        // be sure that the starting data exists
993        err = V_WorkDataExists(type)
994        if(err==1)
995                return(err)
996        endif
997               
998        //copy from current dir (type) to ABS
999        V_CopyHDFToWorkFolder(type,"ABS")       
1000
1001       
1002        w_moncount = V_getMonitorCount(type)            //monitor count in "type"
1003        if(w_moncount == 0)
1004                //zero monitor counts will give divide by zero ---
1005                DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
1006                Return(1)               //report error
1007        Endif
1008       
1009        //calculate scale factor
1010        s1 = defmon/w_moncount          // monitor count (s1 should be 1)
1011        s2 = s_thick/w_thick
1012        s3 = s_trans/w_trans
1013        s4 = s_cross/s_izero
1014        scale = s1*s2*s3*s4
1015
1016        trans_err = V_getSampleTransError(type)
1017       
1018        // kappa comes in as s_izero, so be sure to use 1/kappa_err
1019
1020        // and now loop through all of the detectors
1021        //do the actual absolute scaling here, modifying the data in ABS
1022        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
1023                detStr = StringFromList(ii, ksDetectorListAll, ";")
1024                Wave data = V_getDetectorDataW("ABS",detStr)
1025                Wave data_err = V_getDetectorDataErrW("ABS",detStr)
1026               
1027                data *= s1*s2*s3*s4
1028                data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1029        endfor
1030       
1031        //********* 15APR02
1032        // DO NOT correct for atenuators here - the COR step already does this, putting all of the data one equal
1033        // footing (zero atten) before doing the subtraction.
1034       
1035        Return (0) //no error
1036End
1037
1038
1039//
1040// TODO:
1041//   --         DoAlert 0,"This has not yet been updated for VSANS"
1042//
1043//
1044// match the attenuation of the RAW data to the "type" data
1045// so that they can be properly added
1046//
1047// are the attenuator numbers the same? if so exit
1048//
1049// if not, find the attenuator number for type
1050// - find both attenuation factors
1051//
1052// rescale the raw data to match the ratio of the two attenuation factors
1053// -- adjust the detector count (rw)
1054// -- the linear data
1055//
1056//
1057Function V_Adjust_RAW_Attenuation(type)
1058        String type
1059
1060        DoAlert 0,"This has not yet been updated for VSANS"
1061       
1062        WAVE rw=$("root:Packages:NIST:RAW:realsread")
1063        WAVE linear_data=$("root:Packages:NIST:RAW:linear_data")
1064        WAVE data=$("root:Packages:NIST:RAW:data")
1065        WAVE data_err=$("root:Packages:NIST:RAW:linear_data_error")
1066        WAVE/T tw = $("root:Packages:NIST:RAW:textRead")
1067       
1068        WAVE dest_reals=$("root:Packages:NIST:"+type+":realsread")
1069
1070        Variable dest_atten,raw_atten,tol
1071        Variable lambda,raw_atten_err,raw_AttenFactor,dest_attenFactor,dest_atten_err
1072        String fileStr
1073
1074        dest_atten = dest_reals[3]
1075        raw_atten = rw[3]
1076       
1077        tol = 0.1               // within 0.1 atten units is OK
1078        if(abs(dest_atten - raw_atten) < tol )
1079                return(0)
1080        endif
1081
1082        fileStr = tw[3]
1083        lambda = rw[26]
1084        // TODO access correct values
1085        raw_AttenFactor = 1//AttenuationFactor(fileStr,lambda,raw_atten,raw_atten_err)
1086        dest_AttenFactor = 1//AttenuationFactor(fileStr,lambda,dest_atten,dest_atten_err)
1087               
1088        rw[2] *= dest_AttenFactor/raw_AttenFactor
1089        linear_data *= dest_AttenFactor/raw_AttenFactor
1090       
1091        // to keep "data" and linear_data in sync
1092        data = linear_data
1093       
1094        return(0)
1095End
1096
1097//
1098// TODO:
1099//   --         DoAlert 0,"This has not yet been updated for VSANS"
1100//
1101//************************
1102//unused testing procedure, may not be up-to-date with other procedures
1103//check before re-implementing
1104//
1105Proc V_DIV_a_Workfile(type)
1106        String type
1107        Prompt type,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
1108       
1109        //macro will take whatever is in SELECTED folder and DIVide it by the current
1110        //contents of the DIV folder - the function will check for existence
1111        //before proceeding
1112
1113        Abort "This has not yet been updated for VSANS"
1114       
1115        Variable err
1116        err = V_DIVCorrection(type)             //returns err = 1 if data doesn't exist in specified folders
1117       
1118        if(err)
1119                Abort "error in V_DIVCorrection()"
1120        endif
1121       
1122        //contents are NOT always dumped to CAL, but are in the new type folder
1123       
1124        String newTitle = "WORK_"+type
1125        DoWindow/F VSANS_Data
1126        DoWindow/T VSANS_Data, newTitle
1127        KillStrings/Z newTitle
1128       
1129        //need to update the display with "data" from the correct dataFolder
1130        //reset the current displaytype to "type"
1131        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=Type
1132       
1133        V_UpdateDisplayInformation(type)
1134       
1135End
1136
1137
1138//
1139// TODO:
1140//   --         DoAlert 0,"This has not yet been updated for VSANS"
1141//   -- how is the error propagation handled?
1142//
1143//function will divide the contents of "workType" folder with the contents of
1144//the DIV folder + detStr
1145// all data is linear scale for the calculation
1146//
1147Function V_DIVCorrection(data,data_err,detStr,workType)
1148        Wave data,data_err
1149        String detStr,workType
1150       
1151        //check for existence of data in type and DIV
1152        // if the desired data doesn't exist, let the user know, and abort
1153        String destPath=""
1154
1155        if(WaveExists(data) == 0)
1156                Print "The data wave does not exist in V_DIVCorrection()"
1157                Return(1)               //error condition
1158        Endif
1159       
1160        //check for DIV
1161        // if the DIV workfile doesn't exist, let the user know,and abort
1162
1163        WAVE/Z div_data = $("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")
1164        if(WaveExists(div_data) == 0)
1165                Print "The DIV wave does not exist in V_DIVCorrection()"
1166                Return(1)               //error condition
1167        Endif
1168        //files exist, proceed
1169
1170        data /= div_data
1171       
1172        data_err /= div_data
1173       
1174        Return(0)
1175End
1176
1177
1178//////////////////////////
Note: See TracBrowser for help on using the repository browser.