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

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

proper handling of setback, intent, and purpose. New value of setback is included, and proper definition of how purpose and intent are defined within NICE and the GUI are partly completed now. Searches still need to be updated.

File size: 39.6 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
574
575        sdd = V_getDet_ActualDistance(fname,detStr)             //sdd derived, including setback [cm]
576
577        // this is the ctr in pixels --xx-- (now it is in cm!)
578//      xCtr = V_getDet_beam_center_x(fname,detStr)
579//      yCtr = V_getDet_beam_center_y(fname,detStr)
580        // this is ctr in mm
581        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
582        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
583        lambda = V_getWavelength(fname)
584        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
585        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
586
587// make the new waves
588        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
589        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
590        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
591        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
592        Wave qTot = $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
593        Wave qx = $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
594        Wave qy = $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
595        Wave qz = $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
596
597// calculate all of the q-values
598// sdd is passed in [cm]
599        qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
600        qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
601        qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
602        qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
603       
604       
605        return(0)
606End
607
608
609//function to calculate the overall q-value, given all of the necesary trig inputs
610//
611// TODO:
612// -- verify the calculation (accuracy - in all input conditions)
613// -- verify the units of everything here, it's currently all jumbled and wrong... and repeated
614// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
615//    to each pixel
616//
617//sdd is in [cm]
618// distX and distY are in [mm]
619//wavelength is in Angstroms
620//
621//returned magnitude of Q is in 1/Angstroms
622//
623Function V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
624        Variable xaxval,yaxval,xctr,yctr,sdd,lam
625        Wave distX,distY
626       
627        Variable dx,dy,qval,two_theta,dist
628               
629
630        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
631        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
632        dist = sqrt(dx^2 + dy^2)
633       
634        dist /= 10  // convert mm to cm
635       
636        two_theta = atan(dist/sdd)
637
638        qval = 4*Pi/lam*sin(two_theta/2)
639       
640        return qval
641End
642
643//calculates just the q-value in the x-direction on the detector
644// TODO:
645// -- verify the calculation (accuracy - in all input conditions)
646// -- verify the units of everything here, it's currently all jumbled and wrong... and repeated
647// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
648//    to each pixel
649//
650//
651// this properly accounts for qz
652//
653Function V_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
654        Variable xaxval,yaxval,xctr,yctr,sdd,lam
655        Wave distX,distY
656
657        Variable qx,qval,phi,dx,dy,dist,two_theta
658       
659        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
660       
661
662        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
663        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
664        phi = V_FindPhi(dx,dy)
665       
666        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
667        dist = sqrt(dx^2 + dy^2)
668        dist /= 10  // convert mm to cm
669
670        two_theta = atan(dist/sdd)
671
672        qx = qval*cos(two_theta/2)*cos(phi)
673       
674        return qx
675End
676
677//calculates just the q-value in the y-direction on the detector
678// TODO:
679// -- verify the calculation (accuracy - in all input conditions)
680// -- verify the units of everything here, it's currently all jumbled and wrong... and repeated
681// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
682//    to each pixel
683//
684//
685// this properly accounts for qz
686//
687Function V_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
688        Variable xaxval,yaxval,xctr,yctr,sdd,lam
689        Wave distX,distY
690
691        Variable qy,qval,phi,dx,dy,dist,two_theta
692       
693        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
694       
695
696        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
697        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
698        phi = V_FindPhi(dx,dy)
699       
700        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
701        dist = sqrt(dx^2 + dy^2)
702        dist /= 10  // convert mm to cm
703
704        two_theta = atan(dist/sdd)
705
706        qy = qval*cos(two_theta/2)*sin(phi)
707       
708        return qy
709End
710
711//calculates just the q-value in the z-direction on the detector
712// TODO:
713// -- verify the calculation (accuracy - in all input conditions)
714// -- verify the units of everything here, it's currently all jumbled and wrong... and repeated
715// -- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
716//    to each pixel
717//
718// not actually used for any calculations, but here for completeness if anyone asks, or for 2D data export
719//
720// this properly accounts for qz
721//
722Function V_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
723        Variable xaxval,yaxval,xctr,yctr,sdd,lam
724        Wave distX,distY
725
726        Variable qz,qval,phi,dx,dy,dist,two_theta
727       
728        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
729       
730
731        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
732        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
733       
734        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
735        dist = sqrt(dx^2 + dy^2)
736        dist /= 10  // convert mm to cm
737
738        two_theta = atan(dist/sdd)
739
740        qz = qval*sin(two_theta/2)
741       
742        return qz
743End
744
745
746//
747// TODO -- VERIFY calculations
748// -- This is the actual solid angle per pixel, not a ratio vs. some "unit SA"
749//    Do I just correct for the different area vs. the "nominal" central area?
750// -- decide how to implement - either directly change the data values (as was done in the past)
751//    or use this as a weighting for when the data is binned to I(q). In the second method, 2D data
752//    would need this to be applied before exporting
753// -- do I keep a wave note indicating that this correction has been applied to the data
754//    so that it can be "un-applied"?
755// -- do I calculate theta from geometry directly, or get it from Q (Assuming it's present?)
756//    (probably just from geometry, since I need SDD and dx and dy values...)
757//
758//
759Function V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
760        Wave w,w_err
761        String fname,detStr,destPath
762
763        Variable sdd,xCtr,yCtr,lambda
764
765// get all of the geometry information 
766//      orientation = V_getDet_tubeOrientation(fname,detStr)
767        sdd = V_getDet_ActualDistance(fname,detStr)
768
769
770        // this is ctr in mm
771        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
772        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
773        lambda = V_getWavelength(fname)
774       
775        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
776       
777        Wave data_realDistX = data_realDistX
778        Wave data_realDistY = data_realDistY
779
780        Duplicate/O w solid_angle,tmp_theta,tmp_dist            //in the current df
781
782//// calculate the scattering angle
783//      dx = (distX - xctr)             //delta x in mm
784//      dy = (distY - yctr)             //delta y in mm
785        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
786       
787        tmp_dist /= 10  // convert mm to cm
788        // sdd is in [cm]
789
790        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle
791
792        Variable ii,jj,numx,numy,dx,dy
793        numx = DimSize(tmp_theta,0)
794        numy = DimSize(tmp_theta,1)
795       
796        for(ii=0        ;ii<numx;ii+=1)
797                for(jj=0;jj<numy;jj+=1)
798                       
799                        if(ii==0)               //do a forward difference if ii==0
800                                dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj])        //delta x for the pixel
801                        else
802                                dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj])        //delta x for the pixel
803                        endif
804                       
805                       
806                        if(jj==0)
807                                dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj])        //delta y for the pixel
808                        else
809                                dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1])        //delta y for the pixel
810                        endif
811       
812                        dx /= 10
813                        dy /= 10                // convert mm to cm (since sdd is in cm)
814                        solid_angle[ii][jj] = dx*dy             //this is in cm^2
815                endfor
816        endfor
817       
818        // to cover up any issues w/negative dx or dy
819        solid_angle = abs(solid_angle)
820       
821        // solid_angle correction
822        // == dx*dy*cos^3/sdd^2
823        solid_angle *= (cos(tmp_theta))^3
824        solid_angle /= sdd^2
825       
826        // Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!)
827        w /= solid_angle
828       
829       
830        // TODO:
831        // correctly apply the correction to the error wave (assume a perfect value?)
832        w_err /= solid_angle            //is this correct??
833
834// TODO -- clean up after I'm satisfied computations are correct               
835//      KillWaves/Z tmp_theta,tmp_dist
836       
837        return(0)
838end
839
840
841////////////
842// TODO: all of below is untested code
843//   copied from SANS
844//
845//
846// TODO :
847//   -- DoAlert 0,"This has not yet been updated for VSANS"
848//
849//performs solid angle and non-linear detector corrections to raw data as it is "added" to a work folder
850//function is called by Raw_to_work() and Add_raw_to_work() functions
851//works on the actual data array, assumes that is is already on LINEAR scale
852//
853Function V_DetCorr(data,data_err,realsread,doEfficiency,doTrans)
854        Wave data,data_err,realsread
855        Variable doEfficiency,doTrans
856
857        DoAlert 0,"This has not yet been updated for VSANS"
858       
859        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,xx0,yy0
860        Variable ii,jj,dtdist,dtdis2
861        Variable xi,xd,yd,rad,ratio,domega,xy
862        Variable lambda,trans,trans_err,lat_err,tmp_err,lat_corr
863       
864//      Print "...doing jacobian and non-linear corrections"
865
866        NVAR pixelsX = root:myGlobals:gNPixelsX
867        NVAR pixelsY = root:myGlobals:gNPixelsY
868       
869        //set up values to send to auxiliary trig functions
870        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
871        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
872
873        x0 = realsread[16]
874        y0 = realsread[17]
875        sx = realsread[10]
876        sx3 = realsread[11]
877        sy = realsread[13]
878        sy3 = realsread[14]
879       
880        dtdist = 1000*realsread[18]     //sdd in mm
881        dtdis2 = dtdist^2
882       
883        lambda = realsRead[26]
884        trans = RealsRead[4]
885        trans_err = RealsRead[41]               //new, March 2011
886       
887
888        //waves to contain repeated function calls
889        Make/O/N=(pixelsX) fyy,xx,yy            //Assumes square detector !!!
890        ii=0
891        do
892                xi = ii
893//              fyy[ii] = dc_fy(ii+1,sy,sy3,ycenter)
894//              xx[ii] = dc_fxn(ii+1,sx,sx3,xcenter)
895//              yy[ii] = dc_fym(ii+1,sy,sy3,ycenter)
896                ii+=1
897        while(ii<pixelsX)
898       
899        Make/O/N=(pixelsX,pixelsY) SolidAngle           // testing only
900       
901        ii=0
902        do
903                xi = ii
904//              xd = dc_fx(ii+1,sx,sx3,xcenter)-xx0
905                jj=0
906                do
907                        yd = fyy[jj]-yy0
908                        //rad is the distance of pixel ij from the sample
909                        //domega is the ratio of the solid angle of pixel ij versus center pixel
910                        // product xy = 1 for a detector with a linear spatial response (modern Ordela)
911                        // solid angle calculated, dW^3 >=1, so multiply data to raise measured values to correct values.
912                        rad = sqrt(dtdis2 + xd^2 + yd^2)
913                        domega = rad/dtdist
914                        ratio = domega^3
915                        xy = xx[ii]*yy[jj]
916                       
917                        data[ii][jj] *= xy*ratio
918                       
919                        solidAngle[ii][jj] = xy*ratio           //testing only 
920                        data_err[ii][jj] *= xy*ratio                    //error propagation assumes that SA and Jacobian are exact, so simply scale error
921                       
922                       
923                        // correction factor for detector efficiency JBG memo det_eff_cor2.doc 3/20/07
924                        // correction inserted 11/2007 SRK
925                        // large angle detector efficiency is >= 1 and will "bump up" the measured value at the highest angles
926                        // so divide here to get the correct answer (5/22/08 SRK)
927                        if(doEfficiency)
928//                              data[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
929//                              data_err[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
930//                              solidAngle[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)           //testing only
931                        endif
932                       
933                        // large angle transmission calculation is <= 1 and will "bump down" the measured value at the highest angles
934                        // so divide here to get the correct answer
935                        if(doTrans)
936                       
937                                if(trans<0.1 && ii==0 && jj==0)
938                                        Print "***transmission is less than 0.1*** and is a significant correction"
939                                endif
940                               
941                                if(trans==0)
942                                        if(ii==0 && jj==0)
943                                                Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
944                                        endif
945                                        trans = 1
946                                endif
947                               
948                                // pass in the transmission error, and the error in the correction is returned as the last parameter
949
950//                              lat_corr = V_LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,lat_err)           //moved from 1D avg SRK 11/2007
951
952                                data[ii][jj] /= lat_corr                        //divide by the correction factor
953                                //
954                                //
955                                //
956                                // relative errors add in quadrature
957                                tmp_err = (data_err[ii][jj]/lat_corr)^2 + (lat_err/lat_corr)^2*data[ii][jj]*data[ii][jj]/lat_corr^2
958                                tmp_err = sqrt(tmp_err)
959                               
960                                data_err[ii][jj] = tmp_err
961                               
962//                              solidAngle[ii][jj] = lat_err
963
964                               
965                                //solidAngle[ii][jj] = LargeAngleTransmissionCorr(trans,dtdist,xd,yd)           //testing only
966                        endif
967                       
968                        jj+=1
969                while(jj<pixelsX)
970                ii+=1
971        while(ii<pixelsX)
972       
973        //clean up waves
974       
975        Return(0)
976End
977
978
979
980//
981// DIVIDE the intensity by this correction to get the right answer
982// TODO:
983//   --         DoAlert 0,"This has not yet been updated for VSANS"
984//
985//
986
987// Apply the large angle transmssion correction as the data is converted to WORK
988// so that whether the data is saved as 2D or 1D, the correction has properly been done.
989//
990// This is, however, a SAMPLE dependent calculation, not purely instrument geometry.
991//
992Function V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
993        Wave w,w_err
994        String fname,detStr,destPath
995
996        Variable sdd,xCtr,yCtr,trans,trans_err,uval
997
998// get all of the geometry information 
999//      orientation = V_getDet_tubeOrientation(fname,detStr)
1000        sdd = V_getDet_ActualDistance(fname,detStr)
1001
1002        // this is ctr in mm
1003        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
1004        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
1005        trans = V_getSampleTransmission(fname)
1006        trans_err = V_getSampleTransError(fname)
1007       
1008        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
1009       
1010        Wave data_realDistX = data_realDistX
1011        Wave data_realDistY = data_realDistY
1012
1013        Duplicate/O w lat_corr,tmp_theta,tmp_dist,lat_err,tmp_err               //in the current df
1014
1015//// calculate the scattering angle
1016//      dx = (distX - xctr)             //delta x in mm
1017//      dy = (distY - yctr)             //delta y in mm
1018        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
1019       
1020        tmp_dist /= 10  // convert mm to cm
1021        // sdd is in [cm]
1022
1023        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle
1024
1025        Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp
1026        numx = DimSize(tmp_theta,0)
1027        numy = DimSize(tmp_theta,1)
1028       
1029       
1030        //optical thickness
1031        uval = -ln(trans)               //use natural logarithm
1032       
1033        for(ii=0        ;ii<numx;ii+=1)
1034                for(jj=0;jj<numy;jj+=1)
1035                       
1036                        cos_th = cos(tmp_theta[ii][jj])
1037                        arg = (1-cos_th)/cos_th
1038                       
1039                        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
1040                        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120
1041                        // OR
1042                        if((uval<0.01) || (cos_th>0.99))       
1043                                //small arg, approx correction
1044                                lat_corr[ii][jj] = 1-0.5*uval*arg
1045                        else
1046                                //large arg, exact correction
1047                                lat_corr[ii][jj] = (1-exp(-uval*arg))/(uval*arg)
1048                        endif
1049                         
1050                        // TODO
1051                        // -- properly calculate and apply the 2D error propagation
1052                        if(trans == 1)
1053                                lat_err[ii][jj] = 0             //no correction, no error
1054                        else
1055                                //sigT, calculated from the Taylor expansion
1056                                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3)
1057                                tmp *= tmp
1058                                tmp *= trans_err^2
1059                                tmp = sqrt(tmp)         //sigT
1060                               
1061                                lat_err[ii][jj] = tmp
1062                        endif
1063                         
1064 
1065                endfor
1066        endfor
1067       
1068
1069       
1070        // Here it is! Apply the correction to the intensity (divide -- to get the proper correction)
1071        w /= lat_corr
1072
1073        // relative errors add in quadrature to the current 2D error
1074        tmp_err = (w_err/lat_corr)^2 + (lat_err/lat_corr)^2*w*w/lat_corr^2
1075        tmp_err = sqrt(tmp_err)
1076       
1077        w_err = tmp_err
1078       
1079        // TODO:
1080        // correctly apply the correction to the error wave (assume a perfect value?)
1081        // w_err /= tmp         //is this correct??
1082
1083        // TODO -- clean up after I'm satisfied computations are correct               
1084        KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err
1085       
1086        return(0)
1087end
1088
1089
1090//
1091// TODO:
1092//   --         DoAlert 0,"This has not yet been updated for VSANS"
1093//
1094//test procedure, not called anymore
1095Proc V_AbsoluteScaling(type,c0,c1,c2,c3,c4,c5,I_err)
1096        String type
1097        Variable c0=1,c1=0.1,c2=0.95,c3=0.1,c4=1,c5=32.0,I_err=0.32
1098        Prompt type,"WORK data type",popup,"CAL;COR;SAM"
1099        Prompt c0, "Sample Transmission"
1100        Prompt c1, "Sample Thickness (cm)"
1101        Prompt c2, "Standard Transmission"
1102        Prompt c3, "Standard Thickness (cm)"
1103        Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)"
1104        Prompt c5, "Standard Cross-Section (cm-1)"
1105        Prompt I_err, "error in I(q=0) (one std dev)"
1106
1107        Variable err
1108        //call the function to do the math
1109        //data from "type" will be scaled and deposited in ABS
1110        err = V_Absolute_Scale(type,c0,c1,c2,c3,c4,c5,I_err)
1111       
1112        if(err)
1113                Abort "Error in V_Absolute_Scale()"
1114        endif
1115       
1116        //contents are always dumped to ABS
1117        type = "ABS"
1118       
1119        //need to update the display with "data" from the correct dataFolder
1120        //reset the current display type to "type"
1121        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
1122        gCurDispType = Type     
1123       
1124        V_UpdateDisplayInformation(Type)
1125       
1126End
1127
1128//
1129// TODO:
1130//
1131// kappa comes in as s_izero, so be sure to use 1/kappa_err
1132//
1133//convert the "type" data to absolute scale using the given standard information
1134//s_ is the standard
1135//w_ is the "work" file
1136//both are work files and should already be normalized to 10^8 monitor counts
1137Function V_Absolute_Scale(type,w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err)
1138        String type
1139        Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
1140
1141
1142        Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
1143        Variable scale,trans_err
1144        Variable err,ii
1145        String detStr
1146       
1147        // be sure that the starting data exists
1148        err = V_WorkDataExists(type)
1149        if(err==1)
1150                return(err)
1151        endif
1152               
1153        //copy from current dir (type) to ABS
1154        V_CopyHDFToWorkFolder(type,"ABS")       
1155
1156       
1157        w_moncount = V_getMonitorCount(type)            //monitor count in "type"
1158        if(w_moncount == 0)
1159                //zero monitor counts will give divide by zero ---
1160                DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
1161                Return(1)               //report error
1162        Endif
1163       
1164        //calculate scale factor
1165        s1 = defmon/w_moncount          // monitor count (s1 should be 1)
1166        s2 = s_thick/w_thick
1167        s3 = s_trans/w_trans
1168        s4 = s_cross/s_izero
1169        scale = s1*s2*s3*s4
1170
1171        trans_err = V_getSampleTransError(type)
1172       
1173        // kappa comes in as s_izero, so be sure to use 1/kappa_err
1174
1175        // and now loop through all of the detectors
1176        //do the actual absolute scaling here, modifying the data in ABS
1177        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
1178                detStr = StringFromList(ii, ksDetectorListAll, ";")
1179                Wave data = V_getDetectorDataW("ABS",detStr)
1180                Wave data_err = V_getDetectorDataErrW("ABS",detStr)
1181               
1182                data *= s1*s2*s3*s4
1183                data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1184        endfor
1185       
1186        //********* 15APR02
1187        // DO NOT correct for atenuators here - the COR step already does this, putting all of the data one equal
1188        // footing (zero atten) before doing the subtraction.
1189       
1190        Return (0) //no error
1191End
1192
1193
1194//
1195// TODO:
1196//   --         DoAlert 0,"This has not yet been updated for VSANS"
1197//
1198//
1199// match the attenuation of the RAW data to the "type" data
1200// so that they can be properly added
1201//
1202// are the attenuator numbers the same? if so exit
1203//
1204// if not, find the attenuator number for type
1205// - find both attenuation factors
1206//
1207// rescale the raw data to match the ratio of the two attenuation factors
1208// -- adjust the detector count (rw)
1209// -- the linear data
1210//
1211//
1212Function V_Adjust_RAW_Attenuation(type)
1213        String type
1214
1215        DoAlert 0,"This has not yet been updated for VSANS"
1216       
1217        WAVE rw=$("root:Packages:NIST:RAW:realsread")
1218        WAVE linear_data=$("root:Packages:NIST:RAW:linear_data")
1219        WAVE data=$("root:Packages:NIST:RAW:data")
1220        WAVE data_err=$("root:Packages:NIST:RAW:linear_data_error")
1221        WAVE/T tw = $("root:Packages:NIST:RAW:textRead")
1222       
1223        WAVE dest_reals=$("root:Packages:NIST:"+type+":realsread")
1224
1225        Variable dest_atten,raw_atten,tol
1226        Variable lambda,raw_atten_err,raw_AttenFactor,dest_attenFactor,dest_atten_err
1227        String fileStr
1228
1229        dest_atten = dest_reals[3]
1230        raw_atten = rw[3]
1231       
1232        tol = 0.1               // within 0.1 atten units is OK
1233        if(abs(dest_atten - raw_atten) < tol )
1234                return(0)
1235        endif
1236
1237        fileStr = tw[3]
1238        lambda = rw[26]
1239        // TODO access correct values
1240        raw_AttenFactor = 1//AttenuationFactor(fileStr,lambda,raw_atten,raw_atten_err)
1241        dest_AttenFactor = 1//AttenuationFactor(fileStr,lambda,dest_atten,dest_atten_err)
1242               
1243        rw[2] *= dest_AttenFactor/raw_AttenFactor
1244        linear_data *= dest_AttenFactor/raw_AttenFactor
1245       
1246        // to keep "data" and linear_data in sync
1247        data = linear_data
1248       
1249        return(0)
1250End
1251
1252//
1253// TODO:
1254//   --         DoAlert 0,"This has not yet been updated for VSANS"
1255//
1256//************************
1257//unused testing procedure, may not be up-to-date with other procedures
1258//check before re-implementing
1259//
1260Proc V_DIV_a_Workfile(type)
1261        String type
1262        Prompt type,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
1263       
1264        //macro will take whatever is in SELECTED folder and DIVide it by the current
1265        //contents of the DIV folder - the function will check for existence
1266        //before proceeding
1267
1268        Abort "This has not yet been updated for VSANS"
1269       
1270        Variable err
1271        err = V_DIVCorrection(type)             //returns err = 1 if data doesn't exist in specified folders
1272       
1273        if(err)
1274                Abort "error in V_DIVCorrection()"
1275        endif
1276       
1277        //contents are NOT always dumped to CAL, but are in the new type folder
1278       
1279        String newTitle = "WORK_"+type
1280        DoWindow/F VSANS_Data
1281        DoWindow/T VSANS_Data, newTitle
1282        KillStrings/Z newTitle
1283       
1284        //need to update the display with "data" from the correct dataFolder
1285        //reset the current displaytype to "type"
1286        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=Type
1287       
1288        V_UpdateDisplayInformation(type)
1289       
1290End
1291
1292
1293//
1294// TODO:
1295//   x-         DoAlert 0,"This has not yet been updated for VSANS"
1296//   -- how is the error propagation handled? Be sure it is calculated correctly when DIV is generated
1297//      and is applied correctly here...
1298//
1299//function will divide the contents of "workType" folder with the contents of
1300//the DIV folder + detStr
1301// all data is linear scale for the calculation
1302//
1303Function V_DIVCorrection(data,data_err,detStr,workType)
1304        Wave data,data_err
1305        String detStr,workType
1306       
1307        //check for existence of data in type and DIV
1308        // if the desired data doesn't exist, let the user know, and abort
1309        String destPath=""
1310
1311        if(WaveExists(data) == 0)
1312                Print "The data wave does not exist in V_DIVCorrection()"
1313                Return(1)               //error condition
1314        Endif
1315       
1316        //check for DIV
1317        // if the DIV workfile doesn't exist, let the user know,and abort
1318
1319        WAVE/Z div_data = $("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")
1320        if(WaveExists(div_data) == 0)
1321                Print "The DIV wave does not exist in V_DIVCorrection()"
1322                Return(1)               //error condition
1323        Endif
1324        //files exist, proceed
1325
1326        data /= div_data
1327
1328// TODO: -- correct the error propagation       
1329        data_err /= div_data
1330       
1331        Return(0)
1332End
1333
1334
1335//////////////////////////
Note: See TracBrowser for help on using the repository browser.