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

Last change on this file since 1235 was 1235, checked in by srkline, 3 years ago

Changed the calculation of solid angle to properly reflect the tube geometry. This mode is now the (correct) default calculation. A global switch can be applied for testing to use the old method (with many warning alerts).

File size: 52.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                                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//              "B" is passed, so I need to check for "B" and for panel orientation
954// -if it is detector "B" (not tubes), then the normal solid angle correction applies
955// -if it is a tube panel, then I need to know the orientation, to know which angles
956//    and pixel dimensions to use
957//
958// *** UPDATED 1/2020 SRK
959// -using new calculation since the lateral direction of the tubes does not affect the solid angle
960// projection (see He (2015) and John's memo)
961//
962//
963Function V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
964        Wave w,w_err
965        String fname,detStr,destPath
966
967        Variable sdd,xCtr,yCtr,lambda
968        String orientation
969       
970// get all of the geometry information 
971        orientation = V_getDet_tubeOrientation(fname,detStr)
972        sdd = V_getDet_ActualDistance(fname,detStr)
973
974        // this is ctr in mm
975        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
976        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
977        lambda = V_getWavelength(fname)
978       
979        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
980       
981        Wave data_realDistX = data_realDistX
982        Wave data_realDistY = data_realDistY
983
984        Duplicate/O w solid_angle,tmp_theta,tmp_dist,tmp_theta_i        //in the current df
985
986//// calculate the scattering angle
987//      dx = (distX - xctr)             //delta x in mm
988//      dy = (distY - yctr)             //delta y in mm
989        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
990       
991        tmp_dist /= 10  // convert mm to cm
992        // sdd is in [cm]
993
994        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the (total) scattering angle
995
996        Variable ii,jj,numx,numy,dx,dy
997        numx = DimSize(tmp_theta,0)
998        numy = DimSize(tmp_theta,1)
999
1000        if(cmpstr(detStr,"B")==0)
1001                //detector B is a grid, straightforward cos^3 solid angle
1002                for(ii=0        ;ii<numx;ii+=1)
1003                        for(jj=0;jj<numy;jj+=1)
1004                               
1005                                if(ii==0)               //do a forward difference if ii==0
1006                                        dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj])        //delta x for the pixel
1007                                else
1008                                        dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj])        //delta x for the pixel
1009                                endif
1010                               
1011                               
1012                                if(jj==0)
1013                                        dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj])        //delta y for the pixel
1014                                else
1015                                        dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1])        //delta y for the pixel
1016                                endif
1017               
1018                                dx /= 10
1019                                dy /= 10                // convert mm to cm (since sdd is in cm)
1020                                solid_angle[ii][jj] = dx*dy             //this is in cm^2
1021                        endfor
1022                endfor
1023               
1024                // to cover up any issues w/negative dx or dy
1025                solid_angle = abs(solid_angle)
1026               
1027                // solid_angle correction
1028                // == dx*dy*cos^3/sdd^2
1029                solid_angle *= (cos(tmp_theta))^3
1030                solid_angle /= sdd^2
1031               
1032                // Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!)
1033                w /= solid_angle
1034
1035                // correctly apply the correction to the error wave (assume a perfect value?)
1036                w_err /= solid_angle            //
1037
1038        else
1039                //             
1040                //different calculation for the tubes, different calculation based on XY orientation
1041                //
1042                if(cmpstr(orientation,"vertical")==0)
1043                        // L/R panels, tube axis is y-direction
1044                        // this is now a different tmp_dist
1045                        // convert everything to cm first!
1046                        // sdd is in [cm], everything else is in [mm]
1047                        tmp_dist = (data_realDistY/10 - yctr/10)/sqrt((data_realDistX/10 - xctr/10)^2 + sdd^2)         
1048                        tmp_theta_i = atan(tmp_dist)            //this is theta_y
1049                       
1050                else
1051                        // horizontal orientation (T/B panels)
1052                        // this is now a different tmp_dist
1053                        // convert everything to cm first!
1054                        // sdd is in [cm], everything else is in [mm]
1055                        tmp_dist = (data_realDistX/10 - xctr/10)/sqrt((data_realDistY/10 - yctr/10)^2 + sdd^2)         
1056                        tmp_theta_i = atan(tmp_dist)            //this is theta_x
1057               
1058                endif
1059               
1060                for(ii=0        ;ii<numx;ii+=1)
1061                        for(jj=0;jj<numy;jj+=1)
1062                               
1063                                if(ii==0)               //do a forward difference if ii==0
1064                                        dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj])        //delta x for the pixel
1065                                else
1066                                        dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj])        //delta x for the pixel
1067                                endif
1068                               
1069                               
1070                                if(jj==0)
1071                                        dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj])        //delta y for the pixel
1072                                else
1073                                        dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1])        //delta y for the pixel
1074                                endif
1075               
1076                                dx /= 10
1077                                dy /= 10                // convert mm to cm (since sdd is in cm)
1078                                solid_angle[ii][jj] = dx*dy             //this is in cm^2
1079                        endfor
1080                endfor
1081               
1082                // to cover up any issues w/negative dx or dy
1083                solid_angle = abs(solid_angle)
1084               
1085                // solid_angle correction
1086                // == dx*dy*cos(th)^2*cos(th_i)/sdd^2           using either the theta_x or theta_y value
1087                solid_angle *= (cos(tmp_theta))^2*cos(tmp_theta_i)
1088                solid_angle /= sdd^2
1089               
1090                // Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!)
1091                w /= solid_angle
1092               
1093                //
1094                // correctly apply the correction to the error wave (assume a perfect value?)
1095                w_err /= solid_angle            //
1096       
1097        endif
1098       
1099
1100// DONE x- clean up after I'm satisfied computations are correct               
1101        KillWaves/Z tmp_theta,tmp_dist,tmp_theta_i
1102       
1103        return(0)
1104end
1105
1106// this is the incorrect solid angle correction that does not take into
1107// account the tube geometry. It is correct for the high-res detector (and the 30m Ordela)
1108//
1109// -- only for testing to prove that the cos(th)^2 *cos(th_i) is correct
1110//
1111Function V_SolidAngleCorrection_COS3(w,w_err,fname,detStr,destPath)
1112        Wave w,w_err
1113        String fname,detStr,destPath
1114
1115        Variable sdd,xCtr,yCtr,lambda
1116        String orientation
1117       
1118// get all of the geometry information 
1119        orientation = V_getDet_tubeOrientation(fname,detStr)
1120        sdd = V_getDet_ActualDistance(fname,detStr)
1121
1122        // this is ctr in mm
1123        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
1124        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
1125        lambda = V_getWavelength(fname)
1126       
1127        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
1128       
1129        Wave data_realDistX = data_realDistX
1130        Wave data_realDistY = data_realDistY
1131
1132        Duplicate/O w solid_angle,tmp_theta,tmp_dist    //in the current df
1133
1134//// calculate the scattering angle
1135//      dx = (distX - xctr)             //delta x in mm
1136//      dy = (distY - yctr)             //delta y in mm
1137        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
1138       
1139        tmp_dist /= 10  // convert mm to cm
1140        // sdd is in [cm]
1141
1142        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the (total) scattering angle
1143
1144        Variable ii,jj,numx,numy,dx,dy
1145        numx = DimSize(tmp_theta,0)
1146        numy = DimSize(tmp_theta,1)
1147
1148        for(ii=0        ;ii<numx;ii+=1)
1149                for(jj=0;jj<numy;jj+=1)
1150                       
1151                        if(ii==0)               //do a forward difference if ii==0
1152                                dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj])        //delta x for the pixel
1153                        else
1154                                dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj])        //delta x for the pixel
1155                        endif
1156                       
1157                       
1158                        if(jj==0)
1159                                dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj])        //delta y for the pixel
1160                        else
1161                                dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1])        //delta y for the pixel
1162                        endif
1163       
1164                        dx /= 10
1165                        dy /= 10                // convert mm to cm (since sdd is in cm)
1166                        solid_angle[ii][jj] = dx*dy             //this is in cm^2
1167                endfor
1168        endfor
1169       
1170        // to cover up any issues w/negative dx or dy
1171        solid_angle = abs(solid_angle)
1172       
1173        // solid_angle correction
1174        // == dx*dy*cos^3/sdd^2
1175        solid_angle *= (cos(tmp_theta))^3
1176        solid_angle /= sdd^2
1177       
1178        // Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!)
1179        w /= solid_angle
1180
1181        // correctly apply the correction to the error wave (assume a perfect value?)
1182        w_err /= solid_angle            //
1183       
1184
1185// DONE x- clean up after I'm satisfied computations are correct               
1186        KillWaves/Z tmp_theta,tmp_dist,tmp_theta_i
1187       
1188        return(0)
1189end
1190
1191
1192
1193//
1194// Large angle transmission correction
1195//
1196// DIVIDE the intensity by this correction to get the right answer
1197//
1198//
1199// Apply the large angle transmssion correction as the data is converted to WORK
1200// so that whether the data is saved as 2D or 1D, the correction has properly been done.
1201//
1202// This is, however, a SAMPLE dependent calculation, not purely instrument geometry.
1203//
1204Function V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
1205        Wave w,w_err
1206        String fname,detStr,destPath
1207
1208        Variable sdd,xCtr,yCtr,trans,trans_err,uval
1209
1210// get all of the geometry information 
1211//      orientation = V_getDet_tubeOrientation(fname,detStr)
1212        sdd = V_getDet_ActualDistance(fname,detStr)
1213
1214        // this is ctr in mm
1215        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
1216        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
1217        trans = V_getSampleTransmission(fname)
1218        trans_err = V_getSampleTransError(fname)
1219       
1220        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
1221       
1222        Wave data_realDistX = data_realDistX
1223        Wave data_realDistY = data_realDistY
1224
1225        Duplicate/O w lat_corr,tmp_theta,tmp_dist,lat_err,tmp_err               //in the current df
1226
1227//// calculate the scattering angle
1228//      dx = (distX - xctr)             //delta x in mm
1229//      dy = (distY - yctr)             //delta y in mm
1230        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
1231       
1232        tmp_dist /= 10  // convert mm to cm
1233        // sdd is in [cm]
1234
1235        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle
1236
1237        Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp
1238        numx = DimSize(tmp_theta,0)
1239        numy = DimSize(tmp_theta,1)
1240       
1241       
1242        //optical thickness
1243        uval = -ln(trans)               //use natural logarithm
1244       
1245        for(ii=0        ;ii<numx;ii+=1)
1246                for(jj=0;jj<numy;jj+=1)
1247                       
1248                        cos_th = cos(tmp_theta[ii][jj])
1249                        arg = (1-cos_th)/cos_th
1250                       
1251                        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
1252                        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120
1253                        // OR
1254                        if((uval<0.01) || (cos_th>0.99))       
1255                                //small arg, approx correction
1256                                lat_corr[ii][jj] = 1-0.5*uval*arg
1257                        else
1258                                //large arg, exact correction
1259                                lat_corr[ii][jj] = (1-exp(-uval*arg))/(uval*arg)
1260                        endif
1261                         
1262                        // (DONE)
1263                        // x- properly calculate and apply the 2D error propagation
1264                        if(trans == 1)
1265                                lat_err[ii][jj] = 0             //no correction, no error
1266                        else
1267                                //sigT, calculated from the Taylor expansion
1268                                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3)
1269                                tmp *= tmp
1270                                tmp *= trans_err^2
1271                                tmp = sqrt(tmp)         //sigT
1272                               
1273                                lat_err[ii][jj] = tmp
1274                        endif
1275                         
1276 
1277                endfor
1278        endfor
1279       
1280
1281       
1282        // Here it is! Apply the correction to the intensity (divide -- to get the proper correction)
1283        w /= lat_corr
1284
1285        // relative errors add in quadrature to the current 2D error
1286        tmp_err = (w_err/lat_corr)^2 + (lat_err/lat_corr)^2*w*w/lat_corr^2
1287        tmp_err = sqrt(tmp_err)
1288       
1289        w_err = tmp_err
1290       
1291
1292        // DONE x- clean up after I'm satisfied computations are correct               
1293        KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err
1294       
1295        return(0)
1296end
1297
1298
1299
1300//
1301//test procedure, not called anymore
1302Proc V_AbsoluteScaling(type,c0,c1,c2,c3,c4,c5,I_err)
1303        String type
1304        Variable c0=1,c1=0.1,c2=0.95,c3=0.1,c4=1,c5=32.0,I_err=0.32
1305        Prompt type,"WORK data type",popup,"CAL;COR;SAM"
1306        Prompt c0, "Sample Transmission"
1307        Prompt c1, "Sample Thickness (cm)"
1308        Prompt c2, "Standard Transmission"
1309        Prompt c3, "Standard Thickness (cm)"
1310        Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)"
1311        Prompt c5, "Standard Cross-Section (cm-1)"
1312        Prompt I_err, "error in I(q=0) (one std dev)"
1313
1314        Variable err
1315        //call the function to do the math
1316        //data from "type" will be scaled and deposited in ABS
1317        err = V_Absolute_Scale(type,c0,c1,c2,c3,c4,c5,I_err)
1318       
1319        if(err)
1320                Abort "Error in V_Absolute_Scale()"
1321        endif
1322       
1323        //contents are always dumped to ABS
1324        type = "ABS"
1325       
1326        //need to update the display with "data" from the correct dataFolder
1327        //reset the current display type to "type"
1328        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
1329        gCurDispType = Type     
1330       
1331        V_UpdateDisplayInformation(Type)
1332       
1333End
1334
1335//
1336//
1337// kappa comes in as s_izero, so be sure to use 1/kappa_err
1338//
1339//convert the "type" data to absolute scale using the given standard information
1340//s_ is the standard
1341//w_ is the "work" file
1342//both are work files and should already be normalized to 10^8 monitor counts
1343Function V_Absolute_Scale(type,absStr)
1344        String type,absStr
1345       
1346       
1347        Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
1348
1349        Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
1350        Variable scale,trans_err
1351        Variable err,ii
1352        String detStr
1353       
1354        // be sure that the starting data exists
1355        err = V_WorkDataExists(type)
1356        if(err==1)
1357                return(err)
1358        endif
1359               
1360        //copy from current dir (type) to ABS
1361        V_CopyHDFToWorkFolder(type,"ABS")       
1362
1363// TODO: -- which monitor to use? Here, I think it should already be normalized to 10^8
1364//     
1365//      w_moncount = V_getMonitorCount(type)            //monitor count in "type"
1366       
1367        w_moncount = V_getBeamMonNormData(type)
1368               
1369        if(w_moncount == 0)
1370                //zero monitor counts will give divide by zero ---
1371                DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
1372                Return(1)               //report error
1373        Endif
1374
1375        w_trans = V_getSampleTransmission(type)         //sample transmission
1376        w_thick = V_getSampleThickness(type)            //sample thickness
1377        trans_err = V_getSampleTransError(type)
1378       
1379       
1380        //get the parames from the list
1381        s_trans = NumberByKey("TSTAND", absStr, "=", ";")       //parse the list of values
1382        s_thick = NumberByKey("DSTAND", absStr, "=", ";")
1383        s_izero = NumberByKey("IZERO", absStr, "=", ";")
1384        s_cross = NumberByKey("XSECT", absStr, "=", ";")
1385        kappa_err = NumberByKey("SDEV", absStr, "=", ";")
1386
1387       
1388        //calculate scale factor
1389        s1 = defmon/w_moncount          // monitor count (s1 should be 1)
1390        s2 = s_thick/w_thick
1391        s3 = s_trans/w_trans
1392        s4 = s_cross/s_izero
1393        scale = s1*s2*s3*s4
1394
1395       
1396        // kappa comes in as s_izero, so be sure to use 1/kappa_err
1397
1398        // and now loop through all of the detectors
1399        //do the actual absolute scaling here, modifying the data in ABS
1400        for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
1401                detStr = StringFromList(ii, ksDetectorListNoB, ";")
1402                Wave data = V_getDetectorDataW("ABS",detStr)
1403                Wave data_err = V_getDetectorDataErrW("ABS",detStr)
1404               
1405                data *= scale
1406                data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1407        endfor
1408
1409        // do the back detector separately, if it is set to be used
1410        NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
1411        if(gIgnoreDetB == 0)
1412                detStr = "B"
1413                Wave data = V_getDetectorDataW("ABS",detStr)
1414                Wave data_err = V_getDetectorDataErrW("ABS",detStr)
1415               
1416                //get the parames from the list
1417                s_trans = NumberByKey("TSTAND_B", absStr, "=", ";")     //parse the list of values
1418                s_thick = NumberByKey("DSTAND_B", absStr, "=", ";")
1419                s_izero = NumberByKey("IZERO_B", absStr, "=", ";")
1420                s_cross = NumberByKey("XSECT_B", absStr, "=", ";")
1421                kappa_err = NumberByKey("SDEV_B", absStr, "=", ";")
1422
1423                //calculate scale factor
1424                s1 = defmon/w_moncount          // monitor count (s1 should be 1)
1425                s2 = s_thick/w_thick
1426                s3 = s_trans/w_trans
1427                s4 = s_cross/s_izero
1428                scale = s1*s2*s3*s4
1429               
1430                data *= scale
1431                data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1432        endif
1433       
1434        //********* 15APR02
1435        // DO NOT correct for atenuators here - the COR step already does this, putting all of the data on equal
1436        // footing (zero atten) before doing the subtraction.
1437       
1438        Return (0) //no error
1439End
1440
1441
1442//
1443// TODO:
1444//   --         DoAlert 0,"This has not yet been updated for VSANS"
1445//
1446//
1447// match the attenuation of the RAW data to the "type" data
1448// so that they can be properly added
1449//
1450// are the attenuator numbers the same? if so exit
1451//
1452// if not, find the attenuator number for type
1453// - find both attenuation factors
1454//
1455// rescale the raw data to match the ratio of the two attenuation factors
1456// -- adjust the detector count (rw)
1457// -- the linear data
1458//
1459//
1460Function V_Adjust_RAW_Attenuation(type)
1461        String type
1462
1463        DoAlert 0,"This has not yet been updated for VSANS"
1464       
1465        WAVE rw=$("root:Packages:NIST:RAW:realsread")
1466        WAVE linear_data=$("root:Packages:NIST:RAW:linear_data")
1467        WAVE data=$("root:Packages:NIST:RAW:data")
1468        WAVE data_err=$("root:Packages:NIST:RAW:linear_data_error")
1469        WAVE/T tw = $("root:Packages:NIST:RAW:textRead")
1470       
1471        WAVE dest_reals=$("root:Packages:NIST:"+type+":realsread")
1472
1473        Variable dest_atten,raw_atten,tol
1474        Variable lambda,raw_atten_err,raw_AttenFactor,dest_attenFactor,dest_atten_err
1475        String fileStr
1476
1477        dest_atten = dest_reals[3]
1478        raw_atten = rw[3]
1479       
1480        tol = 0.1               // within 0.1 atten units is OK
1481        if(abs(dest_atten - raw_atten) < tol )
1482                return(0)
1483        endif
1484
1485        fileStr = tw[3]
1486        lambda = rw[26]
1487        // TODO access correct values
1488        raw_AttenFactor = 1//AttenuationFactor(fileStr,lambda,raw_atten,raw_atten_err)
1489        dest_AttenFactor = 1//AttenuationFactor(fileStr,lambda,dest_atten,dest_atten_err)
1490               
1491        rw[2] *= dest_AttenFactor/raw_AttenFactor
1492        linear_data *= dest_AttenFactor/raw_AttenFactor
1493       
1494        // to keep "data" and linear_data in sync
1495        data = linear_data
1496       
1497        return(0)
1498End
1499
1500//
1501// testing procedure, called from a menu selection
1502//
1503Proc V_DIV_a_Workfile(type)
1504        String type
1505        Prompt type,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
1506       
1507        //macro will take whatever is in SELECTED folder and DIVide it by the current
1508        //contents of the DIV folder - the function will check for existence
1509        //before proceeding
1510
1511        Abort "This has not yet been updated for VSANS"
1512       
1513        Variable err
1514        err = V_DIVCorrection(type)             //returns err = 1 if data doesn't exist in specified folders
1515       
1516        if(err)
1517                Abort "error in V_DIVCorrection()"
1518        endif
1519       
1520        //contents are NOT always dumped to CAL, but are in the new type folder
1521       
1522        String newTitle = "WORK_"+type
1523        DoWindow/F VSANS_Data
1524        DoWindow/T VSANS_Data, newTitle
1525        KillStrings/Z newTitle
1526       
1527        //need to update the display with "data" from the correct dataFolder
1528        //reset the current displaytype to "type"
1529        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=Type
1530       
1531        V_UpdateDisplayInformation(type)
1532       
1533End
1534
1535
1536//
1537// TODO:
1538//   x-         DoAlert 0,"This has not yet been updated for VSANS"
1539//   -- how is the error propagation handled? Be sure it is calculated correctly when DIV is generated
1540//      and is applied correctly here...
1541//
1542//function will divide the contents of "workType" folder with the contents of
1543//the DIV folder + detStr
1544// all data is linear scale for the calculation
1545//
1546Function V_DIVCorrection(data,data_err,detStr,workType)
1547        Wave data,data_err
1548        String detStr,workType
1549       
1550        //check for existence of data in type and DIV
1551        // if the desired data doesn't exist, let the user know, and abort
1552        String destPath=""
1553
1554        NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
1555        if(cmpstr(detStr,"B")==0 && gIgnoreDetB)
1556                return(0)
1557        endif
1558
1559
1560        if(WaveExists(data) == 0)
1561                Print "The data wave does not exist in V_DIVCorrection()"
1562                Return(1)               //error condition
1563        Endif
1564       
1565        //check for DIV
1566        // if the DIV workfile doesn't exist, let the user know,and abort
1567        // !! be sure to check first, before trying to access the wave
1568       
1569//      WAVE/Z div_data = $("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")
1570        if(WaveExists($("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")) == 0)
1571                Print "The DIV wave does not exist in V_DIVCorrection()"
1572                Return(1)               //error condition
1573        Endif
1574        if(WaveExists($("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":linear_data_error")) == 0)
1575                Print "The DIV error wave does not exist in V_DIVCorrection()"
1576                Return(1)               //error condition
1577        Endif
1578        //files exist, proceed
1579
1580        WAVE/Z div_data_err = V_getDetectorDataErrW("DIV",detStr)
1581        WAVE/Z div_data = V_getDetectorDataW("DIV",detStr)
1582
1583
1584
1585// do the error propagation first, since data is changed by the correction
1586        data_err = sqrt(data_err^2/div_data^2 + div_data_err^2 * data^2/div_data^4 )
1587
1588// then the correction
1589        data /= div_data
1590
1591       
1592        Return(0)
1593End
1594
1595
1596//////////////////////////
1597// detector corrections to stitch the back detector into one proper image
1598//
1599//
1600//
1601
1602
1603//
1604// to register the image on the back detector panel
1605//
1606// middle portion (552 pix in Y) is held fixed
1607// top portion of image is shifted right and down
1608// bottom portion of image is shifted right and up
1609//
1610// remainder of image is filled with Zero (NaN causes problems converting to WORK)
1611//
1612// currently, data is not added together and averaged, but it could be
1613//
1614Function V_ShiftBackDetImage(w,adjW)
1615        Wave w,adjW
1616
1617        NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
1618
1619// this is necessary for some old data with the 150x150 back (dummy) panel
1620// the proper back detector has an x-dimension of 680 pixels. Don't do the shift
1621// if the dimensions are incorrect.
1622        if(DimSize(w,0) < 680)
1623                adjW=w
1624                return(0)
1625        endif
1626       
1627        adjW=0
1628               
1629        Variable topX,bottomX
1630        Variable topY,bottomY
1631        Variable totalY,ccdX,ccdY
1632       
1633//      topX = 7
1634//      topY = 105
1635       
1636//      bottomX = 5
1637//      bottomY = 35
1638
1639// TODOHIGHRES
1640// the detector pix dimensions are hard-wired, be sure the are correct
1641        switch(gHighResBinning)
1642                case 1:
1643                        topX = kShift_topX_bin1
1644                        topY = kShift_topY_bin1
1645                        bottomX = kShift_bottomX_bin1
1646                        bottomY = kShift_bottomY_bin1
1647                       
1648                        totalY = 6624   // total YDim
1649                        ccdY = 2208             // = YDim/3
1650                        ccdX = 2720             // = xDim
1651                        break
1652                case 4:
1653                        topX = kShift_topX_bin4
1654                        topY = kShift_topY_bin4
1655                        bottomX = kShift_bottomX_bin4
1656                        bottomY = kShift_bottomY_bin4
1657                       
1658                        totalY = 1656   // total YDim
1659                        ccdY = 552              // = YDim/3
1660                        ccdX = 680              // = xDim
1661
1662                       
1663                        break
1664                default:               
1665                        Abort "No binning case matches in V_ShiftBackDetImage"
1666                       
1667        endswitch
1668
1669                // middle
1670                adjW[][ccdY,ccdY+ccdY] = w[p][q]
1671       
1672                //top
1673                adjW[0+topX,ccdX-1][ccdY+ccdY,totalY-1-topY] = w[p-topX][q+topY]
1674               
1675                //bottom
1676                adjW[0+bottomX,ccdX-1][0+bottomY,ccdY-1] = w[p-bottomX][q-bottomY]
1677
1678       
1679        return(0)
1680End
1681
1682
1683Proc pV_MedianFilterBack(folder)
1684        String folder="RAW"
1685       
1686        V_MedianFilterBack(folder)
1687end
1688
1689Function V_MedianFilterBack(folder)
1690        String folder
1691
1692        Wave w = V_getDetectorDataW(folder,"B")
1693
1694        NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
1695        switch(gHighResBinning)
1696                case 1:
1697                        MatrixFilter /N=11 /P=1 median w                        //              /P=n flag sets the number of passes (default is 1 pass)
1698                       
1699                        Print "*** median noise filter 11x11 applied to the back detector (1 pass) ***"
1700                        break
1701                case 4:
1702                        MatrixFilter /N=3 /P=1 median w                 //              /P=n flag sets the number of passes (default is 1 pass)
1703                       
1704                        Print "*** median noise filter 3x3 applied to the back detector (1 pass) ***"
1705                        break
1706                default:
1707                        Abort "No binning case matches in V_MedianFilterBack"
1708        endswitch
1709
1710        return(0)
1711End
1712
1713
1714Proc pV_SubtractReadNoiseBack(folder,ReadNoise)
1715        String folder="RAW"
1716        Variable readNoise=3160
1717       
1718        V_SubtractReadNoiseBack(folder,readNoise)
1719end
1720
1721Function V_SubtractReadNoiseBack(folder,readNoise)
1722        String folder
1723        Variable readNoise
1724
1725                Wave w = V_getDetectorDataW(folder,"B")
1726                w -= readNoise          // a constant value
1727               
1728//              MatrixFilter /N=3 median w
1729//              Print "*** median noise filter applied to the back detector***"
1730       
1731        return(0)
1732End
1733
1734
1735Proc pV_MedianAndReadNoiseBack(folder,ReadNoise)
1736        String folder="RAW"
1737        Variable readNoise=3160
1738       
1739        V_MedianAndReadNoiseBack(folder,readNoise)
1740end
1741
1742Function V_MedianAndReadNoiseBack(folder,readNoise)
1743        String folder
1744        Variable readNoise
1745
1746                Wave w = V_getDetectorDataW(folder,"B")
1747                w -= readNoise          // a constant value
1748
1749                NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
1750                switch(gHighResBinning)
1751                        case 1:
1752                                MatrixFilter /N=11 /P=1 median w                        //              /P=n flag sets the number of passes (default is 1 pass)
1753                               
1754                                Print "*** median noise filter 11x11 applied to the back detector (1 pass) ***"
1755                                break
1756                        case 4:
1757                                MatrixFilter /N=3 /P=1 median w                 //              /P=n flag sets the number of passes (default is 1 pass)
1758                               
1759                                Print "*** median noise filter 3x3 applied to the back detector (1 pass) ***"
1760                                break
1761                        default:
1762                                Abort "No binning case matches in V_MedianAndReadNoiseBack"
1763                endswitch
1764               
1765        return(0)
1766End
1767
Note: See TracBrowser for help on using the repository browser.