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

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

Made preferences a common panel (moved to PlotUtilsMacro?.ipf and globals to root:Packages:NIST:) and added menu items for all packages. Many files had to be modified so that the preferences could be properly accessed

File Open dialog now is set to "All files" so that XML can be selected. I think that all open access that doesn't already have the full path go through this common function.

File size: 42.2 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=6.1
4
5//****************************
6// Vers 1.2 090501
7//
8//****************************
9//
10// Procedures to perform the "Correct" step during data reduction
11//
12// - ther is olny one procedure to perform the subtractions, and a single
13// parameter flags which subtractions are to be done. Different numbers of
14// attenuators during scattering runs are corrected as described in John's memo,
15// with the note that ONLY method (3) is used, which assumes that 'diffuse' scattering
16// is dominant over 'dark current' (note that 'dark current' = shutter CLOSED)
17//
18//
19//do the CORRECT step based on the answers to emp and bkg subtraction
20        //by setting the proper"mode"
21        //1 = both emp and bgd subtraction
22        //2 = only bgd subtraction
23        //3 = only emp subtraction
24        //4 = no subtraction
25        //additional modes 091301
26        //11 = emp, bgd, drk
27        //12 = bgd and drk
28        //13 = emp and drk
29        //14 = no subtractions
30        //
31//********************************
32
33//unused test procedure for Correct() function
34//must be updated to include "mode" parameter before re-use
35//
36Proc CorrectData()
37
38        Variable err
39        String type
40       
41        err = Correct()         
42       
43        if(err)
44                Abort "error in Correct"
45        endif
46       
47        //contents are always dumped to COR
48        type = "COR"
49       
50        //need to update the display with "data" from the correct dataFolder
51        //reset the current displaytype to "type"
52        String/G root:myGlobals:gDataDisplayType=Type
53       
54        fRawWindowHook()
55       
56End
57
58
59//mode describes the type of subtraction that is to be done
60//1 = both emp and bgd subtraction
61//2 = only bgd subtraction
62//3 = only emp subtraction
63//4 = no subtraction
64//
65// + 10 indicates that WORK.DRK is to be used
66//
67//091301 version
68//now simple dispatches to the correct subtraction - logic was too
69//involved to do in one function - unclear and error-prone
70//
71// 081203 version
72// checks for trans==1 in SAM and EMP before dispatching
73// and asks for new value if desired
74//
75Function Correct(mode)
76        Variable mode
77       
78        Variable err=0,trans,newTrans
79       
80        //switch and dispatch based on the required subtractions
81        // always check for SAM data
82        err = WorkDataExists("SAM")
83        if(err==1)
84                return(err)
85        endif
86       
87       
88        //check for trans==1
89        NVAR doCheck=root:Packages:NIST:gDoTransCheck
90        Wave/Z samR=root:Packages:NIST:SAM:RealsRead
91        Wave/Z empR=root:Packages:NIST:EMP:RealsRead
92        if(doCheck)
93                trans = samR[4]
94                newTrans=GetNewTrans(trans,"SAM")               //will change value if necessary
95                if(numtype(newTrans)==0)
96                        samR[4] = newTrans              //avoid user abort assigning NaN
97                endif
98                if(trans != newTrans)
99                        print "Using SAM trans = ",samR[4]
100                endif
101        endif
102       
103        //copy SAM information to COR, wiping out the old contents of the COR folder first
104        //do this even if no correction is dispatched (if incorrect mode)
105        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:Packages:NIST: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.