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

Last change on this file was 1242, checked in by srkline, 3 years ago

updating the IgorVersion? pragma to v7.0 for all files to be consistent.

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