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

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

minor changes to prefix functions with "V_" to avoid conflicts with non-VSANS functions.

File size: 44.1 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_getMonitorCount() is really rescaled to 10^8, and saved is the "true" count
294
295        tmonsam = V_getMonitorCount("SAM")              //monitor count in SAM
296        tsam = V_getSampleTransmission("SAM")           //SAM transmission
297        sam_trans_err = V_getSampleTransError("SAM")
298       
299        tmonemp = V_getMonitorCount("EMP")              //monitor count in EMP
300        temp = V_getSampleTransmission("EMP")                   //trans emp
301        emp_trans_err = V_getSampleTransError("EMP")
302       
303        tmonbgd = V_getMonitorCount("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        //
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_getMonitorCount() is really rescaled to 10^8, and saved is the "true" count
440
441        tmonsam = V_getMonitorCount("SAM")              //monitor count in SAM
442        tsam = V_getSampleTransmission("SAM")           //SAM transmission
443        sam_trans_err = V_getSampleTransError("SAM")
444       
445        tmonbgd = V_getMonitorCount("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_getMonitorCount() is really rescaled to 10^8, and saved is the "true" count
551
552        tmonsam = V_getMonitorCount("SAM")              //monitor count in SAM
553        tsam = V_getSampleTransmission("SAM")           //SAM transmission
554        sam_trans_err = V_getSampleTransError("SAM")
555       
556        tmonemp = V_getMonitorCount("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_getMonitorCount() is really rescaled to 10^8, and saved is the "true" count
710       
711        tmonsam = V_getMonitorCount("SAM")              //monitor count in SAM
712        tsam = V_getSampleTransmission("SAM")           //SAM transmission
713        sam_trans_err = V_getSampleTransError("SAM")
714       
715        tmonemp = V_getMonitorCount("EMP")              //monitor count in EMP
716        temp = V_getSampleTransmission("EMP")                   //trans emp
717        emp_trans_err = V_getSampleTransError("EMP")
718       
719        tmonbgd = V_getMonitorCount("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_getMonitorCount() is really rescaled to 10^8, and saved is the "true" count
873       
874        tmonsam = V_getMonitorCount("SAM")              //monitor count in SAM
875        tsam = V_getSampleTransmission("SAM")           //SAM transmission
876        sam_trans_err = V_getSampleTransError("SAM")
877       
878        tmonbgd = V_getMonitorCount("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_getMonitorCount() is really rescaled to 10^8, and saved is the "true" count
1004       
1005        tmonsam = V_getMonitorCount("SAM")              //monitor count in SAM
1006        tsam = V_getSampleTransmission("SAM")           //SAM transmission
1007        sam_trans_err = V_getSampleTransError("SAM")
1008       
1009        tmonemp = V_getMonitorCount("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_getMonitorCount() is really rescaled to 10^8, and saved is the "true" count
1135       
1136        tmonsam = V_getMonitorCount("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.