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

Last change on this file since 495 was 495, checked in by srkline, 13 years ago

Added a procedure file that does the necessary wrapping for GenCurveFit?, since it's not a perfect drop-in replacement for FuncFit?. Required modification of the FitWrapper? to switch to GenCurveFit? as needed.
Switch between regular L-M and GenOp? using a menu item.a (sets a gobal)
Still ridiculously slow to use.
Can't yet be used with global fitting. Would be a real pain to implement. Can't imagine how slow that would be to use...

Bug fixes in PatchFiles? (default button for filter type) and Correct (use tolerance of +/- 0.01 pixel for determining of what is a "mismatch" of the beam centers. Problem cropped up with ICE, but should be fixed anyways)

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((abs(xshift) >= 0.01) || (abs(yshift) >= 0.01))                      //APR09 - change tolerance to handle ICE "precision"
898                DoAlert 1,"Do you want to ignore the beam center mismatch?"
899                if(V_flag==1)           //yes -> just go on
900                        xshift=0
901                        yshift=0
902                endif
903        else
904                // "mismatch" is simply a python type conversion error
905                xshift=0
906                yshift=0
907        endif
908       
909        If((xshift == 0) && (yshift == 0))
910                data_out=data_in                //no change
911                noadd_mat = 1                   //use all of the data
912                return(0)
913        endif
914       
915        NVAR pixelsX = root:myGlobals:gNPixelsX
916        NVAR pixelsY = root:myGlobals:gNPixelsY
917       
918        Print "beamcenter shift x,y = ",xshift,yshift
919        Make/O/N=1 noadd
920        for(ii=0;ii<pixelsX;ii+=1)
921                for(jj=0;jj<pixelsY;jj+=1)
922                        //get the contribution of the shifted data
923                        data_out[ii][jj] = ShiftSum(data_in,ii,jj,xshift,yshift,noadd)
924                        if(noadd[0])
925                                noadd_mat[ii][jj] = 0   //shift is off the detector
926                        endif
927                endfor
928        endfor
929        return(0)
930End
931
932//utility function that checks if data exists in a data folder
933//checks only for the existence of DATA - no other waves
934//
935Function WorkDataExists(type)
936        String type
937       
938        String destPath=""
939        destPath =  "root:Packages:NIST:"+Type + ":data"
940        if(WaveExists($destpath) == 0)
941                Print "There is no work file in "+type
942                Return(1)               //error condition
943        else
944                //check for log-scaling of the data and adjust if necessary
945                ConvertFolderToLinearScale(type)
946                return(0)
947        Endif
948End
949
950//////////////////
951// bunch of utility junk to catch
952// sample transmission = 1
953// and handle (too many) options
954//
955Function GetNewTrans(oldTrans,type)
956        Variable oldTrans
957        String type
958       
959        Variable newTrans,newCode
960        if (oldTrans!=1)
961                return(oldTrans)                //get out now if trans != 1, don't change anything
962        endif
963        //get input from the user
964        NewDataFolder/O root:myGlobals:tmp_trans
965        Variable/G root:myGlobals:tmp_trans:inputTrans=0.92
966        Variable/G root:myGlobals:tmp_trans:returnCode=0
967        DoTransInput(type)
968        NVAR inputTrans=root:myGlobals:tmp_trans:inputTrans
969        NVAR code=root:myGlobals:tmp_trans:returnCode
970        newTrans=inputTrans             //keep a copy before deleting everything
971        newCode=code
972        if(newCode==4)
973                Abort "Aborting correction. Use the Transmission Panel to calculate transmissions"
974        Endif
975//      printf "You entered %g and the code is %g\r",newTrans,newCode
976//      KillDataFolder root:tmp_trans
977       
978        if(newCode==1)
979                Variable/G root:myGlobals:gDoTransCheck=0       //turn off checking
980        endif
981       
982        if(newcode==2)          //user changed trans value
983                return(newTrans)
984        else
985                return(oldTrans)        //All other cases, user did not change value
986        endif
987end
988
989Function IgnoreNowButton(ctrlName) : ButtonControl
990        String ctrlName
991       
992//      Print "ignore now"
993        NVAR val=root:myGlobals:tmp_trans:returnCode
994        val=0           //code for ignore once
995       
996        DoWindow/K tmp_GetInputPanel            // Kill self
997End
998
999Function DoTransInput(str)
1000        String str
1001       
1002        NewPanel /W=(150,50,361,294)
1003        DoWindow/C tmp_GetInputPanel            // Set to an unlikely name
1004        DrawText 15,23,"The "+str+" Transmission = 1"
1005        DrawText 15,43,"What do you want to do?"
1006        DrawText 15,125,"(Reset this in Preferences)"
1007        SetVariable setvar0,pos={20,170},size={160,17},limits={0,1,0.01}
1008        SetVariable setvar0,value= root:myGlobals:tmp_trans:inputTrans,title="New Transmission"
1009
1010        Button button0,pos={36,56},size={120,20},proc=IgnoreNowButton,title="Ignore This Time"
1011        Button button1,pos={36,86},size={120,20},proc=IgnoreAlwaysButtonProc,title="Ignore Always"
1012        Button button2,pos={36,143},size={120,20},proc=UseNewValueButtonProc,title="Use New Value"
1013        Button button3,pos={36,213},size={120,20},proc=AbortCorrectionButtonProc,title="Abort Correction"
1014        PauseForUser tmp_GetInputPanel
1015End
1016
1017Function IgnoreAlwaysButtonProc(ctrlName) : ButtonControl
1018        String ctrlName
1019
1020//      Print "ignore always"
1021        NVAR val=root:myGlobals:tmp_trans:returnCode
1022        val=1           //code for ignore always
1023        DoWindow/K tmp_GetInputPanel            // Kill self
1024End
1025
1026Function UseNewValueButtonProc(ctrlName) : ButtonControl
1027        String ctrlName
1028
1029//      Print "use new Value"
1030        NVAR val=root:myGlobals:tmp_trans:returnCode
1031        val=2           //code for use new Value
1032        DoWindow/K tmp_GetInputPanel            // Kill self
1033End
1034
1035Function AbortCorrectionButtonProc(ctrlName) : ButtonControl
1036        String ctrlName
1037
1038//      Print "Abort"
1039        NVAR val=root:myGlobals:tmp_trans:returnCode
1040        val=4           //code for abort
1041        DoWindow/K tmp_GetInputPanel            // Kill self
1042End
1043
1044//////////////////////////
1045//**********unused***********
1046//mode describes the type of subtraction that is to be done
1047//1 = both emp and bgd subtraction
1048//2 = only bgd subtraction
1049//3 = only emp subtraction
1050//4 = no subtraction (handled by ExecuteProtocol(), but implemented here as well)
1051//
1052// + 10 indicates that WORK.DRK is to be used
1053//**********unused***********
1054//**worse yet, only partially converted to use DRK files!***********
1055//
1056Function OLD_Correct(mode)
1057        Variable mode
1058       
1059        //Print "mode = ",mode
1060        if(mode==4)
1061                Print "no subtraction required - Correct(mode) should not have been called"
1062                return(1)               //error - correct should not have been called
1063        Endif
1064       
1065        // always check for existence of data in SAM
1066        // if the desired workfile doesn't exist, let the user know, and abort
1067        String destPath
1068        String type = "SAM"
1069        //check for SAM
1070        destPath = "root:Packages:NIST:"+Type + ":data"
1071        if(WaveExists($destpath) == 0)
1072                Print "There is no work file in "+type+"--Aborting"
1073                Return(1)               //error condition
1074        else
1075                //check for log-scaling of the "SAM" data and adjust if necessary
1076                ConvertFolderToLinearScale(type)
1077                Wave sam_data = $"root:Packages:NIST:SAM:data"
1078        Endif
1079       
1080        //check for BGD if mode = 1 or 2 or 11 or 12
1081        if( (mode ==1) || (mode==2) || (mode==11) || (mode==12) )
1082                type = "BGD"
1083                destPath =  "root:Packages:NIST:"+Type + ":data"
1084                if(WaveExists($destpath) == 0)
1085                        Print "There is no work file in "+type+"--Aborting"
1086                        Return(1)               //error condition
1087                else
1088                        //check for log-scaling of the "BGD" data and adjust if necessary
1089                        ConvertFolderToLinearScale(type)
1090                        Wave bgd_data = $"root:Packages:NIST:BGD:data"
1091                Endif
1092        Endif
1093       
1094        // check for EMP data if type 3 or 1 or 13 or 11
1095        if( (mode==1) || (mode==3) || (mode==11) || (mode==13) )
1096                type = "EMP"
1097                destPath =  "root:Packages:NIST:"+Type + ":data"
1098                if(WaveExists($destpath) == 0)
1099                        Print "There is no work file in "+type+"--Aborting"
1100                        Return(1)               //error condition
1101                else
1102                        //check for log-scaling of the "EMP" data and adjust if necessary
1103                        ConvertFolderToLinearScale(type)
1104                        Wave emp_data = $"root:Packages:NIST:EMP:data"
1105                Endif
1106        Endif
1107       
1108        // check for DRK data if type 11,12,13, or 14
1109        if( (mode==11) || (mode==12) || (mode==13) || (mode==14) )
1110                type = "DRK"
1111                destPath =  "root:Packages:NIST:"+Type + ":data"
1112                if(WaveExists($destpath) == 0)
1113                        Print "There is no work file in "+type+"--Aborting"
1114                        Return(1)               //error condition
1115                else
1116                        //check for log-scaling of the "EMP" data and adjust if necessary
1117                        ConvertFolderToLinearScale(type)
1118                        Wave drk_data = $"root:DRK:data"
1119                Endif
1120        Endif
1121       
1122        //necessary files exist, proceed
1123
1124        //make needed wave references to other folders
1125        //NOTE that these references MAY NOT EXIST, depending on the mode
1126        WAVE sam_reals = $"root:Packages:NIST:SAM:realsread"
1127        WAVE sam_ints = $"root:Packages:NIST:SAM:integersread"
1128        WAVE/T sam_text = $"root:Packages:NIST:SAM:textread"
1129        WAVE/Z emp_reals = $"root:Packages:NIST:EMP:realsread"
1130        WAVE/Z emp_ints = $"root:Packages:NIST:EMP:integersread"
1131        WAVE/T/Z emp_text = $"root:Packages:NIST:EMP:textread"
1132        WAVE/Z bgd_reals = $"root:Packages:NIST:BGD:realsread"
1133        WAVE/Z bgd_ints = $"root:Packages:NIST:BGD:integersread"
1134        WAVE/T/Z bgd_text = $"root:Packages:NIST:BGD:textread"
1135       
1136        //find the attenuation of the sample (if any)
1137        Variable SamAttenFactor,lambda,attenNo,err=0
1138        String samfileStr=""
1139        samfileStr = sam_text[3]
1140        lambda = sam_reals[26]
1141        attenNo = sam_reals[3]
1142        SamAttenFactor = AttenuationFactor(samFileStr,lambda,AttenNo)
1143        //if sample trans is zero, do only SAM-BGD subtraction (notify the user)
1144        Variable sam_trans = sam_reals[4]
1145       
1146        //copy SAM information to COR, wiping out the old contents of the COR folder first
1147        err = CopyWorkContents("SAM","COR")     
1148        if(err==1)
1149                Abort "No data in SAM, abort from Correct()"
1150        endif
1151       
1152        //now switch to COR folder
1153        DestPath="root:Packages:NIST:COR"
1154        //make appropriate wave references
1155        WAVE data=$(destPath + ":data")                                 // these wave references point to the SAM data in COR
1156        WAVE/T textread=$(destPath + ":textread")                       //that are to be directly operated on
1157        WAVE integersread=$(destPath + ":integersread")
1158        WAVE realsread=$(destPath + ":realsRead")
1159
1160        NVAR pixelsX = root:myGlobals:gNPixelsX
1161        NVAR pixelsY = root:myGlobals:gNPixelsY
1162       
1163        //Print "done copying data, starting the correct calculations"
1164       
1165        // Start the actual "correct" step here
1166        Variable wcen=0.001,numsam,tmonsam,tsam,rsam,csam,fsam
1167        Variable tmonbgd,fbgd,xshift,yshift,rbgd,cbgd,sh_sum,ii,jj,trans,tmonemp,temp,femp
1168        Variable cemp,remp
1169        //make temporary waves to hold the intermediate results and the shifted arrays
1170        Duplicate/O data cor1,cor2
1171        cor1 = 0                //initialize to zero
1172        cor2 = 0
1173       
1174        //make needed wave references to other folders
1175        Wave sam_reals = $"root:Packages:NIST:SAM:realsread"
1176        Wave bgd_reals = $"root:Packages:NIST:BGD:realsread"
1177        Wave emp_reals = $"root:Packages:NIST:EMP:realsread"
1178       
1179        //get counts, trans, etc. from file headers
1180        numsam = sam_ints[3]            //number of runs in SAM file
1181        tmonsam = sam_reals[0]          //monitor count in SAM
1182        tsam = sam_reals[4]             //SAM transmission
1183        csam = sam_reals[16]            //x center
1184        rsam = sam_reals[17]            //beam (x,y) define center of corrected field
1185        //Print "rsam,csam = ",rsam,csam
1186       
1187        //
1188        //do sam-bgd subtraction if mode (1) or (2)
1189        //else (mode 3), set cor1 = sam_data
1190        if( (mode==1) || (mode==2) )
1191                fsam = 1
1192                tmonbgd = bgd_reals[0]          //monitor count in BGD
1193                fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
1194       
1195                //set up  center shift, relative to SAM
1196                cbgd = bgd_reals[16]
1197                rbgd = bgd_reals[17]
1198                //Print "rbgd,cbgd = ",rbgd,cbgd
1199                xshift = cbgd-csam
1200                yshift = rbgd-rsam
1201                if(abs(xshift) <= wcen)
1202                        xshift = 0
1203                Endif
1204                if(abs(yshift) <= wcen)
1205                        yshift = 0
1206                Endif
1207               
1208                If((xshift != 0) || (yshift != 0))
1209                        DoAlert 1,"Do you want to ignore the beam center mismatch?"
1210                        if(V_flag==1)           //yes -> just go on
1211                                xshift=0
1212                                yshift=0
1213                        endif
1214                endif
1215                //do the sam-bgd subtraction,  deposit result in cor1[][]
1216                If((xshift == 0) && (yshift == 0))
1217                        //great, no shift required
1218                        cor1 = fsam*sam_data - fbgd*bgd_data*SamAttenFactor
1219                else
1220                        //shift required, very time-consuming
1221                        Print "sam-bgd shift x,y = ",xshift,yshift
1222                        Make/O/N=1 noadd                //needed to get noadd condition back from ShiftSum()
1223                        ii=0
1224                        do
1225                                jj=0
1226                                do
1227                                        //get the contribution of shifted data
1228                                        sh_sum = ShiftSum(bgd_data,ii,jj,xshift,yshift,noadd)
1229                                        if(noadd[0])
1230                                                cor1[ii][jj]=0
1231                                        else
1232                                                //add the sam_data + shifted sum
1233                                                cor1[ii][jj] = fsam*sam_data[ii][jj] - fbgd*sh_sum*SamAttenFactor
1234                                        Endif
1235                                        jj+=1
1236                                while(jj<pixelsY)
1237                                ii+=1
1238                        while(ii<pixelsX)
1239                Endif
1240        else                    //switch on mode
1241                cor1 = sam_data         //setup for just EMP subtraction
1242        Endif
1243       
1244        //Print "sam-bgd done"
1245       
1246        if(mode == 2)           //just a BGD subtraction
1247                //we're done, get out w/no error
1248                //set the COR data to the result
1249                data = cor1
1250                //update COR header
1251                textread[1] = date() + " " + time()             //date + time stamp
1252                SetDataFolder root:
1253                KillWaves/Z cor1,cor2
1254                Return(0)
1255        Endif
1256       
1257        //if mode ==1 (ONLY choice left) do the empty-background subtraction
1258        //else mode = 3, set cor2 to emp_data
1259        if(mode==1)             //full subtraction
1260                trans = emp_reals[4]            //EMP transmission
1261                if(trans == 0)
1262                        trans = 1
1263                        DoAlert 0,"Empty cell transmission was zero. It has been reset to one for the calculation"
1264                endif
1265                tmonemp = emp_reals[0]
1266                femp = tmonsam/tmonemp
1267                temp = trans
1268       
1269                //set up center shift, relative to EMP
1270                cemp = emp_reals[16]
1271                remp = emp_reals[17]
1272                //Print "remp,cemp ",remp,cemp
1273                xshift = cbgd - cemp
1274                yshift = rbgd - remp
1275                if(abs(xshift) <= wcen )
1276                        xshift = 0
1277                endif
1278                if(abs(yshift) <= wcen)
1279                        yshift = 0
1280                endif
1281               
1282                If((xshift != 0) || (yshift != 0))
1283                        DoAlert 1,"Do you want to ignore the beam center mismatch?"
1284                        if(V_flag==1)
1285                                xshift=0
1286                                yshift=0
1287                        endif
1288                endif
1289                //do the emp-bgd subtraction,  deposit result in cor2[][]
1290                If((xshift == 0) && (yshift == 0))
1291                        //great, no shift required, DON'T scale this by the attenuation gfactor
1292                        cor2 = femp*emp_data - fbgd*bgd_data
1293                else
1294                        //shift required, very time-consuming
1295                        Print "emp-bgd shift x,y = ",xshift,yshift
1296                        Make/O/N=1 noadd                //needed to get noadd condition back from ShiftSum()
1297                        ii=0
1298                        do
1299                                jj=0
1300                                do
1301                                        //get the contribution of shifted data
1302                                        sh_sum = ShiftSum(bgd_data,ii,jj,xshift,yshift,noadd)
1303                                        if(noadd[0])
1304                                                cor2[ii][jj]=0
1305                                        else
1306                                                //add the sam_data + shifted sum
1307                                                cor2[ii][jj] = femp*emp_data[ii][jj] - fbgd*sh_sum
1308                                        Endif
1309                                        jj+=1
1310                                while(jj<pixelsY)
1311                                ii+=1
1312                        while(ii<pixelsX)
1313                Endif
1314        else            //switch on mode==1 for full subtraction
1315                cor2 = emp_data
1316                //be sure to set the empty center location... for the shift
1317                trans = emp_reals[4]            //EMP transmission
1318                if(trans == 0)
1319                        trans = 1
1320                        DoAlert 0,"Empty cell transmission was zero. It has been reset to one for the calculation"
1321                endif
1322                tmonemp = emp_reals[0]
1323                femp = tmonsam/tmonemp
1324                temp = trans
1325       
1326                //set up center shift, relative to EMP
1327                cemp = emp_reals[16]
1328                remp = emp_reals[17]
1329        Endif
1330       
1331        //Print "emp-bgd done"
1332       
1333        //mode 2 exited, either 1 or 3 apply from here, and are setup properly.
1334       
1335        //set up for final step, data(COR) = cor1 - Tsam/Temp*cor2
1336        //set up shift, relative to SAM
1337        xshift = cemp - csam
1338        yshift = remp - rsam
1339        if(abs(xshift) <= wcen )
1340                xshift = 0
1341        endif
1342        if(abs(yshift) <= wcen)
1343                yshift = 0
1344        endif
1345       
1346        If((xshift != 0) || (yshift != 0))
1347                DoAlert 1,"Do you want to ignore the beam center mismatch?"
1348                if(V_flag==1)
1349                        xshift=0
1350                        yshift=0
1351                endif
1352        endif
1353        //do the cor1-a*cor2 subtraction,  deposit result in data[][] (in the COR folder)
1354        If((xshift == 0) && (yshift == 0))
1355                //great, no shift required
1356                data = cor1 - (tsam/temp)*cor2*SamAttenFactor
1357        else
1358                //shift required, very time-consuming
1359                Print "sam-emp shift x,y = ",xshift,yshift
1360                Make/O/N=1 noadd                //needed to get noadd condition back from ShiftSum()
1361                ii=0
1362                do
1363                        jj=0
1364                        do
1365                                //get the contribution of shifted data
1366                                sh_sum = ShiftSum(cor2,ii,jj,xshift,yshift,noadd)
1367                                if(noadd[0])
1368                                        data[ii][jj]=0
1369                                else
1370                                        //add the sam_data + shifted sum
1371                                        data[ii][jj] = cor1[ii][jj] - (tsam/temp)*sh_sum*SamAttenFactor
1372                                Endif
1373                                jj+=1
1374                        while(jj<pixelsY)
1375                        ii+=1
1376                while(ii<pixelsX)
1377        Endif
1378       
1379        //Print "all done"
1380       
1381        //update COR header
1382        textread[1] = date() + " " + time()             //date + time stamp
1383       
1384        //clean up
1385        SetDataFolder root:Packages:NIST:COR
1386        SetDataFolder root:
1387        KillWaves/Z cor1,cor2,noadd
1388
1389        Return(0)               //all is ok, if you made it to this point
1390End
Note: See TracBrowser for help on using the repository browser.