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

Last change on this file since 1147 was 1147, checked in by srkline, 4 years ago

bug fix for Igor 8 - removing the setIgorOption forceCOLORONCOLOR for windows graphics. No longer needed in Igor 8, and generates an error. removed this line in the initialization step in all packages.

bug fix for the defining of a mask for sector averaging in VSANS. beam center in [cm] was incorrectly converted to [pixels]. No other calculations were affected.

File size: 47.4 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2#pragma version=1.0
3#pragma IgorVersion=6.1
4
5
6
7//
8// functions to apply corrections to the detector panels
9//
10// these are meant to be called by the procedures that convert "raw" data to
11// "adjusted" or corrected data sets
12//
13
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// ** its distance from the nominal beam center of (0,0) **
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 constant values
155        Variable offset,gap
156
157// kPanelTouchingGap is in mm   
158// the gap is split equally between the panel pairs
159// (DONE) -- replace all of this with V_getDet_panel_gap(fname,detStr) once it is added to the file
160
161        gap = V_getDet_panel_gap(fname,detStr)
162
163// TODO:
164// -- once the gap fields have been verified, this check can be removed
165// -- it should only apply to data pre-2018 when the field did not exist in the file
166// -- any VSANS data from 2018+ should read gap from the file.
167
168        if(gap < -100)          //-999999 returned if field is missing from file
169       
170                if(cmpstr(detStr,"FL") == 0 || cmpstr(detStr,"FR") == 0)
171                        gap = 3.5               //mm (measured, JB 1/4/18)
172                endif
173                if(cmpstr(detStr,"FT") == 0 || cmpstr(detStr,"FB") == 0)
174                        gap = 3.3               //mm (measured, JB 2/1/18)
175                endif
176                if(cmpstr(detStr,"ML") == 0 || cmpstr(detStr,"MR") == 0)
177                        gap = 5.9               //mm (measured, JB 1/4/18)
178                endif
179                if(cmpstr(detStr,"MT") == 0 || cmpstr(detStr,"MB") == 0)
180                        gap = 18.3              //mm (measured, JB 2/1/18)
181                endif
182       
183        endif
184
185       
186        if(cmpstr(orientation,"vertical")==0)
187                //      this is data dimensioned as (Ntubes,Npix)
188       
189                // adjust the x postion based on the beam center being nominally (0,0) in units of cm, not pixels
190                if(cmpstr(fname,"VCALC")== 0 )
191                        offset = VCALC_getPanelTranslation(detStr)
192                        offset *= 10                    // convert to units of mm
193//                      if(cmpstr("L",detStr[1]) == 0)
194//                              offset *= -1            //negative value for L
195//                      endif
196                else
197                        //normal case
198                        offset = V_getDet_LateralOffset(fname,detStr)
199                        offset *= 10 //convert cm to mm
200                endif
201               
202        // calculation is in mm, not cm
203        // offset will be a negative value for the L panel, and positive for the R panel
204                if(kBCTR_CM)
205                        if(cmpstr("L",detStr[1]) == 0)
206//                              data_realDistX[][] = offset - (dimX - p)*tube_width                     // TODO should this be dimX-1-p = 47-p?
207                                data_realDistX[][] = offset - (dimX - p - 1/2)*tube_width - gap/2               // TODO should this be dimX-1-p = 47-p?
208                        else
209                        //      right
210//                              data_realDistX[][] = tube_width*(p+1) + offset + gap            //add to the Right det,
211                                data_realDistX[][] = tube_width*(p+1/2) + offset + gap/2                //add to the Right det
212                        endif
213                else
214                        data_realDistX[][] = tube_width*(p)
215                endif
216                data_realDistY[][] = coefW[0][p] + coefW[1][p]*q + coefW[2][p]*q*q
217       
218       
219        elseif(cmpstr(orientation,"horizontal")==0)
220                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
221                data_realDistY[][] = tube_width*q
222
223                if(cmpstr(fname,"VCALC")== 0 )
224                        offset = VCALC_getPanelTranslation(detStr)
225                        offset *= 10                    // convert to units of mm
226//                      if(cmpstr("B",detStr[1]) == 0)
227//                              offset *= -1    // negative value for Bottom det
228//                      endif
229                else
230                        //normal case
231                        offset = V_getDet_VerticalOffset(fname,detStr)
232                        offset *= 10 //convert cm to mm
233                endif
234               
235                if(kBCTR_CM)
236                        if(cmpstr("T",detStr[1]) == 0)
237//                              data_realDistY[][] = tube_width*(q+1) + offset + gap                   
238                                data_realDistY[][] = tube_width*(q+1/2) + offset + gap/2                       
239                        else
240                                // bottom
241//                              data_realDistY[][] = offset - (dimY - q)*tube_width     // TODO should this be dimY-1-q = 47-q?
242                                data_realDistY[][] = offset - (dimY - q - 1/2)*tube_width - gap/2       // TODO should this be dimY-1-q = 47-q?
243                        endif
244                else
245                        data_realDistY[][] = tube_width*(q)
246                endif
247                data_realDistX[][] = coefW[0][q] + coefW[1][q]*p + coefW[2][q]*p*p
248
249        else           
250                DoAlert 0,"Orientation not correctly passed in NonLinearCorrection(). No correction done."
251                return(0)
252        endif
253       
254        return(0)
255end
256
257
258
259
260// TODO:
261// -- the cal_x and y coefficients are totally fake
262// -- the wave assignment may not be correct.. so beware
263//
264//
265Function V_NonLinearCorrection_B(folder,dataW,cal_x,cal_y,detStr,destPath)
266        String folder
267        Wave dataW,cal_x,cal_y
268        String detStr,destPath
269
270        if(cmpstr(detStr,"B") != 0)
271                return(0)
272        endif
273
274Print "***Cal_X and Cal_Y for Back are using file values ***"
275
276//              cal_x[0] = VCALC_getPixSizeX(detStr)*10                 // pixel size in mm  VCALC_getPixSizeX(detStr) is [cm]
277//              cal_x[1] = 1
278//              cal_x[2] = 10000
279//              cal_y[0] = VCALC_getPixSizeY(detStr)*10                 // pixel size in mm  VCALC_getPixSizeX(detStr) is [cm]
280//              cal_y[1] = 1
281//              cal_y[2] = 10000
282
283       
284        // do I count on the orientation as an input, or do I just figure it out on my own?
285        Variable dimX,dimY
286       
287//      Wave dataW = V_getDetectorDataW(folder,detStr)
288       
289        dimX = DimSize(dataW,0)
290        dimY = DimSize(dataW,1)
291
292        // make a wave of the same dimensions, in the same data folder for the distance
293        // ?? or a 3D wave?
294        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
295        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
296        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
297        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
298       
299//      Wave cal_x = V_getDet_cal_x(folder,detStr)
300//      Wave cal_y = V_getDet_cal_y(folder,detStr)
301       
302        data_realDistX[][] = cal_x[0]*p*10              // cal_x and cal_y are in [cm], need mm
303        data_realDistY[][] = cal_y[0]*q*10
304       
305        return(0)
306end
307
308
309//
310//
311// TODO
312// -- VERIFY the calculations
313// -- verify where this needs to be done (if the beam center is changed)
314// -- then the q-calculation needs to be re-done
315// -- the position along the tube length is referenced to tube[0], for no particular reason
316//    It may be better to take an average? but [0] is an ASSUMPTION
317// -- distance along tube is simple interpolation, or do I use the coefficients to
318//    calculate the actual value
319//
320// -- distance in the lateral direction is based on tube width, which is a fixed parameter
321//
322//
323Function V_ConvertBeamCtrPix_to_mm(folder,detStr,destPath)
324        String folder,detStr,destPath
325       
326        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
327        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")     
328
329        String orientation
330        Variable dimX,dimY,xCtr,yCtr
331        dimX = DimSize(data_realDistX,0)
332        dimY = DimSize(data_realDistX,1)
333        if(dimX > dimY)
334                orientation = "horizontal"
335        else
336                orientation = "vertical"
337        endif
338       
339        xCtr = V_getDet_beam_center_x(folder,detStr)
340        yCtr = V_getDet_beam_center_y(folder,detStr)   
341       
342        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
343        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
344        WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
345        WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
346
347        Variable tube_width = V_getDet_tubeWidth(folder,detStr)         //this is also in mm
348
349//
350        strswitch(detStr)       // string switch
351                case "FL":
352                case "ML":
353                        // for Left/Right
354                        // for left
355                        x_mm[0] = data_realDistX[dimX-1][0] + (xCtr-dimX-1)*tube_width
356                        y_mm[0] = data_realDistY[0][yCtr]
357       
358                        break           
359                case "FR":     
360                case "MR":
361                        // for Left/Right
362                        // for right
363                        x_mm[0] = data_realDistX[0][0] + xCtr*tube_width
364                        y_mm[0] = data_realDistY[0][yCtr]
365                       
366                        break
367                case "FT":     
368                case "MT":
369                        // for Top                     
370                        x_mm[0] = data_realDistX[xCtr][0]
371                        y_mm[0] = data_realDistY[0][0] + yCtr*tube_width
372                       
373                        break           
374                case "FB":     
375                case "MB":
376                        // for Bottom                   
377                        x_mm[0] = data_realDistX[xCtr][0]
378                        y_mm[0] = data_realDistY[0][dimY-1] + (yCtr-dimY-1)*tube_width
379                                               
380                        break
381                default:                        // optional default expression executed
382                        Print "No case matched in V_Convert_FittedPix_2_cm"
383                        return(1)
384        endswitch
385               
386        return(0)
387end
388
389//
390//
391// (DONE)
392// x- VERIFY the calculations
393// x- verify where this needs to be done (if the beam center is changed)
394// x- then the q-calculation needs to be re-done
395// x- the position along the tube length is referenced to tube[0], for no particular reason
396//    It may be better to take an average? but [0] is an ASSUMPTION
397// x- distance along tube is simple interpolation
398//
399// x- distance in the lateral direction is based on tube width, which is a fixed parameter
400//
401// the value in pixels is written to the local data folder, NOT to disk (it is recalculated as needed)
402//
403Function V_ConvertBeamCtr_to_pix(folder,detStr,destPath)
404        String folder,detStr,destPath
405       
406        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
407        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")     
408
409        String orientation
410        Variable dimX,dimY,xCtr,yCtr
411        dimX = DimSize(data_realDistX,0)
412        dimY = DimSize(data_realDistX,1)
413        if(dimX > dimY)
414                orientation = "horizontal"
415        else
416                orientation = "vertical"
417        endif
418       
419        xCtr = V_getDet_beam_center_x(folder,detStr)            //these are in cm
420        yCtr = V_getDet_beam_center_y(folder,detStr)   
421       
422        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
423        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
424        WAVE x_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
425        WAVE y_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
426
427        Variable tube_width = V_getDet_tubeWidth(folder,detStr)
428
429        variable edge,delta
430        Variable gap
431
432// kPanelTouchingGap is in mm   
433// the gap is split equally between the panel pairs
434// TODO -- replace all of this with V_getDet_panel_gap(fname,detStr) once it is added to the file
435// these hard-wired values were determined from 6A and WB beam centers. LR values were exactly the same for
436// both beam considitions (+/- 0.0 mm). FTB was +/- 0.8 mm, MTB +/- 2 mm
437
438        gap = V_getDet_panel_gap(folder,detStr)
439
440// TODO:
441// -- once the gap fields have been verified, this check can be removed
442// -- it should only apply to data pre-2018 when the field did not exist in the file
443// -- any VSANS data from 2018+ should read gap from the file.
444
445        if(gap < -100)          //-999999 returned if field is missing from file
446                if(cmpstr(detStr,"FL") == 0 || cmpstr(detStr,"FR") == 0)
447                        gap = 3.5               //mm (measured, JB 1/4/18)
448                endif
449                if(cmpstr(detStr,"FT") == 0 || cmpstr(detStr,"FB") == 0)
450                        gap = 3.3               //mm (measured, JB 2/1/18)
451                endif
452                if(cmpstr(detStr,"ML") == 0 || cmpstr(detStr,"MR") == 0)
453                        gap = 5.9               //mm (measured, JB 1/4/18)
454                endif
455                if(cmpstr(detStr,"MT") == 0 || cmpstr(detStr,"MB") == 0)
456                        gap = 18.3              //mm (measured, JB 2/1/18)
457                endif
458        endif
459
460//
461        if(cmpstr(orientation,"vertical")==0)
462                //      this is data dimensioned as (Ntubes,Npix)
463
464                if(kBCTR_CM)
465                        if(cmpstr("L",detStr[1]) == 0)
466                                Make/O/D/N=(dimX) tmpTube
467                                tmpTube = data_RealDistX[p][0]
468                                FindLevel/P/Q tmpTube xCtr*10
469                                if(V_Flag)
470                                        edge = data_realDistX[47][0]            //tube 47
471                                        delta = abs(xCtr*10 - edge)
472                                        x_pix[0] = dimX-1 + delta/tube_width
473                                else
474                                        // beam center is on the panel, report the pixel value
475                                        x_pix[0] = V_LevelX
476                                endif
477                               
478                        else
479                        // R panel
480                                Make/O/D/N=(dimX) tmpTube
481                                tmpTube = data_RealDistX[p][0]
482                                FindLevel/P/Q tmpTube xCtr*10
483                                if(V_Flag)
484                                        //level not found
485                                        edge = data_realDistX[0][0]
486                                        delta = abs(xCtr*10 - edge + gap)               // how far past the edge of the panel
487                                        x_pix[0] = -delta/tube_width            //since the left edge of the R panel is pixel 0
488                                else
489                                        // beam center is on the panel, report the pixel value
490                                        x_pix[0] = V_LevelX
491                                endif
492                               
493                        endif
494
495                endif
496
497// the y-center will be on the panel in this direction
498                Make/O/D/N=(dimY) tmpTube
499                tmpTube = data_RealDistY[0][p]
500                FindLevel /P/Q tmpTube, yCtr*10
501               
502                y_pix[0] = V_levelX
503                KillWaves/Z tmpTube
504//              Print x_pix[0],y_pix[0]
505               
506        else
507                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
508
509                if(kBCTR_CM)
510                        if(cmpstr("T",detStr[1]) == 0)
511                                Make/O/D/N=(dimY) tmpTube
512                                tmpTube = data_RealDistY[p][0]
513                                FindLevel/P/Q tmpTube yCtr*10
514                                if(V_Flag)
515                                        edge = data_realDistY[0][0]             //tube 0
516                                        delta = abs(yCtr*10 - edge + gap)
517                                        y_pix[0] =  -delta/tube_width           //since the bottom edge of the T panel is pixel 0                               
518                                else
519                                        y_pix[0] = V_LevelX
520                                endif                           
521                               
522                        else
523                        // FM(B) panel
524                                Make/O/D/N=(dimY) tmpTube
525                                tmpTube = data_RealDistY[p][0]
526                                FindLevel/P/Q tmpTube yCtr*10
527                                if(V_Flag)
528                                        edge = data_realDistY[0][47]            //y tube 47
529                                        delta = abs(yCtr*10 - edge)
530                                        y_pix[0] = dimY-1 + delta/tube_width            //since the top edge of the B panels is pixel 47               
531                                else
532                                        y_pix[0] = V_LevelX
533                                endif
534
535                        endif
536                endif
537
538// the x-center will be on the panel in this direction         
539                Make/O/D/N=(dimX) tmpTube
540                tmpTube = data_RealDistX[p][0]
541                FindLevel /P/Q tmpTube, xCtr*10
542               
543                x_pix[0] = V_levelX
544                KillWaves/Z tmpTube
545               
546        endif
547               
548        return(0)
549end
550
551// converts from [cm] beam center to pixels
552//
553// the value in pixels is written to the local data folder, NOT to disk (it is recalculated as needed)
554//
555Function V_ConvertBeamCtr_to_pixB(folder,detStr,destPath)
556        String folder,detStr,destPath
557       
558        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
559        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")     
560
561        Variable dimX,dimY,xCtr,yCtr
562        dimX = DimSize(data_realDistX,0)
563        dimY = DimSize(data_realDistX,1)
564       
565        xCtr = V_getDet_beam_center_x(folder,detStr)                    //these are in cm, *10 to get mm
566        yCtr = V_getDet_beam_center_y(folder,detStr)   
567       
568        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
569        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
570        WAVE x_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
571        WAVE y_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
572
573
574// simple wave lookup
575// can't use x_pix[0] = data_RealDistX(xCtr)[0] since the data has no x-scale and (xCtr) is interpreted
576// as a point value
577
578//
579//xCtr, yCtr are in cm, *10 to get mm to compare to distance array
580
581        Make/O/D/N=(dimX) tmpTube
582        tmpTube = data_RealDistX[p][0]
583        FindLevel /P/Q tmpTube, xCtr*10
584       
585        x_pix[0] = V_levelX
586        KillWaves/Z tmpTube
587       
588       
589        Make/O/D/N=(dimY) tmpTube
590        tmpTube = data_RealDistY[0][p]
591        FindLevel /P/Q tmpTube, yCtr*10
592       
593        y_pix[0] = V_levelX
594        KillWaves/Z tmpTube
595               
596        print "pixel ctr B = ",x_pix[0],y_pix[0]
597               
598        return(0)
599end
600
601//
602//
603// TODO
604// -- VERIFY the calculations
605// -- verify where this needs to be done (if the beam center is changed)
606// -- then the q-calculation needs to be re-done
607//
608// -- not much is known about the "B" detector, so this
609//    all hinges on the non-linear corrections being done correctly for that detector
610//
611//      Variable detCtrX, detCtrY
612//      // get the pixel center of the detector (not the beam center)
613//      detCtrX = trunc( DimSize(dataW,0)/2 )           //
614//      detCtrY = trunc( DimSize(dataW,1)/2 )
615//
616//
617Function V_ConvertBeamCtrPix_to_mmB(folder,detStr,destPath)
618        String folder,detStr,destPath
619       
620       
621//      DoAlert 0,"Error - Beam center is being interpreted as pixels, but needs to be in cm. V_ConvertBeamCtrPix_to_mmB()"
622       
623        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
624        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")     
625       
626        Variable xCtr,yCtr
627        xCtr = V_getDet_beam_center_x(folder,detStr)
628        yCtr = V_getDet_beam_center_y(folder,detStr)   
629       
630        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
631        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
632        WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
633        WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
634
635        x_mm[0] = data_realDistX[xCtr][0]
636        y_mm[0] = data_realDistY[0][yCtr]
637               
638        return(0)
639end
640
641
642
643
644
645
646/////
647//
648// non-linear corrections to the tube pixels
649// - returns the distance in mm (although this may change)
650//
651// c0,c1,c2,pix
652// c0-c2 are the fit coefficients
653// pix is the test pixel
654//
655// returns the distance in mm (relative to ctr pixel)
656// ctr is the center pixel, as defined when fitting to quadratic was done
657//
658Function V_TubePixel_to_mm(c0,c1,c2,pix)
659        Variable c0,c1,c2,pix
660       
661        Variable dist
662        dist = c0 + c1*pix + c2*pix*pix
663       
664        return(dist)
665End
666//
667////
668
669
670//
671// TESTING ONLY
672Proc V_MakeFakeCalibrationWaves()
673        // make these in the RAW data folder, before converting to a work folder
674        // - then they will be "found" by get()
675        // -- only for the tube, not the Back det
676       
677//      DoAlert 0, "re-do this and do a better job of filling the fake calibration data"
678
679        DoAlert 0, "Calibration waves are read in from the data file"
680       
681//      V_fMakeFakeCalibrationWaves()
682End
683
684
685
686//
687// TESTING ONLY
688//
689// orientation does not matter, there are 48 tubes in each bank
690// so dimension (3,48) for everything.
691//
692// -- but the orientation does indicate TB vs LR, which has implications for
693//  the (fictional) dimension of the pixel along the tube axis, at least as far
694// as for making the fake coefficients.
695//
696Function V_fMakeFakeCalibrationWaves()
697
698        Variable ii,pixSize
699        String detStr,fname="RAW",orientation
700       
701        for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
702                detStr = StringFromList(ii, ksDetectorListNoB, ";")
703//              Wave w = V_getDetectorDataW(fname,detStr)
704                Make/O/D/N=(3,48) $("root:Packages:NIST:VSANS:RAW:entry:instrument:detector_"+detStr+":spatial_calibration")
705                Wave calib = $("root:Packages:NIST:VSANS:RAW:entry:instrument:detector_"+detStr+":spatial_calibration")
706                // !!!! this overwrites what is there
707
708                orientation = V_getDet_tubeOrientation(fname,detStr)
709                if(cmpstr(orientation,"vertical")==0)
710                //      this is vertical tube data dimensioned as (Ntubes,Npix)
711                        pixSize = 8.4           //V_getDet_y_pixel_size(fname,detStr)
712                       
713                elseif(cmpstr(orientation,"horizontal")==0)
714                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
715                        pixSize = 4                     //V_getDet_x_pixel_size(fname,detStr)
716                       
717                else           
718                        DoAlert 0,"Orientation not correctly passed in NonLinearCorrection(). No correction done."
719                endif
720               
721                calib[0][] = -(128/2)*pixSize                   //approx (n/2)*pixSixe
722                calib[1][] = pixSize
723                calib[2][] = 2e-4
724               
725        endfor
726       
727        // now fake calibration for "B"
728        Wave cal_x = V_getDet_cal_x("RAW","B")
729        Wave cal_y = V_getDet_cal_y("RAW","B")
730       
731        cal_x = .34             // mm, ignore the other 2 values
732        cal_y = .34             // mm
733        return(0)
734End
735
736//
737// (DONE)
738// x- MUST VERIFY the definition of SDD and how (if) setback is written to the data files
739// x- currently I'm assuming that the SDD is the "nominal" value which is correct for the
740//    L/R panels, but is not correct for the T/B panels (must add in the setback)
741//
742//
743//
744// data_realDistX, Y must be previously generated from running NonLinearCorrection()
745//
746// call with:
747// fname as the WORK folder, "RAW"
748// detStr = detector String, "FL"
749// destPath = path to destination WORK folder ("root:Packages:NIST:VSANS:"+folder)
750//
751Function V_Detector_CalcQVals(fname,detStr,destPath)
752        String fname,detStr,destPath
753
754        String orientation
755        Variable xCtr,yCtr,lambda,sdd
756       
757// get all of the geometry information 
758        orientation = V_getDet_tubeOrientation(fname,detStr)
759
760
761        sdd = V_getDet_ActualDistance(fname,detStr)             //sdd derived, including setback [cm]
762
763        // this is the ctr in pixels --xx-- (now it is in cm!)
764//      xCtr = V_getDet_beam_center_x(fname,detStr)
765//      yCtr = V_getDet_beam_center_y(fname,detStr)
766        // this is ctr in mm
767        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
768        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
769        lambda = V_getWavelength(fname)
770        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
771        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
772
773// make the new waves
774        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
775        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
776        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
777        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
778        Wave qTot = $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
779        Wave qx = $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
780        Wave qy = $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
781        Wave qz = $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
782
783// calculate all of the q-values
784// sdd is passed in [cm]
785
786
787// after adding in the 680x1656 back detector, load time was 7.8s, without multithreading
788// with multithreading, 1.9s
789//       qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
790//              qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
791//              qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
792//              qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)   
793
794        MultiThread qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
795        MultiThread     qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
796        MultiThread     qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
797        MultiThread     qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
798       
799        return(0)
800End
801
802
803//function to calculate the overall q-value, given all of the necesary trig inputs
804//
805// (DONE)
806// x- verify the calculation (accuracy - in all input conditions)
807// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
808// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
809//    to each pixel
810//
811//sdd is in [cm]
812// distX and distY are in [mm]
813//wavelength is in Angstroms
814//
815//returned magnitude of Q is in 1/Angstroms
816//
817ThreadSafe Function V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
818        Variable xaxval,yaxval,xctr,yctr,sdd,lam
819        Wave distX,distY
820       
821        Variable dx,dy,qval,two_theta,dist
822               
823
824        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
825        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
826        dist = sqrt(dx^2 + dy^2)
827       
828        dist /= 10  // convert mm to cm
829       
830        two_theta = atan(dist/sdd)
831
832        qval = 4*Pi/lam*sin(two_theta/2)
833       
834        return qval
835End
836
837//calculates just the q-value in the x-direction on the detector
838// (DONE)
839// x- verify the calculation (accuracy - in all input conditions)
840// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
841// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
842//    to each pixel
843//
844//
845// this properly accounts for qz
846//
847ThreadSafe Function V_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
848        Variable xaxval,yaxval,xctr,yctr,sdd,lam
849        Wave distX,distY
850
851        Variable qx,qval,phi,dx,dy,dist,two_theta
852       
853        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
854       
855
856        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
857        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
858        phi = V_FindPhi(dx,dy)
859       
860        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
861        dist = sqrt(dx^2 + dy^2)
862        dist /= 10  // convert mm to cm
863
864        two_theta = atan(dist/sdd)
865
866        qx = qval*cos(two_theta/2)*cos(phi)
867       
868        return qx
869End
870
871//calculates just the q-value in the y-direction on the detector
872// (DONE)
873// x- verify the calculation (accuracy - in all input conditions)
874// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
875// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
876//    to each pixel
877//
878//
879// this properly accounts for qz
880//
881ThreadSafe Function V_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
882        Variable xaxval,yaxval,xctr,yctr,sdd,lam
883        Wave distX,distY
884
885        Variable qy,qval,phi,dx,dy,dist,two_theta
886       
887        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
888       
889
890        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
891        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
892        phi = V_FindPhi(dx,dy)
893       
894        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
895        dist = sqrt(dx^2 + dy^2)
896        dist /= 10  // convert mm to cm
897
898        two_theta = atan(dist/sdd)
899
900        qy = qval*cos(two_theta/2)*sin(phi)
901       
902        return qy
903End
904
905//calculates just the q-value in the z-direction on the detector
906// (DONE)
907// x- verify the calculation (accuracy - in all input conditions)
908// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
909// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
910//    to each pixel
911//
912// not actually used for any calculations, but here for completeness if anyone asks, or for 2D data export
913//
914// this properly accounts for qz, because it is qz
915//
916ThreadSafe Function V_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
917        Variable xaxval,yaxval,xctr,yctr,sdd,lam
918        Wave distX,distY
919
920        Variable qz,qval,phi,dx,dy,dist,two_theta
921       
922        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
923       
924
925        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
926        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
927       
928        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
929        dist = sqrt(dx^2 + dy^2)
930        dist /= 10  // convert mm to cm
931
932        two_theta = atan(dist/sdd)
933
934        qz = qval*sin(two_theta/2)
935       
936        return qz
937End
938
939
940//
941// (DONE)
942// x- VERIFY calculations
943// x- This is the actual solid angle per pixel, not a ratio vs. some "unit SA"
944//    Do I just correct for the different area vs. the "nominal" central area?
945// x- decide how to implement - YES - directly change the data values (as was done in the past)
946//    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
947//    would need this to be applied before exporting)
948// x- do I keep a wave note indicating that this correction has been applied to the data
949//    so that it can be "un-applied"? NO
950// x- do I calculate theta from geometry directly, or get it from Q (Assuming it's present?)
951//    (YES just from geometry, since I need SDD and dx and dy values...)
952//
953//
954Function V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
955        Wave w,w_err
956        String fname,detStr,destPath
957
958        Variable sdd,xCtr,yCtr,lambda
959
960// get all of the geometry information 
961//      orientation = V_getDet_tubeOrientation(fname,detStr)
962        sdd = V_getDet_ActualDistance(fname,detStr)
963
964
965        // this is ctr in mm
966        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
967        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
968        lambda = V_getWavelength(fname)
969       
970        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
971       
972        Wave data_realDistX = data_realDistX
973        Wave data_realDistY = data_realDistY
974
975        Duplicate/O w solid_angle,tmp_theta,tmp_dist            //in the current df
976
977//// calculate the scattering angle
978//      dx = (distX - xctr)             //delta x in mm
979//      dy = (distY - yctr)             //delta y in mm
980        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
981       
982        tmp_dist /= 10  // convert mm to cm
983        // sdd is in [cm]
984
985        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle
986
987        Variable ii,jj,numx,numy,dx,dy
988        numx = DimSize(tmp_theta,0)
989        numy = DimSize(tmp_theta,1)
990       
991        for(ii=0        ;ii<numx;ii+=1)
992                for(jj=0;jj<numy;jj+=1)
993                       
994                        if(ii==0)               //do a forward difference if ii==0
995                                dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj])        //delta x for the pixel
996                        else
997                                dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj])        //delta x for the pixel
998                        endif
999                       
1000                       
1001                        if(jj==0)
1002                                dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj])        //delta y for the pixel
1003                        else
1004                                dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1])        //delta y for the pixel
1005                        endif
1006       
1007                        dx /= 10
1008                        dy /= 10                // convert mm to cm (since sdd is in cm)
1009                        solid_angle[ii][jj] = dx*dy             //this is in cm^2
1010                endfor
1011        endfor
1012       
1013        // to cover up any issues w/negative dx or dy
1014        solid_angle = abs(solid_angle)
1015       
1016        // solid_angle correction
1017        // == dx*dy*cos^3/sdd^2
1018        solid_angle *= (cos(tmp_theta))^3
1019        solid_angle /= sdd^2
1020       
1021        // Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!)
1022        w /= solid_angle
1023       
1024        //
1025        // correctly apply the correction to the error wave (assume a perfect value?)
1026        w_err /= solid_angle            //
1027
1028// DONE x- clean up after I'm satisfied computations are correct               
1029        KillWaves/Z tmp_theta,tmp_dist
1030       
1031        return(0)
1032end
1033
1034
1035
1036
1037//
1038// Large angle transmission correction
1039//
1040// DIVIDE the intensity by this correction to get the right answer
1041//
1042//
1043// Apply the large angle transmssion correction as the data is converted to WORK
1044// so that whether the data is saved as 2D or 1D, the correction has properly been done.
1045//
1046// This is, however, a SAMPLE dependent calculation, not purely instrument geometry.
1047//
1048Function V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
1049        Wave w,w_err
1050        String fname,detStr,destPath
1051
1052        Variable sdd,xCtr,yCtr,trans,trans_err,uval
1053
1054// get all of the geometry information 
1055//      orientation = V_getDet_tubeOrientation(fname,detStr)
1056        sdd = V_getDet_ActualDistance(fname,detStr)
1057
1058        // this is ctr in mm
1059        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
1060        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
1061        trans = V_getSampleTransmission(fname)
1062        trans_err = V_getSampleTransError(fname)
1063       
1064        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
1065       
1066        Wave data_realDistX = data_realDistX
1067        Wave data_realDistY = data_realDistY
1068
1069        Duplicate/O w lat_corr,tmp_theta,tmp_dist,lat_err,tmp_err               //in the current df
1070
1071//// calculate the scattering angle
1072//      dx = (distX - xctr)             //delta x in mm
1073//      dy = (distY - yctr)             //delta y in mm
1074        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
1075       
1076        tmp_dist /= 10  // convert mm to cm
1077        // sdd is in [cm]
1078
1079        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle
1080
1081        Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp
1082        numx = DimSize(tmp_theta,0)
1083        numy = DimSize(tmp_theta,1)
1084       
1085       
1086        //optical thickness
1087        uval = -ln(trans)               //use natural logarithm
1088       
1089        for(ii=0        ;ii<numx;ii+=1)
1090                for(jj=0;jj<numy;jj+=1)
1091                       
1092                        cos_th = cos(tmp_theta[ii][jj])
1093                        arg = (1-cos_th)/cos_th
1094                       
1095                        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
1096                        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120
1097                        // OR
1098                        if((uval<0.01) || (cos_th>0.99))       
1099                                //small arg, approx correction
1100                                lat_corr[ii][jj] = 1-0.5*uval*arg
1101                        else
1102                                //large arg, exact correction
1103                                lat_corr[ii][jj] = (1-exp(-uval*arg))/(uval*arg)
1104                        endif
1105                         
1106                        // (DONE)
1107                        // x- properly calculate and apply the 2D error propagation
1108                        if(trans == 1)
1109                                lat_err[ii][jj] = 0             //no correction, no error
1110                        else
1111                                //sigT, calculated from the Taylor expansion
1112                                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3)
1113                                tmp *= tmp
1114                                tmp *= trans_err^2
1115                                tmp = sqrt(tmp)         //sigT
1116                               
1117                                lat_err[ii][jj] = tmp
1118                        endif
1119                         
1120 
1121                endfor
1122        endfor
1123       
1124
1125       
1126        // Here it is! Apply the correction to the intensity (divide -- to get the proper correction)
1127        w /= lat_corr
1128
1129        // relative errors add in quadrature to the current 2D error
1130        tmp_err = (w_err/lat_corr)^2 + (lat_err/lat_corr)^2*w*w/lat_corr^2
1131        tmp_err = sqrt(tmp_err)
1132       
1133        w_err = tmp_err
1134       
1135
1136        // DONE x- clean up after I'm satisfied computations are correct               
1137        KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err
1138       
1139        return(0)
1140end
1141
1142
1143
1144//
1145//test procedure, not called anymore
1146Proc V_AbsoluteScaling(type,c0,c1,c2,c3,c4,c5,I_err)
1147        String type
1148        Variable c0=1,c1=0.1,c2=0.95,c3=0.1,c4=1,c5=32.0,I_err=0.32
1149        Prompt type,"WORK data type",popup,"CAL;COR;SAM"
1150        Prompt c0, "Sample Transmission"
1151        Prompt c1, "Sample Thickness (cm)"
1152        Prompt c2, "Standard Transmission"
1153        Prompt c3, "Standard Thickness (cm)"
1154        Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)"
1155        Prompt c5, "Standard Cross-Section (cm-1)"
1156        Prompt I_err, "error in I(q=0) (one std dev)"
1157
1158        Variable err
1159        //call the function to do the math
1160        //data from "type" will be scaled and deposited in ABS
1161        err = V_Absolute_Scale(type,c0,c1,c2,c3,c4,c5,I_err)
1162       
1163        if(err)
1164                Abort "Error in V_Absolute_Scale()"
1165        endif
1166       
1167        //contents are always dumped to ABS
1168        type = "ABS"
1169       
1170        //need to update the display with "data" from the correct dataFolder
1171        //reset the current display type to "type"
1172        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
1173        gCurDispType = Type     
1174       
1175        V_UpdateDisplayInformation(Type)
1176       
1177End
1178
1179//
1180//
1181// kappa comes in as s_izero, so be sure to use 1/kappa_err
1182//
1183//convert the "type" data to absolute scale using the given standard information
1184//s_ is the standard
1185//w_ is the "work" file
1186//both are work files and should already be normalized to 10^8 monitor counts
1187Function V_Absolute_Scale(type,absStr)
1188        String type,absStr
1189       
1190       
1191        Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
1192
1193        Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
1194        Variable scale,trans_err
1195        Variable err,ii
1196        String detStr
1197       
1198        // be sure that the starting data exists
1199        err = V_WorkDataExists(type)
1200        if(err==1)
1201                return(err)
1202        endif
1203               
1204        //copy from current dir (type) to ABS
1205        V_CopyHDFToWorkFolder(type,"ABS")       
1206
1207// TODO: -- which monitor to use? Here, I think it should already be normalized to 10^8
1208//     
1209//      w_moncount = V_getMonitorCount(type)            //monitor count in "type"
1210       
1211        w_moncount = V_getBeamMonNormData(type)
1212               
1213        if(w_moncount == 0)
1214                //zero monitor counts will give divide by zero ---
1215                DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
1216                Return(1)               //report error
1217        Endif
1218
1219        w_trans = V_getSampleTransmission(type)         //sample transmission
1220        w_thick = V_getSampleThickness(type)            //sample thickness
1221        trans_err = V_getSampleTransError(type)
1222       
1223       
1224        //get the parames from the list
1225        s_trans = NumberByKey("TSTAND", absStr, "=", ";")       //parse the list of values
1226        s_thick = NumberByKey("DSTAND", absStr, "=", ";")
1227        s_izero = NumberByKey("IZERO", absStr, "=", ";")
1228        s_cross = NumberByKey("XSECT", absStr, "=", ";")
1229        kappa_err = NumberByKey("SDEV", absStr, "=", ";")
1230
1231       
1232        //calculate scale factor
1233        s1 = defmon/w_moncount          // monitor count (s1 should be 1)
1234        s2 = s_thick/w_thick
1235        s3 = s_trans/w_trans
1236        s4 = s_cross/s_izero
1237        scale = s1*s2*s3*s4
1238
1239       
1240        // kappa comes in as s_izero, so be sure to use 1/kappa_err
1241
1242        // and now loop through all of the detectors
1243        //do the actual absolute scaling here, modifying the data in ABS
1244        for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
1245                detStr = StringFromList(ii, ksDetectorListNoB, ";")
1246                Wave data = V_getDetectorDataW("ABS",detStr)
1247                Wave data_err = V_getDetectorDataErrW("ABS",detStr)
1248               
1249                data *= scale
1250                data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1251        endfor
1252
1253        // do the back detector separately, if it is set to be used
1254        NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
1255        if(gIgnoreDetB == 0)
1256                detStr = "B"
1257                Wave data = V_getDetectorDataW("ABS",detStr)
1258                Wave data_err = V_getDetectorDataErrW("ABS",detStr)
1259               
1260                //get the parames from the list
1261                s_trans = NumberByKey("TSTAND_B", absStr, "=", ";")     //parse the list of values
1262                s_thick = NumberByKey("DSTAND_B", absStr, "=", ";")
1263                s_izero = NumberByKey("IZERO_B", absStr, "=", ";")
1264                s_cross = NumberByKey("XSECT_B", absStr, "=", ";")
1265                kappa_err = NumberByKey("SDEV_B", absStr, "=", ";")
1266
1267                //calculate scale factor
1268                s1 = defmon/w_moncount          // monitor count (s1 should be 1)
1269                s2 = s_thick/w_thick
1270                s3 = s_trans/w_trans
1271                s4 = s_cross/s_izero
1272                scale = s1*s2*s3*s4
1273               
1274                data *= scale
1275                data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1276        endif
1277       
1278        //********* 15APR02
1279        // DO NOT correct for atenuators here - the COR step already does this, putting all of the data on equal
1280        // footing (zero atten) before doing the subtraction.
1281       
1282        Return (0) //no error
1283End
1284
1285
1286//
1287// TODO:
1288//   --         DoAlert 0,"This has not yet been updated for VSANS"
1289//
1290//
1291// match the attenuation of the RAW data to the "type" data
1292// so that they can be properly added
1293//
1294// are the attenuator numbers the same? if so exit
1295//
1296// if not, find the attenuator number for type
1297// - find both attenuation factors
1298//
1299// rescale the raw data to match the ratio of the two attenuation factors
1300// -- adjust the detector count (rw)
1301// -- the linear data
1302//
1303//
1304Function V_Adjust_RAW_Attenuation(type)
1305        String type
1306
1307        DoAlert 0,"This has not yet been updated for VSANS"
1308       
1309        WAVE rw=$("root:Packages:NIST:RAW:realsread")
1310        WAVE linear_data=$("root:Packages:NIST:RAW:linear_data")
1311        WAVE data=$("root:Packages:NIST:RAW:data")
1312        WAVE data_err=$("root:Packages:NIST:RAW:linear_data_error")
1313        WAVE/T tw = $("root:Packages:NIST:RAW:textRead")
1314       
1315        WAVE dest_reals=$("root:Packages:NIST:"+type+":realsread")
1316
1317        Variable dest_atten,raw_atten,tol
1318        Variable lambda,raw_atten_err,raw_AttenFactor,dest_attenFactor,dest_atten_err
1319        String fileStr
1320
1321        dest_atten = dest_reals[3]
1322        raw_atten = rw[3]
1323       
1324        tol = 0.1               // within 0.1 atten units is OK
1325        if(abs(dest_atten - raw_atten) < tol )
1326                return(0)
1327        endif
1328
1329        fileStr = tw[3]
1330        lambda = rw[26]
1331        // TODO access correct values
1332        raw_AttenFactor = 1//AttenuationFactor(fileStr,lambda,raw_atten,raw_atten_err)
1333        dest_AttenFactor = 1//AttenuationFactor(fileStr,lambda,dest_atten,dest_atten_err)
1334               
1335        rw[2] *= dest_AttenFactor/raw_AttenFactor
1336        linear_data *= dest_AttenFactor/raw_AttenFactor
1337       
1338        // to keep "data" and linear_data in sync
1339        data = linear_data
1340       
1341        return(0)
1342End
1343
1344//
1345// testing procedure, called from a menu selection
1346//
1347Proc V_DIV_a_Workfile(type)
1348        String type
1349        Prompt type,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
1350       
1351        //macro will take whatever is in SELECTED folder and DIVide it by the current
1352        //contents of the DIV folder - the function will check for existence
1353        //before proceeding
1354
1355        Abort "This has not yet been updated for VSANS"
1356       
1357        Variable err
1358        err = V_DIVCorrection(type)             //returns err = 1 if data doesn't exist in specified folders
1359       
1360        if(err)
1361                Abort "error in V_DIVCorrection()"
1362        endif
1363       
1364        //contents are NOT always dumped to CAL, but are in the new type folder
1365       
1366        String newTitle = "WORK_"+type
1367        DoWindow/F VSANS_Data
1368        DoWindow/T VSANS_Data, newTitle
1369        KillStrings/Z newTitle
1370       
1371        //need to update the display with "data" from the correct dataFolder
1372        //reset the current displaytype to "type"
1373        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=Type
1374       
1375        V_UpdateDisplayInformation(type)
1376       
1377End
1378
1379
1380//
1381// TODO:
1382//   x-         DoAlert 0,"This has not yet been updated for VSANS"
1383//   -- how is the error propagation handled? Be sure it is calculated correctly when DIV is generated
1384//      and is applied correctly here...
1385//
1386//function will divide the contents of "workType" folder with the contents of
1387//the DIV folder + detStr
1388// all data is linear scale for the calculation
1389//
1390Function V_DIVCorrection(data,data_err,detStr,workType)
1391        Wave data,data_err
1392        String detStr,workType
1393       
1394        //check for existence of data in type and DIV
1395        // if the desired data doesn't exist, let the user know, and abort
1396        String destPath=""
1397
1398        NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
1399        if(cmpstr(detStr,"B")==0 && gIgnoreDetB)
1400                return(0)
1401        endif
1402
1403
1404        if(WaveExists(data) == 0)
1405                Print "The data wave does not exist in V_DIVCorrection()"
1406                Return(1)               //error condition
1407        Endif
1408       
1409        //check for DIV
1410        // if the DIV workfile doesn't exist, let the user know,and abort
1411        // !! be sure to check first, before trying to access the wave
1412       
1413//      WAVE/Z div_data = $("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")
1414        if(WaveExists($("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")) == 0)
1415                Print "The DIV wave does not exist in V_DIVCorrection()"
1416                Return(1)               //error condition
1417        Endif
1418        if(WaveExists($("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":linear_data_error")) == 0)
1419                Print "The DIV error wave does not exist in V_DIVCorrection()"
1420                Return(1)               //error condition
1421        Endif
1422        //files exist, proceed
1423
1424        WAVE/Z div_data_err = V_getDetectorDataErrW("DIV",detStr)
1425        WAVE/Z div_data = V_getDetectorDataW("DIV",detStr)
1426
1427
1428
1429// do the error propagation first, since data is changed by the correction
1430        data_err = sqrt(data_err^2/div_data^2 + div_data_err^2 * data^2/div_data^4 )
1431
1432// then the correction
1433        data /= div_data
1434
1435       
1436        Return(0)
1437End
1438
1439
1440//////////////////////////
1441// detector corrections to stitch the back detector into one proper image
1442//
1443//
1444//
1445
1446
1447//
1448// to register the image on the back detector panel
1449//
1450// middle portion (552 pix in Y) is held fixed
1451// top portion of image is shifted right and down
1452// bottom portion of image is shifted right and up
1453//
1454// remainder of image is filled with Zero (NaN causes problems converting to WORK)
1455//
1456// currently, data is not added together and averaged, but it could be
1457//
1458Function V_ShiftBackDetImage(w,adjW)
1459        Wave w,adjW
1460
1461        NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
1462
1463// this is necessary for some old data with the 150x150 back (dummy) panel
1464        NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
1465        if(gIgnoreDetB == 1)
1466                adjW=w
1467                return(0)
1468        endif
1469       
1470        adjW=0
1471               
1472        Variable topX,bottomX
1473        Variable topY,bottomY
1474        Variable totalY,ccdX,ccdY
1475       
1476//      topX = 7
1477//      topY = 105
1478       
1479//      bottomX = 5
1480//      bottomY = 35
1481
1482// TODOHIGHRES
1483// the detector pix dimensions are hard-wired, be sure the are correct
1484        switch(gHighResBinning)
1485                case 1:
1486                        topX = kShift_topX_bin1
1487                        topY = kShift_topY_bin1
1488                        bottomX = kShift_bottomX_bin1
1489                        bottomY = kShift_bottomY_bin1
1490                       
1491                        totalY = 6624   // total YDim
1492                        ccdY = 2208             // = YDim/3
1493                        ccdX = 2720             // = xDim
1494                        break
1495                case 4:
1496                        topX = kShift_topX_bin4
1497                        topY = kShift_topY_bin4
1498                        bottomX = kShift_bottomX_bin4
1499                        bottomY = kShift_bottomY_bin4
1500                       
1501                        totalY = 1656   // total YDim
1502                        ccdY = 552              // = YDim/3
1503                        ccdX = 680              // = xDim
1504
1505                       
1506                        break
1507                default:               
1508                        Abort "No binning case matches in V_ShiftBackDetImage"
1509                       
1510        endswitch
1511
1512                // middle
1513                adjW[][ccdY,ccdY+ccdY] = w[p][q]
1514       
1515                //top
1516                adjW[0+topX,ccdX-1][ccdY+ccdY,totalY-1-topY] = w[p-topX][q+topY]
1517               
1518                //bottom
1519                adjW[0+bottomX,ccdX-1][0+bottomY,ccdY-1] = w[p-bottomX][q-bottomY]
1520
1521       
1522        return(0)
1523End
1524
1525
1526Proc pV_MedianFilterBack(folder)
1527        String folder="RAW"
1528       
1529        V_MedianFilterBack(folder)
1530end
1531
1532Function V_MedianFilterBack(folder)
1533        String folder
1534
1535        Wave w = V_getDetectorDataW(folder,"B")
1536
1537        NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
1538        switch(gHighResBinning)
1539                case 1:
1540                        MatrixFilter /N=11 /P=1 median w                        //              /P=n flag sets the number of passes (default is 1 pass)
1541                       
1542                        Print "*** median noise filter 11x11 applied to the back detector (1 pass) ***"
1543                        break
1544                case 4:
1545                        MatrixFilter /N=3 /P=1 median w                 //              /P=n flag sets the number of passes (default is 1 pass)
1546                       
1547                        Print "*** median noise filter 3x3 applied to the back detector (1 pass) ***"
1548                        break
1549                default:
1550                        Abort "No binning case matches in V_MedianFilterBack"
1551        endswitch
1552
1553        return(0)
1554End
1555
1556
1557Proc pV_SubtractReadNoiseBack(folder,ReadNoise)
1558        String folder="RAW"
1559        Variable readNoise=3160
1560       
1561        V_SubtractReadNoiseBack(folder,readNoise)
1562end
1563
1564Function V_SubtractReadNoiseBack(folder,readNoise)
1565        String folder
1566        Variable readNoise
1567
1568                Wave w = V_getDetectorDataW(folder,"B")
1569                w -= readNoise          // a constant value
1570               
1571//              MatrixFilter /N=3 median w
1572//              Print "*** median noise filter applied to the back detector***"
1573       
1574        return(0)
1575End
1576
1577
1578Proc pV_MedianAndReadNoiseBack(folder,ReadNoise)
1579        String folder="RAW"
1580        Variable readNoise=3160
1581       
1582        V_MedianAndReadNoiseBack(folder,readNoise)
1583end
1584
1585Function V_MedianAndReadNoiseBack(folder,readNoise)
1586        String folder
1587        Variable readNoise
1588
1589                Wave w = V_getDetectorDataW(folder,"B")
1590                w -= readNoise          // a constant value
1591
1592                NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
1593                switch(gHighResBinning)
1594                        case 1:
1595                                MatrixFilter /N=11 /P=1 median w                        //              /P=n flag sets the number of passes (default is 1 pass)
1596                               
1597                                Print "*** median noise filter 11x11 applied to the back detector (1 pass) ***"
1598                                break
1599                        case 4:
1600                                MatrixFilter /N=3 /P=1 median w                 //              /P=n flag sets the number of passes (default is 1 pass)
1601                               
1602                                Print "*** median noise filter 3x3 applied to the back detector (1 pass) ***"
1603                                break
1604                        default:
1605                                Abort "No binning case matches in V_MedianAndReadNoiseBack"
1606                endswitch
1607               
1608        return(0)
1609End
1610
Note: See TracBrowser for help on using the repository browser.