source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_Correct.ipf @ 1021

Last change on this file since 1021 was 1021, checked in by srkline, 6 years ago

Adjusting testing macros and panels to be listed on the VSANS menu rather than the Macros menu.

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