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

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

Added USANS loader/initializer

Updated units of distance for q-calculation (SDD)

status display bug fixed

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