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

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

adding to the functionality of the data reduction protocol

getting the 1D plotting to work both as a standalone plot and part of the reduction protocol (where the choices are all made and fed in)

added log/lin toggle to 1D display. abc exponent scaling is fundamentally more messy, and has been shelved for now. Revisit if demand is there.

fixed bug in how raw data files are identified, to prevent .abs files from being incorrectly identified as raw (if same name before extension). Works correctly now, but may need tweak in future if naming scheme changes.

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