- Timestamp:
- Jan 19, 2017 7:53:30 AM (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_Correct.ipf
r1019 r1020 1 #pragma rtGlobals=1 // Use modern global access method. 2 #pragma version=5.0 3 #pragma IgorVersion=6.1 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 4 18 5 19 … … 29 43 30 44 // 31 // unusedtest procedure for Correct() function45 // test procedure for Correct() function 32 46 //must be updated to include "mode" parameter before re-use 33 47 // … … 139 153 return(err) 140 154 Endif 141 // TODOerr = V_CorrectMode_1()155 err = V_CorrectMode_1() 142 156 break 143 157 case 2: … … 146 160 return(err) 147 161 Endif 148 // TODOerr = V_CorrectMode_2()162 err = V_CorrectMode_2() 149 163 break 150 164 case 3: … … 165 179 // endif 166 180 167 // TODOerr = V_CorrectMode_3()181 err = V_CorrectMode_3() 168 182 break 169 183 case 4: … … 194 208 return(err) 195 209 Endif 196 // TODOerr = V_CorrectMode_11()210 err = V_CorrectMode_11() 197 211 break 198 212 case 12: … … 205 219 return(err) 206 220 Endif 207 // TODOerr = V_CorrectMode_12()221 err = V_CorrectMode_12() 208 222 break 209 223 case 13: … … 227 241 return(err) 228 242 Endif 229 // TODOerr = V_CorrectMode_13()243 err = V_CorrectMode_13() 230 244 break 231 245 case 14: … … 234 248 return(err) 235 249 Endif 236 // TODOerr = V_CorrectMode_14()250 err = V_CorrectMode_14() 237 251 break 238 252 default: //something wrong … … 253 267 // was freshly loaded. added final copy of cor result to cor:data and cor:linear_data 254 268 // 255 xFunction V_CorrectMode_1() 256 257 //create the necessary wave references 258 WAVE sam_data=$"root:Packages:NIST:SAM:linear_data" 259 WAVE sam_reals=$"root:Packages:NIST:SAM:realsread" 260 WAVE sam_ints=$"root:Packages:NIST:SAM:integersread" 261 WAVE/T sam_text=$"root:Packages:NIST:SAM:textread" 262 WAVE bgd_data=$"root:Packages:NIST:BGD:linear_data" 263 WAVE bgd_reals=$"root:Packages:NIST:BGD:realsread" 264 WAVE bgd_ints=$"root:Packages:NIST:BGD:integersread" 265 WAVE/T bgd_text=$"root:Packages:NIST:BGD:textread" 266 WAVE emp_data=$"root:Packages:NIST:EMP:linear_data" 267 WAVE emp_reals=$"root:Packages:NIST:EMP:realsread" 268 WAVE emp_ints=$"root:Packages:NIST:EMP:integersread" 269 WAVE/T emp_text=$"root:Packages:NIST:EMP:textread" 270 WAVE cor_data=$"root:Packages:NIST:COR:linear_data" 271 WAVE/T cor_text=$"root:Packages:NIST:COR:textread" 272 273 // needed to propagate error 274 WAVE cor_data_display=$"root:Packages:NIST:COR:data" //just for the final copy 275 WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error" 276 WAVE bgd_err =$"root:Packages:NIST:BGD:linear_data_error" 277 WAVE emp_err =$"root:Packages:NIST:EMP:linear_data_error" 278 WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error" 279 280 Variable sam_trans_err,emp_trans_err 281 sam_trans_err = sam_reals[41] 282 emp_trans_err = emp_reals[41] 283 284 285 //get sam and bgd attenuation factors 286 String fileStr="" 287 Variable lambda,attenNo,sam_AttenFactor,bgd_attenFactor,emp_AttenFactor 269 Function V_CorrectMode_1() 270 271 //get SAM, BGD, EMP attenuation factor 272 Variable sam_AttenFactor,sam_atten_err,sam_trans_err 273 Variable bgd_AttenFactor,bgd_atten_err 274 Variable emp_AttenFactor,emp_atten_err,emp_trans_err 275 Variable ii 276 String detStr 288 277 Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd 289 278 Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp 290 Variable sam_atten_err,emp_atten_err,bgd_atten_err 291 fileStr = sam_text[3] 292 lambda = sam_reals[26] 293 attenNo = sam_reals[3] 294 sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err) 295 fileStr = bgd_text[3] 296 lambda = bgd_reals[26] 297 attenNo = bgd_reals[3] 298 bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,bgd_atten_err) 299 fileStr = emp_text[3] 300 lambda = emp_reals[26] 301 attenNo = emp_reals[3] 302 emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,emp_atten_err) 303 304 //get relative monitor counts (should all be 10^8, since normalized in add step) 305 tmonsam = sam_reals[0] //monitor count in SAM 306 tsam = sam_reals[4] //SAM transmission 307 csam = sam_reals[16] //x center 308 rsam = sam_reals[17] //beam (x,y) define center of corrected field 309 tmonbgd = bgd_reals[0] //monitor count in BGD 310 cbgd = bgd_reals[16] 311 rbgd = bgd_reals[17] 312 tmonemp = emp_reals[0] //monitor count in EMP 313 temp = emp_reals[4] //trans emp 314 cemp = emp_reals[16] //beamcenter of EMP 315 remp = emp_reals[17] 316 317 if(temp==0) 318 DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction" 319 temp=1 320 Endif 321 322 NVAR pixelsX = root:myGlobals:gNPixelsX 323 NVAR pixelsY = root:myGlobals:gNPixelsY 324 325 //get the shifted data arrays, EMP and BGD, each relative to SAM 326 Make/D/O/N=(pixelsX,pixelsY) cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp 327 xshift = cbgd-csam 328 yshift = rbgd-rsam 329 if(abs(xshift) <= wcen) 330 xshift = 0 331 Endif 332 if(abs(yshift) <= wcen) 333 yshift = 0 334 Endif 335 GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift) //bgd_temp 336 337 xshift = cemp-csam 338 yshift = remp-rsam 339 if(abs(xshift) <= wcen) 340 xshift = 0 341 Endif 342 if(abs(yshift) <= wcen) 343 yshift = 0 344 Endif 345 GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift) //emp_temp 346 347 //do the subtraction 348 fsam=1 349 femp = tmonsam/tmonemp //this should be ==1 since normalized files 350 fbgd = tmonsam/tmonbgd //this should be ==1 since normalized files 351 cor1 = fsam*sam_data/sam_attenFactor - fbgd*bgd_temp/bgd_attenFactor 352 cor1 -= (tsam/temp)*(femp*emp_temp/emp_attenFactor - fbgd*bgd_temp/bgd_attenFactor) 353 cor1 *= noadd_bgd*noadd_emp //zero out the array mismatch values 354 355 // do the error propagation piecewise 356 Duplicate/O sam_err, tmp_a, tmp_b, tmp_c, tmp_d,c_val,d_val 357 tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2 //sig a ^2 358 359 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 360 361 tmp_c = (sam_trans_err/temp)^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2 362 tmp_c += (tsam/temp^2)^2*emp_trans_err^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2 363 364 tmp_d = (tsam/(temp*emp_attenFactor))^2*(emp_err)^2 + (tsam*emp_data/(temp*emp_attenFactor^2))^2*(emp_atten_err)^2 365 366 cor_err = sqrt(tmp_a + tmp_b + tmp_c + tmp_d) 367 368 //we're done, get out w/no error 369 //set the COR data and linear_data to the result 370 cor_data = cor1 371 cor_data_display = cor1 372 373 //update COR header 374 cor_text[1] = date() + " " + time() //date + time stamp 375 376 KillWaves/Z cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp 377 Killwaves/Z tmp_a,tmp_b,tmp_c,tmp_d,c_val,d_val 378 SetDataFolder root: 379 Return(0) 380 End 381 382 //background only 383 // existence of data checked by dispatching routine 384 // data has already been copied to COR folder 385 xFunction V_CorrectMode_2() 386 387 //create the necessary wave references 388 WAVE sam_data=$"root:Packages:NIST:SAM:linear_data" 389 WAVE sam_reals=$"root:Packages:NIST:SAM:realsread" 390 WAVE sam_ints=$"root:Packages:NIST:SAM:integersread" 391 WAVE/T sam_text=$"root:Packages:NIST:SAM:textread" 392 WAVE bgd_data=$"root:Packages:NIST:BGD:linear_data" 393 WAVE bgd_reals=$"root:Packages:NIST:BGD:realsread" 394 WAVE bgd_ints=$"root:Packages:NIST:BGD:integersread" 395 WAVE/T bgd_text=$"root:Packages:NIST:BGD:textread" 396 WAVE cor_data=$"root:Packages:NIST:COR:linear_data" 397 WAVE/T cor_text=$"root:Packages:NIST:COR:textread" 398 399 // needed to propagate error 400 WAVE cor_data_display=$"root:Packages:NIST:COR:data" //just for the final copy 401 WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error" 402 WAVE bgd_err =$"root:Packages:NIST:BGD:linear_data_error" 403 WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error" 404 405 Variable sam_trans_err 406 sam_trans_err = sam_reals[41] 407 408 409 //get sam and bgd attenuation factors 410 String fileStr="" 411 Variable lambda,attenNo,sam_AttenFactor,bgd_attenFactor 412 Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd 413 Variable wcen=0.001 414 Variable sam_atten_err,bgd_atten_err 415 fileStr = sam_text[3] 416 lambda = sam_reals[26] 417 attenNo = sam_reals[3] 418 sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err) 419 fileStr = bgd_text[3] 420 lambda = bgd_reals[26] 421 attenNo = bgd_reals[3] 422 bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,bgd_atten_err) 423 424 //Print "atten = ",sam_attenFactor,bgd_attenFactor 425 426 //get relative monitor counts (should all be 10^8, since normalized in add step) 427 tmonsam = sam_reals[0] //monitor count in SAM 428 csam = sam_reals[16] //x center 429 rsam = sam_reals[17] //beam (x,y) define center of corrected field 430 tmonbgd = bgd_reals[0] //monitor count in BGD 431 cbgd = bgd_reals[16] 432 rbgd = bgd_reals[17] 433 434 // set up beamcenter shift, relative to SAM 435 xshift = cbgd-csam 436 yshift = rbgd-rsam 437 if(abs(xshift) <= wcen) 438 xshift = 0 439 Endif 440 if(abs(yshift) <= wcen) 441 yshift = 0 442 Endif 443 444 NVAR pixelsX = root:myGlobals:gNPixelsX 445 NVAR pixelsY = root:myGlobals:gNPixelsY 446 //get shifted data arrays, relative to SAM 447 Make/D/O/N=(pixelsX,pixelsY) cor1,bgd_temp,noadd_bgd //temp arrays 448 GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift) //bgd_temp is the BGD 449 450 //do the sam-bgd subtraction, deposit result in cor1 451 fsam = 1 452 fbgd = tmonsam/tmonbgd //this should be ==1 since normalized files 453 454 //print "fsam,fbgd = ",fsam,fbgd 455 456 cor1 = fsam*sam_data/sam_AttenFactor - fbgd*bgd_temp/bgd_AttenFactor 457 cor1 *= noadd_bgd //zeros out regions where arrays do not overlap, one otherwise 458 459 // do the error propagation piecewise 460 Duplicate/O sam_err, tmp_a, tmp_b 461 tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2 //sig a ^2 462 463 tmp_b = (bgd_err/bgd_attenFactor)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2 //sig b ^2 464 465 cor_err = sqrt(tmp_a + tmp_b) 466 467 468 //we're done, get out w/no error 469 //set the COR_data to the result 470 cor_data = cor1 471 cor_data_display = cor1 472 473 //update COR header 474 cor_text[1] = date() + " " + time() //date + time stamp 475 476 KillWaves/Z cor1,bgd_temp,noadd_bgd 477 Killwaves/Z tmp_a,tmp_b 478 479 SetDataFolder root: 480 Return(0) 481 End 482 483 // empty subtraction only 484 // data does exist, checked by dispatch routine 485 // 486 xFunction V_CorrectMode_3() 487 //create the necessary wave references 488 WAVE sam_data=$"root:Packages:NIST:SAM:linear_data" 489 WAVE sam_reals=$"root:Packages:NIST:SAM:realsread" 490 WAVE sam_ints=$"root:Packages:NIST:SAM:integersread" 491 WAVE/T sam_text=$"root:Packages:NIST:SAM:textread" 492 WAVE emp_data=$"root:Packages:NIST:EMP:linear_data" 493 WAVE emp_reals=$"root:Packages:NIST:EMP:realsread" 494 WAVE emp_ints=$"root:Packages:NIST:EMP:integersread" 495 WAVE/T emp_text=$"root:Packages:NIST:EMP:textread" 496 WAVE cor_data=$"root:Packages:NIST:COR:linear_data" 497 WAVE/T cor_text=$"root:Packages:NIST:COR:textread" 498 499 // needed to propagate error 500 WAVE cor_data_display=$"root:Packages:NIST:COR:data" //just for the final copy 501 WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error" 502 WAVE emp_err =$"root:Packages:NIST:EMP:linear_data_error" 503 WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error" 504 505 Variable sam_trans_err,emp_trans_err 506 sam_trans_err = sam_reals[41] 507 emp_trans_err = emp_reals[41] 508 509 //get sam and bgd attenuation factors 510 String fileStr="" 511 Variable lambda,attenNo,sam_AttenFactor,emp_attenFactor 512 Variable tmonsam,fsam,femp,xshift,yshift,rsam,csam,remp,cemp,tmonemp 513 Variable wcen=0.001,tsam,temp 514 Variable sam_atten_err,emp_atten_err 515 fileStr = sam_text[3] 516 lambda = sam_reals[26] 517 attenNo = sam_reals[3] 518 sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err) 519 fileStr = emp_text[3] 520 lambda = emp_reals[26] 521 attenNo = emp_reals[3] 522 emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,emp_atten_err) 523 524 //get relative monitor counts (should all be 10^8, since normalized in add step) 525 tmonsam = sam_reals[0] //monitor count in SAM 526 tsam = sam_reals[4] //SAM transmission 527 csam = sam_reals[16] //x center 528 rsam = sam_reals[17] //beam (x,y) define center of corrected field 529 tmonemp = emp_reals[0] //monitor count in EMP 530 temp = emp_reals[4] //trans emp 531 cemp = emp_reals[16] //beamcenter of EMP 532 remp = emp_reals[17] 533 534 if(temp==0) 535 DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction" 536 temp=1 537 Endif 538 539 //Print "rbgd,cbgd = ",rbgd,cbgd 540 // set up beamcenter shift, relative to SAM 541 xshift = cemp-csam 542 yshift = remp-rsam 543 if(abs(xshift) <= wcen) 544 xshift = 0 545 Endif 546 if(abs(yshift) <= wcen) 547 yshift = 0 548 Endif 549 550 NVAR pixelsX = root:myGlobals:gNPixelsX 551 NVAR pixelsY = root:myGlobals:gNPixelsY 552 //get shifted data arrays, relative to SAM 553 Make/D/O/N=(pixelsX,pixelsY) cor1,emp_temp,noadd_emp //temp arrays 554 GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift) //emp_temp is the EMP 555 556 //do the sam-bgd subtraction, deposit result in cor1 557 fsam = 1 558 femp = tmonsam/tmonemp //this should be ==1 since normalized files 559 560 cor1 = fsam*sam_data/sam_AttenFactor - femp*(tsam/temp)*emp_temp/emp_AttenFactor 561 cor1 *= noadd_emp //zeros out regions where arrays do not overlap, one otherwise 562 563 // do the error propagation piecewise 564 Duplicate/O sam_err, tmp_a, tmp_c ,c_val 565 tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2 //sig a ^2 566 567 tmp_c = (sam_trans_err*emp_data/(temp*emp_attenFactor))^2 + (emp_err*tsam/(temp*emp_attenFactor))^2 568 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 569 570 cor_err = sqrt(tmp_a + tmp_c) 571 572 //we're done, get out w/no error 573 //set the COR data to the result 574 cor_data = cor1 575 cor_data_display = cor1 576 577 //update COR header 578 cor_text[1] = date() + " " + time() //date + time stamp 579 580 KillWaves/Z cor1,emp_temp,noadd_emp 581 Killwaves/Z tmp_a,tmp_c,c_val 582 583 SetDataFolder root: 584 Return(0) 585 End 586 587 // NO subtraction - simply rescales for attenuators 588 // SAM data does exist, checked by dispatch routine 589 // SAM data has already been copied to COR (both are the same at the start of the function) 590 // 591 // 592 // 593 Function V_CorrectMode_4() 594 595 //get SAM attenuation factor 596 Variable sam_AttenFactor,sam_atten_err,ii 597 String detStr 598 279 280 // these values apply to all of the detectors 599 281 sam_AttenFactor = V_getAttenuator_transmission("SAM") 600 282 sam_atten_err = V_getAttenuator_trans_err("SAM") 601 602 283 bgd_AttenFactor = V_getAttenuator_transmission("BGD") 284 bgd_atten_err = V_getAttenuator_trans_err("BGD") 285 emp_AttenFactor = V_getAttenuator_transmission("EMP") 286 emp_atten_err = V_getAttenuator_trans_err("EMP") 287 288 //get relative monitor counts (should all be 10^8, since normalized in add step) 289 // get transmission and trans error for SAM, EMP 290 // TODO -- verify that the call to V_getMonitorCount() is really rescaled to 10^8, and saved is the "true" count 291 292 tmonsam = V_getMonitorCount("SAM") //monitor count in SAM 293 tsam = V_getSampleTransmission("SAM") //SAM transmission 294 sam_trans_err = V_getSampleTransError("SAM") 295 296 tmonemp = V_getMonitorCount("EMP") //monitor count in EMP 297 temp = V_getSampleTransmission("EMP") //trans emp 298 emp_trans_err = V_getSampleTransError("EMP") 299 300 tmonbgd = V_getMonitorCount("BGD") //monitor count in BGD 301 302 303 // and now loop through all of the detectors 603 304 for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1) 604 305 detStr = StringFromList(ii, ksDetectorListAll, ";") … … 607 308 Wave sam_data = V_getDetectorDataW("SAM",detStr) 608 309 Wave sam_err = V_getDetectorDataErrW("SAM",detStr) 310 Wave bgd_data = V_getDetectorDataW("BGD",detStr) 311 Wave bgd_err = V_getDetectorDataErrW("BGD",detStr) 312 Wave emp_data = V_getDetectorDataW("EMP",detStr) 313 Wave emp_err = V_getDetectorDataErrW("EMP",detStr) 314 315 // to check for beam center mismatch -- simply warn, but do no shift 316 // 317 318 csam = V_getDet_beam_center_x("SAM",detStr) //x center 319 rsam = V_getDet_beam_center_y("SAM",detStr) //beam (x,y) define center of corrected field 320 321 cbgd = V_getDet_beam_center_x("BGD",detStr) 322 rbgd = V_getDet_beam_center_y("BGD",detStr) 323 324 cemp = V_getDet_beam_center_x("EMP",detStr) //beamcenter of EMP 325 remp = V_getDet_beam_center_y("EMP",detStr) 326 327 if(temp==0) 328 DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction" 329 temp=1 330 Endif 331 332 333 Duplicate/O cor_data cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp 334 335 // TODO -- document this, make a note, so everyone knows this is not done 336 // skip this part, but duplicate the results of no shift condition 337 // where bgd_temp = input data, and noadd_bgd = 1 (so no data is zeroed out) 338 339 //get the shifted data arrays, EMP and BGD, each relative to SAM 340 341 xshift = cbgd-csam 342 yshift = rbgd-rsam 343 if(abs(xshift) <= wcen) 344 xshift = 0 345 Endif 346 if(abs(yshift) <= wcen) 347 yshift = 0 348 Endif 349 // for the BGD file - alert if needed, generate dummy "pass-through" values 350 // 351 if(xshift != 0 || yshift != 0) 352 DoAlert 0, "Beam center mismatch for BGD file. Data has NOT been corrected." 353 endif 354 bgd_temp = bgd_data // no shift, no effect 355 noadd_bgd = 1 356 //GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift) //bgd_temp 357 358 xshift = cemp-csam 359 yshift = remp-rsam 360 if(abs(xshift) <= wcen) 361 xshift = 0 362 Endif 363 if(abs(yshift) <= wcen) 364 yshift = 0 365 Endif 366 // for the EMP file - alert if needed, generate dummy "pass-through" values 367 // 368 if(xshift != 0 || yshift != 0) 369 DoAlert 0, "Beam center mismatch for EMP file. Data has NOT been corrected." 370 endif 371 emp_temp = emp_data // no shift, no effect 372 noadd_emp = 1 373 //GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift) //emp_temp 374 375 376 // ******* 377 //do the subtraction 378 fsam=1 379 femp = tmonsam/tmonemp //this should be ==1 since normalized files 380 fbgd = tmonsam/tmonbgd //this should be ==1 since normalized files 381 cor1 = fsam*sam_data/sam_attenFactor - fbgd*bgd_temp/bgd_attenFactor 382 cor1 -= (tsam/temp)*(femp*emp_temp/emp_attenFactor - fbgd*bgd_temp/bgd_attenFactor) 383 cor1 *= noadd_bgd*noadd_emp //zero out the array mismatch values 384 385 // do the error propagation piecewise 386 Duplicate/O sam_err, tmp_a, tmp_b, tmp_c, tmp_d,c_val,d_val 387 tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2 //sig a ^2 388 389 tmp_b = (bgd_err/bgd_attenFactor)^2*(tsam/temp - 1)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2*(1-tsam/temp)^2 //sig b ^2 390 391 tmp_c = (sam_trans_err/temp)^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2 392 tmp_c += (tsam/temp^2)^2*emp_trans_err^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2 393 394 tmp_d = (tsam/(temp*emp_attenFactor))^2*(emp_err)^2 + (tsam*emp_data/(temp*emp_attenFactor^2))^2*(emp_atten_err)^2 395 396 cor_err = sqrt(tmp_a + tmp_b + tmp_c + tmp_d) 397 609 398 endfor 610 399 611 cor_data = sam_data/sam_AttenFactor //simply rescale the data 612 613 // do the error propagation piecewise 614 cor_err = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2 //sig a ^2 615 cor_err = sqrt(cor_err) 616 400 //we're done, get out w/no error 401 402 403 //TODO -- do I update COR header? 404 // cor_text[1] = date() + " " + time() //date + time stamp 405 406 KillWaves/Z cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp 407 Killwaves/Z tmp_a,tmp_b,tmp_c,tmp_d,c_val,d_val 408 409 SetDataFolder root: 410 Return(0) 411 End 412 413 //background only 414 // existence of data checked by dispatching routine 415 // data has already been copied to COR folder 416 Function V_CorrectMode_2() 417 418 //get SAM, BGD attenuation factor 419 Variable sam_AttenFactor,sam_atten_err,sam_trans_err 420 Variable bgd_AttenFactor,bgd_atten_err 421 Variable ii 422 String detStr 423 Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd 424 Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp 425 426 // these values apply to all of the detectors 427 sam_AttenFactor = V_getAttenuator_transmission("SAM") 428 sam_atten_err = V_getAttenuator_trans_err("SAM") 429 bgd_AttenFactor = V_getAttenuator_transmission("BGD") 430 bgd_atten_err = V_getAttenuator_trans_err("BGD") 431 432 //get relative monitor counts (should all be 10^8, since normalized in add step) 433 // get transmission and trans error for SAM, EMP 434 // TODO -- verify that the call to V_getMonitorCount() is really rescaled to 10^8, and saved is the "true" count 435 436 tmonsam = V_getMonitorCount("SAM") //monitor count in SAM 437 tsam = V_getSampleTransmission("SAM") //SAM transmission 438 sam_trans_err = V_getSampleTransError("SAM") 439 440 tmonbgd = V_getMonitorCount("BGD") //monitor count in BGD 441 442 443 // and now loop through all of the detectors 444 for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1) 445 detStr = StringFromList(ii, ksDetectorListAll, ";") 446 Wave cor_data = V_getDetectorDataW("COR",detStr) 447 Wave cor_err = V_getDetectorDataErrW("COR",detStr) 448 Wave sam_data = V_getDetectorDataW("SAM",detStr) 449 Wave sam_err = V_getDetectorDataErrW("SAM",detStr) 450 Wave bgd_data = V_getDetectorDataW("BGD",detStr) 451 Wave bgd_err = V_getDetectorDataErrW("BGD",detStr) 452 453 454 // to check for beam center mismatch -- simply warn, but do no shift 455 // 456 457 csam = V_getDet_beam_center_x("SAM",detStr) //x center 458 rsam = V_getDet_beam_center_y("SAM",detStr) //beam (x,y) define center of corrected field 459 460 cbgd = V_getDet_beam_center_x("BGD",detStr) 461 rbgd = V_getDet_beam_center_y("BGD",detStr) 462 463 464 Duplicate/O cor_data cor1,bgd_temp,noadd_bgd 465 466 // TODO -- document this, make a note, so everyone knows this is not done 467 // skip this part, but duplicate the results of no shift condition 468 // where bgd_temp = input data, and noadd_bgd = 1 (so no data is zeroed out) 469 470 //get the shifted data array BGD, relative to SAM 471 472 xshift = cbgd-csam 473 yshift = rbgd-rsam 474 if(abs(xshift) <= wcen) 475 xshift = 0 476 Endif 477 if(abs(yshift) <= wcen) 478 yshift = 0 479 Endif 480 // for the BGD file - alert if needed, generate dummy "pass-through" values 481 // 482 if(xshift != 0 || yshift != 0) 483 DoAlert 0, "Beam center mismatch for BGD file. Data has NOT been corrected." 484 endif 485 bgd_temp = bgd_data // no shift, no effect 486 noadd_bgd = 1 487 //GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift) //bgd_temp 488 489 490 // ********** 491 //do the sam-bgd subtraction, deposit result in cor1 492 fsam = 1 493 fbgd = tmonsam/tmonbgd //this should be ==1 since normalized files 494 495 //print "fsam,fbgd = ",fsam,fbgd 496 497 cor1 = fsam*sam_data/sam_AttenFactor - fbgd*bgd_temp/bgd_AttenFactor 498 cor1 *= noadd_bgd //zeros out regions where arrays do not overlap, one otherwise 499 500 // do the error propagation piecewise 501 Duplicate/O sam_err, tmp_a, tmp_b 502 tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2 //sig a ^2 503 504 tmp_b = (bgd_err/bgd_attenFactor)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2 //sig b ^2 505 506 cor_err = sqrt(tmp_a + tmp_b) 507 508 endfor 509 510 //we're done, get out w/no error 511 512 // TODO -- do I update COR header? 513 //cor_text[1] = date() + " " + time() //date + time stamp 514 515 KillWaves/Z cor1,bgd_temp,noadd_bgd 516 Killwaves/Z tmp_a,tmp_b 517 518 SetDataFolder root: 519 Return(0) 520 End 521 522 // empty subtraction only 523 // data does exist, checked by dispatch routine 524 // 525 Function V_CorrectMode_3() 526 527 //get SAM, EMP attenuation factor 528 Variable sam_AttenFactor,sam_atten_err,sam_trans_err 529 Variable emp_AttenFactor,emp_atten_err,emp_trans_err 530 Variable ii 531 String detStr 532 Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd 533 Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp 534 535 // these values apply to all of the detectors 536 sam_AttenFactor = V_getAttenuator_transmission("SAM") 537 sam_atten_err = V_getAttenuator_trans_err("SAM") 538 emp_AttenFactor = V_getAttenuator_transmission("EMP") 539 emp_atten_err = V_getAttenuator_trans_err("EMP") 540 541 //get relative monitor counts (should all be 10^8, since normalized in add step) 542 // get transmission and trans error for SAM, EMP 543 // TODO -- verify that the call to V_getMonitorCount() is really rescaled to 10^8, and saved is the "true" count 544 545 tmonsam = V_getMonitorCount("SAM") //monitor count in SAM 546 tsam = V_getSampleTransmission("SAM") //SAM transmission 547 sam_trans_err = V_getSampleTransError("SAM") 548 549 tmonemp = V_getMonitorCount("EMP") //monitor count in EMP 550 temp = V_getSampleTransmission("EMP") //trans emp 551 emp_trans_err = V_getSampleTransError("EMP") 552 553 554 // and now loop through all of the detectors 555 for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1) 556 detStr = StringFromList(ii, ksDetectorListAll, ";") 557 Wave cor_data = V_getDetectorDataW("COR",detStr) 558 Wave cor_err = V_getDetectorDataErrW("COR",detStr) 559 Wave sam_data = V_getDetectorDataW("SAM",detStr) 560 Wave sam_err = V_getDetectorDataErrW("SAM",detStr) 561 Wave emp_data = V_getDetectorDataW("EMP",detStr) 562 Wave emp_err = V_getDetectorDataErrW("EMP",detStr) 563 564 // to check for beam center mismatch -- simply warn, but do no shift 565 // 566 567 csam = V_getDet_beam_center_x("SAM",detStr) //x center 568 rsam = V_getDet_beam_center_y("SAM",detStr) //beam (x,y) define center of corrected field 569 570 cemp = V_getDet_beam_center_x("EMP",detStr) //beamcenter of EMP 571 remp = V_getDet_beam_center_y("EMP",detStr) 572 573 if(temp==0) 574 DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction" 575 temp=1 576 Endif 577 578 579 Duplicate/O cor_data cor1,emp_temp,noadd_emp 580 581 // TODO -- document this, make a note, so everyone knows this is not done 582 // skip this part, but duplicate the results of no shift condition 583 // where bgd_temp = input data, and noadd_bgd = 1 (so no data is zeroed out) 584 585 //get the shifted data array EMP, each relative to SAM 586 587 xshift = cemp-csam 588 yshift = remp-rsam 589 if(abs(xshift) <= wcen) 590 xshift = 0 591 Endif 592 if(abs(yshift) <= wcen) 593 yshift = 0 594 Endif 595 // for the EMP file - alert if needed, generate dummy "pass-through" values 596 // 597 if(xshift != 0 || yshift != 0) 598 DoAlert 0, "Beam center mismatch for EMP file. Data has NOT been corrected." 599 endif 600 emp_temp = emp_data // no shift, no effect 601 noadd_emp = 1 602 //GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift) //emp_temp 603 604 // ********** 605 606 //do the sam-bgd subtraction, deposit result in cor1 607 fsam = 1 608 femp = tmonsam/tmonemp //this should be ==1 since normalized files 609 610 cor1 = fsam*sam_data/sam_AttenFactor - femp*(tsam/temp)*emp_temp/emp_AttenFactor 611 cor1 *= noadd_emp //zeros out regions where arrays do not overlap, one otherwise 612 613 // do the error propagation piecewise 614 Duplicate/O sam_err, tmp_a, tmp_c ,c_val 615 tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2 //sig a ^2 616 617 tmp_c = (sam_trans_err*emp_data/(temp*emp_attenFactor))^2 + (emp_err*tsam/(temp*emp_attenFactor))^2 618 tmp_c += (tsam*emp_data*emp_trans_err/(temp*temp*emp_attenFactor))^2 + (tsam*emp_data*emp_atten_err/(temp*emp_attenFactor^2))^2//total of 6 terms 619 620 cor_err = sqrt(tmp_a + tmp_c) 621 622 endfor 623 624 //we're done, get out w/no error 625 626 // TODO -- do I update COR header? 627 //cor_text[1] = date() + " " + time() //date + time stamp 628 629 KillWaves/Z cor1,emp_temp,noadd_emp 630 Killwaves/Z tmp_a,tmp_c,c_val 631 632 SetDataFolder root: 633 Return(0) 634 End 635 636 // NO subtraction - simply rescales for attenuators 637 // SAM data does exist, checked by dispatch routine 638 // SAM data has already been copied to COR (both are the same at the start of the function) 639 // 640 // TODO -- do I need to rescale to sam_trans here ?? 641 // 642 // 643 Function V_CorrectMode_4() 644 645 //get SAM attenuation factor 646 Variable sam_AttenFactor,sam_atten_err,ii 647 String detStr 648 649 sam_AttenFactor = V_getAttenuator_transmission("SAM") 650 sam_atten_err = V_getAttenuator_trans_err("SAM") 651 652 653 for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1) 654 detStr = StringFromList(ii, ksDetectorListAll, ";") 655 Wave cor_data = V_getDetectorDataW("COR",detStr) 656 Wave cor_err = V_getDetectorDataErrW("COR",detStr) 657 Wave sam_data = V_getDetectorDataW("SAM",detStr) 658 Wave sam_err = V_getDetectorDataErrW("SAM",detStr) 659 660 661 cor_data = sam_data/sam_AttenFactor //simply rescale the data 662 663 // do the error propagation piecewise 664 cor_err = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2 //sig a ^2 665 cor_err = sqrt(cor_err) 666 667 endfor 668 617 669 //TODO -- do I want to update COR header? 618 670 // cor_text[1] = date() + " " + time() //date + time stamp … … 622 674 End 623 675 624 xFunction V_CorrectMode_11() 625 //create the necessary wave references 626 WAVE sam_data=$"root:Packages:NIST:SAM:linear_data" 627 WAVE sam_reals=$"root:Packages:NIST:SAM:realsread" 628 WAVE sam_ints=$"root:Packages:NIST:SAM:integersread" 629 WAVE/T sam_text=$"root:Packages:NIST:SAM:textread" 630 WAVE bgd_data=$"root:Packages:NIST:BGD:linear_data" 631 WAVE bgd_reals=$"root:Packages:NIST:BGD:realsread" 632 WAVE bgd_ints=$"root:Packages:NIST:BGD:integersread" 633 WAVE/T bgd_text=$"root:Packages:NIST:BGD:textread" 634 WAVE emp_data=$"root:Packages:NIST:EMP:linear_data" 635 WAVE emp_reals=$"root:Packages:NIST:EMP:realsread" 636 WAVE emp_ints=$"root:Packages:NIST:EMP:integersread" 637 WAVE/T emp_text=$"root:Packages:NIST:EMP:textread" 638 WAVE drk_data=$"root:Packages:NIST:DRK:linear_data" 639 WAVE drk_reals=$"root:Packages:NIST:DRK:realsread" 640 WAVE drk_ints=$"root:Packages:NIST:DRK:integersread" 641 WAVE/T drk_text=$"root:Packages:NIST:DRK:textread" 642 WAVE cor_data=$"root:Packages:NIST:COR:linear_data" 643 WAVE/T cor_text=$"root:Packages:NIST:COR:textread" 644 645 // needed to propagate error 646 WAVE cor_data_display=$"root:Packages:NIST:COR:data" //just for the final copy 647 WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error" 648 WAVE bgd_err =$"root:Packages:NIST:BGD:linear_data_error" 649 WAVE emp_err =$"root:Packages:NIST:EMP:linear_data_error" 650 WAVE drk_err =$"root:Packages:NIST:DRK:linear_data_error" 651 WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error" 652 653 Variable sam_trans_err,emp_trans_err 654 sam_trans_err = sam_reals[41] 655 emp_trans_err = emp_reals[41] 656 657 //get sam and bgd attenuation factors 658 String fileStr="" 659 Variable lambda,attenNo,sam_AttenFactor,bgd_attenFactor,emp_AttenFactor 676 677 678 Function V_CorrectMode_11() 679 680 //get SAM, BGD, EMP attenuation factor 681 Variable sam_AttenFactor,sam_atten_err,sam_trans_err 682 Variable bgd_AttenFactor,bgd_atten_err 683 Variable emp_AttenFactor,emp_atten_err,emp_trans_err 684 Variable ii 685 String detStr 660 686 Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd 661 Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp,time_sam,time_drk,savmon_sam 662 Variable sam_atten_err,bgd_atten_err,emp_atten_err 663 fileStr = sam_text[3] 664 lambda = sam_reals[26] 665 attenNo = sam_reals[3] 666 sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err) 667 fileStr = bgd_text[3] 668 lambda = bgd_reals[26] 669 attenNo = bgd_reals[3] 670 bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,bgd_atten_err) 671 fileStr = emp_text[3] 672 lambda = emp_reals[26] 673 attenNo = emp_reals[3] 674 emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,emp_atten_err) 675 687 Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp 688 Variable savmon_sam,time_sam,time_drk 689 690 // these values apply to all of the detectors 691 sam_AttenFactor = V_getAttenuator_transmission("SAM") 692 sam_atten_err = V_getAttenuator_trans_err("SAM") 693 bgd_AttenFactor = V_getAttenuator_transmission("BGD") 694 bgd_atten_err = V_getAttenuator_trans_err("BGD") 695 emp_AttenFactor = V_getAttenuator_transmission("EMP") 696 emp_atten_err = V_getAttenuator_trans_err("EMP") 697 676 698 //get relative monitor counts (should all be 10^8, since normalized in add step) 677 tmonsam = sam_reals[0] //monitor count in SAM 678 tsam = sam_reals[4] //SAM transmission 679 csam = sam_reals[16] //x center 680 rsam = sam_reals[17] //beam (x,y) define center of corrected field 681 tmonbgd = bgd_reals[0] //monitor count in BGD 682 cbgd = bgd_reals[16] 683 rbgd = bgd_reals[17] 684 tmonemp = emp_reals[0] //monitor count in EMP 685 temp = emp_reals[4] //trans emp 686 cemp = emp_reals[16] //beamcenter of EMP 687 remp = emp_reals[17] 688 savmon_sam=sam_reals[1] //true monitor count in SAM 689 time_sam = sam_ints[2] //count time SAM 690 time_drk = drk_ints[2] //drk count time 691 692 NVAR pixelsX = root:myGlobals:gNPixelsX 693 NVAR pixelsY = root:myGlobals:gNPixelsY 694 //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM 695 Make/D/O/N=(pixelsX,pixelsY) drk_temp, drk_tmp_err 696 drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam) 697 drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam) //temporarily rescale the error of DRK 698 699 if(temp==0) 700 DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction" 701 temp=1 702 Endif 703 704 //get the shifted data arrays, EMP and BGD, each relative to SAM 705 Make/D/O/N=(pixelsX,pixelsY) cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp 706 xshift = cbgd-csam 707 yshift = rbgd-rsam 708 if(abs(xshift) <= wcen) 709 xshift = 0 710 Endif 711 if(abs(yshift) <= wcen) 712 yshift = 0 713 Endif 714 GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift) //bgd_temp 715 716 xshift = cemp-csam 717 yshift = remp-rsam 718 if(abs(xshift) <= wcen) 719 xshift = 0 720 Endif 721 if(abs(yshift) <= wcen) 722 yshift = 0 723 Endif 724 GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift) //emp_temp 725 //always ignore the DRK center shift 726 727 //do the subtraction 728 fsam=1 729 femp = tmonsam/tmonemp //this should be ==1 since normalized files 730 fbgd = tmonsam/tmonbgd //this should be ==1 since normalized files 731 cor1 = fsam*sam_data/sam_attenFactor 732 cor1 -= (tsam/temp)*(femp*emp_temp/emp_attenFactor - fbgd*bgd_temp/bgd_attenFactor) 733 cor1 -= (fbgd*bgd_temp/bgd_attenFactor - drk_temp) 734 cor1 -= drk_temp/sam_attenFactor 735 cor1 *= noadd_bgd*noadd_emp //zero out the array mismatch values 736 737 // do the error propagation piecewise 738 Duplicate/O sam_err, tmp_a, tmp_b, tmp_c, tmp_d,c_val,d_val 739 tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2 //sig a ^2 740 741 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 742 743 tmp_c = (sam_trans_err/temp)^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2 744 tmp_c += (tsam/temp^2)^2*emp_trans_err^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2 745 746 tmp_d = (tsam/(temp*emp_attenFactor))^2*(emp_err)^2 + (tsam*emp_data/(temp*emp_attenFactor^2))^2*(emp_atten_err)^2 747 748 cor_err = sqrt(tmp_a + tmp_b + tmp_c + tmp_d + drk_tmp_err^2) 699 // get transmission and trans error for SAM, EMP 700 // TODO -- verify that the call to V_getMonitorCount() is really rescaled to 10^8, and saved is the "true" count 701 702 tmonsam = V_getMonitorCount("SAM") //monitor count in SAM 703 tsam = V_getSampleTransmission("SAM") //SAM transmission 704 sam_trans_err = V_getSampleTransError("SAM") 705 706 tmonemp = V_getMonitorCount("EMP") //monitor count in EMP 707 temp = V_getSampleTransmission("EMP") //trans emp 708 emp_trans_err = V_getSampleTransError("EMP") 709 710 tmonbgd = V_getMonitorCount("BGD") //monitor count in BGD 711 712 // for proper scaling, get the time and actual monitor counts 713 // TODO -- make sure that these calls are reading the proper values 714 savmon_sam = V_getBeamMonNormSaved_count("SAM") //true monitor count in SAM 715 time_sam = V_getCount_time("SAM") //count time SAM 716 time_drk = V_getCount_time("DRK") //drk count time 717 718 719 // and now loop through all of the detectors 720 for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1) 721 detStr = StringFromList(ii, ksDetectorListAll, ";") 722 Wave cor_data = V_getDetectorDataW("COR",detStr) 723 Wave cor_err = V_getDetectorDataErrW("COR",detStr) 724 Wave sam_data = V_getDetectorDataW("SAM",detStr) 725 Wave sam_err = V_getDetectorDataErrW("SAM",detStr) 726 Wave bgd_data = V_getDetectorDataW("BGD",detStr) 727 Wave bgd_err = V_getDetectorDataErrW("BGD",detStr) 728 Wave emp_data = V_getDetectorDataW("EMP",detStr) 729 Wave emp_err = V_getDetectorDataErrW("EMP",detStr) 730 Wave drk_data = V_getDetectorDataW("DRK",detStr) 731 Wave drk_err = V_getDetectorDataErrW("DRK",detStr) 732 733 // to check for beam center mismatch -- simply warn, but do no shift 734 // 735 736 csam = V_getDet_beam_center_x("SAM",detStr) //x center 737 rsam = V_getDet_beam_center_y("SAM",detStr) //beam (x,y) define center of corrected field 738 739 cbgd = V_getDet_beam_center_x("BGD",detStr) 740 rbgd = V_getDet_beam_center_y("BGD",detStr) 741 742 cemp = V_getDet_beam_center_x("EMP",detStr) //beamcenter of EMP 743 remp = V_getDet_beam_center_y("EMP",detStr) 744 745 //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM 746 Duplicate/O drk_data drk_temp, drk_tmp_err 747 drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam) 748 drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam) //temporarily rescale the error of DRK 749 750 if(temp==0) 751 DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction" 752 temp=1 753 Endif 754 755 Duplicate/O cor_data cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp 756 757 // TODO -- document this, make a note, so everyone knows this is not done 758 // skip this part, but duplicate the results of no shift condition 759 // where bgd_temp = input data, and noadd_bgd = 1 (so no data is zeroed out) 760 761 //get the shifted data arrays, EMP and BGD, each relative to SAM 762 763 xshift = cbgd-csam 764 yshift = rbgd-rsam 765 if(abs(xshift) <= wcen) 766 xshift = 0 767 Endif 768 if(abs(yshift) <= wcen) 769 yshift = 0 770 Endif 771 // for the BGD file - alert if needed, generate dummy "pass-through" values 772 // 773 if(xshift != 0 || yshift != 0) 774 DoAlert 0, "Beam center mismatch for BGD file. Data has NOT been corrected." 775 endif 776 bgd_temp = bgd_data // no shift, no effect 777 noadd_bgd = 1 778 //GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift) //bgd_temp 779 780 xshift = cemp-csam 781 yshift = remp-rsam 782 if(abs(xshift) <= wcen) 783 xshift = 0 784 Endif 785 if(abs(yshift) <= wcen) 786 yshift = 0 787 Endif 788 // for the EMP file - alert if needed, generate dummy "pass-through" values 789 // 790 if(xshift != 0 || yshift != 0) 791 DoAlert 0, "Beam center mismatch for EMP file. Data has NOT been corrected." 792 endif 793 emp_temp = emp_data // no shift, no effect 794 noadd_emp = 1 795 //GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift) //emp_temp 796 797 798 //always ignore the DRK center shift 799 800 // ************ 801 //do the subtraction 802 fsam=1 803 femp = tmonsam/tmonemp //this should be ==1 since normalized files 804 fbgd = tmonsam/tmonbgd //this should be ==1 since normalized files 805 cor1 = fsam*sam_data/sam_attenFactor 806 cor1 -= (tsam/temp)*(femp*emp_temp/emp_attenFactor - fbgd*bgd_temp/bgd_attenFactor) 807 cor1 -= (fbgd*bgd_temp/bgd_attenFactor - drk_temp) 808 cor1 -= drk_temp/sam_attenFactor 809 cor1 *= noadd_bgd*noadd_emp //zero out the array mismatch values 810 811 // do the error propagation piecewise 812 Duplicate/O sam_err, tmp_a, tmp_b, tmp_c, tmp_d,c_val,d_val 813 tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2 //sig a ^2 814 815 tmp_b = (bgd_err/bgd_attenFactor)^2*(tsam/temp - 1)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2*(1-tsam/temp)^2 //sig b ^2 816 817 tmp_c = (sam_trans_err/temp)^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2 818 tmp_c += (tsam/temp^2)^2*emp_trans_err^2*(emp_data/emp_attenFactor-bgd_data/bgd_attenFactor)^2 819 820 tmp_d = (tsam/(temp*emp_attenFactor))^2*(emp_err)^2 + (tsam*emp_data/(temp*emp_attenFactor^2))^2*(emp_atten_err)^2 821 822 cor_err = sqrt(tmp_a + tmp_b + tmp_c + tmp_d + drk_tmp_err^2) 823 824 endfor 749 825 750 826 //we're done, get out w/no error 751 //set the COR data to the result 752 cor_data = cor1 753 cor_data_display = cor1 754 755 //update COR header 756 cor_text[1] = date() + " " + time() //date + time stamp 827 828 829 //TODO -- do I update COR header? 830 // cor_text[1] = date() + " " + time() //date + time stamp 757 831 758 832 KillWaves/Z cor1,bgd_temp,noadd_bgd,emp_temp,noadd_emp,drk_temp … … 765 839 //bgd and drk subtraction 766 840 // 767 xFunction V_CorrectMode_12() 768 //create the necessary wave references 769 WAVE sam_data=$"root:Packages:NIST:SAM:linear_data" 770 WAVE sam_reals=$"root:Packages:NIST:SAM:realsread" 771 WAVE sam_ints=$"root:Packages:NIST:SAM:integersread" 772 WAVE/T sam_text=$"root:Packages:NIST:SAM:textread" 773 WAVE bgd_data=$"root:Packages:NIST:BGD:linear_data" 774 WAVE bgd_reals=$"root:Packages:NIST:BGD:realsread" 775 WAVE bgd_ints=$"root:Packages:NIST:BGD:integersread" 776 WAVE/T bgd_text=$"root:Packages:NIST:BGD:textread" 777 WAVE drk_data=$"root:Packages:NIST:DRK:linear_data" 778 WAVE drk_reals=$"root:Packages:NIST:DRK:realsread" 779 WAVE drk_ints=$"root:Packages:NIST:DRK:integersread" 780 WAVE/T drk_text=$"root:Packages:NIST:DRK:textread" 781 WAVE cor_data=$"root:Packages:NIST:COR:linear_data" 782 WAVE/T cor_text=$"root:Packages:NIST:COR:textread" 783 784 // needed to propagate error 785 WAVE cor_data_display=$"root:Packages:NIST:COR:data" //just for the final copy 786 WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error" 787 WAVE bgd_err =$"root:Packages:NIST:BGD:linear_data_error" 788 WAVE drk_err =$"root:Packages:NIST:DRK:linear_data_error" 789 WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error" 790 791 Variable sam_trans_err 792 sam_trans_err = sam_reals[41] 793 794 795 //get sam and bgd attenuation factors 796 String fileStr="" 797 Variable lambda,attenNo,sam_AttenFactor,bgd_attenFactor 841 Function V_CorrectMode_12() 842 843 //get SAM, BGD, EMP attenuation factor 844 Variable sam_AttenFactor,sam_atten_err,sam_trans_err 845 Variable bgd_AttenFactor,bgd_atten_err 846 Variable emp_AttenFactor,emp_atten_err,emp_trans_err 847 Variable ii 848 String detStr 798 849 Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd 799 Variable wcen=0.001,time_drk,time_sam,savmon_sam,tsam 800 Variable sam_atten_err,bgd_atten_err 801 fileStr = sam_text[3] 802 lambda = sam_reals[26] 803 attenNo = sam_reals[3] 804 sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err) 805 fileStr = bgd_text[3] 806 lambda = bgd_reals[26] 807 attenNo = bgd_reals[3] 808 bgd_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,bgd_atten_err) 809 850 Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp 851 Variable savmon_sam,time_sam,time_drk 852 853 // these values apply to all of the detectors 854 sam_AttenFactor = V_getAttenuator_transmission("SAM") 855 sam_atten_err = V_getAttenuator_trans_err("SAM") 856 bgd_AttenFactor = V_getAttenuator_transmission("BGD") 857 bgd_atten_err = V_getAttenuator_trans_err("BGD") 858 810 859 //get relative monitor counts (should all be 10^8, since normalized in add step) 811 tmonsam = sam_reals[0] //monitor count in SAM 812 tsam = sam_reals[4] //SAM transmission 813 csam = sam_reals[16] //x center 814 rsam = sam_reals[17] //beam (x,y) define center of corrected field 815 tmonbgd = bgd_reals[0] //monitor count in BGD 816 cbgd = bgd_reals[16] 817 rbgd = bgd_reals[17] 818 savmon_sam=sam_reals[1] //true monitor count in SAM 819 time_sam = sam_ints[2] //count time SAM 820 time_drk = drk_ints[2] //drk count time 821 822 NVAR pixelsX = root:myGlobals:gNPixelsX 823 NVAR pixelsY = root:myGlobals:gNPixelsY 824 //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM 825 Make/D/O/N=(pixelsX,pixelsY) drk_temp,drk_tmp_err 826 drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam) 827 drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam) //temporarily rescale the error of DRK 828 829 // set up beamcenter shift, relative to SAM 830 xshift = cbgd-csam 831 yshift = rbgd-rsam 832 if(abs(xshift) <= wcen) 833 xshift = 0 834 Endif 835 if(abs(yshift) <= wcen) 836 yshift = 0 837 Endif 838 //get shifted data arrays, relative to SAM 839 Make/D/O/N=(pixelsX,pixelsY) cor1,bgd_temp,noadd_bgd //temp arrays 840 GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift) //bgd_temp is the BGD 841 //always ignore the DRK center shift 842 843 //do the sam-bgd subtraction, deposit result in cor1 844 fsam = 1 845 fbgd = tmonsam/tmonbgd //this should be ==1 since normalized files 846 847 cor1 = fsam*sam_data/sam_AttenFactor + fbgd*tsam*bgd_temp/bgd_AttenFactor 848 cor1 += -1*(fbgd*bgd_temp/bgd_attenFactor - drk_temp) - drk_temp/sam_attenFactor 849 cor1 *= noadd_bgd //zeros out regions where arrays do not overlap, one otherwise 850 851 // do the error propagation piecewise 852 Duplicate/O sam_err, tmp_a, tmp_b 853 tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2 //sig a ^2 854 855 tmp_b = (bgd_err/bgd_attenFactor)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2 //sig b ^2 856 857 cor_err = sqrt(tmp_a + tmp_b + drk_tmp_err^2) 858 860 // get transmission and trans error for SAM, EMP 861 // TODO -- verify that the call to V_getMonitorCount() is really rescaled to 10^8, and saved is the "true" count 862 863 tmonsam = V_getMonitorCount("SAM") //monitor count in SAM 864 tsam = V_getSampleTransmission("SAM") //SAM transmission 865 sam_trans_err = V_getSampleTransError("SAM") 866 867 tmonbgd = V_getMonitorCount("BGD") //monitor count in BGD 868 869 // for proper scaling, get the time and actual monitor counts 870 // TODO -- make sure that these calls are reading the proper values 871 savmon_sam = V_getBeamMonNormSaved_count("SAM") //true monitor count in SAM 872 time_sam = V_getCount_time("SAM") //count time SAM 873 time_drk = V_getCount_time("DRK") //drk count time 874 875 876 // and now loop through all of the detectors 877 for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1) 878 detStr = StringFromList(ii, ksDetectorListAll, ";") 879 Wave cor_data = V_getDetectorDataW("COR",detStr) 880 Wave cor_err = V_getDetectorDataErrW("COR",detStr) 881 Wave sam_data = V_getDetectorDataW("SAM",detStr) 882 Wave sam_err = V_getDetectorDataErrW("SAM",detStr) 883 Wave bgd_data = V_getDetectorDataW("BGD",detStr) 884 Wave bgd_err = V_getDetectorDataErrW("BGD",detStr) 885 Wave drk_data = V_getDetectorDataW("DRK",detStr) 886 Wave drk_err = V_getDetectorDataErrW("DRK",detStr) 887 888 // to check for beam center mismatch -- simply warn, but do no shift 889 // 890 891 csam = V_getDet_beam_center_x("SAM",detStr) //x center 892 rsam = V_getDet_beam_center_y("SAM",detStr) //beam (x,y) define center of corrected field 893 894 cbgd = V_getDet_beam_center_x("BGD",detStr) 895 rbgd = V_getDet_beam_center_y("BGD",detStr) 896 897 //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM 898 Duplicate/O drk_data drk_temp, drk_tmp_err 899 drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam) 900 drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam) //temporarily rescale the error of DRK 901 902 if(temp==0) 903 DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction" 904 temp=1 905 Endif 906 907 Duplicate/O cor_data cor1,bgd_temp,noadd_bgd 908 909 // TODO -- document this, make a note, so everyone knows this is not done 910 // skip this part, but duplicate the results of no shift condition 911 // where bgd_temp = input data, and noadd_bgd = 1 (so no data is zeroed out) 912 913 //get the shifted data arrays, EMP and BGD, each relative to SAM 914 915 xshift = cbgd-csam 916 yshift = rbgd-rsam 917 if(abs(xshift) <= wcen) 918 xshift = 0 919 Endif 920 if(abs(yshift) <= wcen) 921 yshift = 0 922 Endif 923 // for the BGD file - alert if needed, generate dummy "pass-through" values 924 // 925 if(xshift != 0 || yshift != 0) 926 DoAlert 0, "Beam center mismatch for BGD file. Data has NOT been corrected." 927 endif 928 bgd_temp = bgd_data // no shift, no effect 929 noadd_bgd = 1 930 //GetShiftedArray(bgd_data,bgd_temp,noadd_bgd,xshift,yshift) //bgd_temp 931 932 933 //always ignore the DRK center shift 934 935 // ************ 936 //do the sam-bgd subtraction, deposit result in cor1 937 fsam = 1 938 fbgd = tmonsam/tmonbgd //this should be ==1 since normalized files 939 940 cor1 = fsam*sam_data/sam_AttenFactor + fbgd*tsam*bgd_temp/bgd_AttenFactor 941 cor1 += -1*(fbgd*bgd_temp/bgd_attenFactor - drk_temp) - drk_temp/sam_attenFactor 942 cor1 *= noadd_bgd //zeros out regions where arrays do not overlap, one otherwise 943 944 // do the error propagation piecewise 945 Duplicate/O sam_err, tmp_a, tmp_b 946 tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2 //sig a ^2 947 948 tmp_b = (bgd_err/bgd_attenFactor)^2 + (bgd_atten_err*bgd_data/bgd_attenFactor^2)^2 //sig b ^2 949 950 cor_err = sqrt(tmp_a + tmp_b + drk_tmp_err^2) 951 952 endfor 953 859 954 //we're done, get out w/no error 860 //set the COR_data to the result 861 cor_data = cor1 862 cor_data_display = cor1 863 864 //update COR header 865 cor_text[1] = date() + " " + time() //date + time stamp 955 956 // TODO -- do I update COR header? 957 // cor_text[1] = date() + " " + time() //date + time stamp 866 958 867 959 KillWaves/Z cor1,bgd_temp,noadd_bgd,drk_temp … … 876 968 //scale DRK by monitor count scaling factor and the ratio of couting times 877 969 //to place the DRK file on equal footing 878 xFunction V_CorrectMode_13() 879 //create the necessary wave references 880 WAVE sam_data=$"root:Packages:NIST:SAM:linear_data" 881 WAVE sam_reals=$"root:Packages:NIST:SAM:realsread" 882 WAVE sam_ints=$"root:Packages:NIST:SAM:integersread" 883 WAVE/T sam_text=$"root:Packages:NIST:SAM:textread" 884 WAVE emp_data=$"root:Packages:NIST:EMP:linear_data" 885 WAVE emp_reals=$"root:Packages:NIST:EMP:realsread" 886 WAVE emp_ints=$"root:Packages:NIST:EMP:integersread" 887 WAVE/T emp_text=$"root:Packages:NIST:EMP:textread" 888 WAVE drk_data=$"root:Packages:NIST:DRK:linear_data" 889 WAVE drk_reals=$"root:Packages:NIST:DRK:realsread" 890 WAVE drk_ints=$"root:Packages:NIST:DRK:integersread" 891 WAVE/T drk_text=$"root:Packages:NIST:DRK:textread" 892 WAVE cor_data=$"root:Packages:NIST:COR:linear_data" 893 WAVE/T cor_text=$"root:Packages:NIST:COR:textread" 894 895 // needed to propagate error 896 WAVE cor_data_display=$"root:Packages:NIST:COR:data" //just for the final copy 897 WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error" 898 WAVE emp_err =$"root:Packages:NIST:EMP:linear_data_error" 899 WAVE drk_err =$"root:Packages:NIST:DRK:linear_data_error" 900 WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error" 901 902 Variable sam_trans_err,emp_trans_err 903 sam_trans_err = sam_reals[41] 904 emp_trans_err = emp_reals[41] 905 906 //get sam and bgd attenuation factors (DRK irrelevant) 907 String fileStr="" 908 Variable lambda,attenNo,sam_AttenFactor,emp_attenFactor 909 Variable tmonsam,fsam,femp,xshift,yshift,rsam,csam,remp,cemp,tmonemp 910 Variable wcen=0.001,tsam,temp,savmon_sam,time_sam,time_drk 911 Variable sam_atten_err,emp_atten_err 912 fileStr = sam_text[3] 913 lambda = sam_reals[26] 914 attenNo = sam_reals[3] 915 sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err) 916 fileStr = emp_text[3] 917 lambda = emp_reals[26] 918 attenNo = emp_reals[3] 919 emp_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,emp_atten_err) 920 970 Function V_CorrectMode_13() 971 972 //get SAM, EMP attenuation factor 973 Variable sam_AttenFactor,sam_atten_err,sam_trans_err 974 Variable bgd_AttenFactor,bgd_atten_err 975 Variable emp_AttenFactor,emp_atten_err,emp_trans_err 976 Variable ii 977 String detStr 978 Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd 979 Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp 980 Variable savmon_sam,time_sam,time_drk 981 982 // these values apply to all of the detectors 983 sam_AttenFactor = V_getAttenuator_transmission("SAM") 984 sam_atten_err = V_getAttenuator_trans_err("SAM") 985 emp_AttenFactor = V_getAttenuator_transmission("EMP") 986 emp_atten_err = V_getAttenuator_trans_err("EMP") 987 921 988 //get relative monitor counts (should all be 10^8, since normalized in add step) 922 tmonsam = sam_reals[0] //monitor count in SAM 923 tsam = sam_reals[4] //SAM transmission 924 csam = sam_reals[16] //x center 925 rsam = sam_reals[17] //beam (x,y) define center of corrected field 926 tmonemp = emp_reals[0] //monitor count in EMP 927 temp = emp_reals[4] //trans emp 928 cemp = emp_reals[16] //beamcenter of EMP 929 remp = emp_reals[17] 930 savmon_sam=sam_reals[1] //true monitor count in SAM 931 time_sam = sam_ints[2] //count time SAM 932 time_drk = drk_ints[2] //drk count time 933 934 NVAR pixelsX = root:myGlobals:gNPixelsX 935 NVAR pixelsY = root:myGlobals:gNPixelsY 936 //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM 937 Make/D/O/N=(pixelsX,pixelsY) drk_temp,drk_tmp_err 938 drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam) 939 drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam) //temporarily rescale the error of DRK 940 941 942 if(temp==0) 943 DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction" 944 temp=1 945 Endif 946 947 //Print "rbgd,cbgd = ",rbgd,cbgd 948 // set up beamcenter shift, relative to SAM 949 xshift = cemp-csam 950 yshift = remp-rsam 951 if(abs(xshift) <= wcen) 952 xshift = 0 953 Endif 954 if(abs(yshift) <= wcen) 955 yshift = 0 956 Endif 957 //get shifted data arrays, relative to SAM 958 Make/D/O/N=(pixelsX,pixelsY) cor1,emp_temp,noadd_emp //temp arrays 959 GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift) //emp_temp is the EMP 960 //always ignore beamcenter shift for DRK 961 962 //do the sam-bgd subtraction, deposit result in cor1 963 fsam = 1 964 femp = tmonsam/tmonemp //this should be ==1 since normalized files 965 966 cor1 = fsam*sam_data/sam_AttenFactor - femp*(tsam/temp)*emp_temp/emp_AttenFactor 967 cor1 += drk_temp - drk_temp/sam_attenFactor 968 cor1 *= noadd_emp //zeros out regions where arrays do not overlap, one otherwise 969 970 // do the error propagation piecewise 971 Duplicate/O sam_err, tmp_a, tmp_c, c_val 972 tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2 //sig a ^2 973 974 tmp_c = (sam_trans_err*emp_data/(temp*emp_attenFactor))^2 + (emp_err*tsam/(temp*emp_attenFactor))^2 975 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 976 977 cor_err = sqrt(tmp_a + tmp_c + drk_tmp_err^2) 978 989 // get transmission and trans error for SAM, EMP 990 // TODO -- verify that the call to V_getMonitorCount() is really rescaled to 10^8, and saved is the "true" count 991 992 tmonsam = V_getMonitorCount("SAM") //monitor count in SAM 993 tsam = V_getSampleTransmission("SAM") //SAM transmission 994 sam_trans_err = V_getSampleTransError("SAM") 995 996 tmonemp = V_getMonitorCount("EMP") //monitor count in EMP 997 temp = V_getSampleTransmission("EMP") //trans emp 998 emp_trans_err = V_getSampleTransError("EMP") 999 1000 // for proper scaling, get the time and actual monitor counts 1001 // TODO -- make sure that these calls are reading the proper values 1002 savmon_sam = V_getBeamMonNormSaved_count("SAM") //true monitor count in SAM 1003 time_sam = V_getCount_time("SAM") //count time SAM 1004 time_drk = V_getCount_time("DRK") //drk count time 1005 1006 1007 // and now loop through all of the detectors 1008 for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1) 1009 detStr = StringFromList(ii, ksDetectorListAll, ";") 1010 Wave cor_data = V_getDetectorDataW("COR",detStr) 1011 Wave cor_err = V_getDetectorDataErrW("COR",detStr) 1012 Wave sam_data = V_getDetectorDataW("SAM",detStr) 1013 Wave sam_err = V_getDetectorDataErrW("SAM",detStr) 1014 Wave emp_data = V_getDetectorDataW("EMP",detStr) 1015 Wave emp_err = V_getDetectorDataErrW("EMP",detStr) 1016 Wave drk_data = V_getDetectorDataW("DRK",detStr) 1017 Wave drk_err = V_getDetectorDataErrW("DRK",detStr) 1018 1019 // to check for beam center mismatch -- simply warn, but do no shift 1020 // 1021 1022 csam = V_getDet_beam_center_x("SAM",detStr) //x center 1023 rsam = V_getDet_beam_center_y("SAM",detStr) //beam (x,y) define center of corrected field 1024 1025 cemp = V_getDet_beam_center_x("EMP",detStr) //beamcenter of EMP 1026 remp = V_getDet_beam_center_y("EMP",detStr) 1027 1028 //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM 1029 Duplicate/O drk_data drk_temp, drk_tmp_err 1030 drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam) 1031 drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam) //temporarily rescale the error of DRK 1032 1033 if(temp==0) 1034 DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction" 1035 temp=1 1036 Endif 1037 1038 Duplicate/O cor_data cor1,emp_temp,noadd_emp 1039 1040 // TODO -- document this, make a note, so everyone knows this is not done 1041 // skip this part, but duplicate the results of no shift condition 1042 // where bgd_temp = input data, and noadd_bgd = 1 (so no data is zeroed out) 1043 1044 //get the shifted data arrays, EMP , each relative to SAM 1045 1046 xshift = cemp-csam 1047 yshift = remp-rsam 1048 if(abs(xshift) <= wcen) 1049 xshift = 0 1050 Endif 1051 if(abs(yshift) <= wcen) 1052 yshift = 0 1053 Endif 1054 // for the EMP file - alert if needed, generate dummy "pass-through" values 1055 // 1056 if(xshift != 0 || yshift != 0) 1057 DoAlert 0, "Beam center mismatch for EMP file. Data has NOT been corrected." 1058 endif 1059 emp_temp = emp_data // no shift, no effect 1060 noadd_emp = 1 1061 //GetShiftedArray(emp_data,emp_temp,noadd_emp,xshift,yshift) //emp_temp 1062 1063 1064 //always ignore the DRK center shift 1065 1066 // *************** 1067 //do the sam-bgd subtraction, deposit result in cor1 1068 fsam = 1 1069 femp = tmonsam/tmonemp //this should be ==1 since normalized files 1070 1071 cor1 = fsam*sam_data/sam_AttenFactor - femp*(tsam/temp)*emp_temp/emp_AttenFactor 1072 cor1 += drk_temp - drk_temp/sam_attenFactor 1073 cor1 *= noadd_emp //zeros out regions where arrays do not overlap, one otherwise 1074 1075 // do the error propagation piecewise 1076 Duplicate/O sam_err, tmp_a, tmp_c, c_val 1077 tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2 //sig a ^2 1078 1079 tmp_c = (sam_trans_err*emp_data/(temp*emp_attenFactor))^2 + (emp_err*tsam/(temp*emp_attenFactor))^2 1080 tmp_c += (tsam*emp_data*emp_trans_err/(temp*temp*emp_attenFactor))^2 + (tsam*emp_data*emp_atten_err/(temp*emp_attenFactor^2))^2//total of 6 terms 1081 1082 cor_err = sqrt(tmp_a + tmp_c + drk_tmp_err^2) 1083 1084 endfor 1085 979 1086 //we're done, get out w/no error 980 //set the COR data to the result 981 cor_data = cor1 982 cor_data_display = cor1 983 984 //update COR header 985 cor_text[1] = date() + " " + time() //date + time stamp 1087 1088 // TODO -- do I update COR header? 1089 // cor_text[1] = date() + " " + time() //date + time stamp 986 1090 987 1091 KillWaves/Z cor1,emp_temp,noadd_emp,drk_temp … … 994 1098 // ONLY drk subtraction 995 1099 // 996 xFunction V_CorrectMode_14() 997 //create the necessary wave references 998 WAVE sam_data=$"root:Packages:NIST:SAM:linear_data" 999 WAVE sam_reals=$"root:Packages:NIST:SAM:realsread" 1000 WAVE sam_ints=$"root:Packages:NIST:SAM:integersread" 1001 WAVE/T sam_text=$"root:Packages:NIST:SAM:textread" 1002 WAVE drk_data=$"root:Packages:NIST:DRK:linear_data" 1003 WAVE drk_reals=$"root:Packages:NIST:DRK:realsread" 1004 WAVE drk_ints=$"root:Packages:NIST:DRK:integersread" 1005 WAVE/T drk_text=$"root:Packages:NIST:DRK:textread" 1006 WAVE cor_data=$"root:Packages:NIST:COR:linear_data" 1007 WAVE/T cor_text=$"root:Packages:NIST:COR:textread" 1008 1009 // needed to propagate error 1010 WAVE cor_data_display=$"root:Packages:NIST:COR:data" //just for the final copy 1011 WAVE sam_err =$"root:Packages:NIST:SAM:linear_data_error" 1012 WAVE drk_err =$"root:Packages:NIST:DRK:linear_data_error" 1013 WAVE cor_err =$"root:Packages:NIST:COR:linear_data_error" 1014 1015 Variable sam_trans_err 1016 sam_trans_err = sam_reals[41] 1017 1018 1019 //get sam and bgd attenuation factors 1020 String fileStr="" 1021 Variable lambda,attenNo,sam_AttenFactor,bgd_attenFactor 1100 Function V_CorrectMode_14() 1101 1102 //get SAM, EMP attenuation factor 1103 Variable sam_AttenFactor,sam_atten_err,sam_trans_err 1104 Variable bgd_AttenFactor,bgd_atten_err 1105 Variable emp_AttenFactor,emp_atten_err,emp_trans_err 1106 Variable ii 1107 String detStr 1022 1108 Variable tmonsam,fsam,fbgd,xshift,yshift,rsam,csam,rbgd,cbgd,tmonbgd 1023 Variable wcen=0.001,time_drk,time_sam,savmon_sam,tsam 1024 Variable sam_atten_err 1025 fileStr = sam_text[3] 1026 lambda = sam_reals[26] 1027 attenNo = sam_reals[3] 1028 sam_AttenFactor = AttenuationFactor(fileStr,lambda,AttenNo,sam_atten_err) 1029 1109 Variable wcen=0.001,tsam,temp,remp,cemp,tmonemp,femp 1110 Variable savmon_sam,time_sam,time_drk 1111 1112 // these values apply to all of the detectors 1113 sam_AttenFactor = V_getAttenuator_transmission("SAM") 1114 sam_atten_err = V_getAttenuator_trans_err("SAM") 1115 1116 1030 1117 //get relative monitor counts (should all be 10^8, since normalized in add step) 1031 tmonsam = sam_reals[0] //monitor count in SAM 1032 tsam = sam_reals[4] //SAM transmission 1033 csam = sam_reals[16] //x center 1034 rsam = sam_reals[17] //beam (x,y) define center of corrected field 1035 1036 savmon_sam=sam_reals[1] //true monitor count in SAM 1037 time_sam = sam_ints[2] //count time SAM 1038 time_drk = drk_ints[2] //drk count time 1039 1040 NVAR pixelsX = root:myGlobals:gNPixelsX 1041 NVAR pixelsY = root:myGlobals:gNPixelsY 1042 //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM 1043 Make/D/O/N=(pixelsX,pixelsY) 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 Make/D/O/N=(pixelsX,pixelsY) cor1 //temp arrays 1048 //always ignore the DRK center shift 1049 1050 //do the subtraction, deposit result in cor1 1051 fsam = 1 1052 fbgd = tmonsam/tmonbgd //this should be ==1 since normalized files 1053 1054 //correct sam for attenuators, and do the same to drk, since it was scaled to sam count time 1055 cor1 = fsam*sam_data/sam_AttenFactor - drk_temp/sam_attenFactor 1056 1057 // do the error propagation piecewise 1058 Duplicate/O sam_err, tmp_a 1059 tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2 //sig a ^2 1060 1061 cor_err = sqrt(tmp_a + drk_tmp_err^2) 1062 1118 // get transmission and trans error for SAM, EMP 1119 // TODO -- verify that the call to V_getMonitorCount() is really rescaled to 10^8, and saved is the "true" count 1120 1121 tmonsam = V_getMonitorCount("SAM") //monitor count in SAM 1122 tsam = V_getSampleTransmission("SAM") //SAM transmission 1123 sam_trans_err = V_getSampleTransError("SAM") 1124 1125 1126 // for proper scaling, get the time and actual monitor counts 1127 // TODO -- make sure that these calls are reading the proper values 1128 savmon_sam = V_getBeamMonNormSaved_count("SAM") //true monitor count in SAM 1129 time_sam = V_getCount_time("SAM") //count time SAM 1130 time_drk = V_getCount_time("DRK") //drk count time 1131 1132 1133 // and now loop through all of the detectors 1134 for(ii=0;ii<ItemsInList(ksDetectorListAll);ii+=1) 1135 detStr = StringFromList(ii, ksDetectorListAll, ";") 1136 Wave cor_data = V_getDetectorDataW("COR",detStr) 1137 Wave cor_err = V_getDetectorDataErrW("COR",detStr) 1138 Wave sam_data = V_getDetectorDataW("SAM",detStr) 1139 Wave sam_err = V_getDetectorDataErrW("SAM",detStr) 1140 Wave drk_data = V_getDetectorDataW("DRK",detStr) 1141 Wave drk_err = V_getDetectorDataErrW("DRK",detStr) 1142 1143 1144 //rescale drk to sam cnt time and then multiply by the same monitor scaling as SAM 1145 Duplicate/O drk_data drk_temp, drk_tmp_err 1146 drk_temp = drk_data*(time_sam/time_drk)*(tmonsam/savmon_sam) 1147 drk_tmp_err *= drk_err*(time_sam/time_drk)*(tmonsam/savmon_sam) //temporarily rescale the error of DRK 1148 1149 if(temp==0) 1150 DoAlert 0,"Empty Cell transmission was zero. It has been reset to one for the subtraction" 1151 temp=1 1152 Endif 1153 1154 Duplicate/O cor_data cor1 1155 1156 //always ignore the DRK center shift 1157 1158 // ************ 1159 //do the subtraction, deposit result in cor1 1160 fsam = 1 1161 fbgd = tmonsam/tmonbgd //this should be ==1 since normalized files 1162 1163 //correct sam for attenuators, and do the same to drk, since it was scaled to sam count time 1164 cor1 = fsam*sam_data/sam_AttenFactor - drk_temp/sam_attenFactor 1165 1166 // do the error propagation piecewise 1167 Duplicate/O sam_err, tmp_a 1168 tmp_a = (sam_err/sam_attenFactor)^2 + (sam_atten_err*sam_data/sam_attenFactor^2)^2 //sig a ^2 1169 1170 cor_err = sqrt(tmp_a + drk_tmp_err^2) 1171 1172 endfor 1173 1063 1174 //we're done, get out w/no error 1064 //set the COR_data to the result 1065 cor_data = cor1 1066 cor_data_display = cor1 1067 1068 //update COR header 1069 cor_text[1] = date() + " " + time() //date + time stamp 1175 1176 //TODO -- do I update COR header? 1177 // cor_text[1] = date() + " " + time() //date + time stamp 1070 1178 1071 1179 KillWaves/Z cor1,bgd_temp,noadd_bgd,drk_temp
Note: See TracChangeset
for help on using the changeset viewer.
