source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/Correct.ipf @ 890

Last change on this file since 890 was 794, checked in by srkline, 11 years ago

Lots of changes:
-2D resolution smearing
-2D error propagation

1) 2D resolution smearing has been corrected to use sigma (perp) correctly
rather than phi. This shows up in the quadrature loop in all of the 2D models
and in the Smear_2D "generic" function.

2) 1D resolutionn smearing is now corrected to account for the integration
range of +/- 3 sigma (99.73% of distribution). Now the result is divided by
0.9973 to rescale it to the correct value.

3) Smeared models are now AAO to improve speed and to allow easier use with
functions that are inherently AAO. No changes are needed, since the call is
behind the scenes, replacing Smear_Model_N() with Smear_Model_N_AAO().

4) in PlotUtils_2D added functions to re-bin the QxQy? data into a 1D format
BinQxQy_to_1D(). This also re-bins the errors in two ways, adding the per-pixel
errors in quadrature, or the deviation from the mean of the intensity. Manually
editing the intensity declaration allows 2D->1D binning of smeared models.

5) Per-pixel propagation of errors has been added through the entire data
reduction sequence. Data errors are generated on loading using Poisson
statistics (specifically tailored for accuracy at low counts) and then is
propagated through each manipulation of the data using standard error
propagation. The error matrix is linear_data_error. As a by-product, all of
the math operations on data are explicitly done on linear_data, to avoid
any potential mistakes of log/linear scaling. Results of this propagation
largely match J. Barker's /ERROR routines from the VAX, with some differences
at low data count values (as expected) and at higher count values near the
beam stop (somewhat unexpected). This per-pixel error is ouput in the QxQy_ASCII
data files. NO CHANGE has been made to the 1D data, which uses the deviation from
the mean as the error - since this is correct.

6) Added tables for the uncertainty in attenuator transmission (from JGB)

7) added two more REAL values to the VAX header to store information
necessary for error propagation. These are couting error that are part of
the transmission error and of the absolute scaling error. These add Read/Write?
functions in NCNR_DataReadWrite

The transmission error (one standard deviation) is written at byte 396 (4 bytes)

RealsRead?[41]

The Box Sum error (one standard deviation) is written at byte 400 (4 bytes)

RealsRead?[42]

