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

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

lots of changes to plotting of q-values, generating fake data with non-linear corrections, masking of data, etc.

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