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

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

more changes, bug fixes, detector dead time fix

File size: 33.1 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 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 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        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 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 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 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 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 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//      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 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                     //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// data_realDistX, Y must be previously generated from running NonLinearCorrection()
399//
400// call with:
401// fname as the WORK folder, "RAW"
402// detStr = detector String, "FL"
403// destPath = path to destination WORK folder ("root:Packages:NIST:VSANS:"+folder)
404//
405Function V_Detector_CalcQVals(fname,detStr,destPath)
406        String fname,detStr,destPath
407
408        String orientation
409        Variable xCtr,yCtr,lambda,sdd
410       
411// get all of the geometry information 
412        orientation = V_getDet_tubeOrientation(fname,detStr)
413        sdd = V_getDet_distance(fname,detStr)
414        sdd/=100                // sdd reported in cm, pass in m
415        // this is the ctr in pixels
416//      xCtr = V_getDet_beam_center_x(fname,detStr)
417//      yCtr = V_getDet_beam_center_y(fname,detStr)
418        // this is ctr in mm
419        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
420        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
421        lambda = V_getWavelength(fname)
422        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
423        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
424
425// make the new waves
426        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
427        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
428        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
429        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
430        Wave qTot = $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
431        Wave qx = $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
432        Wave qy = $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
433        Wave qz = $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
434
435// calculate all of the q-values
436        qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
437        qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
438        qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
439        qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
440       
441       
442        return(0)
443End
444
445
446//function to calculate the overall q-value, given all of the necesary trig inputs
447//
448// TODO:
449// -- verify the calculation (accuracy - in all input conditions)
450// -- verify the units of everything here, it's currently all jumbled and wrong... and repeated
451// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
452//    to each pixel
453//
454//sdd is in meters
455//wavelength is in Angstroms
456//
457//returned magnitude of Q is in 1/Angstroms
458//
459Function V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
460        Variable xaxval,yaxval,xctr,yctr,sdd,lam
461        Wave distX,distY
462       
463        Variable dx,dy,qval,two_theta,dist
464               
465        sdd *=100               //convert to cm
466        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
467        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
468        dist = sqrt(dx^2 + dy^2)
469       
470        dist /= 10  // convert mm to cm
471       
472        two_theta = atan(dist/sdd)
473
474        qval = 4*Pi/lam*sin(two_theta/2)
475       
476        return qval
477End
478
479//calculates just the q-value in the x-direction on the detector
480// TODO:
481// -- verify the calculation (accuracy - in all input conditions)
482// -- verify the units of everything here, it's currently all jumbled and wrong... and repeated
483// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
484//    to each pixel
485//
486//
487// this properly accounts for qz
488//
489Function V_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
490        Variable xaxval,yaxval,xctr,yctr,sdd,lam
491        Wave distX,distY
492
493        Variable qx,qval,phi,dx,dy,dist,two_theta
494       
495        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
496       
497        sdd *=100               //convert to cm
498        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
499        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
500        phi = V_FindPhi(dx,dy)
501       
502        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
503        dist = sqrt(dx^2 + dy^2)
504        dist /= 10  // convert mm to cm
505
506        two_theta = atan(dist/sdd)
507
508        qx = qval*cos(two_theta/2)*cos(phi)
509       
510        return qx
511End
512
513//calculates just the q-value in the y-direction on the detector
514// TODO:
515// -- verify the calculation (accuracy - in all input conditions)
516// -- verify the units of everything here, it's currently all jumbled and wrong... and repeated
517// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
518//    to each pixel
519//
520//
521// this properly accounts for qz
522//
523Function V_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
524        Variable xaxval,yaxval,xctr,yctr,sdd,lam
525        Wave distX,distY
526
527        Variable qy,qval,phi,dx,dy,dist,two_theta
528       
529        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
530       
531        sdd *=100               //convert to cm
532        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
533        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
534        phi = V_FindPhi(dx,dy)
535       
536        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
537        dist = sqrt(dx^2 + dy^2)
538        dist /= 10  // convert mm to cm
539
540        two_theta = atan(dist/sdd)
541
542        qy = qval*cos(two_theta/2)*sin(phi)
543       
544        return qy
545End
546
547//calculates just the q-value in the z-direction on the detector
548// TODO:
549// -- verify the calculation (accuracy - in all input conditions)
550// -- verify the units of everything here, it's currently all jumbled and wrong... and repeated
551// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
552//    to each pixel
553//
554// not actually used for anything, but here for completeness if anyone asks
555//
556// this properly accounts for qz
557//
558Function V_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
559        Variable xaxval,yaxval,xctr,yctr,sdd,lam
560        Wave distX,distY
561
562        Variable qz,qval,phi,dx,dy,dist,two_theta
563       
564        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
565       
566        sdd *=100               //convert to cm
567        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
568        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
569       
570        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
571        dist = sqrt(dx^2 + dy^2)
572        dist /= 10  // convert mm to cm
573
574        two_theta = atan(dist/sdd)
575
576        qz = qval*sin(two_theta/2)
577       
578        return qz
579End
580
581
582
583////////////
584// TODO: all of below is untested code
585//   copied from SANS
586//
587//
588// TODO :
589//   -- DoAlert 0,"This has not yet been updated for VSANS"
590//
591//performs solid angle and non-linear detector corrections to raw data as it is "added" to a work folder
592//function is called by Raw_to_work() and Add_raw_to_work() functions
593//works on the actual data array, assumes that is is already on LINEAR scale
594//
595Function DetCorr(data,data_err,realsread,doEfficiency,doTrans)
596        Wave data,data_err,realsread
597        Variable doEfficiency,doTrans
598
599        DoAlert 0,"This has not yet been updated for VSANS"
600       
601        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,xx0,yy0
602        Variable ii,jj,dtdist,dtdis2
603        Variable xi,xd,yd,rad,ratio,domega,xy
604        Variable lambda,trans,trans_err,lat_err,tmp_err,lat_corr
605       
606//      Print "...doing jacobian and non-linear corrections"
607
608        NVAR pixelsX = root:myGlobals:gNPixelsX
609        NVAR pixelsY = root:myGlobals:gNPixelsY
610       
611        //set up values to send to auxiliary trig functions
612        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
613        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
614
615        x0 = realsread[16]
616        y0 = realsread[17]
617        sx = realsread[10]
618        sx3 = realsread[11]
619        sy = realsread[13]
620        sy3 = realsread[14]
621       
622        dtdist = 1000*realsread[18]     //sdd in mm
623        dtdis2 = dtdist^2
624       
625        lambda = realsRead[26]
626        trans = RealsRead[4]
627        trans_err = RealsRead[41]               //new, March 2011
628       
629
630        //waves to contain repeated function calls
631        Make/O/N=(pixelsX) fyy,xx,yy            //Assumes square detector !!!
632        ii=0
633        do
634                xi = ii
635//              fyy[ii] = dc_fy(ii+1,sy,sy3,ycenter)
636//              xx[ii] = dc_fxn(ii+1,sx,sx3,xcenter)
637//              yy[ii] = dc_fym(ii+1,sy,sy3,ycenter)
638                ii+=1
639        while(ii<pixelsX)
640       
641        Make/O/N=(pixelsX,pixelsY) SolidAngle           // testing only
642       
643        ii=0
644        do
645                xi = ii
646//              xd = dc_fx(ii+1,sx,sx3,xcenter)-xx0
647                jj=0
648                do
649                        yd = fyy[jj]-yy0
650                        //rad is the distance of pixel ij from the sample
651                        //domega is the ratio of the solid angle of pixel ij versus center pixel
652                        // product xy = 1 for a detector with a linear spatial response (modern Ordela)
653                        // solid angle calculated, dW^3 >=1, so multiply data to raise measured values to correct values.
654                        rad = sqrt(dtdis2 + xd^2 + yd^2)
655                        domega = rad/dtdist
656                        ratio = domega^3
657                        xy = xx[ii]*yy[jj]
658                       
659                        data[ii][jj] *= xy*ratio
660                       
661                        solidAngle[ii][jj] = xy*ratio           //testing only 
662                        data_err[ii][jj] *= xy*ratio                    //error propagation assumes that SA and Jacobian are exact, so simply scale error
663                       
664                       
665                        // correction factor for detector efficiency JBG memo det_eff_cor2.doc 3/20/07
666                        // correction inserted 11/2007 SRK
667                        // large angle detector efficiency is >= 1 and will "bump up" the measured value at the highest angles
668                        // so divide here to get the correct answer (5/22/08 SRK)
669                        if(doEfficiency)
670                                data[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
671                                data_err[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
672//                              solidAngle[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)           //testing only
673                        endif
674                       
675                        // large angle transmission calculation is <= 1 and will "bump down" the measured value at the highest angles
676                        // so divide here to get the correct answer
677                        if(doTrans)
678                       
679                                if(trans<0.1 && ii==0 && jj==0)
680                                        Print "***transmission is less than 0.1*** and is a significant correction"
681                                endif
682                               
683                                if(trans==0)
684                                        if(ii==0 && jj==0)
685                                                Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
686                                        endif
687                                        trans = 1
688                                endif
689                               
690                                // pass in the transmission error, and the error in the correction is returned as the last parameter
691                                lat_corr = LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,lat_err)             //moved from 1D avg SRK 11/2007
692                                data[ii][jj] /= lat_corr                        //divide by the correction factor
693                                //
694                                //
695                                //
696                                // relative errors add in quadrature
697                                tmp_err = (data_err[ii][jj]/lat_corr)^2 + (lat_err/lat_corr)^2*data[ii][jj]*data[ii][jj]/lat_corr^2
698                                tmp_err = sqrt(tmp_err)
699                               
700                                data_err[ii][jj] = tmp_err
701                               
702//                              solidAngle[ii][jj] = lat_err
703
704                               
705                                //solidAngle[ii][jj] = LargeAngleTransmissionCorr(trans,dtdist,xd,yd)           //testing only
706                        endif
707                       
708                        jj+=1
709                while(jj<pixelsX)
710                ii+=1
711        while(ii<pixelsX)
712       
713        //clean up waves
714       
715        Return(0)
716End
717
718
719
720//distances passed in are in mm
721// dtdist is SDD
722// xd and yd are distances from the beam center to the current pixel
723//
724// TODO:
725//   --         DoAlert 0,"This has not yet been updated for VSANS"
726//
727Function DetEffCorr(lambda,dtdist,xd,yd)
728        Variable lambda,dtdist,xd,yd
729
730        DoAlert 0,"This has not yet been updated for VSANS"
731       
732        Variable theta,cosT,ff,stAl,stHe
733       
734        theta = atan( (sqrt(xd^2 + yd^2))/dtdist )
735        cosT = cos(theta)
736       
737        stAl = 0.00967*lambda*0.8               //dimensionless, constants from JGB memo
738        stHe = 0.146*lambda*2.5
739       
740        ff = exp(-stAl/cosT)*(1-exp(-stHe/cosT)) / ( exp(-stAl)*(1-exp(-stHe)) )
741               
742        return(ff)
743End
744
745// DIVIDE the intensity by this correction to get the right answer
746// TODO:
747//   --         DoAlert 0,"This has not yet been updated for VSANS"
748//
749//
750Function LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,err)
751        Variable trans,dtdist,xd,yd,trans_err,&err
752
753        DoAlert 0,"This has not yet been updated for VSANS"
754       
755        //angle dependent transmission correction
756        Variable uval,arg,cos_th,correction,theta
757       
758        ////this section is the trans_correct() VAX routine
759//      if(trans<0.1)
760//              Print "***transmission is less than 0.1*** and is a significant correction"
761//      endif
762//      if(trans==0)
763//              Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
764//              trans = 1
765//      endif
766       
767        theta = atan( (sqrt(xd^2 + yd^2))/dtdist )              //theta at the input pixel
768       
769        //optical thickness
770        uval = -ln(trans)               //use natural logarithm
771        cos_th = cos(theta)
772        arg = (1-cos_th)/cos_th
773       
774        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
775        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120
776        // OR
777        if((uval<0.01) || (cos_th>0.99))       
778                //small arg, approx correction
779                correction= 1-0.5*uval*arg
780        else
781                //large arg, exact correction
782                correction = (1-exp(-uval*arg))/(uval*arg)
783        endif
784
785        Variable tmp
786       
787        if(trans == 1)
788                err = 0         //no correction, no error
789        else
790                //sigT, calculated from the Taylor expansion
791                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3)
792                tmp *= tmp
793                tmp *= trans_err^2
794                tmp = sqrt(tmp)         //sigT
795               
796                err = tmp
797        endif
798       
799//      Printf "trans error = %g\r",trans_err
800//      Printf "correction = %g +/- %g\r", correction, err
801       
802        //end of transmission/pathlength correction
803
804        return(correction)
805end
806
807
808//
809// TODO:
810//   --         DoAlert 0,"This has not yet been updated for VSANS"
811//
812//test procedure, not called anymore
813Proc AbsoluteScaling(type,c0,c1,c2,c3,c4,c5)
814        String type
815        Variable c0=1,c1=0.1,c2=0.95,c3=0.1,c4=1,c5=32.0
816        Prompt type,"WORK data type",popup,"CAL;COR;SAM"
817        Prompt c0, "Sample Transmission"
818        Prompt c1, "Sample Thickness (cm)"
819        Prompt c2, "Standard Transmission"
820        Prompt c3, "Standard Thickness (cm)"
821        Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)"
822        Prompt c5, "Standard Cross-Section (cm-1)"
823
824        Variable err
825        //call the function to do the math
826        //data from "type" will be scaled and deposited in ABS
827        err = Absolute_Scale(type,c0,c1,c2,c3,c4,c5)
828       
829        if(err)
830                Abort "Error in Absolute_Scale()"
831        endif
832       
833        //contents are always dumped to ABS
834        type = "ABS"
835       
836        String newTitle = "WORK_"+type
837        DoWindow/F SANS_Data
838        DoWindow/T SANS_Data, newTitle
839        KillStrings/Z newTitle
840       
841        //need to update the display with "data" from the correct dataFolder
842        //reset the current display type to "type"
843        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
844        gCurDispType = Type     
845       
846        fRawWindowHook()
847       
848End
849
850//
851// TODO:
852//   --         DoAlert 0,"This has not yet been updated for VSANS"
853//
854//s_ is the standard
855//w_ is the "work" file
856//both are work files and should already be normalized to 10^8 monitor counts
857Function Absolute_Scale(type,w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err)
858        String type
859        Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
860
861        DoAlert 0,"This has not yet been updated for VSANS"
862               
863        //convert the "type" data to absolute scale using the given standard information
864        //copying the "type" waves to ABS
865       
866        //check for existence of data, rescale to linear if needed
867        String destPath
868        //check for "type"
869        if(WaveExists($("root:Packages:NIST:"+Type + ":data")) == 0)
870                Print "There is no work file in "+type+"--Aborting"
871                Return(1)               //error condition
872        Endif
873        //check for log-scaling of the "type" data and adjust if necessary
874        destPath = "root:Packages:NIST:"+Type
875        NVAR gIsLogScale = $(destPath + ":gIsLogScale")
876        if(gIsLogScale)
877                Duplicate/O $(destPath + ":linear_data") $(destPath + ":data")//back to linear scale
878                Variable/G $(destPath + ":gIsLogScale")=0       //the "type" data is not logscale anymore
879        endif
880       
881        //copy "oldtype" information to ABS
882        //overwriting out the old contents of the ABS folder (/O option in Duplicate)
883        //copy over the waves data,vlegend,text,integers,reals(read)
884
885        String oldType= "root:Packages:NIST:"+type              //this is where the data to be absoluted is
886        //copy from current dir (type) to ABS, defined by destPath
887        Duplicate/O $(oldType + ":data"),$"root:Packages:NIST:ABS:data"
888        Duplicate/O $(oldType + ":linear_data"),$"root:Packages:NIST:ABS:linear_data"
889        Duplicate/O $(oldType + ":linear_data_error"),$"root:Packages:NIST:ABS:linear_data_error"
890//      Duplicate/O $(oldType + ":vlegend"),$"root:Packages:NIST:ABS:vlegend"
891        Duplicate/O $(oldType + ":textread"),$"root:Packages:NIST:ABS:textread"
892        Duplicate/O $(oldType + ":integersread"),$"root:Packages:NIST:ABS:integersread"
893        Duplicate/O $(oldType + ":realsread"),$"root:Packages:NIST:ABS:realsread"
894        //need to save a copy of filelist string too (from the current type folder)
895        SVAR oldFileList = $(oldType + ":fileList")
896        //need to copy filelist string too
897        String/G $"root:Packages:NIST:ABS:fileList" = oldFileList
898       
899        //now switch to ABS folder
900        //make appropriate wave references
901        WAVE data=$"root:Packages:NIST:ABS:linear_data"                                 // these wave references point to the "type" data in ABS
902        WAVE data_err=$"root:Packages:NIST:ABS:linear_data_error"                                       // these wave references point to the "type" data in ABS
903        WAVE data_copy=$"root:Packages:NIST:ABS:data"                                   // just for display
904        WAVE/T textread=$"root:Packages:NIST:ABS:textread"                      //that are to be directly operated on
905        WAVE integersread=$"root:Packages:NIST:ABS:integersread"
906        WAVE realsread=$"root:Packages:NIST:ABS:realsread"
907        Variable/G $"root:Packages:NIST:ABS:gIsLogscale"=0                      //make new flag in ABS folder, data is linear scale
908       
909        //do the actual absolute scaling here, modifying the data in ABS
910        Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
911       
912        w_moncount = realsread[0]               //monitor count in "type"
913        if(w_moncount == 0)
914                //zero monitor counts will give divide by zero ---
915                DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
916                Return(1)               //report error
917        Endif
918       
919        //calculate scale factor
920        Variable scale,trans_err
921        s1 = defmon/realsread[0]                //[0] is monitor count (s1 should be 1)
922        s2 = s_thick/w_thick
923        s3 = s_trans/w_trans
924        s4 = s_cross/s_izero
925       
926        // kappa comes in as s_izero, so be sure to use 1/kappa_err
927       
928        data *= s1*s2*s3*s4
929       
930        scale = s1*s2*s3*s4
931        trans_err = realsRead[41]
932       
933//      print scale
934//      print data[0][0]
935       
936        data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
937
938//      print data_err[0][0]
939       
940// keep "data" in sync with linear_data
941        data_copy = data
942       
943        //********* 15APR02
944        // DO NOt correct for atenuators here - the COR step already does this, putting all of the data one equal
945        // footing (zero atten) before doing the subtraction.
946        //
947        //Print "ABS data multiplied by  ",s1*s2*s3*s4/attenFactor
948       
949        //update the ABS header information
950        textread[1] = date() + " " + time()             //date + time stamp
951       
952        Return (0) //no error
953End
954
955
956//
957// TODO:
958//   --         DoAlert 0,"This has not yet been updated for VSANS"
959//
960//
961// match the attenuation of the RAW data to the "type" data
962// so that they can be properly added
963//
964// are the attenuator numbers the same? if so exit
965//
966// if not, find the attenuator number for type
967// - find both attenuation factors
968//
969// rescale the raw data to match the ratio of the two attenuation factors
970// -- adjust the detector count (rw)
971// -- the linear data
972//
973//
974Function Adjust_RAW_Attenuation(type)
975        String type
976
977        DoAlert 0,"This has not yet been updated for VSANS"
978       
979        WAVE rw=$("root:Packages:NIST:RAW:realsread")
980        WAVE linear_data=$("root:Packages:NIST:RAW:linear_data")
981        WAVE data=$("root:Packages:NIST:RAW:data")
982        WAVE data_err=$("root:Packages:NIST:RAW:linear_data_error")
983        WAVE/T tw = $("root:Packages:NIST:RAW:textRead")
984       
985        WAVE dest_reals=$("root:Packages:NIST:"+type+":realsread")
986
987        Variable dest_atten,raw_atten,tol
988        Variable lambda,raw_atten_err,raw_AttenFactor,dest_attenFactor,dest_atten_err
989        String fileStr
990
991        dest_atten = dest_reals[3]
992        raw_atten = rw[3]
993       
994        tol = 0.1               // within 0.1 atten units is OK
995        if(abs(dest_atten - raw_atten) < tol )
996                return(0)
997        endif
998
999        fileStr = tw[3]
1000        lambda = rw[26]
1001        // TODO access correct values
1002        raw_AttenFactor = 1//AttenuationFactor(fileStr,lambda,raw_atten,raw_atten_err)
1003        dest_AttenFactor = 1//AttenuationFactor(fileStr,lambda,dest_atten,dest_atten_err)
1004               
1005        rw[2] *= dest_AttenFactor/raw_AttenFactor
1006        linear_data *= dest_AttenFactor/raw_AttenFactor
1007       
1008        // to keep "data" and linear_data in sync
1009        data = linear_data
1010       
1011        return(0)
1012End
1013
1014//
1015// TODO:
1016//   --         DoAlert 0,"This has not yet been updated for VSANS"
1017//
1018//************************
1019//unused testing procedure, may not be up-to-date with other procedures
1020//check before re-implementing
1021//
1022Proc DIV_a_Workfile(type)
1023        String type
1024        Prompt type,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
1025       
1026        //macro will take whatever is in SELECTED folder and DIVide it by the current
1027        //contents of the DIV folder - the function will check for existence
1028        //before proceeding
1029
1030//DoAlert 0,"This has not yet been updated for VSANS"
1031       
1032        Variable err
1033        err = DIVCorrection(type)               //returns err = 1 if data doesn't exist in specified folders
1034       
1035        if(err)
1036                Abort "error in DIVCorrection()"
1037        endif
1038       
1039        //contents are NOT always dumped to CAL, but are in the new type folder
1040       
1041        String newTitle = "WORK_"+type
1042        DoWindow/F VSANS_Data
1043        DoWindow/T VSANS_Data, newTitle
1044        KillStrings/Z newTitle
1045       
1046        //need to update the display with "data" from the correct dataFolder
1047        //reset the current displaytype to "type"
1048        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=Type
1049       
1050        UpdateDisplayInformation(type)
1051       
1052End
1053
1054
1055//
1056// TODO:
1057//   --         DoAlert 0,"This has not yet been updated for VSANS"
1058//  -- how is the error propagation  handled?
1059//
1060//function will divide the contents of "workType" folder with the contents of
1061//the DIV folder + detStr
1062// all data is linear scale for the calculation
1063//
1064Function DIVCorrection(data,data_err,detStr,workType)
1065        Wave data,data_err
1066        String detStr,workType
1067       
1068        //check for existence of data in type and DIV
1069        // if the desired data doesn't exist, let the user know, and abort
1070        String destPath=""
1071
1072        if(WaveExists(data) == 0)
1073                Print "The data wave does not exist in DIVCorrection()"
1074                Return(1)               //error condition
1075        Endif
1076       
1077        //check for DIV
1078        // if the DIV workfile doesn't exist, let the user know,and abort
1079
1080        WAVE/Z div_data = $("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")
1081        if(WaveExists(div_data) == 0)
1082                Print "The DIV wave does not exist in DIVCorrection()"
1083                Return(1)               //error condition
1084        Endif
1085        //files exist, proceed
1086
1087        data /= div_data
1088       
1089//      data_err /= div_data
1090       
1091        Return(0)
1092End
1093
1094
1095//////////////////////////
Note: See TracBrowser for help on using the repository browser.