source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_Correct.ipf @ 1019

Last change on this file since 1019 was 1019, checked in by srkline, 6 years ago

adding procedures for:

simple save of a DIV file. no functionality to generate a DIV file yet, since I don't know how this will happen.

Simple dump of the file structure in a data "tree"

Verified that the error bars on the I(q) data are correctly calculated as standard error of the mean. There was never an issue with this, or with SANS calculations.

Started filling in "Correct" routines based on the SANS version. Only stubs present currently.

File size: 40.5 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=6.1
4
5
6//
7// Procedures to perform the "Correct" step during data reduction
8//
9// - there is only one procedure to perform the subtractions, and a single
10// parameter flags which subtractions are to be done. Different numbers of
11// attenuators during scattering runs are corrected as described in John's memo,
12// with the note that ONLY method (3) is used, which assumes that 'diffuse' scattering
13// is dominant over 'dark current' (note that 'dark current' = shutter CLOSED)
14//
15//
16//do the CORRECT step based on the answers to emp and bkg subtraction
17//by setting the proper"mode"
18        //1 = both emp and bgd subtraction
19        //2 = only bgd subtraction
20        //3 = only emp subtraction
21        //4 = no subtraction
22        //additional modes 091301
23        //11 = emp, bgd, drk
24        //12 = bgd and drk
25        //13 = emp and drk
26        //14 = no subtractions
27        //
28//********************************
29
30//
31//unused test procedure for Correct() function
32//must be updated to include "mode" parameter before re-use
33//
34Macro V_CorrectData(cor_mode)
35        Variable cor_mode=4
36       
37        Variable err
38        String cor_type
39       
40        err = V_Correct(cor_mode)               
41       
42        if(err)
43                Abort "error in V_Correct"
44        endif
45       
46        //contents are always dumped to COR
47        cor_type = "COR"
48       
49        // set the global string       
50        String/G root:Packages:NIST:VSANS:Globals:gCurDispType=cor_type
51       
52        //need to update the display with "data" from the correct dataFolder
53        UpdateDisplayInformation(cor_type)
54       
55End
56
57//
58// dispatch function to correct data
59//
60// --determines the mode of subtraction
61// --TODO (checks for trans=1 in SAM and EMP) - not currently implemented
62// --checks for existence of data in WORK folders
63// --clears out the COR folder
64// --calls the correct correct
65//
66//
67//mode describes the type of subtraction that is to be done
68//1 = both emp and bgd subtraction
69//2 = only bgd subtraction
70//3 = only emp subtraction
71//4 = no subtraction
72//
73// + 10 indicates that WORK.DRK is to be used
74//
75//091301 version
76//now simply dispatches to the correct subtraction - logic was too
77//involved to do in one function - unclear and error-prone
78//
79// 081203 version
80// checks for trans==1 in SAM and EMP before dispatching
81// and asks for new value if desired
82//
83//
84//
85Function V_Correct(cor_mode)
86        Variable cor_mode
87       
88        Variable err=0,trans,newTrans
89       
90        //switch and dispatch based on the required subtractions
91        // always check for SAM data
92        err = V_WorkDataExists("SAM")
93        if(err==1)
94                return(err)
95        endif
96       
97        // TODO -- this is currently not implemented   
98        //check for trans==1
99//      NVAR doCheck=root:Packages:NIST:gDoTransCheck
100//      Wave/Z samR=root:Packages:NIST:SAM:RealsRead
101//      Wave/Z empR=root:Packages:NIST:EMP:RealsRead
102//      if(doCheck)
103//              trans = samR[4]
104//              newTrans=GetNewTrans(trans,"SAM")               //will change value if necessary
105//              if(numtype(newTrans)==0)
106//                      samR[4] = newTrans              //avoid user abort assigning NaN
107//              endif
108//              if(trans != newTrans)
109//                      print "Using SAM trans = ",samR[4]
110//              endif
111//      endif
112       
113        //copy SAM information to COR, wiping out the old contents of the COR folder first
114        //do this even if no correction is dispatched (if incorrect mode)
115
116        CopyHDFToWorkFolder("SAM","COR")       
117
118       
119//      Print "dispatching to mode = ",mode
120        switch(cor_mode)
121                case 1:
122                        err = V_WorkDataExists("EMP")
123                        if(err==1)
124                                return(err)
125                        Endif
126                        // TODO -- not implemented
127//                      if(doCheck)
128//                              trans = empR[4]
129//                              newTrans=GetNewTrans(trans,"EMP")               //will change value if necessary
130//                              if(numtype(newTrans)==0)
131//                                      empR[4] = newTrans
132//                              endif
133//                              if(trans != newTrans)
134//                                      print "Using EMP trans = ",empR[4]
135//                              endif
136//                      endif
137                        err = V_WorkDataExists("BGD")
138                        if(err==1)
139                                return(err)
140                        Endif
141// TODO                 err = V_CorrectMode_1()
142                        break
143                case 2:
144                        err = V_WorkDataExists("BGD")
145                        if(err==1)
146                                return(err)
147                        Endif
148// TODO                 err = V_CorrectMode_2()
149                        break
150                case 3:
151                        err = V_WorkDataExists("EMP")
152                        if(err==1)
153                                return(err)
154                        Endif
155                        // TODO - not implemented
156//                      if(doCheck)
157//                              trans = empR[4]
158//                              newTrans=GetNewTrans(trans,"EMP")               //will change value if necessary
159//                              if(numtype(newTrans)==0)
160//                                      empR[4] = newTrans
161//                              endif
162//                              if(trans != newTrans)
163//                                      print "Using EMP trans = ",empR[4]
164//                              endif
165//                      endif
166
167// TODO                 err = V_CorrectMode_3()
168                        break
169                case 4:
170                        err = V_CorrectMode_4()
171                        break
172                case 11:
173                        err = V_WorkDataExists("EMP")
174                        if(err==1)
175                                return(err)
176                        Endif
177                        // TODO -- not implemented
178//                      if(doCheck)
179//                              trans = empR[4]
180//                              newTrans=GetNewTrans(trans,"EMP")               //will change value if necessary
181//                              if(numtype(newTrans)==0)
182//                                      empR[4] = newTrans
183//                              endif
184//                              if(trans != newTrans)
185//                                      print "Using EMP trans = ",empR[4]
186//                              endif
187//                      endif
188                        err = V_WorkDataExists("BGD")
189                        if(err==1)
190                                return(err)
191                        Endif
192                        err = V_WorkDataExists("DRK")
193                        if(err==1)
194                                return(err)
195                        Endif
196// TODO                 err = V_CorrectMode_11()
197                        break
198                case 12:
199                        err = V_WorkDataExists("BGD")
200                        if(err==1)
201                                return(err)
202                        Endif
203                        err = V_WorkDataExists("DRK")
204                        if(err==1)
205                                return(err)
206                        Endif
207// TODO                 err = V_CorrectMode_12()
208                        break
209                case 13:
210                        err = V_WorkDataExists("EMP")
211                        if(err==1)
212                                return(err)
213                        Endif
214                        // TODO -- not implemented
215//                      if(doCheck)
216//                              trans = empR[4]
217//                              newTrans=GetNewTrans(trans,"EMP")               //will change value if necessary
218//                              if(numtype(newTrans)==0)
219//                                      empR[4] = newTrans
220//                              endif
221//                              if(trans != newTrans)
222//                                      print "Using EMP trans = ",empR[4]
223//                              endif
224//                      endif
225                        err = V_WorkDataExists("DRK")
226                        if(err==1)
227                                return(err)
228                        Endif
229// TODO                 err = V_CorrectMode_13()
230                        break
231                case 14:
232                        err = V_WorkDataExists("DRK")
233                        if(err==1)
234                                return(err)
235                        Endif
236// TODO                 err = V_CorrectMode_14()
237                        break
238                default:        //something wrong
239                        Print "Incorrect mode in V_Correct()"
240                        return(1)       //error
241        endswitch
242
243        //calculation attempted, return the result
244        return(err)     
245End
246
247// subtraction of both EMP and BGD from SAM
248// data exists, checked by dispatch routine
249//
250// this is the most common use
251// March 2011 added error propagation
252//                                      added explicit reference to use linear_data, instead of trusting that data
253//                                      was freshly loaded. added final copy of cor result to cor:data and cor:linear_data
254//
255xFunction V_CorrectMode_1()
256       
257        //create the necessary wave references
258        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data"
259        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread"
260        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread"
261        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread"
262        WAVE bgd_data=$"root:Packages:NIST:BGD:linear_data"
263        WAVE bgd_reals=$"root:Packages:NIST:BGD:realsread"
264        WAVE bgd_ints=$"root:Packages:NIST:BGD:integersread"
265        WAVE/T bgd_text=$"root:Packages:NIST:BGD:textread"
266        WAVE emp_data=$"root:Packages:NIST:EMP:linear_data"
267        WAVE emp_reals=$"root:Packages:NIST:EMP:realsread"
268        WAVE emp_ints=$"root:Packages:NIST:EMP:integersread"
269        WAVE/T emp_text=$"root:Packages:NIST:EMP:textread"
270        WAVE cor_data=$"root:Packages:NIST:COR:linear_data"
271        WAVE/T cor_text=$"root:Packages:NIST:COR:textread"
272       
273        // needed to propagate error
274        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy
275        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error"
276        WAVE bgd_err =$"root:Packages:NIST:BGD:linear_data_error"
277        WAVE emp_err =$"root:Packages:NIST:EMP:linear_data_error"
278        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error"
279       
280        Variable sam_trans_err,emp_trans_err
281        sam_trans_err = sam_reals[41]
282        emp_trans_err = emp_reals[41]
283       
284       
285        //get sam and bgd attenuation factors
286        String fileStr=""
287        Variable lambda,attenNo,sam_AttenFactor,bgd_attenFactor,emp_AttenFactor
288        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
289        Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp
290        Variable sam_atten_err,emp_atten_err,bgd_atten_err
291        fileStr = sam_text[3]
292        lambda = sam_reals[26]
293        attenNo = sam_reals[3]
294        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err)
295        fileStr = bgd_text[3]
296        lambda = bgd_reals[26]
297        attenNo = bgd_reals[3]
298        bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,bgd_atten_err)
299        fileStr = emp_text[3]
300        lambda = emp_reals[26]
301        attenNo = emp_reals[3]
302        emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,emp_atten_err)
303       
304        //get relative monitor counts (should all be 10^8, since normalized in add step)
305        tmonsam = sam_reals[0]          //monitor count in SAM
306        tsam = sam_reals[4]             //SAM transmission
307        csam = sam_reals[16]            //x center
308        rsam = sam_reals[17]            //beam (x,y) define center of corrected field
309        tmonbgd = bgd_reals[0]          //monitor count in BGD
310        cbgd = bgd_reals[16]
311        rbgd = bgd_reals[17]
312        tmonemp = emp_reals[0]          //monitor count in EMP
313        temp = emp_reals[4]                     //trans emp
314        cemp = emp_reals[16]            //beamcenter of EMP
315        remp = emp_reals[17]
316       
317        if(temp==0)
318                DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction"
319                temp=1
320        Endif
321       
322        NVAR pixelsX = root:myGlobals:gNPixelsX
323        NVAR pixelsY = root:myGlobals:gNPixelsY
324       
325        //get the shifted data arrays, EMP and BGD, each relative to SAM
326        Make/D/O/N=(pixelsX,pixelsY) cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp
327        xshift = cbgd-csam
328        yshift = rbgd-rsam
329        if(abs(xshift) <= wcen)
330                xshift = 0
331        Endif
332        if(abs(yshift) <= wcen)
333                yshift = 0
334        Endif
335        GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift)              //bgd_temp
336       
337        xshift = cemp-csam
338        yshift = remp-rsam
339        if(abs(xshift) <= wcen)
340                xshift = 0
341        Endif
342        if(abs(yshift) <= wcen)
343                yshift = 0
344        Endif
345        GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift)              //emp_temp
346
347        //do the subtraction
348        fsam=1
349        femp = tmonsam/tmonemp          //this should be ==1 since normalized files
350        fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
351        cor1 = fsam*sam_data/sam_attenFactor - fbgd*bgd_temp/bgd_attenFactor
352        cor1 -= (tsam/temp)*(femp*emp_temp/emp_attenFactor - fbgd*bgd_temp/bgd_attenFactor)
353        cor1 *= noadd_bgd*noadd_emp             //zero out the array mismatch values
354
355// do the error propagation piecewise   
356        Duplicate/O sam_err, tmp_a, tmp_b, tmp_c, tmp_d,c_val,d_val
357        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
358       
359        tmp_b = (bgd_err/bgd_attenFactor)^2*(tsam/temp - 1)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2*(1-tsam/temp)^2            //sig b ^2
360
361        tmp_c = (sam_trans_err/temp)^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2
362        tmp_c += (tsam/temp^2)^2*emp_trans_err^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2
363       
364        tmp_d = (tsam/(temp*emp_attenFactor))^2*(emp_err)^2 + (tsam*emp_data/(temp*emp_attenFactor^2))^2*(emp_atten_err)^2
365
366        cor_err = sqrt(tmp_a + tmp_b + tmp_c + tmp_d)
367       
368        //we're done, get out w/no error
369        //set the COR data and linear_data to the result
370        cor_data = cor1
371        cor_data_display = cor1
372       
373        //update COR header
374        cor_text[1] = date() + " " + time()             //date + time stamp
375
376        KillWaves/Z cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp
377        Killwaves/Z tmp_a,tmp_b,tmp_c,tmp_d,c_val,d_val
378        SetDataFolder root:
379        Return(0)
380End
381
382//background only
383// existence of data checked by dispatching routine
384// data has already been copied to COR folder
385xFunction V_CorrectMode_2()
386
387        //create the necessary wave references
388        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data"
389        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread"
390        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread"
391        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread"
392        WAVE bgd_data=$"root:Packages:NIST:BGD:linear_data"
393        WAVE bgd_reals=$"root:Packages:NIST:BGD:realsread"
394        WAVE bgd_ints=$"root:Packages:NIST:BGD:integersread"
395        WAVE/T bgd_text=$"root:Packages:NIST:BGD:textread"
396        WAVE cor_data=$"root:Packages:NIST:COR:linear_data"
397        WAVE/T cor_text=$"root:Packages:NIST:COR:textread"
398
399        // needed to propagate error
400        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy
401        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error"
402        WAVE bgd_err =$"root:Packages:NIST:BGD:linear_data_error"
403        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error"
404       
405        Variable sam_trans_err
406        sam_trans_err = sam_reals[41]
407
408       
409        //get sam and bgd attenuation factors
410        String fileStr=""
411        Variable lambda,attenNo,sam_AttenFactor,bgd_attenFactor
412        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
413        Variable wcen=0.001
414        Variable sam_atten_err,bgd_atten_err
415        fileStr = sam_text[3]
416        lambda = sam_reals[26]
417        attenNo = sam_reals[3]
418        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err)
419        fileStr = bgd_text[3]
420        lambda = bgd_reals[26]
421        attenNo = bgd_reals[3]
422        bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,bgd_atten_err)
423       
424        //Print "atten = ",sam_attenFactor,bgd_attenFactor
425       
426        //get relative monitor counts (should all be 10^8, since normalized in add step)
427        tmonsam = sam_reals[0]          //monitor count in SAM
428        csam = sam_reals[16]            //x center
429        rsam = sam_reals[17]            //beam (x,y) define center of corrected field
430        tmonbgd = bgd_reals[0]          //monitor count in BGD
431        cbgd = bgd_reals[16]
432        rbgd = bgd_reals[17]
433
434        // set up beamcenter shift, relative to SAM
435        xshift = cbgd-csam
436        yshift = rbgd-rsam
437        if(abs(xshift) <= wcen)
438                xshift = 0
439        Endif
440        if(abs(yshift) <= wcen)
441                yshift = 0
442        Endif
443       
444        NVAR pixelsX = root:myGlobals:gNPixelsX
445        NVAR pixelsY = root:myGlobals:gNPixelsY
446        //get shifted data arrays, relative to SAM
447        Make/D/O/N=(pixelsX,pixelsY) cor1,bgd_temp,noadd_bgd            //temp arrays
448        GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift)              //bgd_temp is the BGD
449       
450        //do the sam-bgd subtraction,  deposit result in cor1
451        fsam = 1
452        fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
453       
454        //print "fsam,fbgd = ",fsam,fbgd
455       
456        cor1 = fsam*sam_data/sam_AttenFactor - fbgd*bgd_temp/bgd_AttenFactor
457        cor1 *= noadd_bgd               //zeros out regions where arrays do not overlap, one otherwise
458
459// do the error propagation piecewise   
460        Duplicate/O sam_err, tmp_a, tmp_b
461        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
462       
463        tmp_b = (bgd_err/bgd_attenFactor)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2              //sig b ^2
464
465        cor_err = sqrt(tmp_a + tmp_b)
466
467
468        //we're done, get out w/no error
469        //set the COR_data to the result
470        cor_data = cor1
471        cor_data_display = cor1
472
473        //update COR header
474        cor_text[1] = date() + " " + time()             //date + time stamp
475
476        KillWaves/Z cor1,bgd_temp,noadd_bgd
477        Killwaves/Z tmp_a,tmp_b
478
479        SetDataFolder root:
480        Return(0)
481End
482
483// empty subtraction only
484// data does exist, checked by dispatch routine
485//
486xFunction V_CorrectMode_3()
487        //create the necessary wave references
488        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data"
489        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread"
490        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread"
491        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread"
492        WAVE emp_data=$"root:Packages:NIST:EMP:linear_data"
493        WAVE emp_reals=$"root:Packages:NIST:EMP:realsread"
494        WAVE emp_ints=$"root:Packages:NIST:EMP:integersread"
495        WAVE/T emp_text=$"root:Packages:NIST:EMP:textread"
496        WAVE cor_data=$"root:Packages:NIST:COR:linear_data"
497        WAVE/T cor_text=$"root:Packages:NIST:COR:textread"
498       
499        // needed to propagate error
500        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy
501        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error"
502        WAVE emp_err =$"root:Packages:NIST:EMP:linear_data_error"
503        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error"
504       
505        Variable sam_trans_err,emp_trans_err
506        sam_trans_err = sam_reals[41]
507        emp_trans_err = emp_reals[41]   
508       
509        //get sam and bgd attenuation factors
510        String fileStr=""
511        Variable lambda,attenNo,sam_AttenFactor,emp_attenFactor
512        Variable tmonsam,fsam,femp,xshift,yshift,rsam,csam,remp,cemp,tmonemp
513        Variable wcen=0.001,tsam,temp
514        Variable sam_atten_err,emp_atten_err
515        fileStr = sam_text[3]
516        lambda = sam_reals[26]
517        attenNo = sam_reals[3]
518        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err)
519        fileStr = emp_text[3]
520        lambda = emp_reals[26]
521        attenNo = emp_reals[3]
522        emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,emp_atten_err)
523       
524        //get relative monitor counts (should all be 10^8, since normalized in add step)
525        tmonsam = sam_reals[0]          //monitor count in SAM
526        tsam = sam_reals[4]             //SAM transmission
527        csam = sam_reals[16]            //x center
528        rsam = sam_reals[17]            //beam (x,y) define center of corrected field
529        tmonemp = emp_reals[0]          //monitor count in EMP
530        temp = emp_reals[4]                     //trans emp
531        cemp = emp_reals[16]            //beamcenter of EMP
532        remp = emp_reals[17]
533       
534        if(temp==0)
535                DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction"
536                temp=1
537        Endif
538       
539        //Print "rbgd,cbgd = ",rbgd,cbgd
540        // set up beamcenter shift, relative to SAM
541        xshift = cemp-csam
542        yshift = remp-rsam
543        if(abs(xshift) <= wcen)
544                xshift = 0
545        Endif
546        if(abs(yshift) <= wcen)
547                yshift = 0
548        Endif
549       
550        NVAR pixelsX = root:myGlobals:gNPixelsX
551        NVAR pixelsY = root:myGlobals:gNPixelsY
552        //get shifted data arrays, relative to SAM
553        Make/D/O/N=(pixelsX,pixelsY) cor1,emp_temp,noadd_emp            //temp arrays
554        GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift)              //emp_temp is the EMP
555       
556        //do the sam-bgd subtraction,  deposit result in cor1
557        fsam = 1
558        femp = tmonsam/tmonemp          //this should be ==1 since normalized files
559       
560        cor1 = fsam*sam_data/sam_AttenFactor - femp*(tsam/temp)*emp_temp/emp_AttenFactor
561        cor1 *= noadd_emp               //zeros out regions where arrays do not overlap, one otherwise
562
563// do the error propagation piecewise   
564        Duplicate/O sam_err, tmp_a, tmp_c ,c_val
565        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
566       
567        tmp_c = (sam_trans_err*emp_data/(temp*emp_attenFactor))^2 + (emp_err*tsam/(temp*emp_attenFactor))^2
568        tmp_c += (tsam*emp_data*emp_trans_err/(temp*temp*emp_attenFactor))^2 + (tsam*emp_data*emp_atten_err/(temp*emp_attenFactor^2))^2//total of 6 terms
569
570        cor_err = sqrt(tmp_a + tmp_c)
571       
572        //we're done, get out w/no error
573        //set the COR data to the result
574        cor_data = cor1
575        cor_data_display = cor1
576
577        //update COR header
578        cor_text[1] = date() + " " + time()             //date + time stamp
579
580        KillWaves/Z cor1,emp_temp,noadd_emp
581        Killwaves/Z tmp_a,tmp_c,c_val
582
583        SetDataFolder root:
584        Return(0)
585End
586
587// NO subtraction - simply rescales for attenuators
588// SAM data does exist, checked by dispatch routine
589// SAM data has already been copied to COR (both are the same at the start of the function)
590//
591//
592//
593Function V_CorrectMode_4()
594
595        //get SAM attenuation factor
596        Variable sam_AttenFactor,sam_atten_err,ii
597        String detStr
598       
599        sam_AttenFactor = V_getAttenuator_transmission("SAM")
600        sam_atten_err = V_getAttenuator_trans_err("SAM")
601
602
603        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
604                detStr = StringFromList(ii, ksDetectorListAll, ";")
605                Wave cor_data = V_getDetectorDataW("COR",detStr)
606                Wave cor_err = V_getDetectorDataErrW("COR",detStr)
607                Wave sam_data = V_getDetectorDataW("SAM",detStr)
608                Wave sam_err = V_getDetectorDataErrW("SAM",detStr)
609        endfor
610       
611        cor_data = sam_data/sam_AttenFactor             //simply rescale the data
612
613// do the error propagation piecewise
614        cor_err = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2            //sig a ^2
615        cor_err = sqrt(cor_err)
616
617        //TODO -- do I want to update COR header?
618//      cor_text[1] = date() + " " + time()             //date + time stamp
619
620        SetDataFolder root:
621        Return(0)
622End
623
624xFunction V_CorrectMode_11()
625        //create the necessary wave references
626        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data"
627        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread"
628        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread"
629        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread"
630        WAVE bgd_data=$"root:Packages:NIST:BGD:linear_data"
631        WAVE bgd_reals=$"root:Packages:NIST:BGD:realsread"
632        WAVE bgd_ints=$"root:Packages:NIST:BGD:integersread"
633        WAVE/T bgd_text=$"root:Packages:NIST:BGD:textread"
634        WAVE emp_data=$"root:Packages:NIST:EMP:linear_data"
635        WAVE emp_reals=$"root:Packages:NIST:EMP:realsread"
636        WAVE emp_ints=$"root:Packages:NIST:EMP:integersread"
637        WAVE/T emp_text=$"root:Packages:NIST:EMP:textread"
638        WAVE drk_data=$"root:Packages:NIST:DRK:linear_data"
639        WAVE drk_reals=$"root:Packages:NIST:DRK:realsread"
640        WAVE drk_ints=$"root:Packages:NIST:DRK:integersread"
641        WAVE/T drk_text=$"root:Packages:NIST:DRK:textread"
642        WAVE cor_data=$"root:Packages:NIST:COR:linear_data"
643        WAVE/T cor_text=$"root:Packages:NIST:COR:textread"
644
645        // needed to propagate error
646        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy
647        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error"
648        WAVE bgd_err =$"root:Packages:NIST:BGD:linear_data_error"
649        WAVE emp_err =$"root:Packages:NIST:EMP:linear_data_error"
650        WAVE drk_err =$"root:Packages:NIST:DRK:linear_data_error"
651        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error"
652       
653        Variable sam_trans_err,emp_trans_err
654        sam_trans_err = sam_reals[41]
655        emp_trans_err = emp_reals[41]
656       
657        //get sam and bgd attenuation factors
658        String fileStr=""
659        Variable lambda,attenNo,sam_AttenFactor,bgd_attenFactor,emp_AttenFactor
660        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
661        Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp,time_sam,time_drk,savmon_sam
662        Variable sam_atten_err,bgd_atten_err,emp_atten_err
663        fileStr = sam_text[3]
664        lambda = sam_reals[26]
665        attenNo = sam_reals[3]
666        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err)
667        fileStr = bgd_text[3]
668        lambda = bgd_reals[26]
669        attenNo = bgd_reals[3]
670        bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,bgd_atten_err)
671        fileStr = emp_text[3]
672        lambda = emp_reals[26]
673        attenNo = emp_reals[3]
674        emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,emp_atten_err)
675       
676        //get relative monitor counts (should all be 10^8, since normalized in add step)
677        tmonsam = sam_reals[0]          //monitor count in SAM
678        tsam = sam_reals[4]             //SAM transmission
679        csam = sam_reals[16]            //x center
680        rsam = sam_reals[17]            //beam (x,y) define center of corrected field
681        tmonbgd = bgd_reals[0]          //monitor count in BGD
682        cbgd = bgd_reals[16]
683        rbgd = bgd_reals[17]
684        tmonemp = emp_reals[0]          //monitor count in EMP
685        temp = emp_reals[4]                     //trans emp
686        cemp = emp_reals[16]            //beamcenter of EMP
687        remp = emp_reals[17]
688        savmon_sam=sam_reals[1]         //true monitor count in SAM
689        time_sam = sam_ints[2]          //count time SAM
690        time_drk = drk_ints[2]          //drk count time
691       
692        NVAR pixelsX = root:myGlobals:gNPixelsX
693        NVAR pixelsY = root:myGlobals:gNPixelsY
694        //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM
695        Make/D/O/N=(pixelsX,pixelsY) drk_temp, drk_tmp_err
696        drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam)
697        drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam)                 //temporarily rescale the error of DRK
698       
699        if(temp==0)
700                DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction"
701                temp=1
702        Endif
703       
704        //get the shifted data arrays, EMP and BGD, each relative to SAM
705        Make/D/O/N=(pixelsX,pixelsY) cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp
706        xshift = cbgd-csam
707        yshift = rbgd-rsam
708        if(abs(xshift) <= wcen)
709                xshift = 0
710        Endif
711        if(abs(yshift) <= wcen)
712                yshift = 0
713        Endif
714        GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift)              //bgd_temp
715       
716        xshift = cemp-csam
717        yshift = remp-rsam
718        if(abs(xshift) <= wcen)
719                xshift = 0
720        Endif
721        if(abs(yshift) <= wcen)
722                yshift = 0
723        Endif
724        GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift)              //emp_temp
725        //always ignore the DRK center shift
726       
727        //do the subtraction
728        fsam=1
729        femp = tmonsam/tmonemp          //this should be ==1 since normalized files
730        fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
731        cor1 = fsam*sam_data/sam_attenFactor
732        cor1 -= (tsam/temp)*(femp*emp_temp/emp_attenFactor - fbgd*bgd_temp/bgd_attenFactor)
733        cor1 -= (fbgd*bgd_temp/bgd_attenFactor - drk_temp)
734        cor1 -= drk_temp/sam_attenFactor
735        cor1 *= noadd_bgd*noadd_emp             //zero out the array mismatch values
736       
737// do the error propagation piecewise   
738        Duplicate/O sam_err, tmp_a, tmp_b, tmp_c, tmp_d,c_val,d_val
739        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
740       
741        tmp_b = (bgd_err/bgd_attenFactor)^2*(tsam/temp - 1)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2*(1-tsam/temp)^2            //sig b ^2
742
743        tmp_c = (sam_trans_err/temp)^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2
744        tmp_c += (tsam/temp^2)^2*emp_trans_err^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2
745       
746        tmp_d = (tsam/(temp*emp_attenFactor))^2*(emp_err)^2 + (tsam*emp_data/(temp*emp_attenFactor^2))^2*(emp_atten_err)^2
747
748        cor_err = sqrt(tmp_a + tmp_b + tmp_c + tmp_d + drk_tmp_err^2)
749       
750        //we're done, get out w/no error
751        //set the COR data to the result
752        cor_data = cor1
753        cor_data_display = cor1
754
755        //update COR header
756        cor_text[1] = date() + " " + time()             //date + time stamp
757
758        KillWaves/Z cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp,drk_temp
759        Killwaves/Z tmp_a,tmp_b,tmp_c,tmp_d,c_val,d_val,drk_tmp_err
760
761        SetDataFolder root:
762        Return(0)
763End
764
765//bgd and drk subtraction
766//
767xFunction V_CorrectMode_12()
768        //create the necessary wave references
769        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data"
770        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread"
771        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread"
772        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread"
773        WAVE bgd_data=$"root:Packages:NIST:BGD:linear_data"
774        WAVE bgd_reals=$"root:Packages:NIST:BGD:realsread"
775        WAVE bgd_ints=$"root:Packages:NIST:BGD:integersread"
776        WAVE/T bgd_text=$"root:Packages:NIST:BGD:textread"
777        WAVE drk_data=$"root:Packages:NIST:DRK:linear_data"
778        WAVE drk_reals=$"root:Packages:NIST:DRK:realsread"
779        WAVE drk_ints=$"root:Packages:NIST:DRK:integersread"
780        WAVE/T drk_text=$"root:Packages:NIST:DRK:textread"
781        WAVE cor_data=$"root:Packages:NIST:COR:linear_data"
782        WAVE/T cor_text=$"root:Packages:NIST:COR:textread"
783
784        // needed to propagate error
785        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy
786        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error"
787        WAVE bgd_err =$"root:Packages:NIST:BGD:linear_data_error"
788        WAVE drk_err =$"root:Packages:NIST:DRK:linear_data_error"
789        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error"
790       
791        Variable sam_trans_err
792        sam_trans_err = sam_reals[41]
793       
794       
795        //get sam and bgd attenuation factors
796        String fileStr=""
797        Variable lambda,attenNo,sam_AttenFactor,bgd_attenFactor
798        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
799        Variable wcen=0.001,time_drk,time_sam,savmon_sam,tsam
800        Variable sam_atten_err,bgd_atten_err
801        fileStr = sam_text[3]
802        lambda = sam_reals[26]
803        attenNo = sam_reals[3]
804        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err)
805        fileStr = bgd_text[3]
806        lambda = bgd_reals[26]
807        attenNo = bgd_reals[3]
808        bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,bgd_atten_err)
809       
810        //get relative monitor counts (should all be 10^8, since normalized in add step)
811        tmonsam = sam_reals[0]          //monitor count in SAM
812        tsam = sam_reals[4]             //SAM transmission
813        csam = sam_reals[16]            //x center
814        rsam = sam_reals[17]            //beam (x,y) define center of corrected field
815        tmonbgd = bgd_reals[0]          //monitor count in BGD
816        cbgd = bgd_reals[16]
817        rbgd = bgd_reals[17]
818        savmon_sam=sam_reals[1]         //true monitor count in SAM
819        time_sam = sam_ints[2]          //count time SAM
820        time_drk = drk_ints[2]          //drk count time
821       
822        NVAR pixelsX = root:myGlobals:gNPixelsX
823        NVAR pixelsY = root:myGlobals:gNPixelsY
824        //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM
825        Make/D/O/N=(pixelsX,pixelsY) drk_temp,drk_tmp_err
826        drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam)
827        drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam)                 //temporarily rescale the error of DRK
828
829        // set up beamcenter shift, relative to SAM
830        xshift = cbgd-csam
831        yshift = rbgd-rsam
832        if(abs(xshift) <= wcen)
833                xshift = 0
834        Endif
835        if(abs(yshift) <= wcen)
836                yshift = 0
837        Endif
838        //get shifted data arrays, relative to SAM
839        Make/D/O/N=(pixelsX,pixelsY) cor1,bgd_temp,noadd_bgd            //temp arrays
840        GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift)              //bgd_temp is the BGD
841        //always ignore the DRK center shift
842       
843        //do the sam-bgd subtraction,  deposit result in cor1
844        fsam = 1
845        fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
846       
847        cor1 = fsam*sam_data/sam_AttenFactor + fbgd*tsam*bgd_temp/bgd_AttenFactor
848        cor1 += -1*(fbgd*bgd_temp/bgd_attenFactor - drk_temp) - drk_temp/sam_attenFactor
849        cor1 *= noadd_bgd               //zeros out regions where arrays do not overlap, one otherwise
850
851// do the error propagation piecewise   
852        Duplicate/O sam_err, tmp_a, tmp_b
853        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
854       
855        tmp_b = (bgd_err/bgd_attenFactor)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2              //sig b ^2
856
857        cor_err = sqrt(tmp_a + tmp_b + drk_tmp_err^2)
858
859        //we're done, get out w/no error
860        //set the COR_data to the result
861        cor_data = cor1
862        cor_data_display = cor1
863
864        //update COR header
865        cor_text[1] = date() + " " + time()             //date + time stamp
866
867        KillWaves/Z cor1,bgd_temp,noadd_bgd,drk_temp
868        Killwaves/Z tmp_a,tmp_b,drk_tmp_err
869
870        SetDataFolder root:
871        Return(0)
872End
873
874//EMP and DRK subtractions
875// all data exists, DRK is on a time basis (noNorm)
876//scale DRK by monitor count scaling factor and the ratio of couting times
877//to place the DRK file on equal footing
878xFunction V_CorrectMode_13()
879        //create the necessary wave references
880        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data"
881        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread"
882        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread"
883        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread"
884        WAVE emp_data=$"root:Packages:NIST:EMP:linear_data"
885        WAVE emp_reals=$"root:Packages:NIST:EMP:realsread"
886        WAVE emp_ints=$"root:Packages:NIST:EMP:integersread"
887        WAVE/T emp_text=$"root:Packages:NIST:EMP:textread"
888        WAVE drk_data=$"root:Packages:NIST:DRK:linear_data"
889        WAVE drk_reals=$"root:Packages:NIST:DRK:realsread"
890        WAVE drk_ints=$"root:Packages:NIST:DRK:integersread"
891        WAVE/T drk_text=$"root:Packages:NIST:DRK:textread"
892        WAVE cor_data=$"root:Packages:NIST:COR:linear_data"
893        WAVE/T cor_text=$"root:Packages:NIST:COR:textread"
894
895        // needed to propagate error
896        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy
897        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error"
898        WAVE emp_err =$"root:Packages:NIST:EMP:linear_data_error"
899        WAVE drk_err =$"root:Packages:NIST:DRK:linear_data_error"
900        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error"
901       
902        Variable sam_trans_err,emp_trans_err
903        sam_trans_err = sam_reals[41]
904        emp_trans_err = emp_reals[41]
905       
906        //get sam and bgd attenuation factors (DRK irrelevant)
907        String fileStr=""
908        Variable lambda,attenNo,sam_AttenFactor,emp_attenFactor
909        Variable tmonsam,fsam,femp,xshift,yshift,rsam,csam,remp,cemp,tmonemp
910        Variable wcen=0.001,tsam,temp,savmon_sam,time_sam,time_drk
911        Variable sam_atten_err,emp_atten_err
912        fileStr = sam_text[3]
913        lambda = sam_reals[26]
914        attenNo = sam_reals[3]
915        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err)
916        fileStr = emp_text[3]
917        lambda = emp_reals[26]
918        attenNo = emp_reals[3]
919        emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,emp_atten_err)
920       
921        //get relative monitor counts (should all be 10^8, since normalized in add step)
922        tmonsam = sam_reals[0]          //monitor count in SAM
923        tsam = sam_reals[4]             //SAM transmission
924        csam = sam_reals[16]            //x center
925        rsam = sam_reals[17]            //beam (x,y) define center of corrected field
926        tmonemp = emp_reals[0]          //monitor count in EMP
927        temp = emp_reals[4]                     //trans emp
928        cemp = emp_reals[16]            //beamcenter of EMP
929        remp = emp_reals[17]
930        savmon_sam=sam_reals[1]         //true monitor count in SAM
931        time_sam = sam_ints[2]          //count time SAM
932        time_drk = drk_ints[2]          //drk count time
933       
934        NVAR pixelsX = root:myGlobals:gNPixelsX
935        NVAR pixelsY = root:myGlobals:gNPixelsY
936        //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM
937        Make/D/O/N=(pixelsX,pixelsY) drk_temp,drk_tmp_err
938        drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam)
939        drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam)                 //temporarily rescale the error of DRK
940
941       
942        if(temp==0)
943                DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction"
944                temp=1
945        Endif
946       
947        //Print "rbgd,cbgd = ",rbgd,cbgd
948        // set up beamcenter shift, relative to SAM
949        xshift = cemp-csam
950        yshift = remp-rsam
951        if(abs(xshift) <= wcen)
952                xshift = 0
953        Endif
954        if(abs(yshift) <= wcen)
955                yshift = 0
956        Endif
957        //get shifted data arrays, relative to SAM
958        Make/D/O/N=(pixelsX,pixelsY) cor1,emp_temp,noadd_emp            //temp arrays
959        GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift)              //emp_temp is the EMP
960        //always ignore beamcenter shift for DRK
961       
962        //do the sam-bgd subtraction,  deposit result in cor1
963        fsam = 1
964        femp = tmonsam/tmonemp          //this should be ==1 since normalized files
965       
966        cor1 = fsam*sam_data/sam_AttenFactor - femp*(tsam/temp)*emp_temp/emp_AttenFactor
967        cor1 += drk_temp - drk_temp/sam_attenFactor
968        cor1 *= noadd_emp               //zeros out regions where arrays do not overlap, one otherwise
969
970// do the error propagation piecewise   
971        Duplicate/O sam_err, tmp_a, tmp_c, c_val
972        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
973       
974        tmp_c = (sam_trans_err*emp_data/(temp*emp_attenFactor))^2 + (emp_err*tsam/(temp*emp_attenFactor))^2
975        tmp_c += (tsam*emp_data*emp_trans_err/(temp*temp*emp_attenFactor))^2 + (tsam*emp_data*emp_atten_err/(temp*emp_attenFactor^2))^2//total of 6 terms
976       
977        cor_err = sqrt(tmp_a + tmp_c + drk_tmp_err^2)
978       
979        //we're done, get out w/no error
980        //set the COR data to the result
981        cor_data = cor1
982        cor_data_display = cor1
983
984        //update COR header
985        cor_text[1] = date() + " " + time()             //date + time stamp
986
987        KillWaves/Z cor1,emp_temp,noadd_emp,drk_temp
988        Killwaves/Z tmp_a,tmp_c,c_val,drk_tmp_err
989
990        SetDataFolder root:
991        Return(0)
992End
993
994// ONLY drk subtraction
995//
996xFunction V_CorrectMode_14()
997        //create the necessary wave references
998        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data"
999        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread"
1000        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread"
1001        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread"
1002        WAVE drk_data=$"root:Packages:NIST:DRK:linear_data"
1003        WAVE drk_reals=$"root:Packages:NIST:DRK:realsread"
1004        WAVE drk_ints=$"root:Packages:NIST:DRK:integersread"
1005        WAVE/T drk_text=$"root:Packages:NIST:DRK:textread"
1006        WAVE cor_data=$"root:Packages:NIST:COR:linear_data"
1007        WAVE/T cor_text=$"root:Packages:NIST:COR:textread"
1008
1009        // needed to propagate error
1010        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy
1011        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error"
1012        WAVE drk_err =$"root:Packages:NIST:DRK:linear_data_error"
1013        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error"
1014       
1015        Variable sam_trans_err
1016        sam_trans_err = sam_reals[41]
1017       
1018       
1019        //get sam and bgd attenuation factors
1020        String fileStr=""
1021        Variable lambda,attenNo,sam_AttenFactor,bgd_attenFactor
1022        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
1023        Variable wcen=0.001,time_drk,time_sam,savmon_sam,tsam
1024        Variable sam_atten_err
1025        fileStr = sam_text[3]
1026        lambda = sam_reals[26]
1027        attenNo = sam_reals[3]
1028        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err)
1029       
1030        //get relative monitor counts (should all be 10^8, since normalized in add step)
1031        tmonsam = sam_reals[0]          //monitor count in SAM
1032        tsam = sam_reals[4]             //SAM transmission
1033        csam = sam_reals[16]            //x center
1034        rsam = sam_reals[17]            //beam (x,y) define center of corrected field
1035
1036        savmon_sam=sam_reals[1]         //true monitor count in SAM
1037        time_sam = sam_ints[2]          //count time SAM
1038        time_drk = drk_ints[2]          //drk count time
1039       
1040        NVAR pixelsX = root:myGlobals:gNPixelsX
1041        NVAR pixelsY = root:myGlobals:gNPixelsY
1042        //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM
1043        Make/D/O/N=(pixelsX,pixelsY) drk_temp,drk_tmp_err
1044        drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam)
1045        drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam)                 //temporarily rescale the error of DRK
1046
1047        Make/D/O/N=(pixelsX,pixelsY) cor1       //temp arrays
1048        //always ignore the DRK center shift
1049       
1050        //do the subtraction,  deposit result in cor1
1051        fsam = 1
1052        fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
1053       
1054        //correct sam for attenuators, and do the same to drk, since it was scaled to sam count time
1055        cor1 = fsam*sam_data/sam_AttenFactor  - drk_temp/sam_attenFactor
1056
1057// do the error propagation piecewise   
1058        Duplicate/O sam_err, tmp_a
1059        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
1060
1061        cor_err = sqrt(tmp_a + drk_tmp_err^2)
1062       
1063        //we're done, get out w/no error
1064        //set the COR_data to the result
1065        cor_data = cor1
1066        cor_data_display = cor1
1067
1068        //update COR header
1069        cor_text[1] = date() + " " + time()             //date + time stamp
1070
1071        KillWaves/Z cor1,bgd_temp,noadd_bgd,drk_temp
1072        Killwaves/Z tmp_a,tmp_b,tmp_c,tmp_d,c_val,d_val,drk_tmp_err
1073
1074        SetDataFolder root:
1075        Return(0)
1076End
1077
1078
1079//
1080// For VSANS - this should go away. if there is a mismatch, don't try to fudge it.
1081//
1082//
1083//function to return the shifted contents of a data array for subtraction
1084//(SLOW) if ShiftSum is called
1085//data_in is input
1086//data_out is shifted matrix
1087//noadd_mat =1 if shift matrix is valid, =0 if no data
1088//
1089//if no shift is required, data_in is returned and noadd_mat =1 (all valid)
1090//
1091xFunction V_GetShiftedArray(data_in,data_out,noadd_mat,xshift,yshift)
1092        WAVE data_in,data_out,noadd_mat
1093        Variable xshift,yshift
1094
1095        Variable ii=0,jj=0
1096        noadd_mat = 1           //initialize to 1
1097       
1098        If((xshift != 0) || (yshift != 0))
1099//      If((abs(xshift) >= 0.01) || (abs(yshift) >= 0.01))                      //APR09 - loosen tolerance to handle ICE "precision"
1100                DoAlert 1,"Do you want to ignore the beam center mismatch?"
1101                if(V_flag==1)           //yes -> just go on
1102                        xshift=0
1103                        yshift=0
1104                endif
1105        else
1106                // "mismatch" is simply a python type conversion error
1107                xshift=0
1108                yshift=0
1109        endif
1110       
1111        If((xshift == 0) && (yshift == 0))
1112                data_out=data_in                //no change
1113                noadd_mat = 1                   //use all of the data
1114                return(0)
1115        endif
1116       
1117        NVAR pixelsX = root:myGlobals:gNPixelsX
1118        NVAR pixelsY = root:myGlobals:gNPixelsY
1119       
1120        Print "beamcenter shift x,y = ",xshift,yshift
1121        Make/O/N=1 noadd
1122        for(ii=0;ii<pixelsX;ii+=1)
1123                for(jj=0;jj<pixelsY;jj+=1)
1124                        //get the contribution of the shifted data
1125                        data_out[ii][jj] = ShiftSum(data_in,ii,jj,xshift,yshift,noadd)
1126                        if(noadd[0])
1127                                noadd_mat[ii][jj] = 0   //shift is off the detector
1128                        endif
1129                endfor
1130        endfor
1131        return(0)
1132End
1133
1134//
1135//utility function that checks if data exists in a data folder
1136//checks only for the existence of data in detector FL - no other waves
1137//
1138Function V_WorkDataExists(type)
1139        String type
1140       
1141        String destPath=""
1142        destPath = "root:Packages:NIST:VSANS:"+type+":entry:instrument:detector_FL:data"
1143        if(WaveExists($destpath) == 0)
1144                Print "There is no work file in "+type
1145                Return(1)               //error condition
1146        else
1147                // data exists, assume everything else is OK and proceed.
1148                return(0)
1149        Endif
1150End
1151
1152//////////////////
1153// bunch of utility junk to catch
1154// sample transmission = 1
1155// and handle (too many) options
1156//
1157xFunction V_GetNewTrans(oldTrans,type)
1158        Variable oldTrans
1159        String type
1160       
1161        Variable newTrans,newCode
1162        if (oldTrans!=1)
1163                return(oldTrans)                //get out now if trans != 1, don't change anything
1164        endif
1165        //get input from the user
1166        NewDataFolder/O root:myGlobals:tmp_trans
1167        Variable/G root:myGlobals:tmp_trans:inputTrans=0.92
1168        Variable/G root:myGlobals:tmp_trans:returnCode=0
1169        DoTransInput(type)
1170        NVAR inputTrans=root:myGlobals:tmp_trans:inputTrans
1171        NVAR code=root:myGlobals:tmp_trans:returnCode
1172        newTrans=inputTrans             //keep a copy before deleting everything
1173        newCode=code
1174        if(newCode==4)
1175                Abort "Aborting correction. Use the Transmission Panel to calculate transmissions"
1176        Endif
1177//      printf "You entered %g and the code is %g\r",newTrans,newCode
1178//      KillDataFolder root:tmp_trans
1179       
1180        if(newCode==1)
1181                Variable/G root:Packages:NIST:gDoTransCheck=0   //turn off checking
1182        endif
1183       
1184        if(newcode==2)          //user changed trans value
1185                return(newTrans)
1186        else
1187                return(oldTrans)        //All other cases, user did not change value
1188        endif
1189end
1190
1191xFunction V_IgnoreNowButton(ctrlName) : ButtonControl
1192        String ctrlName
1193       
1194//      Print "ignore now"
1195        NVAR val=root:myGlobals:tmp_trans:returnCode
1196        val=0           //code for ignore once
1197       
1198        DoWindow/K tmp_GetInputPanel            // Kill self
1199End
1200
1201xFunction V_DoTransInput(str)
1202        String str
1203       
1204        NewPanel /W=(150,50,361,294)
1205        DoWindow/C tmp_GetInputPanel            // Set to an unlikely name
1206        DrawText 15,23,"The "+str+" Transmission = 1"
1207        DrawText 15,43,"What do you want to do?"
1208        DrawText 15,125,"(Reset this in Preferences)"
1209        SetVariable setvar0,pos={20,170},size={160,17},limits={0,1,0.01}
1210        SetVariable setvar0,value= root:myGlobals:tmp_trans:inputTrans,title="New Transmission"
1211
1212        Button button0,pos={36,56},size={120,20},proc=IgnoreNowButton,title="Ignore This Time"
1213        Button button1,pos={36,86},size={120,20},proc=IgnoreAlwaysButtonProc,title="Ignore Always"
1214        Button button2,pos={36,143},size={120,20},proc=UseNewValueButtonProc,title="Use New Value"
1215        Button button3,pos={36,213},size={120,20},proc=AbortCorrectionButtonProc,title="Abort Correction"
1216        PauseForUser tmp_GetInputPanel
1217End
1218
1219xFunction V_IgnoreAlwaysButtonProc(ctrlName) : ButtonControl
1220        String ctrlName
1221
1222//      Print "ignore always"
1223        NVAR val=root:myGlobals:tmp_trans:returnCode
1224        val=1           //code for ignore always
1225        DoWindow/K tmp_GetInputPanel            // Kill self
1226End
1227
1228xFunction V_UseNewValueButtonProc(ctrlName) : ButtonControl
1229        String ctrlName
1230
1231//      Print "use new Value"
1232        NVAR val=root:myGlobals:tmp_trans:returnCode
1233        val=2           //code for use new Value
1234        DoWindow/K tmp_GetInputPanel            // Kill self
1235End
1236
1237xFunction V_AbortCorrectionButtonProc(ctrlName) : ButtonControl
1238        String ctrlName
1239
1240//      Print "Abort"
1241        NVAR val=root:myGlobals:tmp_trans:returnCode
1242        val=4           //code for abort
1243        DoWindow/K tmp_GetInputPanel            // Kill self
1244End
1245
Note: See TracBrowser for help on using the repository browser.