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

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

lots of changes to 1D averaging and the plotting routines, detector corrections, and basic reads

File size: 32.5 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:entry:instrument:detector_"+detStr+":data_realDistX")
120        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:entry:instrument:detector_"+detStr+":data_realDistY")
121        Wave data_realDistX = $(destPath + ":entry:entry:instrument:detector_"+detStr+":data_realDistX")
122        Wave data_realDistY = $(destPath + ":entry: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:entry:instrument:detector_"+detStr+":data_realDistX")
167        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:entry:instrument:detector_"+detStr+":data_realDistY")
168        Wave data_realDistX = $(destPath + ":entry:entry:instrument:detector_"+detStr+":data_realDistX")
169        Wave data_realDistY = $(destPath + ":entry: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:entry:instrument:detector_"+detStr+":data_realDistX")
200        Wave data_realDistY = $(destPath + ":entry: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:entry:instrument:detector_"+detStr+":beam_center_x_mm")
216        Make/O/D/N=1 $(destPath + ":entry:entry:instrument:detector_"+detStr+":beam_center_y_mm")
217        WAVE x_mm = $(destPath + ":entry:entry:instrument:detector_"+detStr+":beam_center_x_mm")
218        WAVE y_mm = $(destPath + ":entry: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:entry:instrument:detector_"+detStr+":data_realDistX")
258        Wave data_realDistY = $(destPath + ":entry: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:entry:instrument:detector_"+detStr+":beam_center_x_mm")
265        Make/O/D/N=1 $(destPath + ":entry:entry:instrument:detector_"+detStr+":beam_center_y_mm")
266        WAVE x_mm = $(destPath + ":entry:entry:instrument:detector_"+detStr+":beam_center_x_mm")
267        WAVE y_mm = $(destPath + ":entry: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:entry:instrument:detector_"+detStr+":spatial_calibration")
340                Wave calib = $("root:Packages:NIST:VSANS:RAW:entry: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:entry:instrument:detector_"+detStr+":data_realDistX")
397        Wave data_realDistY = $(destPath + ":entry:entry:instrument:detector_"+detStr+":data_realDistY")
398
399// make the new waves
400        Duplicate/O data_realDistX $(destPath + ":entry:entry:instrument:detector_"+detStr+":qTot_"+detStr)
401        Duplicate/O data_realDistX $(destPath + ":entry:entry:instrument:detector_"+detStr+":qx_"+detStr)
402        Duplicate/O data_realDistX $(destPath + ":entry:entry:instrument:detector_"+detStr+":qy_"+detStr)
403        Duplicate/O data_realDistX $(destPath + ":entry:entry:instrument:detector_"+detStr+":qz_"+detStr)
404        Wave qTot = $(destPath + ":entry:entry:instrument:detector_"+detStr+":qTot_"+detStr)
405        Wave qx = $(destPath + ":entry:entry:instrument:detector_"+detStr+":qx_"+detStr)
406        Wave qy = $(destPath + ":entry:entry:instrument:detector_"+detStr+":qy_"+detStr)
407        Wave qz = $(destPath + ":entry: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 displaytype to "type"
817        String/G root:myGlobals:gDataDisplayType=Type
818       
819        fRawWindowHook()
820       
821End
822
823//
824// TODO:
825//   --         DoAlert 0,"This has not yet been updated for VSANS"
826//
827//s_ is the standard
828//w_ is the "work" file
829//both are work files and should already be normalized to 10^8 monitor counts
830Function Absolute_Scale(type,w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err)
831        String type
832        Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
833
834        DoAlert 0,"This has not yet been updated for VSANS"
835               
836        //convert the "type" data to absolute scale using the given standard information
837        //copying the "type" waves to ABS
838       
839        //check for existence of data, rescale to linear if needed
840        String destPath
841        //check for "type"
842        if(WaveExists($("root:Packages:NIST:"+Type + ":data")) == 0)
843                Print "There is no work file in "+type+"--Aborting"
844                Return(1)               //error condition
845        Endif
846        //check for log-scaling of the "type" data and adjust if necessary
847        destPath = "root:Packages:NIST:"+Type
848        NVAR gIsLogScale = $(destPath + ":gIsLogScale")
849        if(gIsLogScale)
850                Duplicate/O $(destPath + ":linear_data") $(destPath + ":data")//back to linear scale
851                Variable/G $(destPath + ":gIsLogScale")=0       //the "type" data is not logscale anymore
852        endif
853       
854        //copy "oldtype" information to ABS
855        //overwriting out the old contents of the ABS folder (/O option in Duplicate)
856        //copy over the waves data,vlegend,text,integers,reals(read)
857
858        String oldType= "root:Packages:NIST:"+type              //this is where the data to be absoluted is
859        //copy from current dir (type) to ABS, defined by destPath
860        Duplicate/O $(oldType + ":data"),$"root:Packages:NIST:ABS:data"
861        Duplicate/O $(oldType + ":linear_data"),$"root:Packages:NIST:ABS:linear_data"
862        Duplicate/O $(oldType + ":linear_data_error"),$"root:Packages:NIST:ABS:linear_data_error"
863//      Duplicate/O $(oldType + ":vlegend"),$"root:Packages:NIST:ABS:vlegend"
864        Duplicate/O $(oldType + ":textread"),$"root:Packages:NIST:ABS:textread"
865        Duplicate/O $(oldType + ":integersread"),$"root:Packages:NIST:ABS:integersread"
866        Duplicate/O $(oldType + ":realsread"),$"root:Packages:NIST:ABS:realsread"
867        //need to save a copy of filelist string too (from the current type folder)
868        SVAR oldFileList = $(oldType + ":fileList")
869        //need to copy filelist string too
870        String/G $"root:Packages:NIST:ABS:fileList" = oldFileList
871       
872        //now switch to ABS folder
873        //make appropriate wave references
874        WAVE data=$"root:Packages:NIST:ABS:linear_data"                                 // these wave references point to the "type" data in ABS
875        WAVE data_err=$"root:Packages:NIST:ABS:linear_data_error"                                       // these wave references point to the "type" data in ABS
876        WAVE data_copy=$"root:Packages:NIST:ABS:data"                                   // just for display
877        WAVE/T textread=$"root:Packages:NIST:ABS:textread"                      //that are to be directly operated on
878        WAVE integersread=$"root:Packages:NIST:ABS:integersread"
879        WAVE realsread=$"root:Packages:NIST:ABS:realsread"
880        Variable/G $"root:Packages:NIST:ABS:gIsLogscale"=0                      //make new flag in ABS folder, data is linear scale
881       
882        //do the actual absolute scaling here, modifying the data in ABS
883        Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
884       
885        w_moncount = realsread[0]               //monitor count in "type"
886        if(w_moncount == 0)
887                //zero monitor counts will give divide by zero ---
888                DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
889                Return(1)               //report error
890        Endif
891       
892        //calculate scale factor
893        Variable scale,trans_err
894        s1 = defmon/realsread[0]                //[0] is monitor count (s1 should be 1)
895        s2 = s_thick/w_thick
896        s3 = s_trans/w_trans
897        s4 = s_cross/s_izero
898       
899        // kappa comes in as s_izero, so be sure to use 1/kappa_err
900       
901        data *= s1*s2*s3*s4
902       
903        scale = s1*s2*s3*s4
904        trans_err = realsRead[41]
905       
906//      print scale
907//      print data[0][0]
908       
909        data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
910
911//      print data_err[0][0]
912       
913// keep "data" in sync with linear_data
914        data_copy = data
915       
916        //********* 15APR02
917        // DO NOt correct for atenuators here - the COR step already does this, putting all of the data one equal
918        // footing (zero atten) before doing the subtraction.
919        //
920        //Print "ABS data multiplied by  ",s1*s2*s3*s4/attenFactor
921       
922        //update the ABS header information
923        textread[1] = date() + " " + time()             //date + time stamp
924       
925        Return (0) //no error
926End
927
928
929//
930// TODO:
931//   --         DoAlert 0,"This has not yet been updated for VSANS"
932//
933//
934// match the attenuation of the RAW data to the "type" data
935// so that they can be properly added
936//
937// are the attenuator numbers the same? if so exit
938//
939// if not, find the attenuator number for type
940// - find both attenuation factors
941//
942// rescale the raw data to match the ratio of the two attenuation factors
943// -- adjust the detector count (rw)
944// -- the linear data
945//
946//
947Function Adjust_RAW_Attenuation(type)
948        String type
949
950        DoAlert 0,"This has not yet been updated for VSANS"
951       
952        WAVE rw=$("root:Packages:NIST:RAW:realsread")
953        WAVE linear_data=$("root:Packages:NIST:RAW:linear_data")
954        WAVE data=$("root:Packages:NIST:RAW:data")
955        WAVE data_err=$("root:Packages:NIST:RAW:linear_data_error")
956        WAVE/T tw = $("root:Packages:NIST:RAW:textRead")
957       
958        WAVE dest_reals=$("root:Packages:NIST:"+type+":realsread")
959
960        Variable dest_atten,raw_atten,tol
961        Variable lambda,raw_atten_err,raw_AttenFactor,dest_attenFactor,dest_atten_err
962        String fileStr
963
964        dest_atten = dest_reals[3]
965        raw_atten = rw[3]
966       
967        tol = 0.1               // within 0.1 atten units is OK
968        if(abs(dest_atten - raw_atten) < tol )
969                return(0)
970        endif
971
972        fileStr = tw[3]
973        lambda = rw[26]
974        // TODO access correct values
975        raw_AttenFactor = 1//AttenuationFactor(fileStr,lambda,raw_atten,raw_atten_err)
976        dest_AttenFactor = 1//AttenuationFactor(fileStr,lambda,dest_atten,dest_atten_err)
977               
978        rw[2] *= dest_AttenFactor/raw_AttenFactor
979        linear_data *= dest_AttenFactor/raw_AttenFactor
980       
981        // to keep "data" and linear_data in sync
982        data = linear_data
983       
984        return(0)
985End
986
987//
988// TODO:
989//   --         DoAlert 0,"This has not yet been updated for VSANS"
990//
991//************************
992//unused testing procedure, may not be up-to-date with other procedures
993//check before re-implementing
994//
995Macro DIV_a_Workfile(type)
996        String type
997        Prompt type,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
998       
999        //macro will take whatever is in SELECTED folder and DIVide it by the current
1000        //contents of the DIV folder - the function will check for existence
1001        //before proceeding
1002
1003//DoAlert 0,"This has not yet been updated for VSANS"
1004       
1005        Variable err
1006        err = DIVCorrection(type)               //returns err = 1 if data doesn't exist in specified folders
1007       
1008        if(err)
1009                Abort "error in DIVCorrection()"
1010        endif
1011       
1012        //contents are NOT always dumped to CAL, but are in the new type folder
1013       
1014        String newTitle = "WORK_"+type
1015        DoWindow/F VSANS_Data
1016        DoWindow/T VSANS_Data, newTitle
1017        KillStrings/Z newTitle
1018       
1019        //need to update the display with "data" from the correct dataFolder
1020        //reset the current displaytype to "type"
1021        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=Type
1022       
1023        UpdateDisplayInformation(type)
1024       
1025End
1026
1027
1028//
1029// TODO:
1030//   --         DoAlert 0,"This has not yet been updated for VSANS"
1031//
1032//function will divide the contents of "workType" folder with the contents of
1033//the DIV folder + detStr
1034// all data is linear scale for the calculation
1035//
1036Function DIVCorrection(data,data_err,detStr,workType)
1037        Wave data,data_err
1038        String detStr,workType
1039       
1040        //check for existence of data in type and DIV
1041        // if the desired data doesn't exist, let the user know, and abort
1042        String destPath=""
1043
1044        if(WaveExists(data) == 0)
1045                Print "The data wave does not exist in DIVCorrection()"
1046                Return(1)               //error condition
1047        Endif
1048       
1049        //check for DIV
1050        // if the DIV workfile doesn't exist, let the user know,and abort
1051
1052        WAVE/Z div_data = $("root:Packages:NIST:VSANS:DIV:entry:entry:instrument:detector_"+detStr+":data")
1053        if(WaveExists(div_data) == 0)
1054                Print "The DIV wave does not exist in DIVCorrection()"
1055                Return(1)               //error condition
1056        Endif
1057        //files exist, proceed
1058
1059        data /= div_data
1060       
1061//      data_err /= div_data
1062       
1063        Return(0)
1064End
1065
1066
1067//////////////////////////
Note: See TracBrowser for help on using the repository browser.