File size: 50.7 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=6.1
4
5//****************************
6// Vers 1.2 090501
7//
8//****************************
9//
10// Procedures to perform the "Correct" step during data reduction
11//
12// - ther is olny one procedure to perform the subtractions, and a single
13// parameter flags which subtractions are to be done. Different numbers of
14// attenuators during scattering runs are corrected as described in John's memo,
15// with the note that ONLY method (3) is used, which assumes that 'diffuse' scattering
16// is dominant over 'dark current' (note that 'dark current' = shutter CLOSED)
17//
18//
19//do the CORRECT step based on the answers to emp and bkg subtraction
20        //by setting the proper"mode"
21        //1 = both emp and bgd subtraction
22        //2 = only bgd subtraction
23        //3 = only emp subtraction
24        //4 = no subtraction
25        //additional modes 091301
26        //11 = emp, bgd, drk
27        //12 = bgd and drk
28        //13 = emp and drk
29        //14 = no subtractions
30        //
31//********************************
32
33//unused test procedure for Correct() function
34//must be updated to include "mode" parameter before re-use
35//
36Proc CorrectData()
37
38        Variable err
39        String type
40       
41        err = Correct()         
42       
43        if(err)
44                Abort "error in Correct"
45        endif
46       
47        //contents are always dumped to COR
48        type = "COR"
49       
50        //need to update the display with "data" from the correct dataFolder
51        //reset the current displaytype to "type"
52        String/G root:myGlobals:gDataDisplayType=Type
53       
54        fRawWindowHook()
55       
56End
57
58
59//mode describes the type of subtraction that is to be done
60//1 = both emp and bgd subtraction
61//2 = only bgd subtraction
62//3 = only emp subtraction
63//4 = no subtraction
64//
65// + 10 indicates that WORK.DRK is to be used
66//
67//091301 version
68//now simple dispatches to the correct subtraction - logic was too
69//involved to do in one function - unclear and error-prone
70//
71// 081203 version
72// checks for trans==1 in SAM and EMP before dispatching
73// and asks for new value if desired
74//
75Function Correct(mode)
76        Variable mode
77       
78        Variable err=0,trans,newTrans
79       
80        //switch and dispatch based on the required subtractions
81        // always check for SAM data
82        err = WorkDataExists("SAM")
83        if(err==1)
84                return(err)
85        endif
86       
87       
88        //check for trans==1
89        NVAR doCheck=root:Packages:NIST:gDoTransCheck
90        Wave/Z samR=root:Packages:NIST:SAM:RealsRead
91        Wave/Z empR=root:Packages:NIST:EMP:RealsRead
92        if(doCheck)
93                trans = samR[4]
94                newTrans=GetNewTrans(trans,"SAM")               //will change value if necessary
95                if(numtype(newTrans)==0)
96                        samR[4] = newTrans              //avoid user abort assigning NaN
97                endif
98                if(trans != newTrans)
99                        print "Using SAM trans = ",samR[4]
100                endif
101        endif
102       
103        //copy SAM information to COR, wiping out the old contents of the COR folder first
104        //do this even if no correction is dispatched (if incorrect mode)
105        // the Copy converts "SAM" to linear scale, so "COR" is then linear too
106        err = CopyWorkContents("SAM","COR")     
107        if(err==1)
108                Abort "No data in SAM, abort from Correct()"
109        endif
110       
111//      Print "dispatching to mode = ",mode
112        switch(mode)
113                case 1:
114                        err = WorkDataExists("EMP")
115                        if(err==1)
116                                return(err)
117                        Endif
118                        if(doCheck)
119                                trans = empR[4]
120                                newTrans=GetNewTrans(trans,"EMP")               //will change value if necessary
121                                if(numtype(newTrans)==0)
122                                        empR[4] = newTrans
123                                endif
124                                if(trans != newTrans)
125                                        print "Using EMP trans = ",empR[4]
126                                endif
127                        endif
128                        err = WorkDataExists("BGD")
129                        if(err==1)
130                                return(err)
131                        Endif
132                        err = CorrectMode_1()
133                        break
134                case 2:
135                        err = WorkDataExists("BGD")
136                        if(err==1)
137                                return(err)
138                        Endif
139                        err = CorrectMode_2()
140                        break
141                case 3:
142                        err = WorkDataExists("EMP")
143                        if(err==1)
144                                return(err)
145                        Endif
146                        if(doCheck)
147                                trans = empR[4]
148                                newTrans=GetNewTrans(trans,"EMP")               //will change value if necessary
149                                if(numtype(newTrans)==0)
150                                        empR[4] = newTrans
151                                endif
152                                if(trans != newTrans)
153                                        print "Using EMP trans = ",empR[4]
154                                endif
155                        endif
156                        err = CorrectMode_3()
157                        break
158                case 4:
159                        err = CorrectMode_4()
160                        break
161                case 11:
162                        err = WorkDataExists("EMP")
163                        if(err==1)
164                                return(err)
165                        Endif
166                        if(doCheck)
167                                trans = empR[4]
168                                newTrans=GetNewTrans(trans,"EMP")               //will change value if necessary
169                                if(numtype(newTrans)==0)
170                                        empR[4] = newTrans
171                                endif
172                                if(trans != newTrans)
173                                        print "Using EMP trans = ",empR[4]
174                                endif
175                        endif
176                        err = WorkDataExists("BGD")
177                        if(err==1)
178                                return(err)
179                        Endif
180                        err = WorkDataExists("DRK")
181                        if(err==1)
182                                return(err)
183                        Endif
184                        err = CorrectMode_11()
185                        break
186                case 12:
187                        err = WorkDataExists("BGD")
188                        if(err==1)
189                                return(err)
190                        Endif
191                        err = WorkDataExists("DRK")
192                        if(err==1)
193                                return(err)
194                        Endif
195                        err = CorrectMode_12()
196                        break
197                case 13:
198                        err = WorkDataExists("EMP")
199                        if(err==1)
200                                return(err)
201                        Endif
202                        if(doCheck)
203                                trans = empR[4]
204                                newTrans=GetNewTrans(trans,"EMP")               //will change value if necessary
205                                if(numtype(newTrans)==0)
206                                        empR[4] = newTrans
207                                endif
208                                if(trans != newTrans)
209                                        print "Using EMP trans = ",empR[4]
210                                endif
211                        endif
212                        err = WorkDataExists("DRK")
213                        if(err==1)
214                                return(err)
215                        Endif
216                        err = CorrectMode_13()
217                        break
218                case 14:
219                        err = WorkDataExists("DRK")
220                        if(err==1)
221                                return(err)
222                        Endif
223                        err = CorrectMode_14()
224                        break
225                default:        //something wrong
226                        Print "Incorrect mode in Correct()"
227                        return(1)       //error
228        endswitch
229
230        //calculation attempted, return the result
231        return(err)     
232End
233
234// subtraction of bot EMP and BGD from SAM
235// data exists, checked by dispatch routine
236//
237// this is the most common use
238// March 2011 added error propagation
239//                                      added explicit reference to use linear_data, instead of trusting that data
240//                                      was freshly loaded. added final copy of cor result to cor:data and cor:linear_data
241//
242Function CorrectMode_1()
243       
244        //create the necessary wave references
245        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data"
246        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread"
247        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread"
248        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread"
249        WAVE bgd_data=$"root:Packages:NIST:BGD:linear_data"
250        WAVE bgd_reals=$"root:Packages:NIST:BGD:realsread"
251        WAVE bgd_ints=$"root:Packages:NIST:BGD:integersread"
252        WAVE/T bgd_text=$"root:Packages:NIST:BGD:textread"
253        WAVE emp_data=$"root:Packages:NIST:EMP:linear_data"
254        WAVE emp_reals=$"root:Packages:NIST:EMP:realsread"
255        WAVE emp_ints=$"root:Packages:NIST:EMP:integersread"
256        WAVE/T emp_text=$"root:Packages:NIST:EMP:textread"
257        WAVE cor_data=$"root:Packages:NIST:COR:linear_data"
258        WAVE/T cor_text=$"root:Packages:NIST:COR:textread"
259       
260        // needed to propagate error
261        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy
262        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error"
263        WAVE bgd_err =$"root:Packages:NIST:BGD:linear_data_error"
264        WAVE emp_err =$"root:Packages:NIST:EMP:linear_data_error"
265        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error"
266       
267        Variable sam_trans_err,emp_trans_err
268        sam_trans_err = sam_reals[41]
269        emp_trans_err = emp_reals[41]
270       
271       
272        //get sam and bgd attenuation factors
273        String fileStr=""
274        Variable lambda,attenNo,sam_AttenFactor,bgd_attenFactor,emp_AttenFactor
275        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
276        Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp
277        Variable sam_atten_err,emp_atten_err,bgd_atten_err
278        fileStr = sam_text[3]
279        lambda = sam_reals[26]
280        attenNo = sam_reals[3]
281        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err)
282        fileStr = bgd_text[3]
283        lambda = bgd_reals[26]
284        attenNo = bgd_reals[3]
285        bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,bgd_atten_err)
286        fileStr = emp_text[3]
287        lambda = emp_reals[26]
288        attenNo = emp_reals[3]
289        emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,emp_atten_err)
290       
291        //get relative monitor counts (should all be 10^8, since normalized in add step)
292        tmonsam = sam_reals[0]          //monitor count in SAM
293        tsam = sam_reals[4]             //SAM transmission
294        csam = sam_reals[16]            //x center
295        rsam = sam_reals[17]            //beam (x,y) define center of corrected field
296        tmonbgd = bgd_reals[0]          //monitor count in BGD
297        cbgd = bgd_reals[16]
298        rbgd = bgd_reals[17]
299        tmonemp = emp_reals[0]          //monitor count in EMP
300        temp = emp_reals[4]                     //trans emp
301        cemp = emp_reals[16]            //beamcenter of EMP
302        remp = emp_reals[17]
303       
304        if(temp==0)
305                DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction"
306                temp=1
307        Endif
308       
309        NVAR pixelsX = root:myGlobals:gNPixelsX
310        NVAR pixelsY = root:myGlobals:gNPixelsY
311       
312        //get the shifted data arrays, EMP and BGD, each relative to SAM
313        Make/D/O/N=(pixelsX,pixelsY) cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp
314        xshift = cbgd-csam
315        yshift = rbgd-rsam
316        if(abs(xshift) <= wcen)
317                xshift = 0
318        Endif
319        if(abs(yshift) <= wcen)
320                yshift = 0
321        Endif
322        GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift)              //bgd_temp
323       
324        xshift = cemp-csam
325        yshift = remp-rsam
326        if(abs(xshift) <= wcen)
327                xshift = 0
328        Endif
329        if(abs(yshift) <= wcen)
330                yshift = 0
331        Endif
332        GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift)              //emp_temp
333
334        //do the subtraction
335        fsam=1
336        femp = tmonsam/tmonemp          //this should be ==1 since normalized files
337        fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
338        cor1 = fsam*sam_data/sam_attenFactor - fbgd*bgd_temp/bgd_attenFactor
339        cor1 -= (tsam/temp)*(femp*emp_temp/emp_attenFactor - fbgd*bgd_temp/bgd_attenFactor)
340        cor1 *= noadd_bgd*noadd_emp             //zero out the array mismatch values
341
342// do the error propagation piecewise   
343        Duplicate/O sam_err, tmp_a, tmp_b, tmp_c, tmp_d,c_val,d_val
344        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
345       
346        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
347
348        tmp_c = (sam_trans_err/temp)^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2
349        tmp_c += (tsam/temp^2)^2*emp_trans_err^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2
350       
351        tmp_d = (tsam/(temp*emp_attenFactor))^2*(emp_err)^2 + (tsam*emp_data/(temp*emp_attenFactor^2))^2*(emp_atten_err)^2
352
353        cor_err = sqrt(tmp_a + tmp_b + tmp_c + tmp_d)
354       
355        //we're done, get out w/no error
356        //set the COR data and linear_data to the result
357        cor_data = cor1
358        cor_data_display = cor1
359       
360        //update COR header
361        cor_text[1] = date() + " " + time()             //date + time stamp
362
363        KillWaves/Z cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp
364        Killwaves/Z tmp_a,tmp_b,tmp_c,tmp_d,c_val,d_val
365        SetDataFolder root:
366        Return(0)
367End
368
369//background only
370// existence of data checked by dispatching routine
371// data has already been copied to COR folder
372Function CorrectMode_2()
373
374        //create the necessary wave references
375        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data"
376        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread"
377        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread"
378        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread"
379        WAVE bgd_data=$"root:Packages:NIST:BGD:linear_data"
380        WAVE bgd_reals=$"root:Packages:NIST:BGD:realsread"
381        WAVE bgd_ints=$"root:Packages:NIST:BGD:integersread"
382        WAVE/T bgd_text=$"root:Packages:NIST:BGD:textread"
383        WAVE cor_data=$"root:Packages:NIST:COR:linear_data"
384        WAVE/T cor_text=$"root:Packages:NIST:COR:textread"
385
386        // needed to propagate error
387        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy
388        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error"
389        WAVE bgd_err =$"root:Packages:NIST:BGD:linear_data_error"
390        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error"
391       
392        Variable sam_trans_err
393        sam_trans_err = sam_reals[41]
394
395       
396        //get sam and bgd attenuation factors
397        String fileStr=""
398        Variable lambda,attenNo,sam_AttenFactor,bgd_attenFactor
399        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
400        Variable wcen=0.001
401        Variable sam_atten_err,bgd_atten_err
402        fileStr = sam_text[3]
403        lambda = sam_reals[26]
404        attenNo = sam_reals[3]
405        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err)
406        fileStr = bgd_text[3]
407        lambda = bgd_reals[26]
408        attenNo = bgd_reals[3]
409        bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,bgd_atten_err)
410       
411        //Print "atten = ",sam_attenFactor,bgd_attenFactor
412       
413        //get relative monitor counts (should all be 10^8, since normalized in add step)
414        tmonsam = sam_reals[0]          //monitor count in SAM
415        csam = sam_reals[16]            //x center
416        rsam = sam_reals[17]            //beam (x,y) define center of corrected field
417        tmonbgd = bgd_reals[0]          //monitor count in BGD
418        cbgd = bgd_reals[16]
419        rbgd = bgd_reals[17]
420
421        // set up beamcenter shift, relative to SAM
422        xshift = cbgd-csam
423        yshift = rbgd-rsam
424        if(abs(xshift) <= wcen)
425                xshift = 0
426        Endif
427        if(abs(yshift) <= wcen)
428                yshift = 0
429        Endif
430       
431        NVAR pixelsX = root:myGlobals:gNPixelsX
432        NVAR pixelsY = root:myGlobals:gNPixelsY
433        //get shifted data arrays, relative to SAM
434        Make/D/O/N=(pixelsX,pixelsY) cor1,bgd_temp,noadd_bgd            //temp arrays
435        GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift)              //bgd_temp is the BGD
436       
437        //do the sam-bgd subtraction,  deposit result in cor1
438        fsam = 1
439        fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
440       
441        //print "fsam,fbgd = ",fsam,fbgd
442       
443        cor1 = fsam*sam_data/sam_AttenFactor - fbgd*bgd_temp/bgd_AttenFactor
444        cor1 *= noadd_bgd               //zeros out regions where arrays do not overlap, one otherwise
445
446// do the error propagation piecewise   
447        Duplicate/O sam_err, tmp_a, tmp_b
448        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
449       
450        tmp_b = (bgd_err/bgd_attenFactor)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2              //sig b ^2
451
452        cor_err = sqrt(tmp_a + tmp_b)
453
454
455        //we're done, get out w/no error
456        //set the COR_data to the result
457        cor_data = cor1
458        cor_data_display = cor1
459
460        //update COR header
461        cor_text[1] = date() + " " + time()             //date + time stamp
462
463        KillWaves/Z cor1,bgd_temp,noadd_bgd
464        Killwaves/Z tmp_a,tmp_b
465
466        SetDataFolder root:
467        Return(0)
468End
469
470// empty subtraction only
471// data does exist, checked by dispatch routine
472//
473Function CorrectMode_3()
474        //create the necessary wave references
475        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data"
476        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread"
477        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread"
478        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread"
479        WAVE emp_data=$"root:Packages:NIST:EMP:linear_data"
480        WAVE emp_reals=$"root:Packages:NIST:EMP:realsread"
481        WAVE emp_ints=$"root:Packages:NIST:EMP:integersread"
482        WAVE/T emp_text=$"root:Packages:NIST:EMP:textread"
483        WAVE cor_data=$"root:Packages:NIST:COR:linear_data"
484        WAVE/T cor_text=$"root:Packages:NIST:COR:textread"
485       
486        // needed to propagate error
487        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy
488        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error"
489        WAVE emp_err =$"root:Packages:NIST:EMP:linear_data_error"
490        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error"
491       
492        Variable sam_trans_err,emp_trans_err
493        sam_trans_err = sam_reals[41]
494        emp_trans_err = emp_reals[41]   
495       
496        //get sam and bgd attenuation factors
497        String fileStr=""
498        Variable lambda,attenNo,sam_AttenFactor,emp_attenFactor
499        Variable tmonsam,fsam,femp,xshift,yshift,rsam,csam,remp,cemp,tmonemp
500        Variable wcen=0.001,tsam,temp
501        Variable sam_atten_err,emp_atten_err
502        fileStr = sam_text[3]
503        lambda = sam_reals[26]
504        attenNo = sam_reals[3]
505        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err)
506        fileStr = emp_text[3]
507        lambda = emp_reals[26]
508        attenNo = emp_reals[3]
509        emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,emp_atten_err)
510       
511        //get relative monitor counts (should all be 10^8, since normalized in add step)
512        tmonsam = sam_reals[0]          //monitor count in SAM
513        tsam = sam_reals[4]             //SAM transmission
514        csam = sam_reals[16]            //x center
515        rsam = sam_reals[17]            //beam (x,y) define center of corrected field
516        tmonemp = emp_reals[0]          //monitor count in EMP
517        temp = emp_reals[4]                     //trans emp
518        cemp = emp_reals[16]            //beamcenter of EMP
519        remp = emp_reals[17]
520       
521        if(temp==0)
522                DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction"
523                temp=1
524        Endif
525       
526        //Print "rbgd,cbgd = ",rbgd,cbgd
527        // set up beamcenter shift, relative to SAM
528        xshift = cemp-csam
529        yshift = remp-rsam
530        if(abs(xshift) <= wcen)
531                xshift = 0
532        Endif
533        if(abs(yshift) <= wcen)
534                yshift = 0
535        Endif
536       
537        NVAR pixelsX = root:myGlobals:gNPixelsX
538        NVAR pixelsY = root:myGlobals:gNPixelsY
539        //get shifted data arrays, relative to SAM
540        Make/D/O/N=(pixelsX,pixelsY) cor1,emp_temp,noadd_emp            //temp arrays
541        GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift)              //emp_temp is the EMP
542       
543        //do the sam-bgd subtraction,  deposit result in cor1
544        fsam = 1
545        femp = tmonsam/tmonemp          //this should be ==1 since normalized files
546       
547        cor1 = fsam*sam_data/sam_AttenFactor - femp*(tsam/temp)*emp_temp/emp_AttenFactor
548        cor1 *= noadd_emp               //zeros out regions where arrays do not overlap, one otherwise
549
550// do the error propagation piecewise   
551        Duplicate/O sam_err, tmp_a, tmp_c ,c_val
552        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
553       
554        tmp_c = (sam_trans_err*emp_data/(temp*emp_attenFactor))^2 + (emp_err*tsam/(temp*emp_attenFactor))^2
555        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
556
557        cor_err = sqrt(tmp_a + tmp_c)
558       
559        //we're done, get out w/no error
560        //set the COR data to the result
561        cor_data = cor1
562        cor_data_display = cor1
563
564        //update COR header
565        cor_text[1] = date() + " " + time()             //date + time stamp
566
567        KillWaves/Z cor1,emp_temp,noadd_emp
568        Killwaves/Z tmp_a,tmp_c,c_val
569
570        SetDataFolder root:
571        Return(0)
572End
573
574// NO subtraction - simply rescales for attenuators
575// SAM data does exist, checked by dispatch routine
576//
577// !! moves data to COR folder, since the data has been corrected, by rescaling
578//
579Function CorrectMode_4()
580        //create the necessary wave references
581        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data"
582        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread"
583        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread"
584        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread"
585
586        WAVE cor_data=$"root:Packages:NIST:COR:linear_data"
587        WAVE/T cor_text=$"root:Packages:NIST:COR:textread"
588
589        // needed to propagate error
590        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy
591        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error"
592        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error"
593               
594        //get sam and bgd attenuation factors
595        String fileStr=""
596        Variable lambda,attenNo,sam_AttenFactor,sam_atten_err
597        fileStr = sam_text[3]
598        lambda = sam_reals[26]
599        attenNo = sam_reals[3]
600        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err)
601
602        NVAR pixelsX = root:myGlobals:gNPixelsX
603        NVAR pixelsY = root:myGlobals:gNPixelsY
604        Make/D/O/N=(pixelsX,pixelsY) cor1
605       
606        cor1 = sam_data/sam_AttenFactor         //simply rescale the data
607
608// do the error propagation piecewise   
609        Duplicate/O sam_err, tmp_a
610        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
611
612        cor_err = sqrt(tmp_a)
613       
614
615        //we're done, get out w/no error
616        //set the COR data to the result
617        cor_data = cor1
618        cor_data_display = cor1
619
620        //update COR header
621        cor_text[1] = date() + " " + time()             //date + time stamp
622
623        KillWaves/Z cor1
624        Killwaves/Z tmp_a
625
626        SetDataFolder root:
627        Return(0)
628End
629
630Function CorrectMode_11()
631        //create the necessary wave references
632        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data"
633        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread"
634        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread"
635        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread"
636        WAVE bgd_data=$"root:Packages:NIST:BGD:linear_data"
637        WAVE bgd_reals=$"root:Packages:NIST:BGD:realsread"
638        WAVE bgd_ints=$"root:Packages:NIST:BGD:integersread"
639        WAVE/T bgd_text=$"root:Packages:NIST:BGD:textread"
640        WAVE emp_data=$"root:Packages:NIST:EMP:linear_data"
641        WAVE emp_reals=$"root:Packages:NIST:EMP:realsread"
642        WAVE emp_ints=$"root:Packages:NIST:EMP:integersread"
643        WAVE/T emp_text=$"root:Packages:NIST:EMP:textread"
644        WAVE drk_data=$"root:Packages:NIST:DRK:linear_data"
645        WAVE drk_reals=$"root:Packages:NIST:DRK:realsread"
646        WAVE drk_ints=$"root:Packages:NIST:DRK:integersread"
647        WAVE/T drk_text=$"root:Packages:NIST:DRK:textread"
648        WAVE cor_data=$"root:Packages:NIST:COR:linear_data"
649        WAVE/T cor_text=$"root:Packages:NIST:COR:textread"
650
651        // needed to propagate error
652        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy
653        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error"
654        WAVE bgd_err =$"root:Packages:NIST:BGD:linear_data_error"
655        WAVE emp_err =$"root:Packages:NIST:EMP:linear_data_error"
656        WAVE drk_err =$"root:Packages:NIST:DRK:linear_data_error"
657        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error"
658       
659        Variable sam_trans_err,emp_trans_err
660        sam_trans_err = sam_reals[41]
661        emp_trans_err = emp_reals[41]
662       
663        //get sam and bgd attenuation factors
664        String fileStr=""
665        Variable lambda,attenNo,sam_AttenFactor,bgd_attenFactor,emp_AttenFactor
666        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
667        Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp,time_sam,time_drk,savmon_sam
668        Variable sam_atten_err,bgd_atten_err,emp_atten_err
669        fileStr = sam_text[3]
670        lambda = sam_reals[26]
671        attenNo = sam_reals[3]
672        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err)
673        fileStr = bgd_text[3]
674        lambda = bgd_reals[26]
675        attenNo = bgd_reals[3]
676        bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,bgd_atten_err)
677        fileStr = emp_text[3]
678        lambda = emp_reals[26]
679        attenNo = emp_reals[3]
680        emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,emp_atten_err)
681       
682        //get relative monitor counts (should all be 10^8, since normalized in add step)
683        tmonsam = sam_reals[0]          //monitor count in SAM
684        tsam = sam_reals[4]             //SAM transmission
685        csam = sam_reals[16]            //x center
686        rsam = sam_reals[17]            //beam (x,y) define center of corrected field
687        tmonbgd = bgd_reals[0]          //monitor count in BGD
688        cbgd = bgd_reals[16]
689        rbgd = bgd_reals[17]
690        tmonemp = emp_reals[0]          //monitor count in EMP
691        temp = emp_reals[4]                     //trans emp
692        cemp = emp_reals[16]            //beamcenter of EMP
693        remp = emp_reals[17]
694        savmon_sam=sam_reals[1]         //true monitor count in SAM
695        time_sam = sam_ints[2]          //count time SAM
696        time_drk = drk_ints[2]          //drk count time
697       
698        NVAR pixelsX = root:myGlobals:gNPixelsX
699        NVAR pixelsY = root:myGlobals:gNPixelsY
700        //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM
701        Make/D/O/N=(pixelsX,pixelsY) drk_temp, drk_tmp_err
702        drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam)
703        drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam)                 //temporarily rescale the error of DRK
704       
705        if(temp==0)
706                DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction"
707                temp=1
708        Endif
709       
710        //get the shifted data arrays, EMP and BGD, each relative to SAM
711        Make/D/O/N=(pixelsX,pixelsY) cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp
712        xshift = cbgd-csam
713        yshift = rbgd-rsam
714        if(abs(xshift) <= wcen)
715                xshift = 0
716        Endif
717        if(abs(yshift) <= wcen)
718                yshift = 0
719        Endif
720        GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift)              //bgd_temp
721       
722        xshift = cemp-csam
723        yshift = remp-rsam
724        if(abs(xshift) <= wcen)
725                xshift = 0
726        Endif
727        if(abs(yshift) <= wcen)
728                yshift = 0
729        Endif
730        GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift)              //emp_temp
731        //always ignore the DRK center shift
732       
733        //do the subtraction
734        fsam=1
735        femp = tmonsam/tmonemp          //this should be ==1 since normalized files
736        fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
737        cor1 = fsam*sam_data/sam_attenFactor
738        cor1 -= (tsam/temp)*(femp*emp_temp/emp_attenFactor - fbgd*bgd_temp/bgd_attenFactor)
739        cor1 -= (fbgd*bgd_temp/bgd_attenFactor - drk_temp)
740        cor1 -= drk_temp/sam_attenFactor
741        cor1 *= noadd_bgd*noadd_emp             //zero out the array mismatch values
742       
743// do the error propagation piecewise   
744        Duplicate/O sam_err, tmp_a, tmp_b, tmp_c, tmp_d,c_val,d_val
745        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
746       
747        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
748
749        tmp_c = (sam_trans_err/temp)^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2
750        tmp_c += (tsam/temp^2)^2*emp_trans_err^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2
751       
752        tmp_d = (tsam/(temp*emp_attenFactor))^2*(emp_err)^2 + (tsam*emp_data/(temp*emp_attenFactor^2))^2*(emp_atten_err)^2
753
754        cor_err = sqrt(tmp_a + tmp_b + tmp_c + tmp_d + drk_tmp_err^2)
755       
756        //we're done, get out w/no error
757        //set the COR data to the result
758        cor_data = cor1
759        cor_data_display = cor1
760
761        //update COR header
762        cor_text[1] = date() + " " + time()             //date + time stamp
763
764        KillWaves/Z cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp,drk_temp
765        Killwaves/Z tmp_a,tmp_b,tmp_c,tmp_d,c_val,d_val,drk_tmp_err
766
767        SetDataFolder root:
768        Return(0)
769End
770
771//bgd and drk subtraction
772//
773Function CorrectMode_12()
774        //create the necessary wave references
775        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data"
776        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread"
777        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread"
778        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread"
779        WAVE bgd_data=$"root:Packages:NIST:BGD:linear_data"
780        WAVE bgd_reals=$"root:Packages:NIST:BGD:realsread"
781        WAVE bgd_ints=$"root:Packages:NIST:BGD:integersread"
782        WAVE/T bgd_text=$"root:Packages:NIST:BGD:textread"
783        WAVE drk_data=$"root:Packages:NIST:DRK:linear_data"
784        WAVE drk_reals=$"root:Packages:NIST:DRK:realsread"
785        WAVE drk_ints=$"root:Packages:NIST:DRK:integersread"
786        WAVE/T drk_text=$"root:Packages:NIST:DRK:textread"
787        WAVE cor_data=$"root:Packages:NIST:COR:linear_data"
788        WAVE/T cor_text=$"root:Packages:NIST:COR:textread"
789
790        // needed to propagate error
791        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy
792        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error"
793        WAVE bgd_err =$"root:Packages:NIST:BGD:linear_data_error"
794        WAVE drk_err =$"root:Packages:NIST:DRK:linear_data_error"
795        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error"
796       
797        Variable sam_trans_err
798        sam_trans_err = sam_reals[41]
799       
800       
801        //get sam and bgd attenuation factors
802        String fileStr=""
803        Variable lambda,attenNo,sam_AttenFactor,bgd_attenFactor
804        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
805        Variable wcen=0.001,time_drk,time_sam,savmon_sam,tsam
806        Variable sam_atten_err,bgd_atten_err
807        fileStr = sam_text[3]
808        lambda = sam_reals[26]
809        attenNo = sam_reals[3]
810        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err)
811        fileStr = bgd_text[3]
812        lambda = bgd_reals[26]
813        attenNo = bgd_reals[3]
814        bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,bgd_atten_err)
815       
816        //get relative monitor counts (should all be 10^8, since normalized in add step)
817        tmonsam = sam_reals[0]          //monitor count in SAM
818        tsam = sam_reals[4]             //SAM transmission
819        csam = sam_reals[16]            //x center
820        rsam = sam_reals[17]            //beam (x,y) define center of corrected field
821        tmonbgd = bgd_reals[0]          //monitor count in BGD
822        cbgd = bgd_reals[16]
823        rbgd = bgd_reals[17]
824        savmon_sam=sam_reals[1]         //true monitor count in SAM
825        time_sam = sam_ints[2]          //count time SAM
826        time_drk = drk_ints[2]          //drk count time
827       
828        NVAR pixelsX = root:myGlobals:gNPixelsX
829        NVAR pixelsY = root:myGlobals:gNPixelsY
830        //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM
831        Make/D/O/N=(pixelsX,pixelsY) drk_temp,drk_tmp_err
832        drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam)
833        drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam)                 //temporarily rescale the error of DRK
834
835        // set up beamcenter shift, relative to SAM
836        xshift = cbgd-csam
837        yshift = rbgd-rsam
838        if(abs(xshift) <= wcen)
839                xshift = 0
840        Endif
841        if(abs(yshift) <= wcen)
842                yshift = 0
843        Endif
844        //get shifted data arrays, relative to SAM
845        Make/D/O/N=(pixelsX,pixelsY) cor1,bgd_temp,noadd_bgd            //temp arrays
846        GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift)              //bgd_temp is the BGD
847        //always ignore the DRK center shift
848       
849        //do the sam-bgd subtraction,  deposit result in cor1
850        fsam = 1
851        fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
852       
853        cor1 = fsam*sam_data/sam_AttenFactor + fbgd*tsam*bgd_temp/bgd_AttenFactor
854        cor1 += -1*(fbgd*bgd_temp/bgd_attenFactor - drk_temp) - drk_temp/sam_attenFactor
855        cor1 *= noadd_bgd               //zeros out regions where arrays do not overlap, one otherwise
856
857// do the error propagation piecewise   
858        Duplicate/O sam_err, tmp_a, tmp_b
859        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
860       
861        tmp_b = (bgd_err/bgd_attenFactor)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2              //sig b ^2
862
863        cor_err = sqrt(tmp_a + tmp_b + drk_tmp_err^2)
864
865        //we're done, get out w/no error
866        //set the COR_data to the result
867        cor_data = cor1
868        cor_data_display = cor1
869
870        //update COR header
871        cor_text[1] = date() + " " + time()             //date + time stamp
872
873        KillWaves/Z cor1,bgd_temp,noadd_bgd,drk_temp
874        Killwaves/Z tmp_a,tmp_b,drk_tmp_err
875
876        SetDataFolder root:
877        Return(0)
878End
879
880//EMP and DRK subtractions
881// all data exists, DRK is on a time basis (noNorm)
882//scale DRK by monitor count scaling factor and the ratio of couting times
883//to place the DRK file on equal footing
884Function CorrectMode_13()
885        //create the necessary wave references
886        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data"
887        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread"
888        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread"
889        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread"
890        WAVE emp_data=$"root:Packages:NIST:EMP:linear_data"
891        WAVE emp_reals=$"root:Packages:NIST:EMP:realsread"
892        WAVE emp_ints=$"root:Packages:NIST:EMP:integersread"
893        WAVE/T emp_text=$"root:Packages:NIST:EMP:textread"
894        WAVE drk_data=$"root:Packages:NIST:DRK:linear_data"
895        WAVE drk_reals=$"root:Packages:NIST:DRK:realsread"
896        WAVE drk_ints=$"root:Packages:NIST:DRK:integersread"
897        WAVE/T drk_text=$"root:Packages:NIST:DRK:textread"
898        WAVE cor_data=$"root:Packages:NIST:COR:linear_data"
899        WAVE/T cor_text=$"root:Packages:NIST:COR:textread"
900
901        // needed to propagate error
902        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy
903        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error"
904        WAVE emp_err =$"root:Packages:NIST:EMP:linear_data_error"
905        WAVE drk_err =$"root:Packages:NIST:DRK:linear_data_error"
906        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error"
907       
908        Variable sam_trans_err,emp_trans_err
909        sam_trans_err = sam_reals[41]
910        emp_trans_err = emp_reals[41]
911       
912        //get sam and bgd attenuation factors (DRK irrelevant)
913        String fileStr=""
914        Variable lambda,attenNo,sam_AttenFactor,emp_attenFactor
915        Variable tmonsam,fsam,femp,xshift,yshift,rsam,csam,remp,cemp,tmonemp
916        Variable wcen=0.001,tsam,temp,savmon_sam,time_sam,time_drk
917        Variable sam_atten_err,emp_atten_err
918        fileStr = sam_text[3]
919        lambda = sam_reals[26]
920        attenNo = sam_reals[3]
921        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err)
922        fileStr = emp_text[3]
923        lambda = emp_reals[26]
924        attenNo = emp_reals[3]
925        emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,emp_atten_err)
926       
927        //get relative monitor counts (should all be 10^8, since normalized in add step)
928        tmonsam = sam_reals[0]          //monitor count in SAM
929        tsam = sam_reals[4]             //SAM transmission
930        csam = sam_reals[16]            //x center
931        rsam = sam_reals[17]            //beam (x,y) define center of corrected field
932        tmonemp = emp_reals[0]          //monitor count in EMP
933        temp = emp_reals[4]                     //trans emp
934        cemp = emp_reals[16]            //beamcenter of EMP
935        remp = emp_reals[17]
936        savmon_sam=sam_reals[1]         //true monitor count in SAM
937        time_sam = sam_ints[2]          //count time SAM
938        time_drk = drk_ints[2]          //drk count time
939       
940        NVAR pixelsX = root:myGlobals:gNPixelsX
941        NVAR pixelsY = root:myGlobals:gNPixelsY
942        //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM
943        Make/D/O/N=(pixelsX,pixelsY) drk_temp,drk_tmp_err
944        drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam)
945        drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam)                 //temporarily rescale the error of DRK
946
947       
948        if(temp==0)
949                DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction"
950                temp=1
951        Endif
952       
953        //Print "rbgd,cbgd = ",rbgd,cbgd
954        // set up beamcenter shift, relative to SAM
955        xshift = cemp-csam
956        yshift = remp-rsam
957        if(abs(xshift) <= wcen)
958                xshift = 0
959        Endif
960        if(abs(yshift) <= wcen)
961                yshift = 0
962        Endif
963        //get shifted data arrays, relative to SAM
964        Make/D/O/N=(pixelsX,pixelsY) cor1,emp_temp,noadd_emp            //temp arrays
965        GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift)              //emp_temp is the EMP
966        //always ignore beamcenter shift for DRK
967       
968        //do the sam-bgd subtraction,  deposit result in cor1
969        fsam = 1
970        femp = tmonsam/tmonemp          //this should be ==1 since normalized files
971       
972        cor1 = fsam*sam_data/sam_AttenFactor - femp*(tsam/temp)*emp_temp/emp_AttenFactor
973        cor1 += drk_temp - drk_temp/sam_attenFactor
974        cor1 *= noadd_emp               //zeros out regions where arrays do not overlap, one otherwise
975
976// do the error propagation piecewise   
977        Duplicate/O sam_err, tmp_a, tmp_c, c_val
978        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
979       
980        tmp_c = (sam_trans_err*emp_data/(temp*emp_attenFactor))^2 + (emp_err*tsam/(temp*emp_attenFactor))^2
981        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
982       
983        cor_err = sqrt(tmp_a + tmp_c + drk_tmp_err^2)
984       
985        //we're done, get out w/no error
986        //set the COR data to the result
987        cor_data = cor1
988        cor_data_display = cor1
989
990        //update COR header
991        cor_text[1] = date() + " " + time()             //date + time stamp
992
993        KillWaves/Z cor1,emp_temp,noadd_emp,drk_temp
994        Killwaves/Z tmp_a,tmp_c,c_val,drk_tmp_err
995
996        SetDataFolder root:
997        Return(0)
998End
999
1000// ONLY drk subtraction
1001//
1002Function CorrectMode_14()
1003        //create the necessary wave references
1004        WAVE sam_data=$"root:Packages:NIST:SAM:linear_data"
1005        WAVE sam_reals=$"root:Packages:NIST:SAM:realsread"
1006        WAVE sam_ints=$"root:Packages:NIST:SAM:integersread"
1007        WAVE/T sam_text=$"root:Packages:NIST:SAM:textread"
1008        WAVE drk_data=$"root:Packages:NIST:DRK:linear_data"
1009        WAVE drk_reals=$"root:Packages:NIST:DRK:realsread"
1010        WAVE drk_ints=$"root:Packages:NIST:DRK:integersread"
1011        WAVE/T drk_text=$"root:Packages:NIST:DRK:textread"
1012        WAVE cor_data=$"root:Packages:NIST:COR:linear_data"
1013        WAVE/T cor_text=$"root:Packages:NIST:COR:textread"
1014
1015        // needed to propagate error
1016        WAVE cor_data_display=$"root:Packages:NIST:COR:data"            //just for the final copy
1017        WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error"
1018        WAVE drk_err =$"root:Packages:NIST:DRK:linear_data_error"
1019        WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error"
1020       
1021        Variable sam_trans_err
1022        sam_trans_err = sam_reals[41]
1023       
1024       
1025        //get sam and bgd attenuation factors
1026        String fileStr=""
1027        Variable lambda,attenNo,sam_AttenFactor,bgd_attenFactor
1028        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
1029        Variable wcen=0.001,time_drk,time_sam,savmon_sam,tsam
1030        Variable sam_atten_err
1031        fileStr = sam_text[3]
1032        lambda = sam_reals[26]
1033        attenNo = sam_reals[3]
1034        sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err)
1035       
1036        //get relative monitor counts (should all be 10^8, since normalized in add step)
1037        tmonsam = sam_reals[0]          //monitor count in SAM
1038        tsam = sam_reals[4]             //SAM transmission
1039        csam = sam_reals[16]            //x center
1040        rsam = sam_reals[17]            //beam (x,y) define center of corrected field
1041
1042        savmon_sam=sam_reals[1]         //true monitor count in SAM
1043        time_sam = sam_ints[2]          //count time SAM
1044        time_drk = drk_ints[2]          //drk count time
1045       
1046        NVAR pixelsX = root:myGlobals:gNPixelsX
1047        NVAR pixelsY = root:myGlobals:gNPixelsY
1048        //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM
1049        Make/D/O/N=(pixelsX,pixelsY) drk_temp,drk_tmp_err
1050        drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam)
1051        drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam)                 //temporarily rescale the error of DRK
1052
1053        Make/D/O/N=(pixelsX,pixelsY) cor1       //temp arrays
1054        //always ignore the DRK center shift
1055       
1056        //do the subtraction,  deposit result in cor1
1057        fsam = 1
1058        fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
1059       
1060        //correct sam for attenuators, and do the same to drk, since it was scaled to sam count time
1061        cor1 = fsam*sam_data/sam_AttenFactor  - drk_temp/sam_attenFactor
1062
1063// do the error propagation piecewise   
1064        Duplicate/O sam_err, tmp_a
1065        tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
1066
1067        cor_err = sqrt(tmp_a + drk_tmp_err^2)
1068       
1069        //we're done, get out w/no error
1070        //set the COR_data to the result
1071        cor_data = cor1
1072        cor_data_display = cor1
1073
1074        //update COR header
1075        cor_text[1] = date() + " " + time()             //date + time stamp
1076
1077        KillWaves/Z cor1,bgd_temp,noadd_bgd,drk_temp
1078        Killwaves/Z tmp_a,tmp_b,tmp_c,tmp_d,c_val,d_val,drk_tmp_err
1079
1080        SetDataFolder root:
1081        Return(0)
1082End
1083
1084
1085//function to return the shifted contents of a data array for subtraction
1086//(SLOW) if ShiftSum is called
1087//data_in is input
1088//data_out is shifted matrix
1089//noadd_mat =1 if shift matrix is valid, =0 if no data
1090//
1091//if no shift is required, data_in is returned and noadd_mat =1 (all valid)
1092//
1093Function GetShiftedArray(data_in,data_out,noadd_mat,xshift,yshift)
1094        WAVE data_in,data_out,noadd_mat
1095        Variable xshift,yshift
1096
1097        Variable ii=0,jj=0
1098        noadd_mat = 1           //initialize to 1
1099       
1100        If((xshift != 0) || (yshift != 0))
1101//      If((abs(xshift) >= 0.01) || (abs(yshift) >= 0.01))                      //APR09 - loosen tolerance to handle ICE "precision"
1102                DoAlert 1,"Do you want to ignore the beam center mismatch?"
1103                if(V_flag==1)           //yes -> just go on
1104                        xshift=0
1105                        yshift=0
1106                endif
1107        else
1108                // "mismatch" is simply a python type conversion error
1109                xshift=0
1110                yshift=0
1111        endif
1112       
1113        If((xshift == 0) && (yshift == 0))
1114                data_out=data_in                //no change
1115                noadd_mat = 1                   //use all of the data
1116                return(0)
1117        endif
1118       
1119        NVAR pixelsX = root:myGlobals:gNPixelsX
1120        NVAR pixelsY = root:myGlobals:gNPixelsY
1121       
1122        Print "beamcenter shift x,y = ",xshift,yshift
1123        Make/O/N=1 noadd
1124        for(ii=0;ii<pixelsX;ii+=1)
1125                for(jj=0;jj<pixelsY;jj+=1)
1126                        //get the contribution of the shifted data
1127                        data_out[ii][jj] = ShiftSum(data_in,ii,jj,xshift,yshift,noadd)
1128                        if(noadd[0])
1129                                noadd_mat[ii][jj] = 0   //shift is off the detector
1130                        endif
1131                endfor
1132        endfor
1133        return(0)
1134End
1135
1136//utility function that checks if data exists in a data folder
1137//checks only for the existence of DATA - no other waves
1138//
1139Function WorkDataExists(type)
1140        String type
1141       
1142        String destPath=""
1143        destPath =  "root:Packages:NIST:"+Type + ":data"
1144        if(WaveExists($destpath) == 0)
1145                Print "There is no work file in "+type
1146                Return(1)               //error condition
1147        else
1148                //check for log-scaling of the data and adjust if necessary
1149                ConvertFolderToLinearScale(type)
1150                return(0)
1151        Endif
1152End
1153
1154//////////////////
1155// bunch of utility junk to catch
1156// sample transmission = 1
1157// and handle (too many) options
1158//
1159Function GetNewTrans(oldTrans,type)
1160        Variable oldTrans
1161        String type
1162       
1163        Variable newTrans,newCode
1164        if (oldTrans!=1)
1165                return(oldTrans)                //get out now if trans != 1, don't change anything
1166        endif
1167        //get input from the user
1168        NewDataFolder/O root:myGlobals:tmp_trans
1169        Variable/G root:myGlobals:tmp_trans:inputTrans=0.92
1170        Variable/G root:myGlobals:tmp_trans:returnCode=0
1171        DoTransInput(type)
1172        NVAR inputTrans=root:myGlobals:tmp_trans:inputTrans
1173        NVAR code=root:myGlobals:tmp_trans:returnCode
1174        newTrans=inputTrans             //keep a copy before deleting everything
1175        newCode=code
1176        if(newCode==4)
1177                Abort "Aborting correction. Use the Transmission Panel to calculate transmissions"
1178        Endif
1179//      printf "You entered %g and the code is %g\r",newTrans,newCode
1180//      KillDataFolder root:tmp_trans
1181       
1182        if(newCode==1)
1183                Variable/G root:Packages:NIST:gDoTransCheck=0   //turn off checking
1184        endif
1185       
1186        if(newcode==2)          //user changed trans value
1187                return(newTrans)
1188        else
1189                return(oldTrans)        //All other cases, user did not change value
1190        endif
1191end
1192
1193Function IgnoreNowButton(ctrlName) : ButtonControl
1194        String ctrlName
1195       
1196//      Print "ignore now"
1197        NVAR val=root:myGlobals:tmp_trans:returnCode
1198        val=0           //code for ignore once
1199       
1200        DoWindow/K tmp_GetInputPanel            // Kill self
1201End
1202
1203Function DoTransInput(str)
1204        String str
1205       
1206        NewPanel /W=(150,50,361,294)
1207        DoWindow/C tmp_GetInputPanel            // Set to an unlikely name
1208        DrawText 15,23,"The "+str+" Transmission = 1"
1209        DrawText 15,43,"What do you want to do?"
1210        DrawText 15,125,"(Reset this in Preferences)"
1211        SetVariable setvar0,pos={20,170},size={160,17},limits={0,1,0.01}
1212        SetVariable setvar0,value= root:myGlobals:tmp_trans:inputTrans,title="New Transmission"
1213
1214        Button button0,pos={36,56},size={120,20},proc=IgnoreNowButton,title="Ignore This Time"
1215        Button button1,pos={36,86},size={120,20},proc=IgnoreAlwaysButtonProc,title="Ignore Always"
1216        Button button2,pos={36,143},size={120,20},proc=UseNewValueButtonProc,title="Use New Value"
1217        Button button3,pos={36,213},size={120,20},proc=AbortCorrectionButtonProc,title="Abort Correction"
1218        PauseForUser tmp_GetInputPanel
1219End
1220
1221Function IgnoreAlwaysButtonProc(ctrlName) : ButtonControl
1222        String ctrlName
1223
1224//      Print "ignore always"
1225        NVAR val=root:myGlobals:tmp_trans:returnCode
1226        val=1           //code for ignore always
1227        DoWindow/K tmp_GetInputPanel            // Kill self
1228End
1229
1230Function UseNewValueButtonProc(ctrlName) : ButtonControl
1231        String ctrlName
1232
1233//      Print "use new Value"
1234        NVAR val=root:myGlobals:tmp_trans:returnCode
1235        val=2           //code for use new Value
1236        DoWindow/K tmp_GetInputPanel            // Kill self
1237End
1238
1239Function AbortCorrectionButtonProc(ctrlName) : ButtonControl
1240        String ctrlName
1241
1242//      Print "Abort"
1243        NVAR val=root:myGlobals:tmp_trans:returnCode
1244        val=4           //code for abort
1245        DoWindow/K tmp_GetInputPanel            // Kill self
1246End
1247
1248//////////////////////////
1249//**********unused***********
1250//mode describes the type of subtraction that is to be done
1251//1 = both emp and bgd subtraction
1252//2 = only bgd subtraction
1253//3 = only emp subtraction
1254//4 = no subtraction (handled by ExecuteProtocol(), but implemented here as well)
1255//
1256// + 10 indicates that WORK.DRK is to be used
1257//**********unused***********
1258//**worse yet, only partially converted to use DRK files!***********
1259//
1260Function OLD_Correct(mode)
1261        Variable mode
1262       
1263        //Print "mode = ",mode
1264        if(mode==4)
1265                Print "no subtraction required - Correct(mode) should not have been called"
1266                return(1)               //error - correct should not have been called
1267        Endif
1268       
1269        // always check for existence of data in SAM
1270        // if the desired workfile doesn't exist, let the user know, and abort
1271        String destPath
1272        String type = "SAM"
1273        //check for SAM
1274        destPath = "root:Packages:NIST:"+Type + ":data"
1275        if(WaveExists($destpath) == 0)
1276                Print "There is no work file in "+type+"--Aborting"
1277                Return(1)               //error condition
1278        else
1279                //check for log-scaling of the "SAM" data and adjust if necessary
1280                ConvertFolderToLinearScale(type)
1281                Wave sam_data = $"root:Packages:NIST:SAM:data"
1282        Endif
1283       
1284        //check for BGD if mode = 1 or 2 or 11 or 12
1285        if( (mode ==1) || (mode==2) || (mode==11) || (mode==12) )
1286                type = "BGD"
1287                destPath =  "root:Packages:NIST:"+Type + ":data"
1288                if(WaveExists($destpath) == 0)
1289                        Print "There is no work file in "+type+"--Aborting"
1290                        Return(1)               //error condition
1291                else
1292                        //check for log-scaling of the "BGD" data and adjust if necessary
1293                        ConvertFolderToLinearScale(type)
1294                        Wave bgd_data = $"root:Packages:NIST:BGD:data"
1295                Endif
1296        Endif
1297       
1298        // check for EMP data if type 3 or 1 or 13 or 11
1299        if( (mode==1) || (mode==3) || (mode==11) || (mode==13) )
1300                type = "EMP"
1301                destPath =  "root:Packages:NIST:"+Type + ":data"
1302                if(WaveExists($destpath) == 0)
1303                        Print "There is no work file in "+type+"--Aborting"
1304                        Return(1)               //error condition
1305                else
1306                        //check for log-scaling of the "EMP" data and adjust if necessary
1307                        ConvertFolderToLinearScale(type)
1308                        Wave emp_data = $"root:Packages:NIST:EMP:data"
1309                Endif
1310        Endif
1311       
1312        // check for DRK data if type 11,12,13, or 14
1313        if( (mode==11) || (mode==12) || (mode==13) || (mode==14) )
1314                type = "DRK"
1315                destPath =  "root:Packages:NIST:"+Type + ":data"
1316                if(WaveExists($destpath) == 0)
1317                        Print "There is no work file in "+type+"--Aborting"
1318                        Return(1)               //error condition
1319                else
1320                        //check for log-scaling of the "EMP" data and adjust if necessary
1321                        ConvertFolderToLinearScale(type)
1322                        Wave drk_data = $"root:DRK:data"
1323                Endif
1324        Endif
1325       
1326        //necessary files exist, proceed
1327
1328        //make needed wave references to other folders
1329        //NOTE that these references MAY NOT EXIST, depending on the mode
1330        WAVE sam_reals = $"root:Packages:NIST:SAM:realsread"
1331        WAVE sam_ints = $"root:Packages:NIST:SAM:integersread"
1332        WAVE/T sam_text = $"root:Packages:NIST:SAM:textread"
1333        WAVE/Z emp_reals = $"root:Packages:NIST:EMP:realsread"
1334        WAVE/Z emp_ints = $"root:Packages:NIST:EMP:integersread"
1335        WAVE/T/Z emp_text = $"root:Packages:NIST:EMP:textread"
1336        WAVE/Z bgd_reals = $"root:Packages:NIST:BGD:realsread"
1337        WAVE/Z bgd_ints = $"root:Packages:NIST:BGD:integersread"
1338        WAVE/T/Z bgd_text = $"root:Packages:NIST:BGD:textread"
1339       
1340        //find the attenuation of the sample (if any)
1341        Variable SamAttenFactor,lambda,attenNo,err=0
1342        String samfileStr=""
1343        Variable sam_atten_err
1344        samfileStr = sam_text[3]
1345        lambda = sam_reals[26]
1346        attenNo = sam_reals[3]
1347        SamAttenFactor = AttenuationFactor(samFileStr,lambda,AttenNo,sam_atten_err)
1348        //if sample trans is zero, do only SAM-BGD subtraction (notify the user)
1349        Variable sam_trans = sam_reals[4]
1350       
1351        //copy SAM information to COR, wiping out the old contents of the COR folder first
1352        err = CopyWorkContents("SAM","COR")     
1353        if(err==1)
1354                Abort "No data in SAM, abort from Correct()"
1355        endif
1356       
1357        //now switch to COR folder
1358        DestPath="root:Packages:NIST:COR"
1359        //make appropriate wave references
1360        WAVE data=$(destPath + ":data")                                 // these wave references point to the SAM data in COR
1361        WAVE/T textread=$(destPath + ":textread")                       //that are to be directly operated on
1362        WAVE integersread=$(destPath + ":integersread")
1363        WAVE realsread=$(destPath + ":realsRead")
1364
1365        NVAR pixelsX = root:myGlobals:gNPixelsX
1366        NVAR pixelsY = root:myGlobals:gNPixelsY
1367       
1368        //Print "done copying data, starting the correct calculations"
1369       
1370        // Start the actual "correct" step here
1371        Variable wcen=0.001,numsam,tmonsam,tsam,rsam,csam,fsam
1372        Variable tmonbgd,fbgd,xshift,yshift,rbgd,cbgd,sh_sum,ii,jj,trans,tmonemp,temp,femp
1373        Variable cemp,remp
1374        //make temporary waves to hold the intermediate results and the shifted arrays
1375        Duplicate/O data cor1,cor2
1376        cor1 = 0                //initialize to zero
1377        cor2 = 0
1378       
1379        //make needed wave references to other folders
1380        Wave sam_reals = $"root:Packages:NIST:SAM:realsread"
1381        Wave bgd_reals = $"root:Packages:NIST:BGD:realsread"
1382        Wave emp_reals = $"root:Packages:NIST:EMP:realsread"
1383       
1384        //get counts, trans, etc. from file headers
1385        numsam = sam_ints[3]            //number of runs in SAM file
1386        tmonsam = sam_reals[0]          //monitor count in SAM
1387        tsam = sam_reals[4]             //SAM transmission
1388        csam = sam_reals[16]            //x center
1389        rsam = sam_reals[17]            //beam (x,y) define center of corrected field
1390        //Print "rsam,csam = ",rsam,csam
1391       
1392        //
1393        //do sam-bgd subtraction if mode (1) or (2)
1394        //else (mode 3), set cor1 = sam_data
1395        if( (mode==1) || (mode==2) )
1396                fsam = 1
1397                tmonbgd = bgd_reals[0]          //monitor count in BGD
1398                fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
1399       
1400                //set up  center shift, relative to SAM
1401                cbgd = bgd_reals[16]
1402                rbgd = bgd_reals[17]
1403                //Print "rbgd,cbgd = ",rbgd,cbgd
1404                xshift = cbgd-csam
1405                yshift = rbgd-rsam
1406                if(abs(xshift) <= wcen)
1407                        xshift = 0
1408                Endif
1409                if(abs(yshift) <= wcen)
1410                        yshift = 0
1411                Endif
1412               
1413                If((xshift != 0) || (yshift != 0))
1414                        DoAlert 1,"Do you want to ignore the beam center mismatch?"
1415                        if(V_flag==1)           //yes -> just go on
1416                                xshift=0
1417                                yshift=0
1418                        endif
1419                endif
1420                //do the sam-bgd subtraction,  deposit result in cor1[][]
1421                If((xshift == 0) && (yshift == 0))
1422                        //great, no shift required
1423                        cor1 = fsam*sam_data - fbgd*bgd_data*SamAttenFactor
1424                else
1425                        //shift required, very time-consuming
1426                        Print "sam-bgd shift x,y = ",xshift,yshift
1427                        Make/O/N=1 noadd                //needed to get noadd condition back from ShiftSum()
1428                        ii=0
1429                        do
1430                                jj=0
1431                                do
1432                                        //get the contribution of shifted data
1433                                        sh_sum = ShiftSum(bgd_data,ii,jj,xshift,yshift,noadd)
1434                                        if(noadd[0])
1435                                                cor1[ii][jj]=0
1436                                        else
1437                                                //add the sam_data + shifted sum
1438                                                cor1[ii][jj] = fsam*sam_data[ii][jj] - fbgd*sh_sum*SamAttenFactor
1439                                        Endif
1440                                        jj+=1
1441                                while(jj<pixelsY)
1442                                ii+=1
1443                        while(ii<pixelsX)
1444                Endif
1445        else                    //switch on mode
1446                cor1 = sam_data         //setup for just EMP subtraction
1447        Endif
1448       
1449        //Print "sam-bgd done"
1450       
1451        if(mode == 2)           //just a BGD subtraction
1452                //we're done, get out w/no error
1453                //set the COR data to the result
1454                data = cor1
1455                //update COR header
1456                textread[1] = date() + " " + time()             //date + time stamp
1457                SetDataFolder root:
1458                KillWaves/Z cor1,cor2
1459                Return(0)
1460        Endif
1461       
1462        //if mode ==1 (ONLY choice left) do the empty-background subtraction
1463        //else mode = 3, set cor2 to emp_data
1464        if(mode==1)             //full subtraction
1465                trans = emp_reals[4]            //EMP transmission
1466                if(trans == 0)
1467                        trans = 1
1468                        DoAlert 0,"Empty cell transmission was zero. It has been reset to one for the calculation"
1469                endif
1470                tmonemp = emp_reals[0]
1471                femp = tmonsam/tmonemp
1472                temp = trans
1473       
1474                //set up center shift, relative to EMP
1475                cemp = emp_reals[16]
1476                remp = emp_reals[17]
1477                //Print "remp,cemp ",remp,cemp
1478                xshift = cbgd - cemp
1479                yshift = rbgd - remp
1480                if(abs(xshift) <= wcen )
1481                        xshift = 0
1482                endif
1483                if(abs(yshift) <= wcen)
1484                        yshift = 0
1485                endif
1486               
1487                If((xshift != 0) || (yshift != 0))
1488                        DoAlert 1,"Do you want to ignore the beam center mismatch?"
1489                        if(V_flag==1)
1490                                xshift=0
1491                                yshift=0
1492                        endif
1493                endif
1494                //do the emp-bgd subtraction,  deposit result in cor2[][]
1495                If((xshift == 0) && (yshift == 0))
1496                        //great, no shift required, DON'T scale this by the attenuation gfactor
1497                        cor2 = femp*emp_data - fbgd*bgd_data
1498                else
1499                        //shift required, very time-consuming
1500                        Print "emp-bgd shift x,y = ",xshift,yshift
1501                        Make/O/N=1 noadd                //needed to get noadd condition back from ShiftSum()
1502                        ii=0
1503                        do
1504                                jj=0
1505                                do
1506                                        //get the contribution of shifted data
1507                                        sh_sum = ShiftSum(bgd_data,ii,jj,xshift,yshift,noadd)
1508                                        if(noadd[0])
1509                                                cor2[ii][jj]=0
1510                                        else
1511                                                //add the sam_data + shifted sum
1512                                                cor2[ii][jj] = femp*emp_data[ii][jj] - fbgd*sh_sum
1513                                        Endif
1514                                        jj+=1
1515                                while(jj<pixelsY)
1516                                ii+=1
1517                        while(ii<pixelsX)
1518                Endif
1519        else            //switch on mode==1 for full subtraction
1520                cor2 = emp_data
1521                //be sure to set the empty center location... for the shift
1522                trans = emp_reals[4]            //EMP transmission
1523                if(trans == 0)
1524                        trans = 1
1525                        DoAlert 0,"Empty cell transmission was zero. It has been reset to one for the calculation"
1526                endif
1527                tmonemp = emp_reals[0]
1528                femp = tmonsam/tmonemp
1529                temp = trans
1530       
1531                //set up center shift, relative to EMP
1532                cemp = emp_reals[16]
1533                remp = emp_reals[17]
1534        Endif
1535       
1536        //Print "emp-bgd done"
1537       
1538        //mode 2 exited, either 1 or 3 apply from here, and are setup properly.
1539       
1540        //set up for final step, data(COR) = cor1 - Tsam/Temp*cor2
1541        //set up shift, relative to SAM
1542        xshift = cemp - csam
1543        yshift = remp - rsam
1544        if(abs(xshift) <= wcen )
1545                xshift = 0
1546        endif
1547        if(abs(yshift) <= wcen)
1548                yshift = 0
1549        endif
1550       
1551        If((xshift != 0) || (yshift != 0))
1552                DoAlert 1,"Do you want to ignore the beam center mismatch?"
1553                if(V_flag==1)
1554                        xshift=0
1555                        yshift=0
1556                endif
1557        endif
1558        //do the cor1-a*cor2 subtraction,  deposit result in data[][] (in the COR folder)
1559        If((xshift == 0) && (yshift == 0))
1560                //great, no shift required
1561                data = cor1 - (tsam/temp)*cor2*SamAttenFactor
1562        else
1563                //shift required, very time-consuming
1564                Print "sam-emp shift x,y = ",xshift,yshift
1565                Make/O/N=1 noadd                //needed to get noadd condition back from ShiftSum()
1566                ii=0
1567                do
1568                        jj=0
1569                        do
1570                                //get the contribution of shifted data
1571                                sh_sum = ShiftSum(cor2,ii,jj,xshift,yshift,noadd)
1572                                if(noadd[0])
1573                                        data[ii][jj]=0
1574                                else
1575                                        //add the sam_data + shifted sum
1576                                        data[ii][jj] = cor1[ii][jj] - (tsam/temp)*sh_sum*SamAttenFactor
1577                                Endif
1578                                jj+=1
1579                        while(jj<pixelsY)
1580                        ii+=1
1581                while(ii<pixelsX)
1582        Endif
1583       
1584        //Print "all done"
1585       
1586        //update COR header
1587        textread[1] = date() + " " + time()             //date + time stamp
1588       
1589        //clean up
1590        SetDataFolder root:Packages:NIST:COR
1591        SetDataFolder root:
1592        KillWaves/Z cor1,cor2,noadd
1593
1594        Return(0)               //all is ok, if you made it to this point
1595End
Note: See TracBrowser for help on using the repository browser.