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

Last change on this file since 260 was 260, checked in by srkline, 15 years ago

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