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

Last change on this file since 1093 was 1093, checked in by srkline, 5 years ago

New dimensions added for the back detector. many functions neede to be updated to accomodate these changes. Beam center is handled in the same way (in cm, not pixels) as other panels even though this panel is like the 2D detectors on SANS.

Still missing is the real values for caibration, pixel size, dead time, etc. that are yet to be measured.

File size: 44.3 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2#pragma version=1.0
3#pragma IgorVersion=6.1
4
5
6
7//
8// functions to apply corrections to the detector panels
9//
10// these are meant to be called by the procedures that convert "raw" data to
11// "adjusted" or corrected data sets
12//
13
14
15
16
17//
18// detector dead time
19//
20// input is the data array (N tubes x M pixels)
21// input of N x 1 array of dead time values
22//
23// output is the corrected counts in data, overwriting the input data
24//
25// Note that the equation in Roe (eqn 2.15, p. 63) looks different, but it is really the
26// same old equation, just written in a more complex form.
27//
28// (DONE)
29// x- verify the direction of the tubes and indexing
30// x- decide on the appropriate functional form for the tubes
31// x- need count time as input
32// x- be sure I'm working in the right data folder (all waves are passed in)
33// x- clean up when done
34// x- calculate + return the error contribution?
35// x- verify the error propagation
36Function V_DeadTimeCorrectionTubes(dataW,data_errW,dtW,ctTime)
37        Wave dataW,data_errW,dtW
38        Variable ctTime
39       
40        // do I count on the orientation as an input, or do I just figure it out on my own?
41        String orientation
42        Variable dimX,dimY
43        dimX = DimSize(dataW,0)
44        dimY = DimSize(dataw,1)
45        if(dimX > dimY)
46                orientation = "horizontal"
47        else
48                orientation = "vertical"
49        endif
50       
51        // sum the counts in each tube and divide by time for total cr per tube
52        Variable npt
53       
54        if(cmpstr(orientation,"vertical")==0)
55                //      this is data dimensioned as (Ntubes,Npix)
56               
57                MatrixOp/O sumTubes = sumRows(dataW)            // n x 1 result
58                sumTubes /= ctTime              //now count rate per tube
59               
60                dataW[][] = dataW[p][q]/(1-sumTubes[p]*dtW[p])          //correct the data
61                data_errW[][] = data_errW[p][q]/(1-sumTubes[p]*dtW[p])          // propagate the error wave
62
63        elseif(cmpstr(orientation,"horizontal")==0)
64        //      this is data (horizontal) dimensioned as (Npix,Ntubes)
65
66                MatrixOp/O sumTubes = sumCols(dataW)            // 1 x m result
67                sumTubes /= ctTime
68               
69                dataW[][] = dataW[p][q]/(1-sumTubes[q]*dtW[q])
70                data_errW[][] = data_errW[p][q]/(1-sumTubes[q]*dtW[q])
71       
72        else           
73                DoAlert 0,"Orientation not correctly passed in DeadTimeCorrectionTubes(). No correction done."
74        endif
75       
76        return(0)
77end
78
79// test function
80Function V_testDTCor()
81
82        String detStr = ""
83        String fname = "RAW"
84        Variable ctTime
85       
86        detStr = "FR"
87        Wave w = V_getDetectorDataW(fname,detStr)
88        Wave w_err = V_getDetectorDataErrW(fname,detStr)
89        Wave w_dt = V_getDetector_deadtime(fname,detStr)
90
91        ctTime = V_getCount_time(fname)
92       
93//      ctTime = 10
94        V_DeadTimeCorrectionTubes(w,w_err,w_dt,ctTime)
95
96End
97
98
99//
100// Non-linear data correction
101//
102// DOES NOT modify the data, only calculates the spatial relationship
103//
104// input is the data array (N tubes x M pixels)
105// input of N x M array of quadratic coefficients
106//
107// output is wave of corrected real space distance corresponding to each pixel of the data
108// ** its distance from the nominal beam center of (0,0) **
109//
110//
111// (DONE)
112// x- UNITS!!!! currently this is mm, which certainly doesn't match anything else!!!
113//
114// x- verify the direction of the tubes and indexing
115// x- be sure I'm working in the right data folder (it is passed in, and the full path is used)
116// x- clean up when done
117// x- calculate + return the error contribution? (there is none for this operation)
118// x- do I want this to return a wave? (no, default names are generated)
119// x- do I need to write a separate function that returns the distance wave for later calculations?
120// x- do I want to make the distance array 3D to keep the x and y dims together? Calculate them all right now?
121// x- what else do I need to pass to the function? (fname=folder? detStr?)
122// y- (yes,see below) need a separate block or function to handle "B" detector which will be ? different
123//
124//
125Function V_NonLinearCorrection(fname,dataW,coefW,tube_width,detStr,destPath)
126        String fname            //can also be a folder such as "RAW"
127        Wave dataW,coefW
128        Variable tube_width
129        String detStr,destPath
130       
131         
132        // do I count on the orientation as an input, or do I just figure it out on my own?
133        String orientation
134        Variable dimX,dimY
135        dimX = DimSize(dataW,0)
136        dimY = DimSize(dataW,1)
137        if(dimX > dimY)
138                orientation = "horizontal"
139        else
140                orientation = "vertical"
141        endif
142
143        // make a wave of the same dimensions, in the same data folder for the distance
144        // ?? or a 3D wave?
145        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
146        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
147        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
148        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
149       
150        // then per tube, do the quadratic calculation to get the real space distance along the tube
151        // the distance perpendicular to the tube is n*(8.4mm) per tube index
152       
153        // TODO
154        // -- GAP IS HARD-WIRED as constant values
155        Variable offset,gap
156
157// kPanelTouchingGap is in mm   
158// the gap is split equally between the panel pairs
159// (DONE) -- replace all of this with V_getDet_panel_gap(fname,detStr) once it is added to the file
160
161        gap = V_getDet_panel_gap(fname,detStr)
162
163// TODO:
164// -- once the gap fields have been verified, this check can be removed
165// -- it should only apply to data pre-2018 when the field did not exist in the file
166// -- any VSANS data from 2018+ should read gap from the file.
167
168        if(gap < -100)          //-999999 returned if field is missing from file
169       
170                if(cmpstr(detStr,"FL") == 0 || cmpstr(detStr,"FR") == 0)
171                        gap = 3.5               //mm (measured, JB 1/4/18)
172                endif
173                if(cmpstr(detStr,"FT") == 0 || cmpstr(detStr,"FB") == 0)
174                        gap = 3.3               //mm (measured, JB 2/1/18)
175                endif
176                if(cmpstr(detStr,"ML") == 0 || cmpstr(detStr,"MR") == 0)
177                        gap = 5.9               //mm (measured, JB 1/4/18)
178                endif
179                if(cmpstr(detStr,"MT") == 0 || cmpstr(detStr,"MB") == 0)
180                        gap = 18.3              //mm (measured, JB 2/1/18)
181                endif
182       
183        endif
184
185       
186        if(cmpstr(orientation,"vertical")==0)
187                //      this is data dimensioned as (Ntubes,Npix)
188       
189                // adjust the x postion based on the beam center being nominally (0,0) in units of cm, not pixels
190                if(cmpstr(fname,"VCALC")== 0 )
191                        offset = VCALC_getPanelTranslation(detStr)
192                        offset *= 10                    // convert to units of mm
193//                      if(cmpstr("L",detStr[1]) == 0)
194//                              offset *= -1            //negative value for L
195//                      endif
196                else
197                        //normal case
198                        offset = V_getDet_LateralOffset(fname,detStr)
199                        offset *= 10 //convert cm to mm
200                endif
201               
202        // calculation is in mm, not cm
203        // offset will be a negative value for the L panel, and positive for the R panel
204                if(kBCTR_CM)
205                        if(cmpstr("L",detStr[1]) == 0)
206//                              data_realDistX[][] = offset - (dimX - p)*tube_width                     // TODO should this be dimX-1-p = 47-p?
207                                data_realDistX[][] = offset - (dimX - p - 1/2)*tube_width - gap/2               // TODO should this be dimX-1-p = 47-p?
208                        else
209                        //      right
210//                              data_realDistX[][] = tube_width*(p+1) + offset + gap            //add to the Right det,
211                                data_realDistX[][] = tube_width*(p+1/2) + offset + gap/2                //add to the Right det
212                        endif
213                else
214                        data_realDistX[][] = tube_width*(p)
215                endif
216                data_realDistY[][] = coefW[0][p] + coefW[1][p]*q + coefW[2][p]*q*q
217       
218       
219        elseif(cmpstr(orientation,"horizontal")==0)
220                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
221                data_realDistY[][] = tube_width*q
222
223                if(cmpstr(fname,"VCALC")== 0 )
224                        offset = VCALC_getPanelTranslation(detStr)
225                        offset *= 10                    // convert to units of mm
226//                      if(cmpstr("B",detStr[1]) == 0)
227//                              offset *= -1    // negative value for Bottom det
228//                      endif
229                else
230                        //normal case
231                        offset = V_getDet_VerticalOffset(fname,detStr)
232                        offset *= 10 //convert cm to mm
233                endif
234               
235                if(kBCTR_CM)
236                        if(cmpstr("T",detStr[1]) == 0)
237//                              data_realDistY[][] = tube_width*(q+1) + offset + gap                   
238                                data_realDistY[][] = tube_width*(q+1/2) + offset + gap/2                       
239                        else
240                                // bottom
241//                              data_realDistY[][] = offset - (dimY - q)*tube_width     // TODO should this be dimY-1-q = 47-q?
242                                data_realDistY[][] = offset - (dimY - q - 1/2)*tube_width - gap/2       // TODO should this be dimY-1-q = 47-q?
243                        endif
244                else
245                        data_realDistY[][] = tube_width*(q)
246                endif
247                data_realDistX[][] = coefW[0][q] + coefW[1][q]*p + coefW[2][q]*p*p
248
249        else           
250                DoAlert 0,"Orientation not correctly passed in NonLinearCorrection(). No correction done."
251                return(0)
252        endif
253       
254        return(0)
255end
256
257
258
259
260// TODO:
261// -- the cal_x and y coefficients are totally fake
262// -- the wave assignment may not be correct.. so beware
263//
264//
265Function V_NonLinearCorrection_B(folder,dataW,cal_x,cal_y,detStr,destPath)
266        String folder
267        Wave dataW,cal_x,cal_y
268        String detStr,destPath
269
270        if(cmpstr(detStr,"B") != 0)
271                return(0)
272        endif
273
274Print "***Cal_X and Cal_Y for Back are using default values instead of file values ***"
275
276                cal_x[0] = VCALC_getPixSizeX(detStr)*10                 // pixel size in mm  VCALC_getPixSizeX(detStr) is [cm]
277                cal_x[1] = 1
278                cal_x[2] = 10000
279                cal_y[0] = VCALC_getPixSizeY(detStr)*10                 // pixel size in mm  VCALC_getPixSizeX(detStr) is [cm]
280                cal_y[1] = 1
281                cal_y[2] = 10000
282
283       
284        // do I count on the orientation as an input, or do I just figure it out on my own?
285        Variable dimX,dimY
286       
287//      Wave dataW = V_getDetectorDataW(folder,detStr)
288       
289        dimX = DimSize(dataW,0)
290        dimY = DimSize(dataW,1)
291
292        // make a wave of the same dimensions, in the same data folder for the distance
293        // ?? or a 3D wave?
294        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
295        Make/O/D/N=(dimX,dimY) $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
296        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
297        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
298       
299//      Wave cal_x = V_getDet_cal_x(folder,detStr)
300//      Wave cal_y = V_getDet_cal_y(folder,detStr)
301       
302        data_realDistX[][] = cal_x[0]*p
303        data_realDistY[][] = cal_y[0]*q
304       
305        return(0)
306end
307
308
309//
310//
311// TODO
312// -- VERIFY the calculations
313// -- verify where this needs to be done (if the beam center is changed)
314// -- then the q-calculation needs to be re-done
315// -- the position along the tube length is referenced to tube[0], for no particular reason
316//    It may be better to take an average? but [0] is an ASSUMPTION
317// -- distance along tube is simple interpolation, or do I use the coefficients to
318//    calculate the actual value
319//
320// -- distance in the lateral direction is based on tube width, which is a fixed parameter
321//
322//
323Function V_ConvertBeamCtr_to_mm(folder,detStr,destPath)
324        String folder,detStr,destPath
325       
326        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
327        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")     
328
329        String orientation
330        Variable dimX,dimY,xCtr,yCtr
331        dimX = DimSize(data_realDistX,0)
332        dimY = DimSize(data_realDistX,1)
333        if(dimX > dimY)
334                orientation = "horizontal"
335        else
336                orientation = "vertical"
337        endif
338       
339        xCtr = V_getDet_beam_center_x(folder,detStr)
340        yCtr = V_getDet_beam_center_y(folder,detStr)   
341       
342        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
343        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
344        WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
345        WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
346
347        Variable tube_width = V_getDet_tubeWidth(folder,detStr)
348
349//
350        if(cmpstr(orientation,"vertical")==0)
351                //      this is data dimensioned as (Ntubes,Npix)
352//              data_realDistX[][] = tube_width*p
353//              data_realDistY[][] = coefW[0][p] + coefW[1][p]*q + coefW[2][p]*q*q
354                x_mm[0] = tube_width*xCtr
355                y_mm[0] = data_realDistY[0][yCtr]
356        else
357                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
358//              data_realDistX[][] = coefW[0][q] + coefW[1][q]*p + coefW[2][q]*p*p
359//              data_realDistY[][] = tube_width*q
360                x_mm[0] = data_realDistX[xCtr][0]
361                y_mm[0] = tube_width*yCtr
362        endif
363               
364        return(0)
365end
366
367//
368//
369// (DONE)
370// x- VERIFY the calculations
371// x- verify where this needs to be done (if the beam center is changed)
372// x- then the q-calculation needs to be re-done
373// x- the position along the tube length is referenced to tube[0], for no particular reason
374//    It may be better to take an average? but [0] is an ASSUMPTION
375// x- distance along tube is simple interpolation
376//
377// x- distance in the lateral direction is based on tube width, which is a fixed parameter
378//
379// the value in pixels is written to the local data folder, NOT to disk (it is recalculated as needed)
380//
381Function V_ConvertBeamCtr_to_pix(folder,detStr,destPath)
382        String folder,detStr,destPath
383       
384        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
385        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")     
386
387        String orientation
388        Variable dimX,dimY,xCtr,yCtr
389        dimX = DimSize(data_realDistX,0)
390        dimY = DimSize(data_realDistX,1)
391        if(dimX > dimY)
392                orientation = "horizontal"
393        else
394                orientation = "vertical"
395        endif
396       
397        xCtr = V_getDet_beam_center_x(folder,detStr)            //these are in cm
398        yCtr = V_getDet_beam_center_y(folder,detStr)   
399       
400        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
401        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
402        WAVE x_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
403        WAVE y_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
404
405        Variable tube_width = V_getDet_tubeWidth(folder,detStr)
406
407        variable edge,delta
408        Variable gap
409
410// kPanelTouchingGap is in mm   
411// the gap is split equally between the panel pairs
412// TODO -- replace all of this with V_getDet_panel_gap(fname,detStr) once it is added to the file
413// these hard-wired values were determined from 6A and WB beam centers. LR values were exactly the same for
414// both beam considitions (+/- 0.0 mm). FTB was +/- 0.8 mm, MTB +/- 2 mm
415        if(cmpstr(detStr,"FL") == 0 || cmpstr(detStr,"FR") == 0)
416                gap = 3.8               //mm
417        endif
418        if(cmpstr(detStr,"FT") == 0 || cmpstr(detStr,"FB") == 0)
419                gap = 5         //mm
420        endif
421        if(cmpstr(detStr,"ML") == 0 || cmpstr(detStr,"MR") == 0)
422                gap = 5.9               //mm
423        endif
424        if(cmpstr(detStr,"MT") == 0 || cmpstr(detStr,"MB") == 0)
425                gap = 5         //mm
426        endif
427// TODO: this is the line to keep, to replace the hard-wired values
428//      gap = V_getDet_panel_gap(fname,detStr)
429       
430//
431        if(cmpstr(orientation,"vertical")==0)
432                //      this is data dimensioned as (Ntubes,Npix)
433
434                if(kBCTR_CM)
435                        if(cmpstr("L",detStr[1]) == 0)
436                                edge = data_realDistX[47][0]            //tube 47
437                                delta = abs(xCtr*10 - edge)
438                                x_pix[0] = dimX-1 + delta/tube_width
439                        else
440                        // R panel
441                                edge = data_realDistX[0][0]
442                                delta = abs(xCtr*10 - edge + gap)
443                                x_pix[0] = -delta/tube_width            //since the left edge of the R panel is pixel 0
444                        endif
445                endif
446
447                Make/O/D/N=(dimY) tmpTube
448                tmpTube = data_RealDistY[0][p]
449                FindLevel /P/Q tmpTube, yCtr
450               
451                y_pix[0] = V_levelX
452                KillWaves/Z tmpTube
453        else
454                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
455
456                if(kBCTR_CM)
457                        if(cmpstr("T",detStr[1]) == 0)
458                                edge = data_realDistY[0][0]             //tube 0
459                                delta = abs(yCtr*10 - edge + gap)
460                                y_pix[0] =  -delta/tube_width           //since the bottom edge of the T panel is pixel 0
461                        else
462                        // FM(B) panel
463                                edge = data_realDistY[0][47]            //y tube 47
464                                delta = abs(yCtr*10 - edge)
465                                y_pix[0] = dimY-1 + delta/tube_width            //since the top edge of the B panels is pixel 47               
466                        endif
467                endif
468
469               
470                Make/O/D/N=(dimX) tmpTube
471                tmpTube = data_RealDistX[p][0]
472                FindLevel /P/Q tmpTube, xCtr
473               
474                x_pix[0] = V_levelX
475                KillWaves/Z tmpTube
476               
477               
478        endif
479               
480        return(0)
481end
482
483// converts from [cm] beam center to pixels
484//
485// the value in pixels is written to the local data folder, NOT to disk (it is recalculated as needed)
486//
487Function V_ConvertBeamCtr_to_pixB(folder,detStr,destPath)
488        String folder,detStr,destPath
489       
490        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
491        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")     
492
493        Variable dimX,dimY,xCtr,yCtr
494        dimX = DimSize(data_realDistX,0)
495        dimY = DimSize(data_realDistX,1)
496       
497        xCtr = V_getDet_beam_center_x(folder,detStr)                    //these are in cm, *10 to get mm
498        yCtr = V_getDet_beam_center_y(folder,detStr)   
499       
500        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
501        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
502        WAVE x_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_pix")
503        WAVE y_pix = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_pix")
504
505
506// simple wave lookup
507// can't use x_pix[0] = data_RealDistX(xCtr)[0] since the data has no x-scale and (xCtr) is interpreted
508// as a point value
509
510//
511//xCtr, yCtr are in cm, *10 to get mm to compare to distance array
512
513        Make/O/D/N=(dimX) tmpTube
514        tmpTube = data_RealDistX[p][0]
515        FindLevel /P/Q tmpTube, xCtr*10
516       
517        x_pix[0] = V_levelX
518        KillWaves/Z tmpTube
519       
520       
521        Make/O/D/N=(dimY) tmpTube
522        tmpTube = data_RealDistY[0][p]
523        FindLevel /P/Q tmpTube, yCtr*10
524       
525        y_pix[0] = V_levelX
526        KillWaves/Z tmpTube
527               
528        print "pixel ctr B = ",x_pix[0],y_pix[0]
529               
530        return(0)
531end
532
533//
534//
535// TODO
536// -- VERIFY the calculations
537// -- verify where this needs to be done (if the beam center is changed)
538// -- then the q-calculation needs to be re-done
539//
540// -- not much is known about the "B" detector, so this
541//    all hinges on the non-linear corrections being done correctly for that detector
542//
543//      Variable detCtrX, detCtrY
544//      // get the pixel center of the detector (not the beam center)
545//      detCtrX = trunc( DimSize(dataW,0)/2 )           //
546//      detCtrY = trunc( DimSize(dataW,1)/2 )
547//
548//
549Function V_ConvertBeamCtr_to_mmB(folder,detStr,destPath)
550        String folder,detStr,destPath
551       
552        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
553        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")     
554       
555        Variable xCtr,yCtr
556        xCtr = V_getDet_beam_center_x(folder,detStr)
557        yCtr = V_getDet_beam_center_y(folder,detStr)   
558       
559        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
560        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
561        WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
562        WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
563
564        x_mm[0] = data_realDistX[xCtr][0]
565        y_mm[0] = data_realDistY[0][yCtr]
566               
567        return(0)
568end
569
570
571
572
573
574
575/////
576//
577// non-linear corrections to the tube pixels
578// - returns the distance in mm (although this may change)
579//
580// c0,c1,c2,pix
581// c0-c2 are the fit coefficients
582// pix is the test pixel
583//
584// returns the distance in mm (relative to ctr pixel)
585// ctr is the center pixel, as defined when fitting to quadratic was done
586//
587Function V_TubePixel_to_mm(c0,c1,c2,pix)
588        Variable c0,c1,c2,pix
589       
590        Variable dist
591        dist = c0 + c1*pix + c2*pix*pix
592       
593        return(dist)
594End
595//
596////
597
598
599//
600// TESTING ONLY
601Proc V_MakeFakeCalibrationWaves()
602        // make these in the RAW data folder, before converting to a work folder
603        // - then they will be "found" by get()
604        // -- only for the tube, not the Back det
605       
606//      DoAlert 0, "re-do this and do a better job of filling the fake calibration data"
607
608        DoAlert 0, "Calibration waves are read in from the data file"
609       
610//      V_fMakeFakeCalibrationWaves()
611End
612
613
614
615//
616// TESTING ONLY
617//
618// orientation does not matter, there are 48 tubes in each bank
619// so dimension (3,48) for everything.
620//
621// -- but the orientation does indicate TB vs LR, which has implications for
622//  the (fictional) dimension of the pixel along the tube axis, at least as far
623// as for making the fake coefficients.
624//
625Function V_fMakeFakeCalibrationWaves()
626
627        Variable ii,pixSize
628        String detStr,fname="RAW",orientation
629       
630        for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
631                detStr = StringFromList(ii, ksDetectorListNoB, ";")
632//              Wave w = V_getDetectorDataW(fname,detStr)
633                Make/O/D/N=(3,48) $("root:Packages:NIST:VSANS:RAW:entry:instrument:detector_"+detStr+":spatial_calibration")
634                Wave calib = $("root:Packages:NIST:VSANS:RAW:entry:instrument:detector_"+detStr+":spatial_calibration")
635                // !!!! this overwrites what is there
636
637                orientation = V_getDet_tubeOrientation(fname,detStr)
638                if(cmpstr(orientation,"vertical")==0)
639                //      this is vertical tube data dimensioned as (Ntubes,Npix)
640                        pixSize = 8.4           //V_getDet_y_pixel_size(fname,detStr)
641                       
642                elseif(cmpstr(orientation,"horizontal")==0)
643                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
644                        pixSize = 4                     //V_getDet_x_pixel_size(fname,detStr)
645                       
646                else           
647                        DoAlert 0,"Orientation not correctly passed in NonLinearCorrection(). No correction done."
648                endif
649               
650                calib[0][] = -(128/2)*pixSize                   //approx (n/2)*pixSixe
651                calib[1][] = pixSize
652                calib[2][] = 2e-4
653               
654        endfor
655       
656        // now fake calibration for "B"
657        Wave cal_x = V_getDet_cal_x("RAW","B")
658        Wave cal_y = V_getDet_cal_y("RAW","B")
659       
660        cal_x = .34             // mm, ignore the other 2 values
661        cal_y = .34             // mm
662        return(0)
663End
664
665//
666// (DONE)
667// x- MUST VERIFY the definition of SDD and how (if) setback is written to the data files
668// x- currently I'm assuming that the SDD is the "nominal" value which is correct for the
669//    L/R panels, but is not correct for the T/B panels (must add in the setback)
670//
671//
672//
673// data_realDistX, Y must be previously generated from running NonLinearCorrection()
674//
675// call with:
676// fname as the WORK folder, "RAW"
677// detStr = detector String, "FL"
678// destPath = path to destination WORK folder ("root:Packages:NIST:VSANS:"+folder)
679//
680Function V_Detector_CalcQVals(fname,detStr,destPath)
681        String fname,detStr,destPath
682
683        String orientation
684        Variable xCtr,yCtr,lambda,sdd
685       
686// get all of the geometry information 
687        orientation = V_getDet_tubeOrientation(fname,detStr)
688
689
690        sdd = V_getDet_ActualDistance(fname,detStr)             //sdd derived, including setback [cm]
691
692        // this is the ctr in pixels --xx-- (now it is in cm!)
693//      xCtr = V_getDet_beam_center_x(fname,detStr)
694//      yCtr = V_getDet_beam_center_y(fname,detStr)
695        // this is ctr in mm
696        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
697        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
698        lambda = V_getWavelength(fname)
699        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
700        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
701
702// make the new waves
703        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
704        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
705        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
706        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
707        Wave qTot = $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
708        Wave qx = $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
709        Wave qy = $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
710        Wave qz = $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
711
712// calculate all of the q-values
713// sdd is passed in [cm]
714        qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
715        qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
716        qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
717        qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
718       
719       
720        return(0)
721End
722
723
724//function to calculate the overall q-value, given all of the necesary trig inputs
725//
726// (DONE)
727// x- verify the calculation (accuracy - in all input conditions)
728// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
729// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
730//    to each pixel
731//
732//sdd is in [cm]
733// distX and distY are in [mm]
734//wavelength is in Angstroms
735//
736//returned magnitude of Q is in 1/Angstroms
737//
738Function V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
739        Variable xaxval,yaxval,xctr,yctr,sdd,lam
740        Wave distX,distY
741       
742        Variable dx,dy,qval,two_theta,dist
743               
744
745        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
746        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
747        dist = sqrt(dx^2 + dy^2)
748       
749        dist /= 10  // convert mm to cm
750       
751        two_theta = atan(dist/sdd)
752
753        qval = 4*Pi/lam*sin(two_theta/2)
754       
755        return qval
756End
757
758//calculates just the q-value in the x-direction on the detector
759// (DONE)
760// x- verify the calculation (accuracy - in all input conditions)
761// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
762// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
763//    to each pixel
764//
765//
766// this properly accounts for qz
767//
768Function V_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
769        Variable xaxval,yaxval,xctr,yctr,sdd,lam
770        Wave distX,distY
771
772        Variable qx,qval,phi,dx,dy,dist,two_theta
773       
774        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
775       
776
777        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
778        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
779        phi = V_FindPhi(dx,dy)
780       
781        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
782        dist = sqrt(dx^2 + dy^2)
783        dist /= 10  // convert mm to cm
784
785        two_theta = atan(dist/sdd)
786
787        qx = qval*cos(two_theta/2)*cos(phi)
788       
789        return qx
790End
791
792//calculates just the q-value in the y-direction on the detector
793// (DONE)
794// x- verify the calculation (accuracy - in all input conditions)
795// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
796// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
797//    to each pixel
798//
799//
800// this properly accounts for qz
801//
802Function V_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
803        Variable xaxval,yaxval,xctr,yctr,sdd,lam
804        Wave distX,distY
805
806        Variable qy,qval,phi,dx,dy,dist,two_theta
807       
808        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
809       
810
811        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
812        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
813        phi = V_FindPhi(dx,dy)
814       
815        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
816        dist = sqrt(dx^2 + dy^2)
817        dist /= 10  // convert mm to cm
818
819        two_theta = atan(dist/sdd)
820
821        qy = qval*cos(two_theta/2)*sin(phi)
822       
823        return qy
824End
825
826//calculates just the q-value in the z-direction on the detector
827// (DONE)
828// x- verify the calculation (accuracy - in all input conditions)
829// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
830// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
831//    to each pixel
832//
833// not actually used for any calculations, but here for completeness if anyone asks, or for 2D data export
834//
835// this properly accounts for qz, because it is qz
836//
837Function V_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
838        Variable xaxval,yaxval,xctr,yctr,sdd,lam
839        Wave distX,distY
840
841        Variable qz,qval,phi,dx,dy,dist,two_theta
842       
843        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
844       
845
846        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
847        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
848       
849        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
850        dist = sqrt(dx^2 + dy^2)
851        dist /= 10  // convert mm to cm
852
853        two_theta = atan(dist/sdd)
854
855        qz = qval*sin(two_theta/2)
856       
857        return qz
858End
859
860
861//
862// (DONE)
863// x- VERIFY calculations
864// x- This is the actual solid angle per pixel, not a ratio vs. some "unit SA"
865//    Do I just correct for the different area vs. the "nominal" central area?
866// x- decide how to implement - YES - directly change the data values (as was done in the past)
867//    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
868//    would need this to be applied before exporting)
869// x- do I keep a wave note indicating that this correction has been applied to the data
870//    so that it can be "un-applied"? NO
871// x- do I calculate theta from geometry directly, or get it from Q (Assuming it's present?)
872//    (YES just from geometry, since I need SDD and dx and dy values...)
873//
874//
875Function V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
876        Wave w,w_err
877        String fname,detStr,destPath
878
879        Variable sdd,xCtr,yCtr,lambda
880
881// get all of the geometry information 
882//      orientation = V_getDet_tubeOrientation(fname,detStr)
883        sdd = V_getDet_ActualDistance(fname,detStr)
884
885
886        // this is ctr in mm
887        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
888        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
889        lambda = V_getWavelength(fname)
890       
891        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
892       
893        Wave data_realDistX = data_realDistX
894        Wave data_realDistY = data_realDistY
895
896        Duplicate/O w solid_angle,tmp_theta,tmp_dist            //in the current df
897
898//// calculate the scattering angle
899//      dx = (distX - xctr)             //delta x in mm
900//      dy = (distY - yctr)             //delta y in mm
901        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
902       
903        tmp_dist /= 10  // convert mm to cm
904        // sdd is in [cm]
905
906        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle
907
908        Variable ii,jj,numx,numy,dx,dy
909        numx = DimSize(tmp_theta,0)
910        numy = DimSize(tmp_theta,1)
911       
912        for(ii=0        ;ii<numx;ii+=1)
913                for(jj=0;jj<numy;jj+=1)
914                       
915                        if(ii==0)               //do a forward difference if ii==0
916                                dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj])        //delta x for the pixel
917                        else
918                                dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj])        //delta x for the pixel
919                        endif
920                       
921                       
922                        if(jj==0)
923                                dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj])        //delta y for the pixel
924                        else
925                                dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1])        //delta y for the pixel
926                        endif
927       
928                        dx /= 10
929                        dy /= 10                // convert mm to cm (since sdd is in cm)
930                        solid_angle[ii][jj] = dx*dy             //this is in cm^2
931                endfor
932        endfor
933       
934        // to cover up any issues w/negative dx or dy
935        solid_angle = abs(solid_angle)
936       
937        // solid_angle correction
938        // == dx*dy*cos^3/sdd^2
939        solid_angle *= (cos(tmp_theta))^3
940        solid_angle /= sdd^2
941       
942        // Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!)
943        w /= solid_angle
944       
945        //
946        // correctly apply the correction to the error wave (assume a perfect value?)
947        w_err /= solid_angle            //
948
949// DONE x- clean up after I'm satisfied computations are correct               
950        KillWaves/Z tmp_theta,tmp_dist
951       
952        return(0)
953end
954
955
956////////////
957// TODO: all of below is untested code
958//   copied from SANS
959//
960//
961// NOV 2017
962// Currently, this is not called from any VSANS routines. it is only referenced
963// from V_Add_raw_to_work(), which would add two VSANS raw data files together. This has
964// not yet been implemented. I am only keeping this function around to be sure that
965// if/when V_Add_raw_to_work() is implemented, all of the functionality of V_DetCorr() is
966// properly duplicated.
967//
968//
969//
970//performs solid angle and non-linear detector corrections to raw data as it is "added" to a work folder
971//function is called by Raw_to_work() and Add_raw_to_work() functions
972//works on the actual data array, assumes that is is already on LINEAR scale
973//
974Function V_DetCorr(data,data_err,realsread,doEfficiency,doTrans)
975        Wave data,data_err,realsread
976        Variable doEfficiency,doTrans
977
978        DoAlert 0,"This has not yet been updated for VSANS"
979       
980        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,xx0,yy0
981        Variable ii,jj,dtdist,dtdis2
982        Variable xi,xd,yd,rad,ratio,domega,xy
983        Variable lambda,trans,trans_err,lat_err,tmp_err,lat_corr
984       
985//      Print "...doing jacobian and non-linear corrections"
986
987        NVAR pixelsX = root:myGlobals:gNPixelsX
988        NVAR pixelsY = root:myGlobals:gNPixelsY
989       
990        //set up values to send to auxiliary trig functions
991        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
992        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
993
994        x0 = realsread[16]
995        y0 = realsread[17]
996        sx = realsread[10]
997        sx3 = realsread[11]
998        sy = realsread[13]
999        sy3 = realsread[14]
1000       
1001        dtdist = 1000*realsread[18]     //sdd in mm
1002        dtdis2 = dtdist^2
1003       
1004        lambda = realsRead[26]
1005        trans = RealsRead[4]
1006        trans_err = RealsRead[41]               //new, March 2011
1007       
1008
1009        //waves to contain repeated function calls
1010        Make/O/N=(pixelsX) fyy,xx,yy            //Assumes square detector !!!
1011        ii=0
1012        do
1013                xi = ii
1014//              fyy[ii] = dc_fy(ii+1,sy,sy3,ycenter)
1015//              xx[ii] = dc_fxn(ii+1,sx,sx3,xcenter)
1016//              yy[ii] = dc_fym(ii+1,sy,sy3,ycenter)
1017                ii+=1
1018        while(ii<pixelsX)
1019       
1020        Make/O/N=(pixelsX,pixelsY) SolidAngle           // testing only
1021       
1022        ii=0
1023        do
1024                xi = ii
1025//              xd = dc_fx(ii+1,sx,sx3,xcenter)-xx0
1026                jj=0
1027                do
1028                        yd = fyy[jj]-yy0
1029                        //rad is the distance of pixel ij from the sample
1030                        //domega is the ratio of the solid angle of pixel ij versus center pixel
1031                        // product xy = 1 for a detector with a linear spatial response (modern Ordela)
1032                        // solid angle calculated, dW^3 >=1, so multiply data to raise measured values to correct values.
1033                        rad = sqrt(dtdis2 + xd^2 + yd^2)
1034                        domega = rad/dtdist
1035                        ratio = domega^3
1036                        xy = xx[ii]*yy[jj]
1037                       
1038                        data[ii][jj] *= xy*ratio
1039                       
1040                        solidAngle[ii][jj] = xy*ratio           //testing only 
1041                        data_err[ii][jj] *= xy*ratio                    //error propagation assumes that SA and Jacobian are exact, so simply scale error
1042                       
1043                       
1044                        // correction factor for detector efficiency JBG memo det_eff_cor2.doc 3/20/07
1045                        // correction inserted 11/2007 SRK
1046                        // large angle detector efficiency is >= 1 and will "bump up" the measured value at the highest angles
1047                        // so divide here to get the correct answer (5/22/08 SRK)
1048                        if(doEfficiency)
1049//                              data[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
1050//                              data_err[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
1051//                              solidAngle[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)           //testing only
1052                        endif
1053                       
1054                        // large angle transmission calculation is <= 1 and will "bump down" the measured value at the highest angles
1055                        // so divide here to get the correct answer
1056                        if(doTrans)
1057                       
1058                                if(trans<0.1 && ii==0 && jj==0)
1059                                        Print "***transmission is less than 0.1*** and is a significant correction"
1060                                endif
1061                               
1062                                if(trans==0)
1063                                        if(ii==0 && jj==0)
1064                                                Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
1065                                        endif
1066                                        trans = 1
1067                                endif
1068                               
1069                                // pass in the transmission error, and the error in the correction is returned as the last parameter
1070
1071//                              lat_corr = V_LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,lat_err)           //moved from 1D avg SRK 11/2007
1072
1073                                data[ii][jj] /= lat_corr                        //divide by the correction factor
1074                                //
1075                                //
1076                                //
1077                                // relative errors add in quadrature
1078                                tmp_err = (data_err[ii][jj]/lat_corr)^2 + (lat_err/lat_corr)^2*data[ii][jj]*data[ii][jj]/lat_corr^2
1079                                tmp_err = sqrt(tmp_err)
1080                               
1081                                data_err[ii][jj] = tmp_err
1082                               
1083//                              solidAngle[ii][jj] = lat_err
1084
1085                               
1086                                //solidAngle[ii][jj] = LargeAngleTransmissionCorr(trans,dtdist,xd,yd)           //testing only
1087                        endif
1088                       
1089                        jj+=1
1090                while(jj<pixelsX)
1091                ii+=1
1092        while(ii<pixelsX)
1093       
1094        //clean up waves
1095       
1096        Return(0)
1097End
1098
1099
1100//
1101// Large angle transmission correction
1102//
1103// DIVIDE the intensity by this correction to get the right answer
1104//
1105//
1106// Apply the large angle transmssion correction as the data is converted to WORK
1107// so that whether the data is saved as 2D or 1D, the correction has properly been done.
1108//
1109// This is, however, a SAMPLE dependent calculation, not purely instrument geometry.
1110//
1111Function V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
1112        Wave w,w_err
1113        String fname,detStr,destPath
1114
1115        Variable sdd,xCtr,yCtr,trans,trans_err,uval
1116
1117// get all of the geometry information 
1118//      orientation = V_getDet_tubeOrientation(fname,detStr)
1119        sdd = V_getDet_ActualDistance(fname,detStr)
1120
1121        // this is ctr in mm
1122        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
1123        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
1124        trans = V_getSampleTransmission(fname)
1125        trans_err = V_getSampleTransError(fname)
1126       
1127        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
1128       
1129        Wave data_realDistX = data_realDistX
1130        Wave data_realDistY = data_realDistY
1131
1132        Duplicate/O w lat_corr,tmp_theta,tmp_dist,lat_err,tmp_err               //in the current df
1133
1134//// calculate the scattering angle
1135//      dx = (distX - xctr)             //delta x in mm
1136//      dy = (distY - yctr)             //delta y in mm
1137        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
1138       
1139        tmp_dist /= 10  // convert mm to cm
1140        // sdd is in [cm]
1141
1142        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle
1143
1144        Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp
1145        numx = DimSize(tmp_theta,0)
1146        numy = DimSize(tmp_theta,1)
1147       
1148       
1149        //optical thickness
1150        uval = -ln(trans)               //use natural logarithm
1151       
1152        for(ii=0        ;ii<numx;ii+=1)
1153                for(jj=0;jj<numy;jj+=1)
1154                       
1155                        cos_th = cos(tmp_theta[ii][jj])
1156                        arg = (1-cos_th)/cos_th
1157                       
1158                        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
1159                        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120
1160                        // OR
1161                        if((uval<0.01) || (cos_th>0.99))       
1162                                //small arg, approx correction
1163                                lat_corr[ii][jj] = 1-0.5*uval*arg
1164                        else
1165                                //large arg, exact correction
1166                                lat_corr[ii][jj] = (1-exp(-uval*arg))/(uval*arg)
1167                        endif
1168                         
1169                        // (DONE)
1170                        // x- properly calculate and apply the 2D error propagation
1171                        if(trans == 1)
1172                                lat_err[ii][jj] = 0             //no correction, no error
1173                        else
1174                                //sigT, calculated from the Taylor expansion
1175                                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3)
1176                                tmp *= tmp
1177                                tmp *= trans_err^2
1178                                tmp = sqrt(tmp)         //sigT
1179                               
1180                                lat_err[ii][jj] = tmp
1181                        endif
1182                         
1183 
1184                endfor
1185        endfor
1186       
1187
1188       
1189        // Here it is! Apply the correction to the intensity (divide -- to get the proper correction)
1190        w /= lat_corr
1191
1192        // relative errors add in quadrature to the current 2D error
1193        tmp_err = (w_err/lat_corr)^2 + (lat_err/lat_corr)^2*w*w/lat_corr^2
1194        tmp_err = sqrt(tmp_err)
1195       
1196        w_err = tmp_err
1197       
1198
1199        // DONE x- clean up after I'm satisfied computations are correct               
1200        KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err
1201       
1202        return(0)
1203end
1204
1205
1206
1207//
1208//test procedure, not called anymore
1209Proc V_AbsoluteScaling(type,c0,c1,c2,c3,c4,c5,I_err)
1210        String type
1211        Variable c0=1,c1=0.1,c2=0.95,c3=0.1,c4=1,c5=32.0,I_err=0.32
1212        Prompt type,"WORK data type",popup,"CAL;COR;SAM"
1213        Prompt c0, "Sample Transmission"
1214        Prompt c1, "Sample Thickness (cm)"
1215        Prompt c2, "Standard Transmission"
1216        Prompt c3, "Standard Thickness (cm)"
1217        Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)"
1218        Prompt c5, "Standard Cross-Section (cm-1)"
1219        Prompt I_err, "error in I(q=0) (one std dev)"
1220
1221        Variable err
1222        //call the function to do the math
1223        //data from "type" will be scaled and deposited in ABS
1224        err = V_Absolute_Scale(type,c0,c1,c2,c3,c4,c5,I_err)
1225       
1226        if(err)
1227                Abort "Error in V_Absolute_Scale()"
1228        endif
1229       
1230        //contents are always dumped to ABS
1231        type = "ABS"
1232       
1233        //need to update the display with "data" from the correct dataFolder
1234        //reset the current display type to "type"
1235        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
1236        gCurDispType = Type     
1237       
1238        V_UpdateDisplayInformation(Type)
1239       
1240End
1241
1242//
1243//
1244// kappa comes in as s_izero, so be sure to use 1/kappa_err
1245//
1246//convert the "type" data to absolute scale using the given standard information
1247//s_ is the standard
1248//w_ is the "work" file
1249//both are work files and should already be normalized to 10^8 monitor counts
1250Function V_Absolute_Scale(type,w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err)
1251        String type
1252        Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
1253
1254
1255        Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
1256        Variable scale,trans_err
1257        Variable err,ii
1258        String detStr
1259       
1260        // be sure that the starting data exists
1261        err = V_WorkDataExists(type)
1262        if(err==1)
1263                return(err)
1264        endif
1265               
1266        //copy from current dir (type) to ABS
1267        V_CopyHDFToWorkFolder(type,"ABS")       
1268
1269// TODO: -- which monitor to use? Here, I think it should already be normalized to 10^8
1270//     
1271//      w_moncount = V_getMonitorCount(type)            //monitor count in "type"
1272       
1273        w_moncount = V_getBeamMonNormData(type)
1274       
1275       
1276        if(w_moncount == 0)
1277                //zero monitor counts will give divide by zero ---
1278                DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
1279                Return(1)               //report error
1280        Endif
1281       
1282        //calculate scale factor
1283        s1 = defmon/w_moncount          // monitor count (s1 should be 1)
1284        s2 = s_thick/w_thick
1285        s3 = s_trans/w_trans
1286        s4 = s_cross/s_izero
1287        scale = s1*s2*s3*s4
1288
1289        trans_err = V_getSampleTransError(type)
1290       
1291        // kappa comes in as s_izero, so be sure to use 1/kappa_err
1292
1293        // and now loop through all of the detectors
1294        //do the actual absolute scaling here, modifying the data in ABS
1295        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
1296                detStr = StringFromList(ii, ksDetectorListAll, ";")
1297                Wave data = V_getDetectorDataW("ABS",detStr)
1298                Wave data_err = V_getDetectorDataErrW("ABS",detStr)
1299               
1300                data *= scale
1301                data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1302        endfor
1303       
1304        //********* 15APR02
1305        // DO NOT correct for atenuators here - the COR step already does this, putting all of the data on equal
1306        // footing (zero atten) before doing the subtraction.
1307       
1308        Return (0) //no error
1309End
1310
1311
1312//
1313// TODO:
1314//   --         DoAlert 0,"This has not yet been updated for VSANS"
1315//
1316//
1317// match the attenuation of the RAW data to the "type" data
1318// so that they can be properly added
1319//
1320// are the attenuator numbers the same? if so exit
1321//
1322// if not, find the attenuator number for type
1323// - find both attenuation factors
1324//
1325// rescale the raw data to match the ratio of the two attenuation factors
1326// -- adjust the detector count (rw)
1327// -- the linear data
1328//
1329//
1330Function V_Adjust_RAW_Attenuation(type)
1331        String type
1332
1333        DoAlert 0,"This has not yet been updated for VSANS"
1334       
1335        WAVE rw=$("root:Packages:NIST:RAW:realsread")
1336        WAVE linear_data=$("root:Packages:NIST:RAW:linear_data")
1337        WAVE data=$("root:Packages:NIST:RAW:data")
1338        WAVE data_err=$("root:Packages:NIST:RAW:linear_data_error")
1339        WAVE/T tw = $("root:Packages:NIST:RAW:textRead")
1340       
1341        WAVE dest_reals=$("root:Packages:NIST:"+type+":realsread")
1342
1343        Variable dest_atten,raw_atten,tol
1344        Variable lambda,raw_atten_err,raw_AttenFactor,dest_attenFactor,dest_atten_err
1345        String fileStr
1346
1347        dest_atten = dest_reals[3]
1348        raw_atten = rw[3]
1349       
1350        tol = 0.1               // within 0.1 atten units is OK
1351        if(abs(dest_atten - raw_atten) < tol )
1352                return(0)
1353        endif
1354
1355        fileStr = tw[3]
1356        lambda = rw[26]
1357        // TODO access correct values
1358        raw_AttenFactor = 1//AttenuationFactor(fileStr,lambda,raw_atten,raw_atten_err)
1359        dest_AttenFactor = 1//AttenuationFactor(fileStr,lambda,dest_atten,dest_atten_err)
1360               
1361        rw[2] *= dest_AttenFactor/raw_AttenFactor
1362        linear_data *= dest_AttenFactor/raw_AttenFactor
1363       
1364        // to keep "data" and linear_data in sync
1365        data = linear_data
1366       
1367        return(0)
1368End
1369
1370//
1371// testing procedure, called from a menu selection
1372//
1373Proc V_DIV_a_Workfile(type)
1374        String type
1375        Prompt type,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
1376       
1377        //macro will take whatever is in SELECTED folder and DIVide it by the current
1378        //contents of the DIV folder - the function will check for existence
1379        //before proceeding
1380
1381        Abort "This has not yet been updated for VSANS"
1382       
1383        Variable err
1384        err = V_DIVCorrection(type)             //returns err = 1 if data doesn't exist in specified folders
1385       
1386        if(err)
1387                Abort "error in V_DIVCorrection()"
1388        endif
1389       
1390        //contents are NOT always dumped to CAL, but are in the new type folder
1391       
1392        String newTitle = "WORK_"+type
1393        DoWindow/F VSANS_Data
1394        DoWindow/T VSANS_Data, newTitle
1395        KillStrings/Z newTitle
1396       
1397        //need to update the display with "data" from the correct dataFolder
1398        //reset the current displaytype to "type"
1399        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=Type
1400       
1401        V_UpdateDisplayInformation(type)
1402       
1403End
1404
1405
1406//
1407// TODO:
1408//   x-         DoAlert 0,"This has not yet been updated for VSANS"
1409//   -- how is the error propagation handled? Be sure it is calculated correctly when DIV is generated
1410//      and is applied correctly here...
1411//
1412//function will divide the contents of "workType" folder with the contents of
1413//the DIV folder + detStr
1414// all data is linear scale for the calculation
1415//
1416Function V_DIVCorrection(data,data_err,detStr,workType)
1417        Wave data,data_err
1418        String detStr,workType
1419       
1420        //check for existence of data in type and DIV
1421        // if the desired data doesn't exist, let the user know, and abort
1422        String destPath=""
1423
1424        if(WaveExists(data) == 0)
1425                Print "The data wave does not exist in V_DIVCorrection()"
1426                Return(1)               //error condition
1427        Endif
1428       
1429        //check for DIV
1430        // if the DIV workfile doesn't exist, let the user know,and abort
1431        // !! be sure to check first, before trying to access the wave
1432       
1433//      WAVE/Z div_data = $("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")
1434        if(WaveExists($("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")) == 0)
1435                Print "The DIV wave does not exist in V_DIVCorrection()"
1436                Return(1)               //error condition
1437        Endif
1438        if(WaveExists($("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":linear_data_error")) == 0)
1439                Print "The DIV error wave does not exist in V_DIVCorrection()"
1440                Return(1)               //error condition
1441        Endif
1442        //files exist, proceed
1443
1444        WAVE/Z div_data_err = V_getDetectorDataErrW("DIV",detStr)
1445        WAVE/Z div_data = V_getDetectorDataW("DIV",detStr)
1446
1447
1448
1449// do the error propagation first, since data is changed by the correction
1450        data_err = sqrt(data_err^2/div_data^2 + div_data_err^2 * data^2/div_data^4 )
1451
1452// then the correction
1453        data /= div_data
1454
1455       
1456        Return(0)
1457End
1458
1459
1460//////////////////////////
Note: See TracBrowser for help on using the repository browser.