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

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

a few changes to the VSANS menu

updated VSANS reduction documentation for posting on the trac page

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