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

Last change on this file since 992 was 992, checked in by srkline, 7 years ago

removed the doubled "entry" field from the VSANS file load.

appears now to work fine with R/W routines and with VCALC.

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