source: sans/SANSReduction/trunk/Put in User Procedures/SANS_Reduction_v5.00/Correct.ipf @ 167

Last change on this file since 167 was 167, checked in by srkline, 16 years ago

Changed NSORT to display rescaling and overlap constants behind the OS file save dialog, allowing the user to cancel the save it something is amiss.

No significant change in trans=1 panel behavior in Correct.ipf

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