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

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

Added two model functions for white beam smearing.

Many other small changes for processing of the back detector, shuffling of VSANS menu items, and consistent naming of V_ procedures.

File size: 44.9 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       
553        DoAlert 0,"Error - Beam center is being interpreted as pixels, but needs to be in cm. V_ConvertBeamCtr_to_mmB()"
554       
555        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
556        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")     
557       
558        Variable xCtr,yCtr
559        xCtr = V_getDet_beam_center_x(folder,detStr)
560        yCtr = V_getDet_beam_center_y(folder,detStr)   
561       
562        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
563        Make/O/D/N=1 $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
564        WAVE x_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_x_mm")
565        WAVE y_mm = $(destPath + ":entry:instrument:detector_"+detStr+":beam_center_y_mm")
566
567        x_mm[0] = data_realDistX[xCtr][0]
568        y_mm[0] = data_realDistY[0][yCtr]
569               
570        return(0)
571end
572
573
574
575
576
577
578/////
579//
580// non-linear corrections to the tube pixels
581// - returns the distance in mm (although this may change)
582//
583// c0,c1,c2,pix
584// c0-c2 are the fit coefficients
585// pix is the test pixel
586//
587// returns the distance in mm (relative to ctr pixel)
588// ctr is the center pixel, as defined when fitting to quadratic was done
589//
590Function V_TubePixel_to_mm(c0,c1,c2,pix)
591        Variable c0,c1,c2,pix
592       
593        Variable dist
594        dist = c0 + c1*pix + c2*pix*pix
595       
596        return(dist)
597End
598//
599////
600
601
602//
603// TESTING ONLY
604Proc V_MakeFakeCalibrationWaves()
605        // make these in the RAW data folder, before converting to a work folder
606        // - then they will be "found" by get()
607        // -- only for the tube, not the Back det
608       
609//      DoAlert 0, "re-do this and do a better job of filling the fake calibration data"
610
611        DoAlert 0, "Calibration waves are read in from the data file"
612       
613//      V_fMakeFakeCalibrationWaves()
614End
615
616
617
618//
619// TESTING ONLY
620//
621// orientation does not matter, there are 48 tubes in each bank
622// so dimension (3,48) for everything.
623//
624// -- but the orientation does indicate TB vs LR, which has implications for
625//  the (fictional) dimension of the pixel along the tube axis, at least as far
626// as for making the fake coefficients.
627//
628Function V_fMakeFakeCalibrationWaves()
629
630        Variable ii,pixSize
631        String detStr,fname="RAW",orientation
632       
633        for(ii=0;ii<ItemsInList(ksDetectorListNoB);ii+=1)
634                detStr = StringFromList(ii, ksDetectorListNoB, ";")
635//              Wave w = V_getDetectorDataW(fname,detStr)
636                Make/O/D/N=(3,48) $("root:Packages:NIST:VSANS:RAW:entry:instrument:detector_"+detStr+":spatial_calibration")
637                Wave calib = $("root:Packages:NIST:VSANS:RAW:entry:instrument:detector_"+detStr+":spatial_calibration")
638                // !!!! this overwrites what is there
639
640                orientation = V_getDet_tubeOrientation(fname,detStr)
641                if(cmpstr(orientation,"vertical")==0)
642                //      this is vertical tube data dimensioned as (Ntubes,Npix)
643                        pixSize = 8.4           //V_getDet_y_pixel_size(fname,detStr)
644                       
645                elseif(cmpstr(orientation,"horizontal")==0)
646                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
647                        pixSize = 4                     //V_getDet_x_pixel_size(fname,detStr)
648                       
649                else           
650                        DoAlert 0,"Orientation not correctly passed in NonLinearCorrection(). No correction done."
651                endif
652               
653                calib[0][] = -(128/2)*pixSize                   //approx (n/2)*pixSixe
654                calib[1][] = pixSize
655                calib[2][] = 2e-4
656               
657        endfor
658       
659        // now fake calibration for "B"
660        Wave cal_x = V_getDet_cal_x("RAW","B")
661        Wave cal_y = V_getDet_cal_y("RAW","B")
662       
663        cal_x = .34             // mm, ignore the other 2 values
664        cal_y = .34             // mm
665        return(0)
666End
667
668//
669// (DONE)
670// x- MUST VERIFY the definition of SDD and how (if) setback is written to the data files
671// x- currently I'm assuming that the SDD is the "nominal" value which is correct for the
672//    L/R panels, but is not correct for the T/B panels (must add in the setback)
673//
674//
675//
676// data_realDistX, Y must be previously generated from running NonLinearCorrection()
677//
678// call with:
679// fname as the WORK folder, "RAW"
680// detStr = detector String, "FL"
681// destPath = path to destination WORK folder ("root:Packages:NIST:VSANS:"+folder)
682//
683Function V_Detector_CalcQVals(fname,detStr,destPath)
684        String fname,detStr,destPath
685
686        String orientation
687        Variable xCtr,yCtr,lambda,sdd
688       
689// get all of the geometry information 
690        orientation = V_getDet_tubeOrientation(fname,detStr)
691
692
693        sdd = V_getDet_ActualDistance(fname,detStr)             //sdd derived, including setback [cm]
694
695        // this is the ctr in pixels --xx-- (now it is in cm!)
696//      xCtr = V_getDet_beam_center_x(fname,detStr)
697//      yCtr = V_getDet_beam_center_y(fname,detStr)
698        // this is ctr in mm
699        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
700        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
701        lambda = V_getWavelength(fname)
702        Wave data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
703        Wave data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
704
705// make the new waves
706        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
707        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
708        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
709        Duplicate/O data_realDistX $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
710        Wave qTot = $(destPath + ":entry:instrument:detector_"+detStr+":qTot_"+detStr)
711        Wave qx = $(destPath + ":entry:instrument:detector_"+detStr+":qx_"+detStr)
712        Wave qy = $(destPath + ":entry:instrument:detector_"+detStr+":qy_"+detStr)
713        Wave qz = $(destPath + ":entry:instrument:detector_"+detStr+":qz_"+detStr)
714
715// calculate all of the q-values
716// sdd is passed in [cm]
717
718
719// after adding in the 680x1656 back detector, load time was 7.8s, without multithreading
720// with multithreading, 1.9s
721//       qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
722//              qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
723//              qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
724//              qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)   
725
726        MultiThread qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
727        MultiThread     qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
728        MultiThread     qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
729        MultiThread     qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lambda,data_realDistX,data_realDistY)
730       
731        return(0)
732End
733
734
735//function to calculate the overall q-value, given all of the necesary trig inputs
736//
737// (DONE)
738// x- verify the calculation (accuracy - in all input conditions)
739// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
740// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
741//    to each pixel
742//
743//sdd is in [cm]
744// distX and distY are in [mm]
745//wavelength is in Angstroms
746//
747//returned magnitude of Q is in 1/Angstroms
748//
749ThreadSafe Function V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
750        Variable xaxval,yaxval,xctr,yctr,sdd,lam
751        Wave distX,distY
752       
753        Variable dx,dy,qval,two_theta,dist
754               
755
756        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
757        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
758        dist = sqrt(dx^2 + dy^2)
759       
760        dist /= 10  // convert mm to cm
761       
762        two_theta = atan(dist/sdd)
763
764        qval = 4*Pi/lam*sin(two_theta/2)
765       
766        return qval
767End
768
769//calculates just the q-value in the x-direction on the detector
770// (DONE)
771// x- verify the calculation (accuracy - in all input conditions)
772// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
773// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
774//    to each pixel
775//
776//
777// this properly accounts for qz
778//
779ThreadSafe Function V_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
780        Variable xaxval,yaxval,xctr,yctr,sdd,lam
781        Wave distX,distY
782
783        Variable qx,qval,phi,dx,dy,dist,two_theta
784       
785        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
786       
787
788        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
789        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
790        phi = V_FindPhi(dx,dy)
791       
792        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
793        dist = sqrt(dx^2 + dy^2)
794        dist /= 10  // convert mm to cm
795
796        two_theta = atan(dist/sdd)
797
798        qx = qval*cos(two_theta/2)*cos(phi)
799       
800        return qx
801End
802
803//calculates just the q-value in the y-direction on the detector
804// (DONE)
805// x- verify the calculation (accuracy - in all input conditions)
806// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
807// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
808//    to each pixel
809//
810//
811// this properly accounts for qz
812//
813ThreadSafe Function V_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
814        Variable xaxval,yaxval,xctr,yctr,sdd,lam
815        Wave distX,distY
816
817        Variable qy,qval,phi,dx,dy,dist,two_theta
818       
819        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
820       
821
822        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
823        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
824        phi = V_FindPhi(dx,dy)
825       
826        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
827        dist = sqrt(dx^2 + dy^2)
828        dist /= 10  // convert mm to cm
829
830        two_theta = atan(dist/sdd)
831
832        qy = qval*cos(two_theta/2)*sin(phi)
833       
834        return qy
835End
836
837//calculates just the q-value in the z-direction on the detector
838// (DONE)
839// x- verify the calculation (accuracy - in all input conditions)
840// x- verify the units of everything here, it's currently all jumbled and wrong... and repeated
841// x- the input data_realDistX and Y are essentially lookup tables of the real space distance corresponding
842//    to each pixel
843//
844// not actually used for any calculations, but here for completeness if anyone asks, or for 2D data export
845//
846// this properly accounts for qz, because it is qz
847//
848ThreadSafe Function V_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
849        Variable xaxval,yaxval,xctr,yctr,sdd,lam
850        Wave distX,distY
851
852        Variable qz,qval,phi,dx,dy,dist,two_theta
853       
854        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,distX,distY)
855       
856
857        dx = (distX[xaxval][yaxval] - xctr)             //delta x in mm
858        dy = (distY[xaxval][yaxval] - yctr)             //delta y in mm
859       
860        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
861        dist = sqrt(dx^2 + dy^2)
862        dist /= 10  // convert mm to cm
863
864        two_theta = atan(dist/sdd)
865
866        qz = qval*sin(two_theta/2)
867       
868        return qz
869End
870
871
872//
873// (DONE)
874// x- VERIFY calculations
875// x- This is the actual solid angle per pixel, not a ratio vs. some "unit SA"
876//    Do I just correct for the different area vs. the "nominal" central area?
877// x- decide how to implement - YES - directly change the data values (as was done in the past)
878//    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
879//    would need this to be applied before exporting)
880// x- do I keep a wave note indicating that this correction has been applied to the data
881//    so that it can be "un-applied"? NO
882// x- do I calculate theta from geometry directly, or get it from Q (Assuming it's present?)
883//    (YES just from geometry, since I need SDD and dx and dy values...)
884//
885//
886Function V_SolidAngleCorrection(w,w_err,fname,detStr,destPath)
887        Wave w,w_err
888        String fname,detStr,destPath
889
890        Variable sdd,xCtr,yCtr,lambda
891
892// get all of the geometry information 
893//      orientation = V_getDet_tubeOrientation(fname,detStr)
894        sdd = V_getDet_ActualDistance(fname,detStr)
895
896
897        // this is ctr in mm
898        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
899        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
900        lambda = V_getWavelength(fname)
901       
902        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
903       
904        Wave data_realDistX = data_realDistX
905        Wave data_realDistY = data_realDistY
906
907        Duplicate/O w solid_angle,tmp_theta,tmp_dist            //in the current df
908
909//// calculate the scattering angle
910//      dx = (distX - xctr)             //delta x in mm
911//      dy = (distY - yctr)             //delta y in mm
912        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
913       
914        tmp_dist /= 10  // convert mm to cm
915        // sdd is in [cm]
916
917        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle
918
919        Variable ii,jj,numx,numy,dx,dy
920        numx = DimSize(tmp_theta,0)
921        numy = DimSize(tmp_theta,1)
922       
923        for(ii=0        ;ii<numx;ii+=1)
924                for(jj=0;jj<numy;jj+=1)
925                       
926                        if(ii==0)               //do a forward difference if ii==0
927                                dx = (data_realDistX[ii+1][jj] - data_realDistX[ii][jj])        //delta x for the pixel
928                        else
929                                dx = (data_realDistX[ii][jj] - data_realDistX[ii-1][jj])        //delta x for the pixel
930                        endif
931                       
932                       
933                        if(jj==0)
934                                dy = (data_realDistY[ii][jj+1] - data_realDistY[ii][jj])        //delta y for the pixel
935                        else
936                                dy = (data_realDistY[ii][jj] - data_realDistY[ii][jj-1])        //delta y for the pixel
937                        endif
938       
939                        dx /= 10
940                        dy /= 10                // convert mm to cm (since sdd is in cm)
941                        solid_angle[ii][jj] = dx*dy             //this is in cm^2
942                endfor
943        endfor
944       
945        // to cover up any issues w/negative dx or dy
946        solid_angle = abs(solid_angle)
947       
948        // solid_angle correction
949        // == dx*dy*cos^3/sdd^2
950        solid_angle *= (cos(tmp_theta))^3
951        solid_angle /= sdd^2
952       
953        // Here it is! Apply the correction to the intensity (I divide -- to get the counts per solid angle!!)
954        w /= solid_angle
955       
956        //
957        // correctly apply the correction to the error wave (assume a perfect value?)
958        w_err /= solid_angle            //
959
960// DONE x- clean up after I'm satisfied computations are correct               
961        KillWaves/Z tmp_theta,tmp_dist
962       
963        return(0)
964end
965
966
967////////////
968// TODO: all of below is untested code
969//   copied from SANS
970//
971//
972// NOV 2017
973// Currently, this is not called from any VSANS routines. it is only referenced
974// from V_Add_raw_to_work(), which would add two VSANS raw data files together. This has
975// not yet been implemented. I am only keeping this function around to be sure that
976// if/when V_Add_raw_to_work() is implemented, all of the functionality of V_DetCorr() is
977// properly duplicated.
978//
979//
980//
981//performs solid angle and non-linear detector corrections to raw data as it is "added" to a work folder
982//function is called by Raw_to_work() and Add_raw_to_work() functions
983//works on the actual data array, assumes that is is already on LINEAR scale
984//
985Function V_DetCorr(data,data_err,realsread,doEfficiency,doTrans)
986        Wave data,data_err,realsread
987        Variable doEfficiency,doTrans
988
989        DoAlert 0,"This has not yet been updated for VSANS"
990       
991        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,xx0,yy0
992        Variable ii,jj,dtdist,dtdis2
993        Variable xi,xd,yd,rad,ratio,domega,xy
994        Variable lambda,trans,trans_err,lat_err,tmp_err,lat_corr
995       
996//      Print "...doing jacobian and non-linear corrections"
997
998        NVAR pixelsX = root:myGlobals:gNPixelsX
999        NVAR pixelsY = root:myGlobals:gNPixelsY
1000       
1001        //set up values to send to auxiliary trig functions
1002        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
1003        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
1004
1005        x0 = realsread[16]
1006        y0 = realsread[17]
1007        sx = realsread[10]
1008        sx3 = realsread[11]
1009        sy = realsread[13]
1010        sy3 = realsread[14]
1011       
1012        dtdist = 1000*realsread[18]     //sdd in mm
1013        dtdis2 = dtdist^2
1014       
1015        lambda = realsRead[26]
1016        trans = RealsRead[4]
1017        trans_err = RealsRead[41]               //new, March 2011
1018       
1019
1020        //waves to contain repeated function calls
1021        Make/O/N=(pixelsX) fyy,xx,yy            //Assumes square detector !!!
1022        ii=0
1023        do
1024                xi = ii
1025//              fyy[ii] = dc_fy(ii+1,sy,sy3,ycenter)
1026//              xx[ii] = dc_fxn(ii+1,sx,sx3,xcenter)
1027//              yy[ii] = dc_fym(ii+1,sy,sy3,ycenter)
1028                ii+=1
1029        while(ii<pixelsX)
1030       
1031        Make/O/N=(pixelsX,pixelsY) SolidAngle           // testing only
1032       
1033        ii=0
1034        do
1035                xi = ii
1036//              xd = dc_fx(ii+1,sx,sx3,xcenter)-xx0
1037                jj=0
1038                do
1039                        yd = fyy[jj]-yy0
1040                        //rad is the distance of pixel ij from the sample
1041                        //domega is the ratio of the solid angle of pixel ij versus center pixel
1042                        // product xy = 1 for a detector with a linear spatial response (modern Ordela)
1043                        // solid angle calculated, dW^3 >=1, so multiply data to raise measured values to correct values.
1044                        rad = sqrt(dtdis2 + xd^2 + yd^2)
1045                        domega = rad/dtdist
1046                        ratio = domega^3
1047                        xy = xx[ii]*yy[jj]
1048                       
1049                        data[ii][jj] *= xy*ratio
1050                       
1051                        solidAngle[ii][jj] = xy*ratio           //testing only 
1052                        data_err[ii][jj] *= xy*ratio                    //error propagation assumes that SA and Jacobian are exact, so simply scale error
1053                       
1054                       
1055                        // correction factor for detector efficiency JBG memo det_eff_cor2.doc 3/20/07
1056                        // correction inserted 11/2007 SRK
1057                        // large angle detector efficiency is >= 1 and will "bump up" the measured value at the highest angles
1058                        // so divide here to get the correct answer (5/22/08 SRK)
1059                        if(doEfficiency)
1060//                              data[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
1061//                              data_err[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
1062//                              solidAngle[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)           //testing only
1063                        endif
1064                       
1065                        // large angle transmission calculation is <= 1 and will "bump down" the measured value at the highest angles
1066                        // so divide here to get the correct answer
1067                        if(doTrans)
1068                       
1069                                if(trans<0.1 && ii==0 && jj==0)
1070                                        Print "***transmission is less than 0.1*** and is a significant correction"
1071                                endif
1072                               
1073                                if(trans==0)
1074                                        if(ii==0 && jj==0)
1075                                                Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
1076                                        endif
1077                                        trans = 1
1078                                endif
1079                               
1080                                // pass in the transmission error, and the error in the correction is returned as the last parameter
1081
1082//                              lat_corr = V_LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,lat_err)           //moved from 1D avg SRK 11/2007
1083
1084                                data[ii][jj] /= lat_corr                        //divide by the correction factor
1085                                //
1086                                //
1087                                //
1088                                // relative errors add in quadrature
1089                                tmp_err = (data_err[ii][jj]/lat_corr)^2 + (lat_err/lat_corr)^2*data[ii][jj]*data[ii][jj]/lat_corr^2
1090                                tmp_err = sqrt(tmp_err)
1091                               
1092                                data_err[ii][jj] = tmp_err
1093                               
1094//                              solidAngle[ii][jj] = lat_err
1095
1096                               
1097                                //solidAngle[ii][jj] = LargeAngleTransmissionCorr(trans,dtdist,xd,yd)           //testing only
1098                        endif
1099                       
1100                        jj+=1
1101                while(jj<pixelsX)
1102                ii+=1
1103        while(ii<pixelsX)
1104       
1105        //clean up waves
1106       
1107        Return(0)
1108End
1109
1110
1111//
1112// Large angle transmission correction
1113//
1114// DIVIDE the intensity by this correction to get the right answer
1115//
1116//
1117// Apply the large angle transmssion correction as the data is converted to WORK
1118// so that whether the data is saved as 2D or 1D, the correction has properly been done.
1119//
1120// This is, however, a SAMPLE dependent calculation, not purely instrument geometry.
1121//
1122Function V_LargeAngleTransmissionCorr(w,w_err,fname,detStr,destPath)
1123        Wave w,w_err
1124        String fname,detStr,destPath
1125
1126        Variable sdd,xCtr,yCtr,trans,trans_err,uval
1127
1128// get all of the geometry information 
1129//      orientation = V_getDet_tubeOrientation(fname,detStr)
1130        sdd = V_getDet_ActualDistance(fname,detStr)
1131
1132        // this is ctr in mm
1133        xCtr = V_getDet_beam_center_x_mm(fname,detStr)
1134        yCtr = V_getDet_beam_center_y_mm(fname,detStr)
1135        trans = V_getSampleTransmission(fname)
1136        trans_err = V_getSampleTransError(fname)
1137       
1138        SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr)
1139       
1140        Wave data_realDistX = data_realDistX
1141        Wave data_realDistY = data_realDistY
1142
1143        Duplicate/O w lat_corr,tmp_theta,tmp_dist,lat_err,tmp_err               //in the current df
1144
1145//// calculate the scattering angle
1146//      dx = (distX - xctr)             //delta x in mm
1147//      dy = (distY - yctr)             //delta y in mm
1148        tmp_dist = sqrt((data_realDistX - xctr)^2 + (data_realDistY - yctr)^2)
1149       
1150        tmp_dist /= 10  // convert mm to cm
1151        // sdd is in [cm]
1152
1153        tmp_theta = atan(tmp_dist/sdd)          //this is two_theta, the scattering angle
1154
1155        Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp
1156        numx = DimSize(tmp_theta,0)
1157        numy = DimSize(tmp_theta,1)
1158       
1159       
1160        //optical thickness
1161        uval = -ln(trans)               //use natural logarithm
1162       
1163        for(ii=0        ;ii<numx;ii+=1)
1164                for(jj=0;jj<numy;jj+=1)
1165                       
1166                        cos_th = cos(tmp_theta[ii][jj])
1167                        arg = (1-cos_th)/cos_th
1168                       
1169                        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
1170                        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120
1171                        // OR
1172                        if((uval<0.01) || (cos_th>0.99))       
1173                                //small arg, approx correction
1174                                lat_corr[ii][jj] = 1-0.5*uval*arg
1175                        else
1176                                //large arg, exact correction
1177                                lat_corr[ii][jj] = (1-exp(-uval*arg))/(uval*arg)
1178                        endif
1179                         
1180                        // (DONE)
1181                        // x- properly calculate and apply the 2D error propagation
1182                        if(trans == 1)
1183                                lat_err[ii][jj] = 0             //no correction, no error
1184                        else
1185                                //sigT, calculated from the Taylor expansion
1186                                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3)
1187                                tmp *= tmp
1188                                tmp *= trans_err^2
1189                                tmp = sqrt(tmp)         //sigT
1190                               
1191                                lat_err[ii][jj] = tmp
1192                        endif
1193                         
1194 
1195                endfor
1196        endfor
1197       
1198
1199       
1200        // Here it is! Apply the correction to the intensity (divide -- to get the proper correction)
1201        w /= lat_corr
1202
1203        // relative errors add in quadrature to the current 2D error
1204        tmp_err = (w_err/lat_corr)^2 + (lat_err/lat_corr)^2*w*w/lat_corr^2
1205        tmp_err = sqrt(tmp_err)
1206       
1207        w_err = tmp_err
1208       
1209
1210        // DONE x- clean up after I'm satisfied computations are correct               
1211        KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err
1212       
1213        return(0)
1214end
1215
1216
1217
1218//
1219//test procedure, not called anymore
1220Proc V_AbsoluteScaling(type,c0,c1,c2,c3,c4,c5,I_err)
1221        String type
1222        Variable c0=1,c1=0.1,c2=0.95,c3=0.1,c4=1,c5=32.0,I_err=0.32
1223        Prompt type,"WORK data type",popup,"CAL;COR;SAM"
1224        Prompt c0, "Sample Transmission"
1225        Prompt c1, "Sample Thickness (cm)"
1226        Prompt c2, "Standard Transmission"
1227        Prompt c3, "Standard Thickness (cm)"
1228        Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)"
1229        Prompt c5, "Standard Cross-Section (cm-1)"
1230        Prompt I_err, "error in I(q=0) (one std dev)"
1231
1232        Variable err
1233        //call the function to do the math
1234        //data from "type" will be scaled and deposited in ABS
1235        err = V_Absolute_Scale(type,c0,c1,c2,c3,c4,c5,I_err)
1236       
1237        if(err)
1238                Abort "Error in V_Absolute_Scale()"
1239        endif
1240       
1241        //contents are always dumped to ABS
1242        type = "ABS"
1243       
1244        //need to update the display with "data" from the correct dataFolder
1245        //reset the current display type to "type"
1246        SVAR gCurDispType = root:Packages:NIST:VSANS:Globals:gCurDispType
1247        gCurDispType = Type     
1248       
1249        V_UpdateDisplayInformation(Type)
1250       
1251End
1252
1253//
1254//
1255// kappa comes in as s_izero, so be sure to use 1/kappa_err
1256//
1257//convert the "type" data to absolute scale using the given standard information
1258//s_ is the standard
1259//w_ is the "work" file
1260//both are work files and should already be normalized to 10^8 monitor counts
1261Function V_Absolute_Scale(type,w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err)
1262        String type
1263        Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
1264
1265
1266        Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
1267        Variable scale,trans_err
1268        Variable err,ii
1269        String detStr
1270       
1271        // be sure that the starting data exists
1272        err = V_WorkDataExists(type)
1273        if(err==1)
1274                return(err)
1275        endif
1276               
1277        //copy from current dir (type) to ABS
1278        V_CopyHDFToWorkFolder(type,"ABS")       
1279
1280// TODO: -- which monitor to use? Here, I think it should already be normalized to 10^8
1281//     
1282//      w_moncount = V_getMonitorCount(type)            //monitor count in "type"
1283       
1284        w_moncount = V_getBeamMonNormData(type)
1285       
1286       
1287        if(w_moncount == 0)
1288                //zero monitor counts will give divide by zero ---
1289                DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
1290                Return(1)               //report error
1291        Endif
1292       
1293        //calculate scale factor
1294        s1 = defmon/w_moncount          // monitor count (s1 should be 1)
1295        s2 = s_thick/w_thick
1296        s3 = s_trans/w_trans
1297        s4 = s_cross/s_izero
1298        scale = s1*s2*s3*s4
1299
1300        trans_err = V_getSampleTransError(type)
1301       
1302        // kappa comes in as s_izero, so be sure to use 1/kappa_err
1303
1304        // and now loop through all of the detectors
1305        //do the actual absolute scaling here, modifying the data in ABS
1306        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
1307                detStr = StringFromList(ii, ksDetectorListAll, ";")
1308                Wave data = V_getDetectorDataW("ABS",detStr)
1309                Wave data_err = V_getDetectorDataErrW("ABS",detStr)
1310               
1311                data *= scale
1312                data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1313        endfor
1314       
1315        //********* 15APR02
1316        // DO NOT correct for atenuators here - the COR step already does this, putting all of the data on equal
1317        // footing (zero atten) before doing the subtraction.
1318       
1319        Return (0) //no error
1320End
1321
1322
1323//
1324// TODO:
1325//   --         DoAlert 0,"This has not yet been updated for VSANS"
1326//
1327//
1328// match the attenuation of the RAW data to the "type" data
1329// so that they can be properly added
1330//
1331// are the attenuator numbers the same? if so exit
1332//
1333// if not, find the attenuator number for type
1334// - find both attenuation factors
1335//
1336// rescale the raw data to match the ratio of the two attenuation factors
1337// -- adjust the detector count (rw)
1338// -- the linear data
1339//
1340//
1341Function V_Adjust_RAW_Attenuation(type)
1342        String type
1343
1344        DoAlert 0,"This has not yet been updated for VSANS"
1345       
1346        WAVE rw=$("root:Packages:NIST:RAW:realsread")
1347        WAVE linear_data=$("root:Packages:NIST:RAW:linear_data")
1348        WAVE data=$("root:Packages:NIST:RAW:data")
1349        WAVE data_err=$("root:Packages:NIST:RAW:linear_data_error")
1350        WAVE/T tw = $("root:Packages:NIST:RAW:textRead")
1351       
1352        WAVE dest_reals=$("root:Packages:NIST:"+type+":realsread")
1353
1354        Variable dest_atten,raw_atten,tol
1355        Variable lambda,raw_atten_err,raw_AttenFactor,dest_attenFactor,dest_atten_err
1356        String fileStr
1357
1358        dest_atten = dest_reals[3]
1359        raw_atten = rw[3]
1360       
1361        tol = 0.1               // within 0.1 atten units is OK
1362        if(abs(dest_atten - raw_atten) < tol )
1363                return(0)
1364        endif
1365
1366        fileStr = tw[3]
1367        lambda = rw[26]
1368        // TODO access correct values
1369        raw_AttenFactor = 1//AttenuationFactor(fileStr,lambda,raw_atten,raw_atten_err)
1370        dest_AttenFactor = 1//AttenuationFactor(fileStr,lambda,dest_atten,dest_atten_err)
1371               
1372        rw[2] *= dest_AttenFactor/raw_AttenFactor
1373        linear_data *= dest_AttenFactor/raw_AttenFactor
1374       
1375        // to keep "data" and linear_data in sync
1376        data = linear_data
1377       
1378        return(0)
1379End
1380
1381//
1382// testing procedure, called from a menu selection
1383//
1384Proc V_DIV_a_Workfile(type)
1385        String type
1386        Prompt type,"WORK data type",popup,"SAM;EMP;BGD;ADJ;"
1387       
1388        //macro will take whatever is in SELECTED folder and DIVide it by the current
1389        //contents of the DIV folder - the function will check for existence
1390        //before proceeding
1391
1392        Abort "This has not yet been updated for VSANS"
1393       
1394        Variable err
1395        err = V_DIVCorrection(type)             //returns err = 1 if data doesn't exist in specified folders
1396       
1397        if(err)
1398                Abort "error in V_DIVCorrection()"
1399        endif
1400       
1401        //contents are NOT always dumped to CAL, but are in the new type folder
1402       
1403        String newTitle = "WORK_"+type
1404        DoWindow/F VSANS_Data
1405        DoWindow/T VSANS_Data, newTitle
1406        KillStrings/Z newTitle
1407       
1408        //need to update the display with "data" from the correct dataFolder
1409        //reset the current displaytype to "type"
1410        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=Type
1411       
1412        V_UpdateDisplayInformation(type)
1413       
1414End
1415
1416
1417//
1418// TODO:
1419//   x-         DoAlert 0,"This has not yet been updated for VSANS"
1420//   -- how is the error propagation handled? Be sure it is calculated correctly when DIV is generated
1421//      and is applied correctly here...
1422//
1423//function will divide the contents of "workType" folder with the contents of
1424//the DIV folder + detStr
1425// all data is linear scale for the calculation
1426//
1427Function V_DIVCorrection(data,data_err,detStr,workType)
1428        Wave data,data_err
1429        String detStr,workType
1430       
1431        //check for existence of data in type and DIV
1432        // if the desired data doesn't exist, let the user know, and abort
1433        String destPath=""
1434
1435        if(WaveExists(data) == 0)
1436                Print "The data wave does not exist in V_DIVCorrection()"
1437                Return(1)               //error condition
1438        Endif
1439       
1440        //check for DIV
1441        // if the DIV workfile doesn't exist, let the user know,and abort
1442        // !! be sure to check first, before trying to access the wave
1443       
1444//      WAVE/Z div_data = $("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")
1445        if(WaveExists($("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data")) == 0)
1446                Print "The DIV wave does not exist in V_DIVCorrection()"
1447                Return(1)               //error condition
1448        Endif
1449        if(WaveExists($("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":linear_data_error")) == 0)
1450                Print "The DIV error wave does not exist in V_DIVCorrection()"
1451                Return(1)               //error condition
1452        Endif
1453        //files exist, proceed
1454
1455        WAVE/Z div_data_err = V_getDetectorDataErrW("DIV",detStr)
1456        WAVE/Z div_data = V_getDetectorDataW("DIV",detStr)
1457
1458
1459
1460// do the error propagation first, since data is changed by the correction
1461        data_err = sqrt(data_err^2/div_data^2 + div_data_err^2 * data^2/div_data^4 )
1462
1463// then the correction
1464        data /= div_data
1465
1466       
1467        Return(0)
1468End
1469
1470
1471//////////////////////////
Note: See TracBrowser for help on using the repository browser.