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

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

additions to VCALC procedures to correctly account for panel motion (individual, not symmetric). Updated the plotting routines to all (mostly) pass through the same subroutines so that additional averaging modes will be easier to add.

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