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

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

changes here are mostly to facilitate the generation and display of DIV files. Set up a simple sequence to normalize DIV data sets and copy the normalized panels to the temporary folder for saving. Made a simple panel to display the 4 panels from a carriage, and to do simple subtract or divide comparison to track changes in DIV over time.

For mask drawing, made the overlay semi-transparent so that the original data can be seen through it - but kept the toggle. log/lin scaling is now based on the VSANS preference.

Adjusted the mouse tracking of the pixel on the RAW data display to properly track the T/B panels with their half-size pixels.

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