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

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

LOTS of changes to accommodate the beam center being reported in cm rather than pixels. Required a lot of changes to VCALC (to fill in simulated data), and in the reading and interpreting of data for display, and most importantly, the calculation of q.

There may still be a few residual bugs with this. I am still re-testing with new sample data sets.

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