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

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

updating the IgorVersion? pragma to v7.0 for all files to be consistent.

File size: 68.4 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2#pragma version=1.0
3#pragma IgorVersion = 7.00
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// Constants for detector efficiency and shadowing
16//
17// V_TubeEfficiencyShadowCorr()
18//
19// JAN 2020
20///////////////
21Constant kTube_ri = 0.372               // inner radius of tube [cm]
22Constant kTube_cc = 0.84                        // center to center spacing [cm]
23Constant kTube_ss = 0.025               // stainless steel shell thickness [cm]
24
25Constant kSig_2b_He = 0.146             // abs xs for 2 bar He(3) [cm-1 A-1] (multiply this by wavelength)
26Constant kSig_8b_He = 0.593             // abs xs for 8 bar He(3) [cm-1 A-1] (multiply this by wavelength)
27Constant kSig_Al = 0.00967              // abs xs for Al [cm-1 A-1] (multiply this by wavelength)
28Constant kSig_ss = 0.146                // abs xs for 304 SS [cm-1 A-1] (multiply this by wavelength)
29
30
31
32
33//
34// detector dead time
35//
36// input is the data array (N tubes x M pixels)
37// input of N x 1 array of dead time values
38//
39// output is the corrected counts in data, overwriting the input data
40//
41// Note that the equation in Roe (eqn 2.15, p. 63) looks different, but it is really the
42// same old equation, just written in a more complex form.
43//
44// (DONE)
45// x- verify the direction of the tubes and indexing
46// x- decide on the appropriate functional form for the tubes
47// x- need count time as input
48// x- be sure I'm working in the right data folder (all waves are passed in)
49// x- clean up when done
50// x- calculate + return the error contribution?
51// x- verify the error propagation
52Function V_DeadTimeCorrectionTubes(dataW,data_errW,dtW,ctTime)
53        Wave dataW,data_errW,dtW
54        Variable ctTime
55       
56        // do I count on the orientation as an input, or do I just figure it out on my own?
57        String orientation
58        Variable dimX,dimY
59        dimX = DimSize(dataW,0)
60        dimY = DimSize(dataw,1)
61        if(dimX > dimY)
62                orientation = "horizontal"
63        else
64                orientation = "vertical"
65        endif
66       
67        // sum the counts in each tube and divide by time for total cr per tube
68        Variable npt
69       
70        if(cmpstr(orientation,"vertical")==0)
71                //      this is data dimensioned as (Ntubes,Npix)
72               
73                MatrixOp/O sumTubes = sumRows(dataW)            // n x 1 result
74                sumTubes /= ctTime              //now count rate per tube
75               
76                dataW[][] = dataW[p][q]/(1-sumTubes[p]*dtW[p])          //correct the data
77                data_errW[][] = data_errW[p][q]/(1-sumTubes[p]*dtW[p])          // propagate the error wave
78
79        elseif(cmpstr(orientation,"horizontal")==0)
80        //      this is data (horizontal) dimensioned as (Npix,Ntubes)
81
82                MatrixOp/O sumTubes = sumCols(dataW)            // 1 x m result
83                sumTubes /= ctTime
84               
85                dataW[][] = dataW[p][q]/(1-sumTubes[q]*dtW[q])
86                data_errW[][] = data_errW[p][q]/(1-sumTubes[q]*dtW[q])
87       
88        else           
89                DoAlert 0,"Orientation not correctly passed in DeadTimeCorrectionTubes(). No correction done."
90        endif
91       
92        return(0)
93end
94
95// test function
96Function V_testDTCor()
97
98        String detStr = ""
99        String fname = "RAW"
100        Variable ctTime
101       
102        detStr = "FR"
103        Wave w = V_getDetectorDataW(fname,detStr)
104        Wave w_err = V_getDetectorDataErrW(fname,detStr)
105        Wave w_dt = V_getDetector_deadtime(fname,detStr)
106
107        ctTime = V_getCount_time(fname)
108       
109//      ctTime = 10
110        V_DeadTimeCorrectionTubes(w,w_err,w_dt,ctTime)
111
112End
113
114
115//
116// Non-linear data correction
117//
118// DOES NOT modify the data, only calculates the spatial relationship
119//
120// input is the data array (N tubes x M pixels)
121// input of N x M array of quadratic coefficients
122//
123// output is wave of corrected real space distance corresponding to each pixel of the data
124// ** its distance from the nominal beam center of (0,0) **
125//
126//
127// (DONE)
128// x- UNITS!!!! currently this is mm, which certainly doesn't match anything else!!!
129//
130// x- verify the direction of the tubes and indexing
131// x- be sure I'm working in the right data folder (it is passed in, and the full path is used)
132// x- clean up when done
133// x- calculate + return the error contribution? (there is none for this operation)
134// x- do I want this to return a wave? (no, default names are generated)
135// x- do I need to write a separate function that returns the distance wave for later calculations?
136// x- do I want to make the distance array 3D to keep the x and y dims together? Calculate them all right now?
137// x- what else do I need to pass to the function? (fname=folder? detStr?)
138// y- (yes,see below) need a separate block or function to handle "B" detector which will be ? different
139//
140//
141Function V_NonLinearCorrection(fname,dataW,coefW,tube_width,detStr,destPath)
142        String fname            //can also be a folder such as "RAW"
143        Wave dataW,coefW
144        Variable tube_width
145        String detStr,destPath
146       
147         
148        // do I count on the orientation as an input, or do I just figure it out on my own?
149        String orientation
150        Variable dimX,dimY
151        dimX = DimSize(dataW,0)
152        dimY = DimSize(dataW,1)
153        if(dimX > dimY)
154                orientation = "horizontal"
155        else
156                orientation = "vertical"
157        endif
158
159        // make a wave of the same dimensions, in the same data folder for the distance
160        // ?? or a 3D wave?
161        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
162        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
163        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
164        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
165       
166        // then per tube, do the quadratic calculation to get the real space distance along the tube
167        // the distance perpendicular to the tube is n*(8.4mm) per tube index
168       
169        // DONE
170        // -- GAP was hard-wired, but in 2018 proper values for all 4 gaps were measured
171        // and added to the file header for each detector panel. there is now a read from the
172        // header to get the gap value
173        Variable offset,gap
174
175        gap = V_getDet_panel_gap(fname,detStr)
176
177// DONE:
178// -- in case of error, V_getDet_panel_gap() will return -999999
179// -- it should only apply to data pre-2018 when the field did not exist in the file
180// -- any VSANS data from 2018+ should read gap from the file and bypass the if()
181
182        if(gap < -100)          //-999999 returned if field is missing from file
183       
184                if(cmpstr(detStr,"FL") == 0 || cmpstr(detStr,"FR") == 0)
185                        gap = 3.5               //mm (measured, JB 1/4/18)
186                endif
187                if(cmpstr(detStr,"FT") == 0 || cmpstr(detStr,"FB") == 0)
188                        gap = 3.3               //mm (measured, JB 2/1/18)
189                endif
190                if(cmpstr(detStr,"ML") == 0 || cmpstr(detStr,"MR") == 0)
191                        gap = 5.9               //mm (measured, JB 1/4/18)
192                endif
193                if(cmpstr(detStr,"MT") == 0 || cmpstr(detStr,"MB") == 0)
194                        gap = 18.3              //mm (measured, JB 2/1/18)
195                endif
196       
197        endif
198
199       
200        if(cmpstr(orientation,"vertical")==0)
201                //      this is data dimensioned as (Ntubes,Npix)
202       
203                // adjust the x postion based on the beam center being nominally (0,0) in units of cm, not pixels
204                if(cmpstr(fname,"VCALC")== 0 )
205                        offset = VCALC_getPanelTranslation(detStr)
206                        offset *= 10                    // convert to units of mm
207//                      if(cmpstr("L",detStr[1]) == 0)
208//                              offset *= -1            //negative value for L
209//                      endif
210                else
211                        //normal case
212                        offset = V_getDet_LateralOffset(fname,detStr)
213                        offset *= 10 //convert cm to mm
214                endif
215               
216        // calculation is in mm, not cm
217        // offset will be a negative value for the L panel, and positive for the R panel
218                if(kBCTR_CM)
219                        if(cmpstr("L",detStr[1]) == 0)
220//                              data_realDistX[][] = offset - (dimX - p)*tube_width                     // TODO should this be dimX-1-p = 47-p?
221                                data_realDistX[][] = offset - (dimX - p - 1/2)*tube_width - gap/2               // TODO should this be dimX-1-p = 47-p?
222                        else
223                        //      right
224//                              data_realDistX[][] = tube_width*(p+1) + offset + gap            //add to the Right det,
225                                data_realDistX[][] = tube_width*(p+1/2) + offset + gap/2                //add to the Right det
226                        endif
227                else
228                        data_realDistX[][] = tube_width*(p)
229                endif
230                data_realDistY[][] = coefW[0][p] + coefW[1][p]*q + coefW[2][p]*q*q
231       
232       
233        elseif(cmpstr(orientation,"horizontal")==0)
234                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
235                data_realDistY[][] = tube_width*q
236
237                if(cmpstr(fname,"VCALC")== 0 )
238                        offset = VCALC_getPanelTranslation(detStr)
239                        offset *= 10                    // convert to units of mm
240//                      if(cmpstr("B",detStr[1]) == 0)
241//                              offset *= -1    // negative value for Bottom det
242//                      endif
243                else
244                        //normal case
245                        offset = V_getDet_VerticalOffset(fname,detStr)
246                        offset *= 10 //convert cm to mm
247                endif
248               
249                if(kBCTR_CM)
250                        if(cmpstr("T",detStr[1]) == 0)
251//                              data_realDistY[][] = tube_width*(q+1) + offset + gap                   
252                                data_realDistY[][] = tube_width*(q+1/2) + offset + gap/2                       
253                        else
254                                // bottom
255//                              data_realDistY[][] = offset - (dimY - q)*tube_width     // TODO should this be dimY-1-q = 47-q?
256                                data_realDistY[][] = offset - (dimY - q - 1/2)*tube_width - gap/2       // TODO should this be dimY-1-q = 47-q?
257                        endif
258                else
259                        data_realDistY[][] = tube_width*(q)
260                endif
261                data_realDistX[][] = coefW[0][q] + coefW[1][q]*p + coefW[2][q]*p*p
262
263        else           
264                DoAlert 0,"Orientation not correctly passed in NonLinearCorrection(). No correction done."
265                return(0)
266        endif
267       
268        return(0)
269end
270
271
272
273
274// TODO:
275// -- the cal_x and y coefficients are totally fake
276// -- the wave assignment may not be correct.. so beware
277//
278//
279Function V_NonLinearCorrection_B(folder,dataW,cal_x,cal_y,detStr,destPath)
280        String folder
281        Wave dataW,cal_x,cal_y
282        String detStr,destPath
283
284        if(cmpstr(detStr,"B") != 0)
285                return(0)
286        endif
287
288Print "***Cal_X and Cal_Y for Back are using file values ***"
289
290//              cal_x[0] = VCALC_getPixSizeX(detStr)*10                 // pixel size in mm  VCALC_getPixSizeX(detStr) is [cm]
291//              cal_x[1] = 1
292//              cal_x[2] = 10000
293//              cal_y[0] = VCALC_getPixSizeY(detStr)*10                 // pixel size in mm  VCALC_getPixSizeX(detStr) is [cm]
294//              cal_y[1] = 1
295//              cal_y[2] = 10000
296
297       
298        // do I count on the orientation as an input, or do I just figure it out on my own?
299        Variable dimX,dimY
300       
301//      Wave dataW = V_getDetectorDataW(folder,detStr)
302       
303        dimX = DimSize(dataW,0)
304        dimY = DimSize(dataW,1)
305
306        // make a wave of the same dimensions, in the same data folder for the distance
307        // ?? or a 3D wave?
308        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
309        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
310        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
311        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
312       
313//      Wave cal_x = V_getDet_cal_x(folder,detStr)
314//      Wave cal_y = V_getDet_cal_y(folder,detStr)
315       
316        data_realDistX[][] = cal_x[0]*p*10              // cal_x and cal_y are in [cm], need mm
317        data_realDistY[][] = cal_y[0]*q*10
318       
319        return(0)
320end
321
322
323//
324//
325// TODO
326// -- VERIFY the calculations
327// -- verify where this needs to be done (if the beam center is changed)
328// -- then the q-calculation needs to be re-done
329// -- the position along the tube length is referenced to tube[0], for no particular reason
330//    It may be better to take an average? but [0] is an ASSUMPTION
331// -- distance along tube is simple interpolation, or do I use the coefficients to
332//    calculate the actual value
333//
334// -- distance in the lateral direction is based on tube width, which is a fixed parameter
335//
336//
337Function V_ConvertBeamCtrPix_to_mm(folder,detStr,destPath)
338        String folder,detStr,destPath
339       
340        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
341        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")     
342
343        String orientation
344        Variable dimX,dimY,xCtr,yCtr
345        dimX = DimSize(data_realDistX,0)
346        dimY = DimSize(data_realDistX,1)
347        if(dimX > dimY)
348                orientation = "horizontal"
349        else
350                orientation = "vertical"
351        endif
352       
353        xCtr = V_getDet_beam_center_x(folder,detStr)
354        yCtr = V_getDet_beam_center_y(folder,detStr)   
355       
356        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
357        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
358        WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
359        WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
360
361        Variable tube_width = V_getDet_tubeWidth(folder,detStr)         //this is also in mm
362
363//
364        strswitch(detStr)       // string switch
365                case "FL":
366                case "ML":
367                        // for Left/Right
368                        // for left
369                        x_mm[0] = data_realDistX[dimX-1][0] + (xCtr-dimX-1)*tube_width
370                        y_mm[0] = data_realDistY[0][yCtr]
371       
372                        break           
373                case "FR":     
374                case "MR":
375                        // for Left/Right
376                        // for right
377                        x_mm[0] = data_realDistX[0][0] + xCtr*tube_width
378                        y_mm[0] = data_realDistY[0][yCtr]
379                       
380                        break
381                case "FT":     
382                case "MT":
383                        // for Top                     
384                        x_mm[0] = data_realDistX[xCtr][0]
385                        y_mm[0] = data_realDistY[0][0] + yCtr*tube_width
386                       
387                        break           
388                case "FB":     
389                case "MB":
390                        // for Bottom                   
391                        x_mm[0] = data_realDistX[xCtr][0]
392                        y_mm[0] = data_realDistY[0][dimY-1] + (yCtr-dimY-1)*tube_width
393                                               
394                        break
395                default:                        // optional default expression executed
396                        Print "No case matched in V_Convert_FittedPix_2_cm"
397                        return(1)
398        endswitch
399               
400        return(0)
401end
402
403//
404//
405// (DONE)
406// x- VERIFY the calculations
407// x- verify where this needs to be done (if the beam center is changed)
408// x- then the q-calculation needs to be re-done
409// x- the position along the tube length is referenced to tube[0], for no particular reason
410//    It may be better to take an average? but [0] is an ASSUMPTION
411// x- distance along tube is simple interpolation
412//
413// x- distance in the lateral direction is based on tube width, which is a fixed parameter
414//
415// the value in pixels is written to the local data folder, NOT to disk (it is recalculated as needed)
416//
417Function V_ConvertBeamCtr_to_pix(folder,detStr,destPath)
418        String folder,detStr,destPath
419       
420        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
421        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")     
422
423        String orientation
424        Variable dimX,dimY,xCtr,yCtr
425        dimX = DimSize(data_realDistX,0)
426        dimY = DimSize(data_realDistX,1)
427        if(dimX > dimY)
428                orientation = "horizontal"
429        else
430                orientation = "vertical"
431        endif
432       
433        xCtr = V_getDet_beam_center_x(folder,detStr)            //these are in cm
434        yCtr = V_getDet_beam_center_y(folder,detStr)   
435       
436        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
437        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
438        WAVE x_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
439        WAVE y_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
440
441        Variable tube_width = V_getDet_tubeWidth(folder,detStr)
442
443        variable edge,delta
444        Variable gap
445
446// the gap is split equally between the panel pairs
447// DONE -- replace hard-wired values with V_getDet_panel_gap(fname,detStr) once it is added to the file
448
449        gap = V_getDet_panel_gap(folder,detStr)
450
451// DONE:
452// -- check in case of error, value should be read from header
453// -- it should only apply to data pre-2018 when the field did not exist in the file
454// -- any VSANS data from 2018+ should read gap from the file.
455
456        if(gap < -100)          //-999999 returned if field is missing from file
457                if(cmpstr(detStr,"FL") == 0 || cmpstr(detStr,"FR") == 0)
458                        gap = 3.5               //mm (measured, JB 1/4/18)
459                endif
460                if(cmpstr(detStr,"FT") == 0 || cmpstr(detStr,"FB") == 0)
461                        gap = 3.3               //mm (measured, JB 2/1/18)
462                endif
463                if(cmpstr(detStr,"ML") == 0 || cmpstr(detStr,"MR") == 0)
464                        gap = 5.9               //mm (measured, JB 1/4/18)
465                endif
466                if(cmpstr(detStr,"MT") == 0 || cmpstr(detStr,"MB") == 0)
467                        gap = 18.3              //mm (measured, JB 2/1/18)
468                endif
469        endif
470
471//
472        if(cmpstr(orientation,"vertical")==0)
473                //      this is data dimensioned as (Ntubes,Npix)
474
475                if(kBCTR_CM)
476                        if(cmpstr("L",detStr[1]) == 0)
477                                Make/O/D/N=(dimX) tmpTube
478                                tmpTube = data_RealDistX[p][0]
479                                FindLevel/P/Q tmpTube xCtr*10
480                                if(V_Flag)
481                                        edge = data_realDistX[47][0]            //tube 47
482                                        delta = abs(xCtr*10 - edge)
483                                        x_pix[0] = dimX-1 + delta/tube_width
484                                else
485                                        // beam center is on the panel, report the pixel value
486                                        x_pix[0] = V_LevelX
487                                endif
488                               
489                        else
490                        // R panel
491                                Make/O/D/N=(dimX) tmpTube
492                                tmpTube = data_RealDistX[p][0]
493                                FindLevel/P/Q tmpTube xCtr*10
494                                if(V_Flag)
495                                        //level not found
496                                        edge = data_realDistX[0][0]
497                                        delta = abs(xCtr*10 - edge + gap)               // how far past the edge of the panel
498                                        x_pix[0] = -delta/tube_width            //since the left edge of the R panel is pixel 0
499                                else
500                                        // beam center is on the panel, report the pixel value
501                                        x_pix[0] = V_LevelX
502                                endif
503                               
504                        endif
505
506                endif
507
508// the y-center will be on the panel in this direction
509                Make/O/D/N=(dimY) tmpTube
510                tmpTube = data_RealDistY[0][p]
511                FindLevel /P/Q tmpTube, yCtr*10
512               
513                y_pix[0] = V_levelX
514                KillWaves/Z tmpTube
515//              Print x_pix[0],y_pix[0]
516               
517        else
518                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
519
520                if(kBCTR_CM)
521                        if(cmpstr("T",detStr[1]) == 0)
522                                Make/O/D/N=(dimY) tmpTube
523                                tmpTube = data_RealDistY[p][0]
524                                FindLevel/P/Q tmpTube yCtr*10
525                                if(V_Flag)
526                                        edge = data_realDistY[0][0]             //tube 0
527                                        delta = abs(yCtr*10 - edge + gap)
528                                        y_pix[0] =  -delta/tube_width           //since the bottom edge of the T panel is pixel 0                               
529                                else
530                                        y_pix[0] = V_LevelX
531                                endif                           
532                               
533                        else
534                        // FM(B) panel
535                                Make/O/D/N=(dimY) tmpTube
536                                tmpTube = data_RealDistY[p][0]
537                                FindLevel/P/Q tmpTube yCtr*10
538                                if(V_Flag)
539                                        edge = data_realDistY[0][47]            //y tube 47
540                                        delta = abs(yCtr*10 - edge)
541                                        y_pix[0] = dimY-1 + delta/tube_width            //since the top edge of the B panels is pixel 47               
542                                else
543                                        y_pix[0] = V_LevelX
544                                endif
545
546                        endif
547                endif
548
549// the x-center will be on the panel in this direction         
550                Make/O/D/N=(dimX) tmpTube
551                tmpTube = data_RealDistX[p][0]
552                FindLevel /P/Q tmpTube, xCtr*10
553               
554                x_pix[0] = V_levelX
555                KillWaves/Z tmpTube
556               
557        endif
558               
559        return(0)
560end
561
562// converts from [cm] beam center to pixels
563//
564// the value in pixels is written to the local data folder, NOT to disk (it is recalculated as needed)
565//
566Function V_ConvertBeamCtr_to_pixB(folder,detStr,destPath)
567        String folder,detStr,destPath
568       
569        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
570        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")     
571
572        Variable dimX,dimY,xCtr,yCtr
573        dimX = DimSize(data_realDistX,0)
574        dimY = DimSize(data_realDistX,1)
575       
576        xCtr = V_getDet_beam_center_x(folder,detStr)                    //these are in cm, *10 to get mm
577        yCtr = V_getDet_beam_center_y(folder,detStr)   
578       
579        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
580        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
581        WAVE x_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
582        WAVE y_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
583
584
585// simple wave lookup
586// can't use x_pix[0] = data_RealDistX(xCtr)[0] since the data has no x-scale and (xCtr) is interpreted
587// as a point value
588
589//
590//xCtr, yCtr are in cm, *10 to get mm to compare to distance array
591
592        Make/O/D/N=(dimX) tmpTube
593        tmpTube = data_RealDistX[p][0]
594        FindLevel /P/Q tmpTube, xCtr*10
595       
596        x_pix[0] = V_levelX
597        KillWaves/Z tmpTube
598       
599       
600        Make/O/D/N=(dimY) tmpTube
601        tmpTube = data_RealDistY[0][p]
602        FindLevel /P/Q tmpTube, yCtr*10
603       
604        y_pix[0] = V_levelX
605        KillWaves/Z tmpTube
606               
607        print "pixel ctr B = ",x_pix[0],y_pix[0]
608               
609        return(0)
610end
611
612//
613//
614// TODO
615// -- VERIFY the calculations
616// -- verify where this needs to be done (if the beam center is changed)
617// -- then the q-calculation needs to be re-done
618//
619// -- not much is known about the "B" detector, so this
620//    all hinges on the non-linear corrections being done correctly for that detector
621//
622//      Variable detCtrX, detCtrY
623//      // get the pixel center of the detector (not the beam center)
624//      detCtrX = trunc( DimSize(dataW,0)/2 )           //
625//      detCtrY = trunc( DimSize(dataW,1)/2 )
626//
627//
628Function V_ConvertBeamCtrPix_to_mmB(folder,detStr,destPath)
629        String folder,detStr,destPath
630       
631       
632//      DoAlert 0,"Error - Beam center is being interpreted as pixels, but needs to be in cm. V_ConvertBeamCtrPix_to_mmB()"
633       
634        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
635        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")     
636       
637        Variable xCtr,yCtr
638        xCtr = V_getDet_beam_center_x(folder,detStr)
639        yCtr = V_getDet_beam_center_y(folder,detStr)   
640       
641        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
642        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
643        WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
644        WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
645
646        x_mm[0] = data_realDistX[xCtr][0]
647        y_mm[0] = data_realDistY[0][yCtr]
648               
649        return(0)
650end
651
652
653
654
655
656
657/////
658//
659// non-linear corrections to the tube pixels
660// - returns the distance in mm (although this may change)
661//
662// c0,c1,c2,pix
663// c0-c2 are the fit coefficients
664// pix is the test pixel
665//
666// returns the distance in mm (relative to ctr pixel)
667// ctr is the center pixel, as defined when fitting to quadratic was done
668//
669Function V_TubePixel_to_mm(c0,c1,c2,pix)
670        Variable c0,c1,c2,pix
671       
672        Variable dist
673        dist = c0 + c1*pix + c2*pix*pix
674       
675        return(dist)
676End
677//
678////
679
680
681//
682// TESTING ONLY
683Proc V_MakeFakeCalibrationWaves()
684        // make these in the RAW data folder, before converting to a work folder
685        // - then they will be "found" by get()
686        // -- only for the tube, not the Back det
687       
688//      DoAlert 0, "re-do this and do a better job of filling the fake calibration data"
689
690        DoAlert 0, "Calibration waves are read in from the data file"
691       
692//      V_fMakeFakeCalibrationWaves()
693End
694
695
696
697//
698// TESTING ONLY
699//
700// orientation does not matter, there are 48 tubes in each bank
701// so dimension (3,48) for everything.
702//
703// -- but the orientation does indicate TB vs LR, which has implications for
704//  the (fictional) dimension of the pixel along the tube axis, at least as far
705// as for making the fake coefficients.
706//
707Function V_fMakeFakeCalibrationWaves()
708
709        Variable ii,pixSize
710        String detStr,fname="RAW",orientation
711       
712        for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
713                detStr = StringFromList(ii, ksDetectorListNoB, ";")
714//              Wave w = V_getDetectorDataW(fname,detStr)
715                Make/O/D/N=(3,48) $("root:Packages:NIST:VSANS:RAW:entry:instrument:detector_"+detStr+":spatial_calibration")
716                Wave calib = $("root:Packages:NIST:VSANS:RAW:entry:instrument:detector_"+detStr+":spatial_calibration")
717                // !!!! this overwrites what is there
718
719                orientation = V_getDet_tubeOrientation(fname,detStr)
720                if(cmpstr(orientation,"vertical")==0)
721                //      this is vertical tube data dimensioned as (Ntubes,Npix)
722                        pixSize = 8.4           //V_getDet_y_pixel_size(fname,detStr)
723                       
724                elseif(cmpstr(orientation,"horizontal")==0)
725                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
726                        pixSize = 4                     //V_getDet_x_pixel_size(fname,detStr)
727                       
728                else           
729                        DoAlert 0,"Orientation not correctly passed in NonLinearCorrection(). No correction done."
730                endif
731               
732                calib[0][] = -(128/2)*pixSize                   //approx (n/2)*pixSixe
733                calib[1][] = pixSize
734                calib[2][] = 2e-4
735               
736        endfor
737       
738        // now fake calibration for "B"
739        Wave cal_x = V_getDet_cal_x("RAW","B")
740        Wave cal_y = V_getDet_cal_y("RAW","B")
741       
742        cal_x = .34             // mm, ignore the other 2 values
743        cal_y = .34             // mm
744        return(0)
745End
746
747//
748// (DONE)
749// x- MUST VERIFY the definition of SDD and how (if) setback is written to the data files
750// x- currently I'm assuming that the SDD is the "nominal" value which is correct for the
751//    L/R panels, but is not correct for the T/B panels (must add in the setback)
752//
753//
754//
755// data_realDistX, Y must be previously generated from running NonLinearCorrection()
756//
757// call with:
758// fname as the WORK folder, "RAW"
759// detStr = detector String, "FL"
760// destPath = path to destination WORK folder ("root:Packages:NIST:VSANS:"+folder)
761//
762Function V_Detector_CalcQVals(fname,detStr,destPath)
763        String fname,detStr,destPath
764
765        String orientation
766        Variable xCtr,yCtr,lambda,sdd
767       
768// get all of the geometry information 
769        orientation = V_getDet_tubeOrientation(fname,detStr)
770
771
772        sdd = V_getDet_ActualDistance(fname,detStr)             //sdd derived, including setback [cm]
773
774        // this is the ctr in pixels --xx-- (now it is in cm!)
775//      xCtr = V_getDet_beam_center_x(fname,detStr)
776//      yCtr = V_getDet_beam_center_y(fname,detStr)
777        // this is ctr in mm
778        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
779        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
780        lambda = V_getWavelength(fname)
781        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
782        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
783
784// make the new waves
785        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
786        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
787        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
788        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
789        Wave qTot = $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
790        Wave qx = $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
791        Wave qy = $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
792        Wave qz = $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
793
794// calculate all of the q-values
795// sdd is passed in [cm]
796
797
798// after adding in the 680x1656 back detector, load time was 7.8s, without multithreading
799// with multithreading, 1.9s
800//       qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
801//              qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
802//              qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
803//              qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)   
804
805        MultiThread qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
806        MultiThread     qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
807        MultiThread     qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
808        MultiThread     qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
809       
810        return(0)
811End
812
813
814//function to calculate the overall q-value, given all of the necesary trig inputs
815//
816// (DONE)
817// x- verify the calculation (accuracy - in all input conditions)
818// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
819// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
820//    to each pixel
821//
822//sdd is in [cm]
823// distX and distY are in [mm]
824//wavelength is in Angstroms
825//
826//returned magnitude of Q is in 1/Angstroms
827//
828ThreadSafe Function V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
829        Variable xaxval,yaxval,xctr,yctr,sdd,lam
830        Wave distX,distY
831       
832        Variable dx,dy,qval,two_theta,dist
833               
834
835        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
836        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
837        dist = sqrt(dx^2 + dy^2)
838       
839        dist /= 10  // convert mm to cm
840       
841        two_theta = atan(dist/sdd)
842
843        qval = 4*Pi/lam*sin(two_theta/2)
844       
845        return qval
846End
847
848//calculates just the q-value in the x-direction on the detector
849// (DONE)
850// x- verify the calculation (accuracy - in all input conditions)
851// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
852// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
853//    to each pixel
854//
855//
856// this properly accounts for qz
857//
858ThreadSafe Function V_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
859        Variable xaxval,yaxval,xctr,yctr,sdd,lam
860        Wave distX,distY
861
862        Variable qx,qval,phi,dx,dy,dist,two_theta
863       
864        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
865       
866
867        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
868        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
869        phi = V_FindPhi(dx,dy)
870       
871        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
872        dist = sqrt(dx^2 + dy^2)
873        dist /= 10  // convert mm to cm
874
875        two_theta = atan(dist/sdd)
876
877        qx = qval*cos(two_theta/2)*cos(phi)
878       
879        return qx
880End
881
882//calculates just the q-value in the y-direction on the detector
883// (DONE)
884// x- verify the calculation (accuracy - in all input conditions)
885// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
886// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
887//    to each pixel
888//
889//
890// this properly accounts for qz
891//
892ThreadSafe Function V_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
893        Variable xaxval,yaxval,xctr,yctr,sdd,lam
894        Wave distX,distY
895
896        Variable qy,qval,phi,dx,dy,dist,two_theta
897       
898        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
899       
900
901        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
902        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
903        phi = V_FindPhi(dx,dy)
904       
905        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
906        dist = sqrt(dx^2 + dy^2)
907        dist /= 10  // convert mm to cm
908
909        two_theta = atan(dist/sdd)
910
911        qy = qval*cos(two_theta/2)*sin(phi)
912       
913        return qy
914End
915
916//calculates just the q-value in the z-direction on the detector
917// (DONE)
918// x- verify the calculation (accuracy - in all input conditions)
919// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
920// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
921//    to each pixel
922//
923// not actually used for any calculations, but here for completeness if anyone asks, or for 2D data export
924//
925// this properly accounts for qz, because it is qz
926//
927ThreadSafe Function V_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
928        Variable xaxval,yaxval,xctr,yctr,sdd,lam
929        Wave distX,distY
930
931        Variable qz,qval,phi,dx,dy,dist,two_theta
932       
933        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
934       
935
936        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
937        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
938       
939        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
940        dist = sqrt(dx^2 + dy^2)
941        dist /= 10  // convert mm to cm
942
943        two_theta = atan(dist/sdd)
944
945        qz = qval*sin(two_theta/2)
946       
947        return qz
948End
949
950
951//
952// (DONE)
953// x- VERIFY calculations
954// x- This is the actual solid angle per pixel, not a ratio vs. some "unit SA"
955//    Do I just correct for the different area vs. the "nominal" central area?
956// x- decide how to implement - YES - directly change the data values (as was done in the past)
957//    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
958//    would need this to be applied before exporting)
959// x- do I keep a wave note indicating that this correction has been applied to the data
960//    so that it can be "un-applied"? NO
961// x- do I calculate theta from geometry directly, or get it from Q (Assuming it's present?)
962//    (YES just from geometry, since I need SDD and dx and dy values...)
963//
964//              "B" is passed, so I need to check for "B" and for panel orientation
965// -if it is detector "B" (not tubes), then the normal solid angle correction applies
966// -if it is a tube panel, then I need to know the orientation, to know which angles
967//    and pixel dimensions to use
968//
969// *** UPDATED 1/2020 SRK
970// -using new calculation since the lateral direction of the tubes does not affect the solid angle
971// projection (see He (2015) and John's memo)
972//
973//
974Function V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
975        Wave w,w_err
976        String fname,detStr,destPath
977
978        Variable sdd,xCtr,yCtr,lambda
979        String orientation
980       
981// get all of the geometry information 
982        orientation = V_getDet_tubeOrientation(fname,detStr)
983        sdd = V_getDet_ActualDistance(fname,detStr)
984
985        // this is ctr in mm
986        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
987        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
988        lambda = V_getWavelength(fname)
989       
990        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
991       
992        Wave data_realDistX = data_realDistX
993        Wave data_realDistY = data_realDistY
994
995        Duplicate/O w solid_angle,tmp_theta,tmp_dist,tmp_theta_i        //in the current df
996
997//// calculate the scattering angle
998//      dx = (distX - xctr)             //delta x in mm
999//      dy = (distY - yctr)             //delta y in mm
1000        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
1001       
1002        tmp_dist /= 10  // convert mm to cm
1003        // sdd is in [cm]
1004
1005        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the (total) scattering angle
1006
1007        Variable ii,jj,numx,numy,dx,dy
1008        numx = DimSize(tmp_theta,0)
1009        numy = DimSize(tmp_theta,1)
1010
1011        if(cmpstr(detStr,"B")==0)
1012                //detector B is a grid, straightforward cos^3 solid angle
1013                for(ii=0        ;ii<numx;ii+=1)
1014                        for(jj=0;jj<numy;jj+=1)
1015                               
1016                                if(ii==0)               //do a forward difference if ii==0
1017                                        dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj])        //delta x for the pixel
1018                                else
1019                                        dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj])        //delta x for the pixel
1020                                endif
1021                               
1022                               
1023                                if(jj==0)
1024                                        dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj])        //delta y for the pixel
1025                                else
1026                                        dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1])        //delta y for the pixel
1027                                endif
1028               
1029                                dx /= 10
1030                                dy /= 10                // convert mm to cm (since sdd is in cm)
1031                                solid_angle[ii][jj] = dx*dy             //this is in cm^2
1032                        endfor
1033                endfor
1034               
1035                // to cover up any issues w/negative dx or dy
1036                solid_angle = abs(solid_angle)
1037               
1038                // solid_angle correction
1039                // == dx*dy*cos^3/sdd^2
1040                solid_angle *= (cos(tmp_theta))^3
1041                solid_angle /= sdd^2
1042               
1043                // Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!)
1044                w /= solid_angle
1045
1046                // correctly apply the correction to the error wave (assume a perfect value?)
1047                w_err /= solid_angle            //
1048
1049        else
1050                //             
1051                //different calculation for the tubes, different calculation based on XY orientation
1052                //
1053                if(cmpstr(orientation,"vertical")==0)
1054                        // L/R panels, tube axis is y-direction
1055                        // this is now a different tmp_dist
1056                        // convert everything to cm first!
1057                        // sdd is in [cm], everything else is in [mm]
1058                        tmp_dist = (data_realDistY/10 - yctr/10)/sqrt((data_realDistX/10 - xctr/10)^2 + sdd^2)         
1059                        tmp_theta_i = atan(tmp_dist)            //this is theta_y
1060                       
1061                else
1062                        // horizontal orientation (T/B panels)
1063                        // this is now a different tmp_dist
1064                        // convert everything to cm first!
1065                        // sdd is in [cm], everything else is in [mm]
1066                        tmp_dist = (data_realDistX/10 - xctr/10)/sqrt((data_realDistY/10 - yctr/10)^2 + sdd^2)         
1067                        tmp_theta_i = atan(tmp_dist)            //this is theta_x
1068               
1069                endif
1070               
1071                for(ii=0        ;ii<numx;ii+=1)
1072                        for(jj=0;jj<numy;jj+=1)
1073                               
1074                                if(ii==0)               //do a forward difference if ii==0
1075                                        dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj])        //delta x for the pixel
1076                                else
1077                                        dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj])        //delta x for the pixel
1078                                endif
1079                               
1080                               
1081                                if(jj==0)
1082                                        dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj])        //delta y for the pixel
1083                                else
1084                                        dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1])        //delta y for the pixel
1085                                endif
1086               
1087                                dx /= 10
1088                                dy /= 10                // convert mm to cm (since sdd is in cm)
1089                                solid_angle[ii][jj] = dx*dy             //this is in cm^2
1090                        endfor
1091                endfor
1092               
1093                // to cover up any issues w/negative dx or dy
1094                solid_angle = abs(solid_angle)
1095               
1096                // solid_angle correction
1097                // == dx*dy*cos(th)^2*cos(th_i)/sdd^2           using either the theta_x or theta_y value
1098                solid_angle *= (cos(tmp_theta))^2*cos(tmp_theta_i)
1099                solid_angle /= sdd^2
1100               
1101                // Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!)
1102                w /= solid_angle
1103               
1104                //
1105                // correctly apply the correction to the error wave (assume a perfect value?)
1106                w_err /= solid_angle            //
1107       
1108        endif
1109       
1110
1111// DONE x- clean up after I'm satisfied computations are correct               
1112        KillWaves/Z tmp_theta,tmp_dist,tmp_theta_i
1113       
1114        return(0)
1115end
1116
1117// this is the incorrect solid angle correction that does not take into
1118// account the tube geometry. It is correct for the high-res detector (and the 30m Ordela)
1119//
1120// -- only for testing to prove that the cos(th)^2 *cos(th_i) is correct
1121//
1122Function V_SolidAngleCorrection_COS3(w,w_err,fname,detStr,destPath)
1123        Wave w,w_err
1124        String fname,detStr,destPath
1125
1126        Variable sdd,xCtr,yCtr,lambda
1127        String orientation
1128       
1129// get all of the geometry information 
1130        orientation = V_getDet_tubeOrientation(fname,detStr)
1131        sdd = V_getDet_ActualDistance(fname,detStr)
1132
1133        // this is ctr in mm
1134        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
1135        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
1136        lambda = V_getWavelength(fname)
1137       
1138        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
1139       
1140        Wave data_realDistX = data_realDistX
1141        Wave data_realDistY = data_realDistY
1142
1143        Duplicate/O w solid_angle,tmp_theta,tmp_dist    //in the current df
1144
1145//// calculate the scattering angle
1146//      dx = (distX - xctr)             //delta x in mm
1147//      dy = (distY - yctr)             //delta y in mm
1148        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
1149       
1150        tmp_dist /= 10  // convert mm to cm
1151        // sdd is in [cm]
1152
1153        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the (total) scattering angle
1154
1155        Variable ii,jj,numx,numy,dx,dy
1156        numx = DimSize(tmp_theta,0)
1157        numy = DimSize(tmp_theta,1)
1158
1159        for(ii=0        ;ii<numx;ii+=1)
1160                for(jj=0;jj<numy;jj+=1)
1161                       
1162                        if(ii==0)               //do a forward difference if ii==0
1163                                dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj])        //delta x for the pixel
1164                        else
1165                                dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj])        //delta x for the pixel
1166                        endif
1167                       
1168                       
1169                        if(jj==0)
1170                                dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj])        //delta y for the pixel
1171                        else
1172                                dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1])        //delta y for the pixel
1173                        endif
1174       
1175                        dx /= 10
1176                        dy /= 10                // convert mm to cm (since sdd is in cm)
1177                        solid_angle[ii][jj] = dx*dy             //this is in cm^2
1178                endfor
1179        endfor
1180       
1181        // to cover up any issues w/negative dx or dy
1182        solid_angle = abs(solid_angle)
1183       
1184        // solid_angle correction
1185        // == dx*dy*cos^3/sdd^2
1186        solid_angle *= (cos(tmp_theta))^3
1187        solid_angle /= sdd^2
1188       
1189        // Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!)
1190        w /= solid_angle
1191
1192        // correctly apply the correction to the error wave (assume a perfect value?)
1193        w_err /= solid_angle            //
1194       
1195
1196// DONE x- clean up after I'm satisfied computations are correct               
1197        KillWaves/Z tmp_theta,tmp_dist,tmp_theta_i
1198       
1199        return(0)
1200end
1201
1202
1203
1204//
1205// Large angle transmission correction
1206//
1207// DIVIDE the intensity by this correction to get the right answer
1208//
1209//
1210// Apply the large angle transmssion correction as the data is converted to WORK
1211// so that whether the data is saved as 2D or 1D, the correction has properly been done.
1212//
1213// This is, however, a SAMPLE dependent calculation, not purely instrument geometry.
1214//
1215Function V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
1216        Wave w,w_err
1217        String fname,detStr,destPath
1218
1219        Variable sdd,xCtr,yCtr,trans,trans_err,uval
1220
1221// get all of the geometry information 
1222//      orientation = V_getDet_tubeOrientation(fname,detStr)
1223        sdd = V_getDet_ActualDistance(fname,detStr)
1224
1225        // this is ctr in mm
1226        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
1227        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
1228        trans = V_getSampleTransmission(fname)
1229        trans_err = V_getSampleTransError(fname)
1230       
1231        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
1232       
1233        Wave data_realDistX = data_realDistX
1234        Wave data_realDistY = data_realDistY
1235
1236        Duplicate/O w lat_corr,tmp_theta,tmp_dist,lat_err,tmp_err               //in the current df
1237
1238//// calculate the scattering angle
1239//      dx = (distX - xctr)             //delta x in mm
1240//      dy = (distY - yctr)             //delta y in mm
1241        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
1242       
1243        tmp_dist /= 10  // convert mm to cm
1244        // sdd is in [cm]
1245
1246        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle
1247
1248        Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp
1249        numx = DimSize(tmp_theta,0)
1250        numy = DimSize(tmp_theta,1)
1251       
1252       
1253        //optical thickness
1254        uval = -ln(trans)               //use natural logarithm
1255       
1256        for(ii=0        ;ii<numx;ii+=1)
1257                for(jj=0;jj<numy;jj+=1)
1258                       
1259                        cos_th = cos(tmp_theta[ii][jj])
1260                        arg = (1-cos_th)/cos_th
1261                       
1262                        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
1263                        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120
1264                        // OR
1265                        if((uval<0.01) || (cos_th>0.99))       
1266                                //small arg, approx correction
1267                                lat_corr[ii][jj] = 1-0.5*uval*arg
1268                        else
1269                                //large arg, exact correction
1270                                lat_corr[ii][jj] = (1-exp(-uval*arg))/(uval*arg)
1271                        endif
1272                         
1273                        // (DONE)
1274                        // x- properly calculate and apply the 2D error propagation
1275                        if(trans == 1)
1276                                lat_err[ii][jj] = 0             //no correction, no error
1277                        else
1278                                //sigT, calculated from the Taylor expansion
1279                                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3)
1280                                tmp *= tmp
1281                                tmp *= trans_err^2
1282                                tmp = sqrt(tmp)         //sigT
1283                               
1284                                lat_err[ii][jj] = tmp
1285                        endif
1286                         
1287 
1288                endfor
1289        endfor
1290       
1291
1292       
1293        // Here it is! Apply the correction to the intensity (divide -- to get the proper correction)
1294        w /= lat_corr
1295
1296        // relative errors add in quadrature to the current 2D error
1297        tmp_err = (w_err/lat_corr)^2 + (lat_err/lat_corr)^2*w*w/lat_corr^2
1298        tmp_err = sqrt(tmp_err)
1299       
1300        w_err = tmp_err
1301       
1302
1303        // DONE x- clean up after I'm satisfied computations are correct               
1304        KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err
1305       
1306        return(0)
1307end
1308
1309
1310
1311//
1312//test procedure, not called anymore
1313Proc V_AbsoluteScaling(type,c0,c1,c2,c3,c4,c5,I_err)
1314        String type
1315        Variable c0=1,c1=0.1,c2=0.95,c3=0.1,c4=1,c5=32.0,I_err=0.32
1316        Prompt type,"WORK data type",popup,"CAL;COR;SAM"
1317        Prompt c0, "Sample Transmission"
1318        Prompt c1, "Sample Thickness (cm)"
1319        Prompt c2, "Standard Transmission"
1320        Prompt c3, "Standard Thickness (cm)"
1321        Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)"
1322        Prompt c5, "Standard Cross-Section (cm-1)"
1323        Prompt I_err, "error in I(q=0) (one std dev)"
1324
1325        Variable err
1326        //call the function to do the math
1327        //data from "type" will be scaled and deposited in ABS
1328        err = V_Absolute_Scale(type,c0,c1,c2,c3,c4,c5,I_err)
1329       
1330        if(err)
1331                Abort "Error in V_Absolute_Scale()"
1332        endif
1333       
1334        //contents are always dumped to ABS
1335        type = "ABS"
1336       
1337        //need to update the display with "data" from the correct dataFolder
1338        //reset the current display type to "type"
1339        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
1340        gCurDispType = Type     
1341       
1342        V_UpdateDisplayInformation(Type)
1343       
1344End
1345
1346//
1347//
1348// kappa comes in as s_izero, so be sure to use 1/kappa_err
1349//
1350//convert the "type" data to absolute scale using the given standard information
1351//s_ is the standard
1352//w_ is the "work" file
1353//both are work files and should already be normalized to 10^8 monitor counts
1354Function V_Absolute_Scale(type,absStr)
1355        String type,absStr
1356       
1357       
1358        Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
1359
1360        Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
1361        Variable scale,trans_err
1362        Variable err,ii
1363        String detStr
1364       
1365        // be sure that the starting data exists
1366        err = V_WorkDataExists(type)
1367        if(err==1)
1368                return(err)
1369        endif
1370               
1371        //copy from current dir (type) to ABS
1372        V_CopyHDFToWorkFolder(type,"ABS")       
1373
1374// TODO: -- which monitor to use? Here, I think it should already be normalized to 10^8
1375//     
1376//      w_moncount = V_getMonitorCount(type)            //monitor count in "type"
1377       
1378        w_moncount = V_getBeamMonNormData(type)
1379               
1380        if(w_moncount == 0)
1381                //zero monitor counts will give divide by zero ---
1382                DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
1383                Return(1)               //report error
1384        Endif
1385
1386        w_trans = V_getSampleTransmission(type)         //sample transmission
1387        w_thick = V_getSampleThickness(type)            //sample thickness
1388        trans_err = V_getSampleTransError(type)
1389       
1390       
1391        //get the parames from the list
1392        s_trans = NumberByKey("TSTAND", absStr, "=", ";")       //parse the list of values
1393        s_thick = NumberByKey("DSTAND", absStr, "=", ";")
1394        s_izero = NumberByKey("IZERO", absStr, "=", ";")
1395        s_cross = NumberByKey("XSECT", absStr, "=", ";")
1396        kappa_err = NumberByKey("SDEV", absStr, "=", ";")
1397
1398       
1399        //calculate scale factor
1400        s1 = defmon/w_moncount          // monitor count (s1 should be 1)
1401        s2 = s_thick/w_thick
1402        s3 = s_trans/w_trans
1403        s4 = s_cross/s_izero
1404        scale = s1*s2*s3*s4
1405
1406       
1407        // kappa comes in as s_izero, so be sure to use 1/kappa_err
1408
1409        // and now loop through all of the detectors
1410        //do the actual absolute scaling here, modifying the data in ABS
1411        for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
1412                detStr = StringFromList(ii, ksDetectorListNoB, ";")
1413                Wave data = V_getDetectorDataW("ABS",detStr)
1414                Wave data_err = V_getDetectorDataErrW("ABS",detStr)
1415               
1416                data *= scale
1417                data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1418        endfor
1419
1420        // do the back detector separately, if it is set to be used
1421        NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
1422        if(gIgnoreDetB == 0)
1423                detStr = "B"
1424                Wave data = V_getDetectorDataW("ABS",detStr)
1425                Wave data_err = V_getDetectorDataErrW("ABS",detStr)
1426               
1427                //get the parames from the list
1428                s_trans = NumberByKey("TSTAND_B", absStr, "=", ";")     //parse the list of values
1429                s_thick = NumberByKey("DSTAND_B", absStr, "=", ";")
1430                s_izero = NumberByKey("IZERO_B", absStr, "=", ";")
1431                s_cross = NumberByKey("XSECT_B", absStr, "=", ";")
1432                kappa_err = NumberByKey("SDEV_B", absStr, "=", ";")
1433
1434                //calculate scale factor
1435                s1 = defmon/w_moncount          // monitor count (s1 should be 1)
1436                s2 = s_thick/w_thick
1437                s3 = s_trans/w_trans
1438                s4 = s_cross/s_izero
1439                scale = s1*s2*s3*s4
1440               
1441                data *= scale
1442                data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1443        endif
1444       
1445        //********* 15APR02
1446        // DO NOT correct for atenuators here - the COR step already does this, putting all of the data on equal
1447        // footing (zero atten) before doing the subtraction.
1448       
1449        Return (0) //no error
1450End
1451
1452
1453//
1454// TODO:
1455//   --         DoAlert 0,"This has not yet been updated for VSANS"
1456//
1457//
1458// match the attenuation of the RAW data to the "type" data
1459// so that they can be properly added
1460//
1461// are the attenuator numbers the same? if so exit
1462//
1463// if not, find the attenuator number for type
1464// - find both attenuation factors
1465//
1466// rescale the raw data to match the ratio of the two attenuation factors
1467// -- adjust the detector count (rw)
1468// -- the linear data
1469//
1470//
1471Function V_Adjust_RAW_Attenuation(type)
1472        String type
1473
1474        DoAlert 0,"This has not yet been updated for VSANS"
1475       
1476        WAVE rw=$("root:Packages:NIST:RAW:realsread")
1477        WAVE linear_data=$("root:Packages:NIST:RAW:linear_data")
1478        WAVE data=$("root:Packages:NIST:RAW:data")
1479        WAVE data_err=$("root:Packages:NIST:RAW:linear_data_error")
1480        WAVE/T tw = $("root:Packages:NIST:RAW:textRead")
1481       
1482        WAVE dest_reals=$("root:Packages:NIST:"+type+":realsread")
1483
1484        Variable dest_atten,raw_atten,tol
1485        Variable lambda,raw_atten_err,raw_AttenFactor,dest_attenFactor,dest_atten_err
1486        String fileStr
1487
1488        dest_atten = dest_reals[3]
1489        raw_atten = rw[3]
1490       
1491        tol = 0.1               // within 0.1 atten units is OK
1492        if(abs(dest_atten - raw_atten) < tol )
1493                return(0)
1494        endif
1495
1496        fileStr = tw[3]
1497        lambda = rw[26]
1498        // TODO access correct values
1499        raw_AttenFactor = 1//AttenuationFactor(fileStr,lambda,raw_atten,raw_atten_err)
1500        dest_AttenFactor = 1//AttenuationFactor(fileStr,lambda,dest_atten,dest_atten_err)
1501               
1502        rw[2] *= dest_AttenFactor/raw_AttenFactor
1503        linear_data *= dest_AttenFactor/raw_AttenFactor
1504       
1505        // to keep "data" and linear_data in sync
1506        data = linear_data
1507       
1508        return(0)
1509End
1510
1511//
1512// testing procedure, called from a menu selection
1513//
1514Proc V_DIV_a_Workfile(type)
1515        String type
1516        Prompt type,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
1517       
1518        //macro will take whatever is in SELECTED folder and DIVide it by the current
1519        //contents of the DIV folder - the function will check for existence
1520        //before proceeding
1521
1522        Abort "This has not yet been updated for VSANS"
1523       
1524        Variable err
1525        err = V_DIVCorrection(type)             //returns err = 1 if data doesn't exist in specified folders
1526       
1527        if(err)
1528                Abort "error in V_DIVCorrection()"
1529        endif
1530       
1531        //contents are NOT always dumped to CAL, but are in the new type folder
1532       
1533        String newTitle = "WORK_"+type
1534        DoWindow/F VSANS_Data
1535        DoWindow/T VSANS_Data, newTitle
1536        KillStrings/Z newTitle
1537       
1538        //need to update the display with "data" from the correct dataFolder
1539        //reset the current displaytype to "type"
1540        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=Type
1541       
1542        V_UpdateDisplayInformation(type)
1543       
1544End
1545
1546
1547//
1548// TODO:
1549//   x-         DoAlert 0,"This has not yet been updated for VSANS"
1550//   -- how is the error propagation handled? Be sure it is calculated correctly when DIV is generated
1551//      and is applied correctly here...
1552//
1553//function will divide the contents of "workType" folder with the contents of
1554//the DIV folder + detStr
1555// all data is linear scale for the calculation
1556//
1557Function V_DIVCorrection(data,data_err,detStr,workType)
1558        Wave data,data_err
1559        String detStr,workType
1560       
1561        //check for existence of data in type and DIV
1562        // if the desired data doesn't exist, let the user know, and abort
1563        String destPath=""
1564
1565        NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
1566        if(cmpstr(detStr,"B")==0 && gIgnoreDetB)
1567                return(0)
1568        endif
1569
1570
1571        if(WaveExists(data) == 0)
1572                Print "The data wave does not exist in V_DIVCorrection()"
1573                Return(1)               //error condition
1574        Endif
1575       
1576        //check for DIV
1577        // if the DIV workfile doesn't exist, let the user know,and abort
1578        // !! be sure to check first, before trying to access the wave
1579       
1580//      WAVE/Z div_data = $("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")
1581        if(WaveExists($("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")) == 0)
1582                Print "The DIV wave does not exist in V_DIVCorrection()"
1583                Return(1)               //error condition
1584        Endif
1585        if(WaveExists($("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":linear_data_error")) == 0)
1586                Print "The DIV error wave does not exist in V_DIVCorrection()"
1587                Return(1)               //error condition
1588        Endif
1589        //files exist, proceed
1590
1591        WAVE/Z div_data_err = V_getDetectorDataErrW("DIV",detStr)
1592        WAVE/Z div_data = V_getDetectorDataW("DIV",detStr)
1593
1594
1595
1596// do the error propagation first, since data is changed by the correction
1597        data_err = sqrt(data_err^2/div_data^2 + div_data_err^2 * data^2/div_data^4 )
1598
1599// then the correction
1600        data /= div_data
1601
1602       
1603        Return(0)
1604End
1605
1606
1607//////////////////////////
1608// detector corrections to stitch the back detector into one proper image
1609//
1610//
1611//
1612
1613
1614//
1615// to register the image on the back detector panel
1616//
1617// middle portion (552 pix in Y) is held fixed
1618// top portion of image is shifted right and down
1619// bottom portion of image is shifted right and up
1620//
1621// remainder of image is filled with Zero (NaN causes problems converting to WORK)
1622//
1623// currently, data is not added together and averaged, but it could be
1624//
1625Function V_ShiftBackDetImage(w,adjW)
1626        Wave w,adjW
1627
1628        NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
1629
1630// this is necessary for some old data with the 150x150 back (dummy) panel
1631// the proper back detector has an x-dimension of 680 pixels. Don't do the shift
1632// if the dimensions are incorrect.
1633        if(DimSize(w,0) < 680)
1634                adjW=w
1635                return(0)
1636        endif
1637       
1638        adjW=0
1639               
1640        Variable topX,bottomX
1641        Variable topY,bottomY
1642        Variable totalY,ccdX,ccdY
1643       
1644//      topX = 7
1645//      topY = 105
1646       
1647//      bottomX = 5
1648//      bottomY = 35
1649
1650// TODOHIGHRES
1651// the detector pix dimensions are hard-wired, be sure the are correct
1652        switch(gHighResBinning)
1653                case 1:
1654                        topX = kShift_topX_bin1
1655                        topY = kShift_topY_bin1
1656                        bottomX = kShift_bottomX_bin1
1657                        bottomY = kShift_bottomY_bin1
1658                       
1659                        totalY = 6624   // total YDim
1660                        ccdY = 2208             // = YDim/3
1661                        ccdX = 2720             // = xDim
1662                        break
1663                case 4:
1664                        topX = kShift_topX_bin4
1665                        topY = kShift_topY_bin4
1666                        bottomX = kShift_bottomX_bin4
1667                        bottomY = kShift_bottomY_bin4
1668                       
1669                        totalY = 1656   // total YDim
1670                        ccdY = 552              // = YDim/3
1671                        ccdX = 680              // = xDim
1672
1673                       
1674                        break
1675                default:               
1676                        Abort "No binning case matches in V_ShiftBackDetImage"
1677                       
1678        endswitch
1679
1680                // middle
1681                adjW[][ccdY,ccdY+ccdY] = w[p][q]
1682       
1683                //top
1684                adjW[0+topX,ccdX-1][ccdY+ccdY,totalY-1-topY] = w[p-topX][q+topY]
1685               
1686                //bottom
1687                adjW[0+bottomX,ccdX-1][0+bottomY,ccdY-1] = w[p-bottomX][q-bottomY]
1688
1689       
1690        return(0)
1691End
1692
1693
1694Proc pV_MedianFilterBack(folder)
1695        String folder="RAW"
1696       
1697        V_MedianFilterBack(folder)
1698end
1699
1700Function V_MedianFilterBack(folder)
1701        String folder
1702
1703        Wave w = V_getDetectorDataW(folder,"B")
1704
1705        NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
1706        switch(gHighResBinning)
1707                case 1:
1708                        MatrixFilter /N=11 /P=1 median w                        //              /P=n flag sets the number of passes (default is 1 pass)
1709                       
1710                        Print "*** median noise filter 11x11 applied to the back detector (1 pass) ***"
1711                        break
1712                case 4:
1713                        MatrixFilter /N=3 /P=1 median w                 //              /P=n flag sets the number of passes (default is 1 pass)
1714                       
1715                        Print "*** median noise filter 3x3 applied to the back detector (1 pass) ***"
1716                        break
1717                default:
1718                        Abort "No binning case matches in V_MedianFilterBack"
1719        endswitch
1720
1721        return(0)
1722End
1723
1724
1725Proc pV_SubtractReadNoiseBack(folder,ReadNoise)
1726        String folder="RAW"
1727        Variable readNoise=3160
1728       
1729        V_SubtractReadNoiseBack(folder,readNoise)
1730end
1731
1732Function V_SubtractReadNoiseBack(folder,readNoise)
1733        String folder
1734        Variable readNoise
1735
1736                Wave w = V_getDetectorDataW(folder,"B")
1737                w -= readNoise          // a constant value
1738               
1739//              MatrixFilter /N=3 median w
1740//              Print "*** median noise filter applied to the back detector***"
1741       
1742        return(0)
1743End
1744
1745
1746Proc pV_MedianAndReadNoiseBack(folder,ReadNoise)
1747        String folder="RAW"
1748        Variable readNoise=3160
1749       
1750        V_MedianAndReadNoiseBack(folder,readNoise)
1751end
1752
1753Function V_MedianAndReadNoiseBack(folder,readNoise)
1754        String folder
1755        Variable readNoise
1756
1757                Wave w = V_getDetectorDataW(folder,"B")
1758                w -= readNoise          // a constant value
1759
1760                NVAR gHighResBinning = root:Packages:NIST:VSANS:Globals:gHighResBinning
1761                switch(gHighResBinning)
1762                        case 1:
1763                                MatrixFilter /N=11 /P=1 median w                        //              /P=n flag sets the number of passes (default is 1 pass)
1764                               
1765                                Print "*** median noise filter 11x11 applied to the back detector (1 pass) ***"
1766                                break
1767                        case 4:
1768                                MatrixFilter /N=3 /P=1 median w                 //              /P=n flag sets the number of passes (default is 1 pass)
1769                               
1770                                Print "*** median noise filter 3x3 applied to the back detector (1 pass) ***"
1771                                break
1772                        default:
1773                                Abort "No binning case matches in V_MedianAndReadNoiseBack"
1774                endswitch
1775               
1776        return(0)
1777End
1778
1779
1780
1781////////////////
1782// Detector efficiency and shadowing
1783///////////////
1784
1785//
1786// Tube efficiency + shadowing
1787//
1788//
1789// -- check for the existence of the proper tables (correct wavelength)
1790//  -- generate tables if needed (one-time calculation)
1791//
1792// interpolate the table for the correction - to avoid repeated integration
1793//
1794// store the tables in: root:Packages:NIST:VSANS:Globals:Efficiency:
1795//
1796Function V_TubeEfficiencyShadowCorr(w,w_err,fname,detStr,destPath)
1797        Wave w,w_err
1798        String fname,detStr,destPath
1799
1800        Variable sdd,xCtr,yCtr,lambda
1801        String orientation
1802
1803// if the panel is "B", exit - since it is not tubes, and this should not be called
1804        if(cmpstr(detStr,"B")==0)
1805                return(1)
1806        endif
1807
1808// get all of the geometry information 
1809        orientation = V_getDet_tubeOrientation(fname,detStr)
1810        sdd = V_getDet_ActualDistance(fname,detStr)
1811
1812        // this is ctr in mm
1813        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
1814        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
1815        lambda = V_getWavelength(fname)
1816       
1817        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
1818       
1819        Wave data_realDistX = data_realDistX
1820        Wave data_realDistY = data_realDistY
1821
1822        Duplicate/O w tmp_theta_x,tmp_theta_y,tmp_dist,tmp_corr         //in the current df
1823
1824//// calculate the scattering angles theta_x and theta_y
1825
1826// flip the definitions of x and y for the T/B panels so that x is always lateral WRT the tubes
1827// and y is always along the length of the tubes
1828
1829        if(cmpstr(orientation,"vertical")==0)
1830                // L/R panels, tube axis is y-direction
1831                // this is now a different tmp_dist
1832                // convert everything to cm first!
1833                // sdd is in [cm], everything else is in [mm]
1834                tmp_dist = (data_realDistY/10 - yctr/10)/sqrt((data_realDistX/10 - xctr/10)^2 + sdd^2)         
1835                tmp_theta_y = atan(tmp_dist)            //this is theta_y
1836                tmp_theta_x = atan( (data_realDistX/10 - xctr/10)/sdd )
1837               
1838        else
1839                // horizontal orientation (T/B panels)
1840                // this is now a different tmp_dist
1841                // convert everything to cm first!
1842                // sdd is in [cm], everything else is in [mm]
1843                tmp_dist = (data_realDistX/10 - xctr/10)/sqrt((data_realDistY/10 - yctr/10)^2 + sdd^2)         
1844                tmp_theta_y = atan(tmp_dist)            //this is theta_y, along tube direction
1845                tmp_theta_x = atan( (data_realDistY/10 - yctr/10)/sdd )         // this is laterally across tubes
1846        endif
1847
1848
1849// identify if the 2D efficiency wave has been generated for the data wavelength
1850//
1851// if so, declare
1852// if not, generate
1853
1854        if(WaveExists($"root:Packages:NIST:VSANS:Globals:Efficiency:eff") == 0)
1855                // generate the proper efficiency wave, at lambda
1856                NewDataFolder/O root:Packages:NIST:VSANS:Globals:Efficiency
1857                Print "recalculating efficiency table ..."
1858                V_TubeShadowEfficiencyTables_oneLam(lambda)
1859                // declare the wave
1860                Wave/Z effW = root:Packages:NIST:VSANS:Globals:Efficiency:eff
1861        else
1862                Wave/Z effW = root:Packages:NIST:VSANS:Globals:Efficiency:eff
1863                //is the efficiency at the correct wavelength?
1864                string str=note(effW)
1865//              Print "Note = ",str
1866               
1867                if(V_CloseEnough(lambda,NumberByKey("LAMBDA", str,"="),0.1))            //absolute difference of < 0.1 A
1868                                // yes, proceed, no need to do anything
1869                else
1870                        // no, regenerate the efficiency and then proceed (wave already declared)
1871                        Print "recalculating efficiency table ..."
1872                        V_TubeShadowEfficiencyTables_oneLam(lambda)
1873                endif
1874        endif
1875       
1876       
1877        Variable ii,jj,numx,numy,xAngle,yAngle
1878        numx = DimSize(w,0)
1879        numy = DimSize(w,1)
1880
1881// loop over all of the pixels of the panel and find the interpolated correction (save as a wave)
1882//
1883        for(ii=0        ;ii<numx;ii+=1)
1884                for(jj=0;jj<numy;jj+=1)
1885
1886                        // from the angles, find the (x,y) point to interpolate to get the efficiency
1887               
1888                        xAngle = tmp_theta_x[ii][jj]
1889                        yAngle = tmp_theta_y[ii][jj]
1890
1891                        xAngle = abs(xAngle)
1892                        yAngle = abs(yAngle)
1893                       
1894//                      the x and y scaling of the eff wave (2D) was set when it was generated (in radians)
1895//               simply reading the scaled xy value does not interpolate!!
1896//                      tmp_corr[ii][jj] = effW(xAngle)(yAngle)         // NO, returns "stepped" values
1897                        tmp_corr[ii][jj] = Interp2D(effW,xAngle,yAngle)
1898
1899                endfor
1900        endfor
1901//     
1902//     
1903// apply the correction and calculate the error
1904//     
1905// Here it is! Apply the correction to the intensity (divide -- to get the proper correction)
1906        w /= tmp_corr
1907//
1908// relative errors add in quadrature to the current 2D error
1909// assume that this numerical calculation of efficiency is exact
1910//
1911//      tmp_err = (w_err/tmp_corr)^2 + (lat_err/lat_corr)^2*w*w/lat_corr^2
1912//      tmp_err = sqrt(tmp_err)
1913//     
1914//      w_err = tmp_err
1915//     
1916
1917        // TODO
1918        // - clean up after I'm satisfied computations are correct             
1919//      KillWaves/Z tmp_theta_x,tmp_theta_y,tmp_dist,tmp_err,tmp_corr
1920       
1921        return(0)
1922end
1923
1924
1925
1926// the actual integration of the efficiency for an individual pixel
1927Function V_Efficiency_Integral(pWave,in_u)
1928        Wave pWave
1929        Variable in_u
1930       
1931        Variable lambda,th_x,th_y,u_p,integrand,T_sh,max_x,d_ss,d_He
1932       
1933        lambda = pWave[0]
1934        th_x = pWave[1]
1935        th_y = pWave[2]
1936       
1937        u_p = in_u + kTube_cc * cos(th_x)
1938       
1939        // calculate shadow if th_x > 23.727 deg. th_x is input in radians
1940        max_x = 23.727 / 360 * 2*pi
1941        if(th_x < max_x)
1942                T_sh = 1
1943        else
1944               
1945                // get d_ss
1946                if(abs(u_p) < kTube_ri)
1947                        d_ss = sqrt( (kTube_ri + kTube_ss)^2 - in_u^2 ) - sqrt( kTube_ri^2 - in_u^2     )
1948                elseif (abs(u_p) < (kTube_ri + kTube_ss))
1949                        d_ss = sqrt( (kTube_ri + kTube_ss)^2 - in_u^2 )
1950                else
1951                        d_ss = 0
1952                endif
1953               
1954                // get d_He
1955                if(abs(u_p) < kTube_ri)
1956                        d_He = 2 * sqrt(        kTube_ri^2 - in_u^2     )
1957                else
1958                        d_He = 0
1959                endif
1960               
1961                //calculate T_sh               
1962                T_sh = exp(-2*kSig_ss*lambda*d_ss/cos(th_y)) * exp(-kSig_8b_He*lambda*d_He/cos(th_y))
1963               
1964        endif
1965       
1966       
1967        // calculate the integrand
1968       
1969        //note that the in_u value is used here to find d_ss and d_he (not u_p)
1970        // get d_ss
1971        if(abs(in_u) < kTube_ri)
1972                d_ss = sqrt( (kTube_ri + kTube_ss)^2 - in_u^2 ) - sqrt( kTube_ri^2 - in_u^2     )
1973        elseif (abs(in_u) < (kTube_ri + kTube_ss))
1974                d_ss = sqrt( (kTube_ri + kTube_ss)^2 - in_u^2 )
1975        else
1976                d_ss = 0
1977        endif
1978       
1979        // get d_He
1980        if(abs(in_u) < kTube_ri)
1981                d_He = 2 * sqrt(        kTube_ri^2 - in_u^2     )
1982        else
1983                d_He = 0
1984        endif
1985       
1986        integrand = T_sh*exp(-kSig_ss*lambda*d_ss/cos(th_y))*( 1-exp(-kSig_8b_He*lambda*d_He/cos(th_y)) )
1987
1988        return(integrand)
1989end
1990
1991//
1992// Tube efficiency + shadowing
1993//
1994// function to generate the table for interpolation
1995//
1996// table is generated for a specific wavelength and normalized to eff(lam,0,0)
1997//
1998// below 24 deg (theta_x), there is no shadowing, so the table rows are all identical
1999//
2000// Only one table is stored, and the wavelength of that table is stored in the wave note
2001// -- detector correction checks the note, and recalculates the table if needed
2002// (calculation takes approx 5 seconds)
2003//
2004Function V_TubeShadowEfficiencyTables_oneLam(lambda)
2005        Variable lambda
2006               
2007// storage location for tables
2008        SetDataFolder root:Packages:NIST:VSANS:Globals:Efficiency
2009
2010//make waves that will be filed with the scattering angles and the result of the calculation
2011//
2012
2013//// fill arrays with the scattering angles theta_x and theta_y
2014// 0 < x < 50
2015// 0 < y < 50
2016
2017// *** the definitions of x and y for the T/B panels is flipped so that x is always lateral WRT the tubes
2018// and y is always along the length of the tubes
2019
2020        Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp,normVal
2021        numx = 25
2022        numy = numx
2023
2024        Make/O/D/N=(numx,numy) eff
2025        Make/O/D/N=(numx) theta_x, theta_y,eff_with_shadow,lam_cos
2026       
2027        SetScale x 0,(numx*2)/360*2*pi,"", eff
2028        SetScale y 0,(numy*2)/360*2*pi,"", eff
2029
2030        Note/K eff              // clear the note
2031        Note eff "LAMBDA="+num2str(lambda)
2032       
2033//      theta_x = p*2
2034        theta_y = p     *2      // value range from 0->45, changes if you change numx
2035       
2036        //convert degrees to radians
2037//      theta_x = theta_x/360*2*pi
2038        theta_y = theta_y/360*2*pi
2039
2040//      Make/O/D/N=12 lam_wave
2041//      lam_wave = {0.5,0.7,1,1.5,2,3,4,6,8,10,15,20}
2042       
2043//      Make/O/D/N=(12*numx) eff_withX_to_interp,lam_cos_theta_y
2044//      eff_withX_to_interp=0
2045//      lam_cos_theta_y=0
2046       
2047        Make/O/D/N=3 pWave
2048        pWave[0] = lambda
2049       
2050
2051        for(ii=0        ;ii<numx;ii+=1)
2052
2053                for(jj=0;jj<numx;jj+=1)
2054                               
2055                                pWave[1] = indexToScale(eff,ii,0)               //set theta x
2056                                pWave[2] = indexToScale(eff,jj,1)               //set theta y
2057       
2058                                eff_with_shadow[jj] = Integrate1D(V_Efficiency_Integral,-kTube_ri,kTube_ri,2,0,pWave)           // adaptive Gaussian quadrature
2059                                eff_with_shadow[jj] /= (2*kTube_ri)
2060                               
2061                                eff[ii][jj] = eff_with_shadow[jj]
2062                endfor
2063               
2064                //eff[ii][] = eff_with_shadow[q]
2065        endfor
2066       
2067        lam_cos = lambda/cos(theta_y)
2068       
2069        Sort lam_cos,eff_with_shadow,lam_cos   
2070       
2071//     
2072//////  // value for normalization at current wavelength
2073        pWave[0] = lambda
2074        pWave[1] = 0
2075        pWave[2] = 0
2076////   
2077        normVal = Integrate1D(V_Efficiency_Integral,-kTube_ri,kTube_ri,2,0,pWave)
2078        normVal /= (2*kTube_ri)
2079//     
2080//      print normVal
2081//     
2082        eff_with_shadow /= normVal              // eff(lam,th_x,th_y) / eff(lam,0,0)
2083
2084        eff /= normVal
2085       
2086        // TODO
2087        // - clean up after I'm satisfied computations are correct             
2088//      KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err
2089       
2090        SetDataFolder root:
2091        return(0)
2092end
2093
2094
2095
2096//
2097// Tube efficiency + shadowing
2098//
2099//
2100// TESTING function to generate the tables for interpolation
2101// and various combinations of the corrections for plotting
2102//
2103Function V_TubeShadowEfficiencyTables_withX()
2104
2105
2106        Variable lambda
2107        lambda = 6
2108       
2109        Variable theta_val=3                    //the single theta_x value that is used
2110       
2111// TODO
2112// -- better storage location for tables
2113// bad place for now...
2114        SetDataFolder root:
2115
2116//make waves that will be filed with the scattering angles and the result of the calculation
2117//
2118
2119//// fill arrays with the scattering angles theta_x and theta_y
2120// 0 < x < 50
2121// 0 < y < 50
2122
2123// *** the definitions of x and y for the T/B panels is flipped so that x is always lateral WRT the tubes
2124// and y is always along the length of the tubes
2125
2126        Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp,normVal
2127        numx = 10
2128        numy = 10
2129
2130//      Make/O/D/N=(numx,numy) eff
2131        Make/O/D/N=(numx) theta_x, theta_y,eff_with_shadow,lam_cos
2132       
2133        theta_x = p*5
2134        theta_y = p*5           // value range from 0->45, changes if you change numx
2135       
2136        //convert degrees to radians
2137        theta_x = theta_x/360*2*pi
2138        theta_y = theta_y/360*2*pi
2139
2140        Make/O/D/N=12 lam_wave
2141        lam_wave = {0.5,0.7,1,1.5,2,3,4,6,8,10,15,20}
2142       
2143        Make/O/D/N=(12*numx) eff_withX_to_interp,lam_cos_theta_y
2144        eff_withX_to_interp=0
2145        lam_cos_theta_y=0
2146       
2147        Make/O/D/N=3 pWave
2148
2149        for(jj=0;jj<12;jj+=1)
2150
2151                pWave[0] = lam_wave[jj]
2152
2153                for(ii=0        ;ii<numx;ii+=1)
2154                       
2155                                pWave[1] = theta_val/360*2*pi           //set theta x to any value
2156                                pWave[2] = theta_y[ii]
2157       
2158                                eff_with_shadow[ii] = Integrate1D(V_Efficiency_Integral,-kTube_ri,kTube_ri,2,0,pWave)           // adaptive Gaussian quadrature
2159                                eff_with_shadow[ii] /= (2*kTube_ri)
2160                               
2161                endfor
2162               
2163                lam_cos = lam_wave[jj]/cos(theta_y)
2164
2165// messy indexing for the concatentation               
2166                lam_cos_theta_y[jj*numx,(jj+1)*numx-1] = lam_cos[p-jj*numx]
2167                eff_withX_to_interp[jj*numx,(jj+1)*numx-1] = eff_with_shadow[p-jj*numx]
2168               
2169        endfor
2170       
2171        Sort lam_cos_theta_y,eff_withX_to_interp,lam_cos_theta_y       
2172       
2173//     
2174////////        // value for normalization at what wavelength???
2175//      pWave[0] = 6
2176//      pWave[1] = 0
2177//      pWave[2] = 0
2178////// 
2179//      normVal = Integrate1D(V_Efficiency_Integral,-kTube_ri,kTube_ri,2,0,pWave)
2180//      normVal /= (2*kTube_ri)
2181////   
2182//      print normVal
2183////   
2184//      eff_withX_to_interp /= normVal          // eff(lam,th_x,th_y) / eff(lam,0,0)
2185
2186        // TODO
2187        // - clean up after I'm satisfied computations are correct             
2188//      KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err
2189       
2190        return(0)
2191end
2192
2193
2194
2195/////////////
2196
2197//
2198//
2199// testing function to calculate the correction for the attenuation
2200// of the scattered beam by windows downstream of the sample
2201// (the back window of the sample block, the Si window)
2202//
2203// For implementation, this function could be made identical
2204// to the large angle transmission correction, since the math is
2205// identical - only the Tw value is different (and should ideally be
2206// quite close to 1). With Tw near 1, this would be a few percent correction
2207// at the largest scattering angles.
2208//
2209//
2210Function V_WindowTransmission(tw)
2211        Variable tw
2212       
2213        Make/O/D/N=100 theta,method1,method2,arg
2214       
2215        theta = p/2
2216        theta = theta/360*2*pi          //convert to radians
2217       
2218//      method1 = exp( -ln(tw)/cos(theta) )/tw
2219       
2220        Variable tau
2221        tau = -ln(tw)
2222        arg = (1-cos(theta))/cos(theta)
2223       
2224        if(tau < 0.01)
2225                method2 = 1 - 0.5*tau*arg       
2226        else
2227                method2 = ( 1 - exp(-tau*arg) )/(tau*arg)
2228        endif
2229       
2230        return(0)
2231end
2232
2233
2234//
2235// Large angle transmission correction for the downstream window
2236//
2237// DIVIDE the intensity by this correction to get the right answer
2238//
2239// -- this is a duplication of the math for the large angle
2240// sample tranmission correction. Same situation, but now the
2241// scattered neutrons are attenuated by whatever windows are
2242// downstream of the scattering event.
2243// (the back window of the sample block, the Si window)
2244//
2245// For implementation, this function is made identical
2246// to the large angle transmission correction, since the math is
2247// identical - only the Tw value is different (and should ideally be
2248// quite close to 1). With Tw near 1, this would be a few percent correction
2249// at the largest scattering angles.
2250//
2251//
2252Function V_DownstreamWindowTransmission(w,w_err,fname,detStr,destPath)
2253        Wave w,w_err
2254        String fname,detStr,destPath
2255
2256        Variable sdd,xCtr,yCtr,uval
2257
2258// get all of the geometry information 
2259//      orientation = V_getDet_tubeOrientation(fname,detStr)
2260        sdd = V_getDet_ActualDistance(fname,detStr)
2261
2262        // this is ctr in mm
2263        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
2264        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
2265
2266// get the value of the overall transmission of the downstream components
2267// + error if available.
2268//      trans = V_getSampleTransmission(fname)
2269//      trans_err = V_getSampleTransError(fname)
2270// TODO -- HARD WIRED values, need to set a global or find a place in the header (instrument block?) (reduction?)
2271// currently globals are forced to one in WorkFolderUtils.ipf as the correction is done
2272        NVAR trans = root:Packages:NIST:VSANS:Globals:gDownstreamWinTrans
2273        NVAR trans_err = root:Packages:NIST:VSANS:Globals:gDownstreamWinTransErr       
2274
2275        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
2276       
2277        Wave data_realDistX = data_realDistX
2278        Wave data_realDistY = data_realDistY
2279
2280        Duplicate/O w dwt_corr,tmp_theta,tmp_dist,dwt_err,tmp_err               //in the current df
2281
2282//// calculate the scattering angle
2283//      dx = (distX - xctr)             //delta x in mm
2284//      dy = (distY - yctr)             //delta y in mm
2285        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
2286       
2287        tmp_dist /= 10  // convert mm to cm
2288        // sdd is in [cm]
2289
2290        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle
2291
2292        Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp
2293        numx = DimSize(tmp_theta,0)
2294        numy = DimSize(tmp_theta,1)
2295       
2296       
2297        //optical thickness
2298        uval = -ln(trans)               //use natural logarithm
2299       
2300        for(ii=0        ;ii<numx;ii+=1)
2301                for(jj=0;jj<numy;jj+=1)
2302                       
2303                        cos_th = cos(tmp_theta[ii][jj])
2304                        arg = (1-cos_th)/cos_th
2305                       
2306                        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
2307                        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120
2308                        // OR
2309                        if((uval<0.01) || (cos_th>0.99))       
2310                                //small arg, approx correction
2311                                dwt_corr[ii][jj] = 1-0.5*uval*arg
2312                        else
2313                                //large arg, exact correction
2314                                dwt_corr[ii][jj] = (1-exp(-uval*arg))/(uval*arg)
2315                        endif
2316                         
2317                        // (DONE)
2318                        // x- properly calculate and apply the 2D error propagation
2319                        if(trans == 1)
2320                                dwt_err[ii][jj] = 0             //no correction, no error
2321                        else
2322                                //sigT, calculated from the Taylor expansion
2323                                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3)
2324                                tmp *= tmp
2325                                tmp *= trans_err^2
2326                                tmp = sqrt(tmp)         //sigT
2327                               
2328                                dwt_err[ii][jj] = tmp
2329                        endif
2330                         
2331                endfor
2332        endfor
2333       
2334        // Here it is! Apply the correction to the intensity (divide -- to get the proper correction)
2335        w /= dwt_corr
2336
2337        // relative errors add in quadrature to the current 2D error
2338        tmp_err = (w_err/dwt_corr)^2 + (dwt_err/dwt_corr)^2*w*w/dwt_corr^2
2339        tmp_err = sqrt(tmp_err)
2340       
2341        w_err = tmp_err
2342       
2343        // DONE x- clean up after I'm satisfied computations are correct               
2344        KillWaves/Z tmp_theta,tmp_dist,tmp_err,dwt_err
2345       
2346        return(0)
2347end
2348
Note: See TracBrowser for help on using the repository browser.