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

Last change on this file since 1062 was 1062, checked in by srkline, 5 years ago

many changes to the VCALC code and a few changes to the main code to get the units consistent, and centimeters everywhere possible. The real space distance array and the non-linear calibrations are still defined and calculated in mm. This can hopefully be changed in the future. Some constants in the data file will need to be updated to cm, such as the T/B setback, which has been confirmed to be 41.0 cm

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