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

Last change on this file since 502 was 502, checked in by srkline, 14 years ago

List of small changes:
Lin_Fits - added I(0) error to one calculation, Need to propagate this to all logical places. (and submit a ticket)
MainPanel? - correctly handle cancel of raw data display
Transmssion - add a "guess" option of a user-supplied "N" characters
Correct - revert back to zero error for beam center mismatch
NCNR_DataRW - add utility to patch monitor counts (comment out of menu later)

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