 Timestamp:
 May 20, 2016 4:28:56 PM (7 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_DetectorCorrections.ipf
r995 r999 370 370 if(cmpstr(orientation,"vertical")==0) 371 371 // this is vertical tube data dimensioned as (Ntubes,Npix) 372 pixSize = 8 372 pixSize = 8.4 //V_getDet_y_pixel_size(fname,detStr) 373 373 374 374 elseif(cmpstr(orientation,"horizontal")==0) … … 396 396 397 397 // 398 // TODO: 399 //  MUST VERIFY the definition of SDD and how (if) setback is written to the data files 400 //  currently I'm assuming that the SDD is the "nominal" value which is correct for the 401 // L/R panels, but is not correct for the T/B panels (must add in the setback) 402 // 403 // 404 // 398 405 // data_realDistX, Y must be previously generated from running NonLinearCorrection() 399 406 // … … 411 418 // get all of the geometry information 412 419 orientation = V_getDet_tubeOrientation(fname,detStr) 413 sdd = V_getDet_distance(fname,detStr) 420 // sdd = V_getDet_distance(fname,detStr) //[cm] 421 // sdd += V_getDet_TBSetback(fname,detStr)/10 // written in [mm], convert to [cm], returns 0 for L/R/B panels 422 423 sdd = V_getDet_ActualDistance(fname,detStr) //sdd derived, including setback [cm] 414 424 sdd/=100 // sdd reported in cm, pass in m 415 425 // this is the ctr in pixels … … 552 562 // to each pixel 553 563 // 554 // not actually used for any thing, but here for completeness if anyone asks564 // not actually used for any calculations, but here for completeness if anyone asks, or for 2D data export 555 565 // 556 566 // this properly accounts for qz … … 579 589 End 580 590 591 592 // 593 // TODO  VERIFY calculations 594 //  This is the actual solid angle per pixel, not a ratio vs. some "unit SA" 595 // Do I just correct for the different area vs. the "nominal" central area? 596 //  decide how to implement  either directly change the data values (as was done in the past) 597 // or use this as a weighting for when the data is binned to I(q). In the second method, 2D data 598 // would need this to be applied before exporting 599 //  do I keep a wave note indicating that this correction has been applied to the data 600 // so that it can be "unapplied"? 601 //  do I calculate theta from geometry directly, or get it from Q (Assuming it's present?) 602 // (probably just from geometry, since I need SDD and dx and dy values...) 603 // 604 // 605 Function SolidAngleCorrection(w,w_err,fname,detStr,destPath) 606 Wave w,w_err 607 String fname,detStr,destPath 608 609 Variable sdd,xCtr,yCtr,lambda 610 611 // get all of the geometry information 612 // orientation = V_getDet_tubeOrientation(fname,detStr) 613 sdd = V_getDet_ActualDistance(fname,detStr) 614 sdd/=100 // sdd in cm, pass in m 615 616 // this is ctr in mm 617 xCtr = V_getDet_beam_center_x_mm(fname,detStr) 618 yCtr = V_getDet_beam_center_y_mm(fname,detStr) 619 lambda = V_getWavelength(fname) 620 621 SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr) 622 623 Wave data_realDistX = data_realDistX 624 Wave data_realDistY = data_realDistY 625 626 Duplicate/O w solid_angle,tmp_theta,tmp_dist //in the current df 627 628 //// calculate the scattering angle 629 // dx = (distX  xctr) //delta x in mm 630 // dy = (distY  yctr) //delta y in mm 631 tmp_dist = sqrt((data_realDistX  xctr)^2 + (data_realDistY  yctr)^2) 632 633 tmp_dist /= 10 // convert mm to cm 634 sdd *=100 //convert to cm 635 636 tmp_theta = atan(tmp_dist/sdd) //this is two_theta, the scattering angle 637 638 Variable ii,jj,numx,numy,dx,dy 639 numx = DimSize(tmp_theta,0) 640 numy = DimSize(tmp_theta,1) 641 642 for(ii=0 ;ii<numx;ii+=1) 643 for(jj=0;jj<numy;jj+=1) 644 645 if(ii==0) //do a forward difference if ii==0 646 dx = (data_realDistX[ii+1][jj]  data_realDistX[ii][jj]) //delta x for the pixel 647 else 648 dx = (data_realDistX[ii][jj]  data_realDistX[ii1][jj]) //delta x for the pixel 649 endif 650 651 652 if(jj==0) 653 dy = (data_realDistY[ii][jj+1]  data_realDistY[ii][jj]) //delta y for the pixel 654 else 655 dy = (data_realDistY[ii][jj]  data_realDistY[ii][jj1]) //delta y for the pixel 656 endif 657 658 dx /= 10 659 dy /= 10 // convert mm to cm (since sdd is in cm) 660 solid_angle[ii][jj] = dx*dy //this is in cm^2 661 endfor 662 endfor 663 664 // to cover up any issues w/negative dx or dy 665 solid_angle = abs(solid_angle) 666 667 // solid_angle correction 668 // == dx*dy*cos^3/sdd^2 669 solid_angle *= (cos(tmp_theta))^3 670 solid_angle /= sdd^2 671 672 // Here it is! Apply the correction to the intensity (I divide  to get the counts per solid angle!!) 673 w /= solid_angle 674 675 676 // TODO: 677 // correctly apply the correction to the error wave (assume a perfect value?) 678 // w_err /= solid_angle //is this correct?? 679 680 // TODO  clean up after I'm satisfied computations are correct 681 // KillWaves/Z tmp_theta,tmp_dist 682 683 return(0) 684 end 581 685 582 686 … … 668 772 // so divide here to get the correct answer (5/22/08 SRK) 669 773 if(doEfficiency) 670 data[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)671 data_err[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)774 // data[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd) 775 // data_err[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd) 672 776 // solidAngle[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd) //testing only 673 777 endif … … 718 822 719 823 720 //distances passed in are in mm721 // dtdist is SDD722 // xd and yd are distances from the beam center to the current pixel723 //724 // TODO:725 //  DoAlert 0,"This has not yet been updated for VSANS"726 //727 Function DetEffCorr(lambda,dtdist,xd,yd)728 Variable lambda,dtdist,xd,yd729 730 DoAlert 0,"This has not yet been updated for VSANS"731 732 Variable theta,cosT,ff,stAl,stHe733 734 theta = atan( (sqrt(xd^2 + yd^2))/dtdist )735 cosT = cos(theta)736 737 stAl = 0.00967*lambda*0.8 //dimensionless, constants from JGB memo738 stHe = 0.146*lambda*2.5739 740 ff = exp(stAl/cosT)*(1exp(stHe/cosT)) / ( exp(stAl)*(1exp(stHe)) )741 742 return(ff)743 End744 824 745 825 // DIVIDE the intensity by this correction to get the right answer … … 1031 1111 1032 1112 Variable err 1033 err = DIVCorrection(type) //returns err = 1 if data doesn't exist in specified folders1113 err = V_DIVCorrection(type) //returns err = 1 if data doesn't exist in specified folders 1034 1114 1035 1115 if(err) 1036 Abort "error in DIVCorrection()"1116 Abort "error in V_DIVCorrection()" 1037 1117 endif 1038 1118 … … 1056 1136 // TODO: 1057 1137 //  DoAlert 0,"This has not yet been updated for VSANS" 1058 //  how is the error propagationhandled?1138 //  how is the error propagation handled? 1059 1139 // 1060 1140 //function will divide the contents of "workType" folder with the contents of … … 1062 1142 // all data is linear scale for the calculation 1063 1143 // 1064 Function DIVCorrection(data,data_err,detStr,workType)1144 Function V_DIVCorrection(data,data_err,detStr,workType) 1065 1145 Wave data,data_err 1066 1146 String detStr,workType … … 1071 1151 1072 1152 if(WaveExists(data) == 0) 1073 Print "The data wave does not exist in DIVCorrection()"1153 Print "The data wave does not exist in V_DIVCorrection()" 1074 1154 Return(1) //error condition 1075 1155 Endif … … 1080 1160 WAVE/Z div_data = $("root:Packages:NIST:VSANS:DIV:entry:instrument:detector_"+detStr+":data") 1081 1161 if(WaveExists(div_data) == 0) 1082 Print "The DIV wave does not exist in DIVCorrection()"1162 Print "The DIV wave does not exist in V_DIVCorrection()" 1083 1163 Return(1) //error condition 1084 1164 Endif
Note: See TracChangeset
for help on using the changeset viewer.