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

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

extensive changes to accomodate 1x1 binning of the HighRes? detector. It is implemented as a global flag. Currently only 4x4 and 1x1 are allowed. 1x1 has never been tested in reality, only simulated data - so my assumed dimensions may not be correct. look for TODOHIGHRES in the file for places that may need to be updated for different file dimensions. Testing of the simulated data is proving to be excruciatingly slow, but passable for a test. Speed optimization will be needed if this is the final solution. Memory management will also be an issue since every "copy" of the highRes matrix is enormous. Carry as few of these around as possible in the future to keep the experiment size to something reasonable.

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