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

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

bug fixes and performance tweaks

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