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

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

Important change -- re-worked the calculation of how the detector panel "gap" is treated. It is now split equally between the panels rather than assigned to one panel. Makes no difference in the q calculation, simply a philosophical change that should make the measured beam centers make a little more sense, as they will be symmetric.

Changed the appearance and function of the panel to select the "trim" values for the combined 1D data sets. Data no longer needs to be saved as ITX and re-read into the panel. Simply uses the data binning popup and the current data.

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