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

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

lots of changes here:
many little fixes to clean up TODO items and marke them DONE

changed the handling of the panel "gap" to split the gap evenly. Q-calculations have been re-verified with this change.

re-named the list of "bin Type" values, and added a few more choices. Streamlined how the averaging and plotting works with this list so that it can be more easily modified as different combinations of binning are envisioned. This resulted in a lot of excess code being cut out and replaced with cleaner logic. This change has also been verified to work as intended.

Attenuation is now always calculated from the table. The table also by (NEW) definition has values for the white beam (one waelength) and graphite (multiple possible wavelengths) where the wavelengths are artificially scaled (*1000) or *1e6) so that the interpolations can be done internally without the need for multiple attenuator tables.

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