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

Last change on this file since 1129 was 1106, checked in by srkline, 4 years ago

Added functions to check that files in a protocol, including the sample data are all from the same configuration. In SANS, the only flag was a beam center mismatch. In VSANS, a more complete check of the configuration in necessary.

Changed all instances of calls to the read and normalize the monitor count to use the "norm"al monitor, not the value listed in the Control block.

File size: 44.2 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        V_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//
269// TODO -- verify the operation
270//  -- check that all "get" reads are returning proper values (+units)
271//
272Function V_CorrectMode_1()
273
274        //get SAM, BGD, EMP attenuation factor
275        Variable sam_AttenFactor,sam_atten_err,sam_trans_err
276        Variable bgd_AttenFactor,bgd_atten_err
277        Variable emp_AttenFactor,emp_atten_err,emp_trans_err
278        Variable ii
279        String detStr
280        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
281        Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp
282       
283        // these values apply to all of the detectors
284        sam_AttenFactor = V_getAttenuator_transmission("SAM")
285        sam_atten_err = V_getAttenuator_trans_err("SAM")
286        bgd_AttenFactor = V_getAttenuator_transmission("BGD")
287        bgd_atten_err = V_getAttenuator_trans_err("BGD")
288        emp_AttenFactor = V_getAttenuator_transmission("EMP")
289        emp_atten_err = V_getAttenuator_trans_err("EMP")
290
291        //get relative monitor counts (should all be 10^8, since normalized in add step)
292        // get transmission and trans error for SAM, EMP
293        // TODO -- verify that the  call to V_getBeamMonNormData() is really rescaled to 10^8, and saved is the "true" count
294
295        tmonsam = V_getBeamMonNormData("SAM")           //monitor count in SAM
296        tsam = V_getSampleTransmission("SAM")           //SAM transmission
297        sam_trans_err = V_getSampleTransError("SAM")
298       
299        tmonemp = V_getBeamMonNormData("EMP")           //monitor count in EMP
300        temp = V_getSampleTransmission("EMP")                   //trans emp
301        emp_trans_err = V_getSampleTransError("EMP")
302       
303        tmonbgd = V_getBeamMonNormData("BGD")           //monitor count in BGD
304
305
306        // and now loop through all of the detectors
307        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
308                detStr = StringFromList(ii, ksDetectorListAll, ";")
309                Wave cor_data = V_getDetectorDataW("COR",detStr)
310                Wave cor_err = V_getDetectorDataErrW("COR",detStr)
311                Wave sam_data = V_getDetectorDataW("SAM",detStr)
312                Wave sam_err = V_getDetectorDataErrW("SAM",detStr)
313                Wave bgd_data = V_getDetectorDataW("BGD",detStr)
314                Wave bgd_err = V_getDetectorDataErrW("BGD",detStr)
315                Wave emp_data = V_getDetectorDataW("EMP",detStr)
316                Wave emp_err = V_getDetectorDataErrW("EMP",detStr)
317               
318        // to check for beam center mismatch -- simply warn, but do no shift
319        // NOTE - these are now reading the beam center in cm, not pixels
320
321                csam = V_getDet_beam_center_x("SAM",detStr)             //x center
322                rsam = V_getDet_beam_center_y("SAM",detStr)             //beam (x,y) define center of corrected field
323       
324                cbgd = V_getDet_beam_center_x("BGD",detStr)
325                rbgd = V_getDet_beam_center_y("BGD",detStr)
326       
327                cemp = V_getDet_beam_center_x("EMP",detStr)             //beamcenter of EMP
328                remp = V_getDet_beam_center_y("EMP",detStr)
329               
330                if(temp==0)
331                        DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction"
332                        temp=1
333                Endif
334       
335       
336                Duplicate/O cor_data cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp
337
338                // TODO -- document this, make a note, so everyone knows this is not done
339                // skip this part, but duplicate the results of no shift condition
340                //  where bgd_temp = input data, and noadd_bgd = 1 (so no data is zeroed out)
341               
342                        //get the shifted data arrays, EMP and BGD, each relative to SAM
343       
344                xshift = cbgd-csam
345                yshift = rbgd-rsam
346                if(abs(xshift) <= wcen)
347                        xshift = 0
348                Endif
349                if(abs(yshift) <= wcen)
350                        yshift = 0
351                Endif
352                // for the BGD file - alert if needed, generate dummy "pass-through" values
353                //
354                if(xshift != 0 || yshift != 0)
355                        DoAlert 0, "Beam center mismatch for BGD file. Data has NOT been corrected."
356                endif
357                bgd_temp = bgd_data             // no shift, no effect
358                noadd_bgd = 1
359                //GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift)            //bgd_temp
360       
361                xshift = cemp-csam
362                yshift = remp-rsam
363                if(abs(xshift) <= wcen)
364                        xshift = 0
365                Endif
366                if(abs(yshift) <= wcen)
367                        yshift = 0
368                Endif
369                // for the EMP file - alert if needed, generate dummy "pass-through" values
370                //
371                if(xshift != 0 || yshift != 0)
372                        DoAlert 0, "Beam center mismatch for EMP file. Data has NOT been corrected."
373                endif
374                emp_temp = emp_data // no shift, no effect
375                noadd_emp = 1
376                //GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift)            //emp_temp
377
378
379                // *******
380                //do the subtraction
381                fsam=1
382                femp = tmonsam/tmonemp          //this should be ==1 since normalized files
383                fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
384                cor1 = fsam*sam_data/sam_attenFactor - fbgd*bgd_temp/bgd_attenFactor
385                cor1 -= (tsam/temp)*(femp*emp_temp/emp_attenFactor - fbgd*bgd_temp/bgd_attenFactor)
386                cor1 *= noadd_bgd*noadd_emp             //zero out the array mismatch values
387       
388                cor_data = cor1         //the final result
389               
390                // do the error propagation piecewise   
391                Duplicate/O sam_err, tmp_a, tmp_b, tmp_c, tmp_d,c_val,d_val
392                tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
393               
394                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
395       
396                tmp_c = (sam_trans_err/temp)^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2
397                tmp_c += (tsam/temp^2)^2*emp_trans_err^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2
398               
399                tmp_d = (tsam/(temp*emp_attenFactor))^2*(emp_err)^2 + (tsam*emp_data/(temp*emp_attenFactor^2))^2*(emp_atten_err)^2
400       
401                cor_err = sqrt(tmp_a + tmp_b + tmp_c + tmp_d)
402       
403        endfor
404       
405        //we're done, get out w/no error
406
407       
408        //TODO -- do I update COR header?
409//      cor_text[1] = date() + " " + time()             //date + time stamp
410
411        KillWaves/Z cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp
412        Killwaves/Z tmp_a,tmp_b,tmp_c,tmp_d,c_val,d_val
413       
414        SetDataFolder root:
415        Return(0)
416End
417
418//background only
419// existence of data checked by dispatching routine
420// data has already been copied to COR folder
421Function V_CorrectMode_2()
422
423        //get SAM, BGD attenuation factor
424        Variable sam_AttenFactor,sam_atten_err,sam_trans_err
425        Variable bgd_AttenFactor,bgd_atten_err
426        Variable ii
427        String detStr
428        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
429        Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp
430       
431        // these values apply to all of the detectors
432        sam_AttenFactor = V_getAttenuator_transmission("SAM")
433        sam_atten_err = V_getAttenuator_trans_err("SAM")
434        bgd_AttenFactor = V_getAttenuator_transmission("BGD")
435        bgd_atten_err = V_getAttenuator_trans_err("BGD")
436
437        //get relative monitor counts (should all be 10^8, since normalized in add step)
438        // get transmission and trans error for SAM, EMP
439        // TODO -- verify that the  call to V_getBeamMonNormData() is really rescaled to 10^8, and saved is the "true" count
440
441        tmonsam = V_getBeamMonNormData("SAM")           //monitor count in SAM
442        tsam = V_getSampleTransmission("SAM")           //SAM transmission
443        sam_trans_err = V_getSampleTransError("SAM")
444       
445        tmonbgd = V_getBeamMonNormData("BGD")           //monitor count in BGD
446
447
448        // and now loop through all of the detectors
449        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
450                detStr = StringFromList(ii, ksDetectorListAll, ";")
451                Wave cor_data = V_getDetectorDataW("COR",detStr)
452                Wave cor_err = V_getDetectorDataErrW("COR",detStr)
453                Wave sam_data = V_getDetectorDataW("SAM",detStr)
454                Wave sam_err = V_getDetectorDataErrW("SAM",detStr)
455                Wave bgd_data = V_getDetectorDataW("BGD",detStr)
456                Wave bgd_err = V_getDetectorDataErrW("BGD",detStr)
457
458               
459        // to check for beam center mismatch -- simply warn, but do no shift
460        //
461
462                csam = V_getDet_beam_center_x("SAM",detStr)             //x center
463                rsam = V_getDet_beam_center_y("SAM",detStr)             //beam (x,y) define center of corrected field
464       
465                cbgd = V_getDet_beam_center_x("BGD",detStr)
466                rbgd = V_getDet_beam_center_y("BGD",detStr)
467
468       
469                Duplicate/O cor_data cor1,bgd_temp,noadd_bgd
470
471                // TODO -- document this, make a note, so everyone knows this is not done
472                // skip this part, but duplicate the results of no shift condition
473                //  where bgd_temp = input data, and noadd_bgd = 1 (so no data is zeroed out)
474               
475                        //get the shifted data array BGD, relative to SAM
476       
477                xshift = cbgd-csam
478                yshift = rbgd-rsam
479                if(abs(xshift) <= wcen)
480                        xshift = 0
481                Endif
482                if(abs(yshift) <= wcen)
483                        yshift = 0
484                Endif
485                // for the BGD file - alert if needed, generate dummy "pass-through" values
486                //
487                if(xshift != 0 || yshift != 0)
488                        DoAlert 0, "Beam center mismatch for BGD file. Data has NOT been corrected."
489                endif
490                bgd_temp = bgd_data             // no shift, no effect
491                noadd_bgd = 1
492                //GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift)            //bgd_temp
493       
494
495                // **********   
496                //do the sam-bgd subtraction,  deposit result in cor1
497                fsam = 1
498                fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
499               
500                //print "fsam,fbgd = ",fsam,fbgd
501               
502                cor1 = fsam*sam_data/sam_AttenFactor - fbgd*bgd_temp/bgd_AttenFactor
503                cor1 *= noadd_bgd               //zeros out regions where arrays do not overlap, one otherwise
504       
505                cor_data = cor1         //the final result
506
507        // do the error propagation piecewise   
508                Duplicate/O sam_err, tmp_a, tmp_b
509                tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
510               
511                tmp_b = (bgd_err/bgd_attenFactor)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2              //sig b ^2
512       
513                cor_err = sqrt(tmp_a + tmp_b)
514
515        endfor
516
517        //we're done, get out w/no error
518
519        // TODO -- do I update COR header?
520        //cor_text[1] = date() + " " + time()           //date + time stamp
521
522        KillWaves/Z cor1,bgd_temp,noadd_bgd
523        Killwaves/Z tmp_a,tmp_b
524
525        SetDataFolder root:
526        Return(0)
527End
528
529// empty subtraction only
530// data does exist, checked by dispatch routine
531//
532Function V_CorrectMode_3()
533
534        //get SAM, EMP attenuation factor
535        Variable sam_AttenFactor,sam_atten_err,sam_trans_err
536        Variable emp_AttenFactor,emp_atten_err,emp_trans_err
537        Variable ii
538        String detStr
539        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
540        Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp
541       
542        // these values apply to all of the detectors
543        sam_AttenFactor = V_getAttenuator_transmission("SAM")
544        sam_atten_err = V_getAttenuator_trans_err("SAM")
545        emp_AttenFactor = V_getAttenuator_transmission("EMP")
546        emp_atten_err = V_getAttenuator_trans_err("EMP")
547
548        //get relative monitor counts (should all be 10^8, since normalized in add step)
549        // get transmission and trans error for SAM, EMP
550        // TODO -- verify that the  call to V_getBeamMonNormData() is really rescaled to 10^8, and saved is the "true" count
551
552        tmonsam = V_getBeamMonNormData("SAM")           //monitor count in SAM
553        tsam = V_getSampleTransmission("SAM")           //SAM transmission
554        sam_trans_err = V_getSampleTransError("SAM")
555       
556        tmonemp = V_getBeamMonNormData("EMP")           //monitor count in EMP
557        temp = V_getSampleTransmission("EMP")                   //trans emp
558        emp_trans_err = V_getSampleTransError("EMP")
559       
560
561        // and now loop through all of the detectors
562        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
563                detStr = StringFromList(ii, ksDetectorListAll, ";")
564                Wave cor_data = V_getDetectorDataW("COR",detStr)
565                Wave cor_err = V_getDetectorDataErrW("COR",detStr)
566                Wave sam_data = V_getDetectorDataW("SAM",detStr)
567                Wave sam_err = V_getDetectorDataErrW("SAM",detStr)
568                Wave emp_data = V_getDetectorDataW("EMP",detStr)
569                Wave emp_err = V_getDetectorDataErrW("EMP",detStr)
570               
571        // to check for beam center mismatch -- simply warn, but do no shift
572        //
573
574                csam = V_getDet_beam_center_x("SAM",detStr)             //x center
575                rsam = V_getDet_beam_center_y("SAM",detStr)             //beam (x,y) define center of corrected field
576       
577                cemp = V_getDet_beam_center_x("EMP",detStr)             //beamcenter of EMP
578                remp = V_getDet_beam_center_y("EMP",detStr)
579               
580                if(temp==0)
581                        DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction"
582                        temp=1
583                Endif
584       
585       
586                Duplicate/O cor_data cor1,emp_temp,noadd_emp
587
588                // TODO -- document this, make a note, so everyone knows this is not done
589                // skip this part, but duplicate the results of no shift condition
590                //  where bgd_temp = input data, and noadd_bgd = 1 (so no data is zeroed out)
591               
592                        //get the shifted data array EMP, each relative to SAM
593       
594                xshift = cemp-csam
595                yshift = remp-rsam
596                if(abs(xshift) <= wcen)
597                        xshift = 0
598                Endif
599                if(abs(yshift) <= wcen)
600                        yshift = 0
601                Endif
602                // for the EMP file - alert if needed, generate dummy "pass-through" values
603                //
604                if(xshift != 0 || yshift != 0)
605                        DoAlert 0, "Beam center mismatch for EMP file. Data has NOT been corrected."
606                endif
607                emp_temp = emp_data // no shift, no effect
608                noadd_emp = 1
609                //GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift)            //emp_temp
610
611                // **********
612       
613                //do the sam-bgd subtraction,  deposit result in cor1
614                fsam = 1
615                femp = tmonsam/tmonemp          //this should be ==1 since normalized files
616               
617                cor1 = fsam*sam_data/sam_AttenFactor - femp*(tsam/temp)*emp_temp/emp_AttenFactor
618                cor1 *= noadd_emp               //zeros out regions where arrays do not overlap, one otherwise
619
620                cor_data = cor1         //the final result
621       
622        // do the error propagation piecewise   
623                Duplicate/O sam_err, tmp_a, tmp_c ,c_val
624                tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
625               
626                tmp_c = (sam_trans_err*emp_data/(temp*emp_attenFactor))^2 + (emp_err*tsam/(temp*emp_attenFactor))^2
627                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
628       
629                cor_err = sqrt(tmp_a + tmp_c)
630       
631        endfor
632       
633        //we're done, get out w/no error
634
635        // TODO -- do I update COR header?
636        //cor_text[1] = date() + " " + time()           //date + time stamp
637
638        KillWaves/Z cor1,emp_temp,noadd_emp
639        Killwaves/Z tmp_a,tmp_c,c_val
640
641        SetDataFolder root:
642        Return(0)
643End
644
645// NO subtraction - simply rescales for attenuators
646// SAM data does exist, checked by dispatch routine
647// SAM data has already been copied to COR (both are the same at the start of the function)
648//
649//  TODO -- do I need to rescale to sam_trans here ??
650//
651//
652Function V_CorrectMode_4()
653
654        //get SAM attenuation factor
655        Variable sam_AttenFactor,sam_atten_err,ii
656        String detStr
657       
658        sam_AttenFactor = V_getAttenuator_transmission("SAM")
659        sam_atten_err = V_getAttenuator_trans_err("SAM")
660
661
662        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
663                detStr = StringFromList(ii, ksDetectorListAll, ";")
664                Wave cor_data = V_getDetectorDataW("COR",detStr)
665                Wave cor_err = V_getDetectorDataErrW("COR",detStr)
666                Wave sam_data = V_getDetectorDataW("SAM",detStr)
667                Wave sam_err = V_getDetectorDataErrW("SAM",detStr)
668       
669       
670                cor_data = sam_data/sam_AttenFactor             //simply rescale the data
671       
672        // do the error propagation piecewise
673                cor_err = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2            //sig a ^2
674                cor_err = sqrt(cor_err)
675
676        endfor
677       
678        //TODO -- do I want to update COR header?
679//      cor_text[1] = date() + " " + time()             //date + time stamp
680
681        SetDataFolder root:
682        Return(0)
683End
684
685
686
687Function V_CorrectMode_11()
688
689        //get SAM, BGD, EMP attenuation factor
690        Variable sam_AttenFactor,sam_atten_err,sam_trans_err
691        Variable bgd_AttenFactor,bgd_atten_err
692        Variable emp_AttenFactor,emp_atten_err,emp_trans_err
693        Variable ii
694        String detStr
695        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
696        Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp
697        Variable savmon_sam,time_sam,time_drk
698       
699        // these values apply to all of the detectors
700        sam_AttenFactor = V_getAttenuator_transmission("SAM")
701        sam_atten_err = V_getAttenuator_trans_err("SAM")
702        bgd_AttenFactor = V_getAttenuator_transmission("BGD")
703        bgd_atten_err = V_getAttenuator_trans_err("BGD")
704        emp_AttenFactor = V_getAttenuator_transmission("EMP")
705        emp_atten_err = V_getAttenuator_trans_err("EMP")
706
707        //get relative monitor counts (should all be 10^8, since normalized in add step)
708        // get transmission and trans error for SAM, EMP
709        // TODO -- verify that the  call to V_getBeamMonNormData() is really rescaled to 10^8, and saved is the "true" count
710       
711        tmonsam = V_getBeamMonNormData("SAM")           //monitor count in SAM
712        tsam = V_getSampleTransmission("SAM")           //SAM transmission
713        sam_trans_err = V_getSampleTransError("SAM")
714       
715        tmonemp = V_getBeamMonNormData("EMP")           //monitor count in EMP
716        temp = V_getSampleTransmission("EMP")                   //trans emp
717        emp_trans_err = V_getSampleTransError("EMP")
718       
719        tmonbgd = V_getBeamMonNormData("BGD")           //monitor count in BGD
720
721        // for proper scaling, get the time and actual monitor counts
722        // TODO -- make sure that these calls are reading the proper values
723        savmon_sam = V_getBeamMonNormSaved_count("SAM")         //true monitor count in SAM
724        time_sam = V_getCount_time("SAM")               //count time SAM
725        time_drk = V_getCount_time("DRK")       //drk count time
726
727
728        // and now loop through all of the detectors
729        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
730                detStr = StringFromList(ii, ksDetectorListAll, ";")
731                Wave cor_data = V_getDetectorDataW("COR",detStr)
732                Wave cor_err = V_getDetectorDataErrW("COR",detStr)
733                Wave sam_data = V_getDetectorDataW("SAM",detStr)
734                Wave sam_err = V_getDetectorDataErrW("SAM",detStr)
735                Wave bgd_data = V_getDetectorDataW("BGD",detStr)
736                Wave bgd_err = V_getDetectorDataErrW("BGD",detStr)
737                Wave emp_data = V_getDetectorDataW("EMP",detStr)
738                Wave emp_err = V_getDetectorDataErrW("EMP",detStr)
739                Wave drk_data = V_getDetectorDataW("DRK",detStr)
740                Wave drk_err = V_getDetectorDataErrW("DRK",detStr)
741               
742        // to check for beam center mismatch -- simply warn, but do no shift
743        //
744
745                csam = V_getDet_beam_center_x("SAM",detStr)             //x center
746                rsam = V_getDet_beam_center_y("SAM",detStr)             //beam (x,y) define center of corrected field
747       
748                cbgd = V_getDet_beam_center_x("BGD",detStr)
749                rbgd = V_getDet_beam_center_y("BGD",detStr)
750       
751                cemp = V_getDet_beam_center_x("EMP",detStr)             //beamcenter of EMP
752                remp = V_getDet_beam_center_y("EMP",detStr)
753               
754                //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM
755                Duplicate/O drk_data drk_temp, drk_tmp_err
756                drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam)
757                drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam)                 //temporarily rescale the error of DRK 
758               
759                if(temp==0)
760                        DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction"
761                        temp=1
762                Endif
763       
764                Duplicate/O cor_data cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp
765
766                // TODO -- document this, make a note, so everyone knows this is not done
767                // skip this part, but duplicate the results of no shift condition
768                //  where bgd_temp = input data, and noadd_bgd = 1 (so no data is zeroed out)
769               
770                        //get the shifted data arrays, EMP and BGD, each relative to SAM
771       
772                xshift = cbgd-csam
773                yshift = rbgd-rsam
774                if(abs(xshift) <= wcen)
775                        xshift = 0
776                Endif
777                if(abs(yshift) <= wcen)
778                        yshift = 0
779                Endif
780                // for the BGD file - alert if needed, generate dummy "pass-through" values
781                //
782                if(xshift != 0 || yshift != 0)
783                        DoAlert 0, "Beam center mismatch for BGD file. Data has NOT been corrected."
784                endif
785                bgd_temp = bgd_data             // no shift, no effect
786                noadd_bgd = 1
787                //GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift)            //bgd_temp
788       
789                xshift = cemp-csam
790                yshift = remp-rsam
791                if(abs(xshift) <= wcen)
792                        xshift = 0
793                Endif
794                if(abs(yshift) <= wcen)
795                        yshift = 0
796                Endif
797                // for the EMP file - alert if needed, generate dummy "pass-through" values
798                //
799                if(xshift != 0 || yshift != 0)
800                        DoAlert 0, "Beam center mismatch for EMP file. Data has NOT been corrected."
801                endif
802                emp_temp = emp_data // no shift, no effect
803                noadd_emp = 1
804                //GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift)            //emp_temp
805
806
807                //always ignore the DRK center shift
808       
809                // ************
810                //do the subtraction
811                fsam=1
812                femp = tmonsam/tmonemp          //this should be ==1 since normalized files
813                fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
814                cor1 = fsam*sam_data/sam_attenFactor
815                cor1 -= (tsam/temp)*(femp*emp_temp/emp_attenFactor - fbgd*bgd_temp/bgd_attenFactor)
816                cor1 -= (fbgd*bgd_temp/bgd_attenFactor - drk_temp)
817                cor1 -= drk_temp/sam_attenFactor
818                cor1 *= noadd_bgd*noadd_emp             //zero out the array mismatch values
819
820                cor_data = cor1         //the final result
821               
822        // do the error propagation piecewise   
823                Duplicate/O sam_err, tmp_a, tmp_b, tmp_c, tmp_d,c_val,d_val
824                tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
825               
826                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
827       
828                tmp_c = (sam_trans_err/temp)^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2
829                tmp_c += (tsam/temp^2)^2*emp_trans_err^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2
830               
831                tmp_d = (tsam/(temp*emp_attenFactor))^2*(emp_err)^2 + (tsam*emp_data/(temp*emp_attenFactor^2))^2*(emp_atten_err)^2
832       
833                cor_err = sqrt(tmp_a + tmp_b + tmp_c + tmp_d + drk_tmp_err^2)
834       
835        endfor
836       
837        //we're done, get out w/no error
838
839
840        //TODO -- do I update COR header?
841//      cor_text[1] = date() + " " + time()             //date + time stamp
842
843        KillWaves/Z cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp,drk_temp
844        Killwaves/Z tmp_a,tmp_b,tmp_c,tmp_d,c_val,d_val,drk_tmp_err
845
846        SetDataFolder root:
847        Return(0)
848End
849
850//bgd and drk subtraction
851//
852Function V_CorrectMode_12()
853
854        //get SAM, BGD, EMP attenuation factor
855        Variable sam_AttenFactor,sam_atten_err,sam_trans_err
856        Variable bgd_AttenFactor,bgd_atten_err
857        Variable emp_AttenFactor,emp_atten_err,emp_trans_err
858        Variable ii
859        String detStr
860        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
861        Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp
862        Variable savmon_sam,time_sam,time_drk
863       
864        // these values apply to all of the detectors
865        sam_AttenFactor = V_getAttenuator_transmission("SAM")
866        sam_atten_err = V_getAttenuator_trans_err("SAM")
867        bgd_AttenFactor = V_getAttenuator_transmission("BGD")
868        bgd_atten_err = V_getAttenuator_trans_err("BGD")
869
870        //get relative monitor counts (should all be 10^8, since normalized in add step)
871        // get transmission and trans error for SAM, EMP
872        // TODO -- verify that the  call to V_getBeamMonNormData() is really rescaled to 10^8, and saved is the "true" count
873       
874        tmonsam = V_getBeamMonNormData("SAM")           //monitor count in SAM
875        tsam = V_getSampleTransmission("SAM")           //SAM transmission
876        sam_trans_err = V_getSampleTransError("SAM")
877       
878        tmonbgd = V_getBeamMonNormData("BGD")           //monitor count in BGD
879
880        // for proper scaling, get the time and actual monitor counts
881        // TODO -- make sure that these calls are reading the proper values
882        savmon_sam = V_getBeamMonNormSaved_count("SAM")         //true monitor count in SAM
883        time_sam = V_getCount_time("SAM")               //count time SAM
884        time_drk = V_getCount_time("DRK")       //drk count time
885
886
887        // and now loop through all of the detectors
888        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
889                detStr = StringFromList(ii, ksDetectorListAll, ";")
890                Wave cor_data = V_getDetectorDataW("COR",detStr)
891                Wave cor_err = V_getDetectorDataErrW("COR",detStr)
892                Wave sam_data = V_getDetectorDataW("SAM",detStr)
893                Wave sam_err = V_getDetectorDataErrW("SAM",detStr)
894                Wave bgd_data = V_getDetectorDataW("BGD",detStr)
895                Wave bgd_err = V_getDetectorDataErrW("BGD",detStr)
896                Wave drk_data = V_getDetectorDataW("DRK",detStr)
897                Wave drk_err = V_getDetectorDataErrW("DRK",detStr)
898               
899        // to check for beam center mismatch -- simply warn, but do no shift
900        //
901
902                csam = V_getDet_beam_center_x("SAM",detStr)             //x center
903                rsam = V_getDet_beam_center_y("SAM",detStr)             //beam (x,y) define center of corrected field
904       
905                cbgd = V_getDet_beam_center_x("BGD",detStr)
906                rbgd = V_getDet_beam_center_y("BGD",detStr)
907               
908                //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM
909                Duplicate/O drk_data drk_temp, drk_tmp_err
910                drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam)
911                drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam)                 //temporarily rescale the error of DRK 
912               
913                if(temp==0)
914                        DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction"
915                        temp=1
916                Endif
917       
918                Duplicate/O cor_data cor1,bgd_temp,noadd_bgd
919
920                // TODO -- document this, make a note, so everyone knows this is not done
921                // skip this part, but duplicate the results of no shift condition
922                //  where bgd_temp = input data, and noadd_bgd = 1 (so no data is zeroed out)
923               
924                        //get the shifted data arrays, EMP and BGD, each relative to SAM
925       
926                xshift = cbgd-csam
927                yshift = rbgd-rsam
928                if(abs(xshift) <= wcen)
929                        xshift = 0
930                Endif
931                if(abs(yshift) <= wcen)
932                        yshift = 0
933                Endif
934                // for the BGD file - alert if needed, generate dummy "pass-through" values
935                //
936                if(xshift != 0 || yshift != 0)
937                        DoAlert 0, "Beam center mismatch for BGD file. Data has NOT been corrected."
938                endif
939                bgd_temp = bgd_data             // no shift, no effect
940                noadd_bgd = 1
941                //GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift)            //bgd_temp
942
943
944                //always ignore the DRK center shift
945
946                // ************
947                //do the sam-bgd subtraction,  deposit result in cor1
948                fsam = 1
949                fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
950               
951                cor1 = fsam*sam_data/sam_AttenFactor + fbgd*tsam*bgd_temp/bgd_AttenFactor
952                cor1 += -1*(fbgd*bgd_temp/bgd_attenFactor - drk_temp) - drk_temp/sam_attenFactor
953                cor1 *= noadd_bgd               //zeros out regions where arrays do not overlap, one otherwise
954
955                cor_data = cor1         //the final result
956       
957        // do the error propagation piecewise   
958                Duplicate/O sam_err, tmp_a, tmp_b
959                tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
960               
961                tmp_b = (bgd_err/bgd_attenFactor)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2              //sig b ^2
962       
963                cor_err = sqrt(tmp_a + tmp_b + drk_tmp_err^2)
964
965        endfor
966       
967        //we're done, get out w/no error
968       
969        // TODO -- do I update COR header?
970//      cor_text[1] = date() + " " + time()             //date + time stamp
971
972        KillWaves/Z cor1,bgd_temp,noadd_bgd,drk_temp
973        Killwaves/Z tmp_a,tmp_b,drk_tmp_err
974
975        SetDataFolder root:
976        Return(0)
977End
978
979//EMP and DRK subtractions
980// all data exists, DRK is on a time basis (noNorm)
981//scale DRK by monitor count scaling factor and the ratio of couting times
982//to place the DRK file on equal footing
983Function V_CorrectMode_13()
984
985        //get SAM, EMP attenuation factor
986        Variable sam_AttenFactor,sam_atten_err,sam_trans_err
987        Variable bgd_AttenFactor,bgd_atten_err
988        Variable emp_AttenFactor,emp_atten_err,emp_trans_err
989        Variable ii
990        String detStr
991        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
992        Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp
993        Variable savmon_sam,time_sam,time_drk
994       
995        // these values apply to all of the detectors
996        sam_AttenFactor = V_getAttenuator_transmission("SAM")
997        sam_atten_err = V_getAttenuator_trans_err("SAM")
998        emp_AttenFactor = V_getAttenuator_transmission("EMP")
999        emp_atten_err = V_getAttenuator_trans_err("EMP")
1000
1001        //get relative monitor counts (should all be 10^8, since normalized in add step)
1002        // get transmission and trans error for SAM, EMP
1003        // TODO -- verify that the  call to V_getBeamMonNormData() is really rescaled to 10^8, and saved is the "true" count
1004       
1005        tmonsam = V_getBeamMonNormData("SAM")           //monitor count in SAM
1006        tsam = V_getSampleTransmission("SAM")           //SAM transmission
1007        sam_trans_err = V_getSampleTransError("SAM")
1008       
1009        tmonemp = V_getBeamMonNormData("EMP")           //monitor count in EMP
1010        temp = V_getSampleTransmission("EMP")                   //trans emp
1011        emp_trans_err = V_getSampleTransError("EMP")
1012
1013        // for proper scaling, get the time and actual monitor counts
1014        // TODO -- make sure that these calls are reading the proper values
1015        savmon_sam = V_getBeamMonNormSaved_count("SAM")         //true monitor count in SAM
1016        time_sam = V_getCount_time("SAM")               //count time SAM
1017        time_drk = V_getCount_time("DRK")       //drk count time
1018
1019
1020        // and now loop through all of the detectors
1021        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
1022                detStr = StringFromList(ii, ksDetectorListAll, ";")
1023                Wave cor_data = V_getDetectorDataW("COR",detStr)
1024                Wave cor_err = V_getDetectorDataErrW("COR",detStr)
1025                Wave sam_data = V_getDetectorDataW("SAM",detStr)
1026                Wave sam_err = V_getDetectorDataErrW("SAM",detStr)
1027                Wave emp_data = V_getDetectorDataW("EMP",detStr)
1028                Wave emp_err = V_getDetectorDataErrW("EMP",detStr)
1029                Wave drk_data = V_getDetectorDataW("DRK",detStr)
1030                Wave drk_err = V_getDetectorDataErrW("DRK",detStr)
1031               
1032        // to check for beam center mismatch -- simply warn, but do no shift
1033        //
1034
1035                csam = V_getDet_beam_center_x("SAM",detStr)             //x center
1036                rsam = V_getDet_beam_center_y("SAM",detStr)             //beam (x,y) define center of corrected field
1037       
1038                cemp = V_getDet_beam_center_x("EMP",detStr)             //beamcenter of EMP
1039                remp = V_getDet_beam_center_y("EMP",detStr)
1040               
1041                //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM
1042                Duplicate/O drk_data drk_temp, drk_tmp_err
1043                drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam)
1044                drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam)                 //temporarily rescale the error of DRK 
1045               
1046                if(temp==0)
1047                        DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction"
1048                        temp=1
1049                Endif
1050       
1051                Duplicate/O cor_data cor1,emp_temp,noadd_emp
1052
1053                // TODO -- document this, make a note, so everyone knows this is not done
1054                // skip this part, but duplicate the results of no shift condition
1055                //  where bgd_temp = input data, and noadd_bgd = 1 (so no data is zeroed out)
1056               
1057                        //get the shifted data arrays, EMP , each relative to SAM
1058       
1059                xshift = cemp-csam
1060                yshift = remp-rsam
1061                if(abs(xshift) <= wcen)
1062                        xshift = 0
1063                Endif
1064                if(abs(yshift) <= wcen)
1065                        yshift = 0
1066                Endif
1067                // for the EMP file - alert if needed, generate dummy "pass-through" values
1068                //
1069                if(xshift != 0 || yshift != 0)
1070                        DoAlert 0, "Beam center mismatch for EMP file. Data has NOT been corrected."
1071                endif
1072                emp_temp = emp_data // no shift, no effect
1073                noadd_emp = 1
1074                //GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift)            //emp_temp
1075
1076
1077                //always ignore the DRK center shift
1078               
1079                // ***************     
1080                //do the sam-bgd subtraction,  deposit result in cor1
1081                fsam = 1
1082                femp = tmonsam/tmonemp          //this should be ==1 since normalized files
1083               
1084                cor1 = fsam*sam_data/sam_AttenFactor - femp*(tsam/temp)*emp_temp/emp_AttenFactor
1085                cor1 += drk_temp - drk_temp/sam_attenFactor
1086                cor1 *= noadd_emp               //zeros out regions where arrays do not overlap, one otherwise
1087
1088                cor_data = cor1         //the final result
1089       
1090        // do the error propagation piecewise   
1091                Duplicate/O sam_err, tmp_a, tmp_c, c_val
1092                tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
1093               
1094                tmp_c = (sam_trans_err*emp_data/(temp*emp_attenFactor))^2 + (emp_err*tsam/(temp*emp_attenFactor))^2
1095                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
1096               
1097                cor_err = sqrt(tmp_a + tmp_c + drk_tmp_err^2)
1098
1099        endfor
1100               
1101        //we're done, get out w/no error
1102
1103        // TODO -- do I update COR header?
1104//      cor_text[1] = date() + " " + time()             //date + time stamp
1105
1106        KillWaves/Z cor1,emp_temp,noadd_emp,drk_temp
1107        Killwaves/Z tmp_a,tmp_c,c_val,drk_tmp_err
1108
1109        SetDataFolder root:
1110        Return(0)
1111End
1112
1113// ONLY drk subtraction
1114//
1115Function V_CorrectMode_14()
1116
1117        //get SAM, EMP attenuation factor
1118        Variable sam_AttenFactor,sam_atten_err,sam_trans_err
1119        Variable bgd_AttenFactor,bgd_atten_err
1120        Variable emp_AttenFactor,emp_atten_err,emp_trans_err
1121        Variable ii
1122        String detStr
1123        Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd
1124        Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp
1125        Variable savmon_sam,time_sam,time_drk
1126       
1127        // these values apply to all of the detectors
1128        sam_AttenFactor = V_getAttenuator_transmission("SAM")
1129        sam_atten_err = V_getAttenuator_trans_err("SAM")
1130
1131
1132        //get relative monitor counts (should all be 10^8, since normalized in add step)
1133        // get transmission and trans error for SAM, EMP
1134        // TODO -- verify that the  call to V_getBeamMonNormData() is really rescaled to 10^8, and saved is the "true" count
1135       
1136        tmonsam = V_getBeamMonNormData("SAM")           //monitor count in SAM
1137        tsam = V_getSampleTransmission("SAM")           //SAM transmission
1138        sam_trans_err = V_getSampleTransError("SAM")
1139       
1140
1141        // for proper scaling, get the time and actual monitor counts
1142        // TODO -- make sure that these calls are reading the proper values
1143        savmon_sam = V_getBeamMonNormSaved_count("SAM")         //true monitor count in SAM
1144        time_sam = V_getCount_time("SAM")               //count time SAM
1145        time_drk = V_getCount_time("DRK")       //drk count time
1146
1147
1148        // and now loop through all of the detectors
1149        for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1)
1150                detStr = StringFromList(ii, ksDetectorListAll, ";")
1151                Wave cor_data = V_getDetectorDataW("COR",detStr)
1152                Wave cor_err = V_getDetectorDataErrW("COR",detStr)
1153                Wave sam_data = V_getDetectorDataW("SAM",detStr)
1154                Wave sam_err = V_getDetectorDataErrW("SAM",detStr)
1155                Wave drk_data = V_getDetectorDataW("DRK",detStr)
1156                Wave drk_err = V_getDetectorDataErrW("DRK",detStr)
1157               
1158               
1159                //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM
1160                Duplicate/O drk_data drk_temp, drk_tmp_err
1161                drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam)
1162                drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam)                 //temporarily rescale the error of DRK 
1163               
1164                if(temp==0)
1165                        DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction"
1166                        temp=1
1167                Endif
1168       
1169                Duplicate/O cor_data cor1
1170
1171                //always ignore the DRK center shift
1172
1173                // ************
1174                //do the subtraction,  deposit result in cor1
1175                fsam = 1
1176                fbgd = tmonsam/tmonbgd  //this should be ==1 since normalized files
1177               
1178                //correct sam for attenuators, and do the same to drk, since it was scaled to sam count time
1179                cor1 = fsam*sam_data/sam_AttenFactor  - drk_temp/sam_attenFactor
1180
1181                cor_data = cor1         //the final result
1182       
1183        // do the error propagation piecewise   
1184                Duplicate/O sam_err, tmp_a
1185                tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2              //sig a ^2
1186       
1187                cor_err = sqrt(tmp_a + drk_tmp_err^2)
1188
1189        endfor
1190               
1191        //we're done, get out w/no error
1192       
1193        //TODO -- do I update COR header?
1194//      cor_text[1] = date() + " " + time()             //date + time stamp
1195
1196        KillWaves/Z cor1,bgd_temp,noadd_bgd,drk_temp
1197        Killwaves/Z tmp_a,tmp_b,tmp_c,tmp_d,c_val,d_val,drk_tmp_err
1198
1199        SetDataFolder root:
1200        Return(0)
1201End
1202
1203
1204//
1205// For VSANS - this should go away. if there is a mismatch, don't try to fudge it.
1206//
1207//
1208//function to return the shifted contents of a data array for subtraction
1209//(SLOW) if ShiftSum is called
1210//data_in is input
1211//data_out is shifted matrix
1212//noadd_mat =1 if shift matrix is valid, =0 if no data
1213//
1214//if no shift is required, data_in is returned and noadd_mat =1 (all valid)
1215//
1216xFunction V_GetShiftedArray(data_in,data_out,noadd_mat,xshift,yshift)
1217        WAVE data_in,data_out,noadd_mat
1218        Variable xshift,yshift
1219
1220        Variable ii=0,jj=0
1221        noadd_mat = 1           //initialize to 1
1222       
1223        If((xshift != 0) || (yshift != 0))
1224//      If((abs(xshift) >= 0.01) || (abs(yshift) >= 0.01))                      //APR09 - loosen tolerance to handle ICE "precision"
1225                DoAlert 1,"Do you want to ignore the beam center mismatch?"
1226                if(V_flag==1)           //yes -> just go on
1227                        xshift=0
1228                        yshift=0
1229                endif
1230        else
1231                // "mismatch" is simply a python type conversion error
1232                xshift=0
1233                yshift=0
1234        endif
1235       
1236        If((xshift == 0) && (yshift == 0))
1237                data_out=data_in                //no change
1238                noadd_mat = 1                   //use all of the data
1239                return(0)
1240        endif
1241       
1242        NVAR pixelsX = root:myGlobals:gNPixelsX
1243        NVAR pixelsY = root:myGlobals:gNPixelsY
1244       
1245        Print "beamcenter shift x,y = ",xshift,yshift
1246        Make/O/N=1 noadd
1247        for(ii=0;ii<pixelsX;ii+=1)
1248                for(jj=0;jj<pixelsY;jj+=1)
1249                        //get the contribution of the shifted data
1250                        data_out[ii][jj] = ShiftSum(data_in,ii,jj,xshift,yshift,noadd)
1251                        if(noadd[0])
1252                                noadd_mat[ii][jj] = 0   //shift is off the detector
1253                        endif
1254                endfor
1255        endfor
1256        return(0)
1257End
1258
1259//
1260//utility function that checks if data exists in a data folder
1261//checks only for the existence of data in detector FL - no other waves
1262//
1263Function V_WorkDataExists(type)
1264        String type
1265       
1266        String destPath=""
1267        destPath = "root:Packages:NIST:VSANS:"+type+":entry:instrument:detector_FL:data"
1268        if(WaveExists($destpath) == 0)
1269                Print "There is no work file in "+type
1270                Return(1)               //error condition
1271        else
1272                // data exists, assume everything else is OK and proceed.
1273                return(0)
1274        Endif
1275End
1276
1277//////////////////
1278// bunch of utility junk to catch
1279// sample transmission = 1
1280// and handle (too many) options
1281//
1282xFunction V_GetNewTrans(oldTrans,type)
1283        Variable oldTrans
1284        String type
1285       
1286        Variable newTrans,newCode
1287        if (oldTrans!=1)
1288                return(oldTrans)                //get out now if trans != 1, don't change anything
1289        endif
1290        //get input from the user
1291        NewDataFolder/O root:myGlobals:tmp_trans
1292        Variable/G root:myGlobals:tmp_trans:inputTrans=0.92
1293        Variable/G root:myGlobals:tmp_trans:returnCode=0
1294        DoTransInput(type)
1295        NVAR inputTrans=root:myGlobals:tmp_trans:inputTrans
1296        NVAR code=root:myGlobals:tmp_trans:returnCode
1297        newTrans=inputTrans             //keep a copy before deleting everything
1298        newCode=code
1299        if(newCode==4)
1300                Abort "Aborting correction. Use the Transmission Panel to calculate transmissions"
1301        Endif
1302//      printf "You entered %g and the code is %g\r",newTrans,newCode
1303//      KillDataFolder root:tmp_trans
1304       
1305        if(newCode==1)
1306                Variable/G root:Packages:NIST:gDoTransCheck=0   //turn off checking
1307        endif
1308       
1309        if(newcode==2)          //user changed trans value
1310                return(newTrans)
1311        else
1312                return(oldTrans)        //All other cases, user did not change value
1313        endif
1314end
1315
1316xFunction V_IgnoreNowButton(ctrlName) : ButtonControl
1317        String ctrlName
1318       
1319//      Print "ignore now"
1320        NVAR val=root:myGlobals:tmp_trans:returnCode
1321        val=0           //code for ignore once
1322       
1323        DoWindow/K tmp_GetInputPanel            // Kill self
1324End
1325
1326xFunction V_DoTransInput(str)
1327        String str
1328       
1329        NewPanel /W=(150,50,361,294)
1330        DoWindow/C tmp_GetInputPanel            // Set to an unlikely name
1331        DrawText 15,23,"The "+str+" Transmission = 1"
1332        DrawText 15,43,"What do you want to do?"
1333        DrawText 15,125,"(Reset this in Preferences)"
1334        SetVariable setvar0,pos={20,170},size={160,17},limits={0,1,0.01}
1335        SetVariable setvar0,value= root:myGlobals:tmp_trans:inputTrans,title="New Transmission"
1336
1337        Button button0,pos={36,56},size={120,20},proc=IgnoreNowButton,title="Ignore This Time"
1338        Button button1,pos={36,86},size={120,20},proc=IgnoreAlwaysButtonProc,title="Ignore Always"
1339        Button button2,pos={36,143},size={120,20},proc=UseNewValueButtonProc,title="Use New Value"
1340        Button button3,pos={36,213},size={120,20},proc=AbortCorrectionButtonProc,title="Abort Correction"
1341        PauseForUser tmp_GetInputPanel
1342End
1343
1344xFunction V_IgnoreAlwaysButtonProc(ctrlName) : ButtonControl
1345        String ctrlName
1346
1347//      Print "ignore always"
1348        NVAR val=root:myGlobals:tmp_trans:returnCode
1349        val=1           //code for ignore always
1350        DoWindow/K tmp_GetInputPanel            // Kill self
1351End
1352
1353xFunction V_UseNewValueButtonProc(ctrlName) : ButtonControl
1354        String ctrlName
1355
1356//      Print "use new Value"
1357        NVAR val=root:myGlobals:tmp_trans:returnCode
1358        val=2           //code for use new Value
1359        DoWindow/K tmp_GetInputPanel            // Kill self
1360End
1361
1362xFunction V_AbortCorrectionButtonProc(ctrlName) : ButtonControl
1363        String ctrlName
1364
1365//      Print "Abort"
1366        NVAR val=root:myGlobals:tmp_trans:returnCode
1367        val=4           //code for abort
1368        DoWindow/K tmp_GetInputPanel            // Kill self
1369End
1370
Note: See TracBrowser for help on using the repository browser.