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

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

minor changes to allow patching of more incorrect metadata

added simple way to estimate beam centers from a single measurement (stiil to be thoroughly verified!)

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