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

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

additions to the reduction for new corrections to the tube detectors for the angle dependent efficiency and the angle dependent tube shadowing. Calculations have been added, new folder for the efficiency has been added, and preference checkboxes have been updated.

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