 Timestamp:
 Jan 24, 2020 1:59:51 PM (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_DetectorCorrections.ipf
r1235 r1236 11 11 // "adjusted" or corrected data sets 12 12 // 13 14 //////////////// 15 // Constants for detector efficiency and shadowing 16 // 17 // V_TubeEfficiencyShadowCorr() 18 // 19 // JAN 2020 20 /////////////// 21 Constant kTube_ri = 0.372 // inner radius of tube [cm] 22 Constant kTube_cc = 0.84 // center to center spacing [cm] 23 Constant kTube_ss = 0.025 // stainless steel shell thickness [cm] 24 25 Constant kSig_2b_He = 0.146 // abs xs for 2 bar He(3) [cm1 A1] (multiply this by wavelength) 26 Constant kSig_8b_He = 0.593 // abs xs for 8 bar He(3) [cm1 A1] (multiply this by wavelength) 27 Constant kSig_Al = 0.00967 // abs xs for Al [cm1 A1] (multiply this by wavelength) 28 Constant kSig_ss = 0.146 // abs xs for 304 SS [cm1 A1] (multiply this by wavelength) 13 29 14 30 … … 1766 1782 End 1767 1783 1784 1785 1786 //////////////// 1787 // Detector efficiency and shadowing 1788 /////////////// 1789 1790 // 1791 // Tube efficiency + shadowing 1792 // 1793 // 1794 //  check for the existence of the proper tables (correct wavelength) 1795 //  generate tables if needed (onetime calculation) 1796 // 1797 // interpolate the table for the correction  to avoid repeated integration 1798 // 1799 // store the tables in: root:Packages:NIST:VSANS:Globals:Efficiency: 1800 // 1801 Function V_TubeEfficiencyShadowCorr(w,w_err,fname,detStr,destPath) 1802 Wave w,w_err 1803 String fname,detStr,destPath 1804 1805 Variable sdd,xCtr,yCtr,lambda 1806 String orientation 1807 1808 // if the panel is "B", exit  since it is not tubes, and this should not be called 1809 if(cmpstr(detStr,"B")==0) 1810 return(1) 1811 endif 1812 1813 // get all of the geometry information 1814 orientation = V_getDet_tubeOrientation(fname,detStr) 1815 sdd = V_getDet_ActualDistance(fname,detStr) 1816 1817 // this is ctr in mm 1818 xCtr = V_getDet_beam_center_x_mm(fname,detStr) 1819 yCtr = V_getDet_beam_center_y_mm(fname,detStr) 1820 lambda = V_getWavelength(fname) 1821 1822 SetDataFolder $(destPath + ":entry:instrument:detector_"+detStr) 1823 1824 Wave data_realDistX = data_realDistX 1825 Wave data_realDistY = data_realDistY 1826 1827 Duplicate/O w tmp_theta_x,tmp_theta_y,tmp_dist,tmp_corr //in the current df 1828 1829 //// calculate the scattering angles theta_x and theta_y 1830 1831 // flip the definitions of x and y for the T/B panels so that x is always lateral WRT the tubes 1832 // and y is always along the length of the tubes 1833 1834 if(cmpstr(orientation,"vertical")==0) 1835 // L/R panels, tube axis is ydirection 1836 // this is now a different tmp_dist 1837 // convert everything to cm first! 1838 // sdd is in [cm], everything else is in [mm] 1839 tmp_dist = (data_realDistY/10  yctr/10)/sqrt((data_realDistX/10  xctr/10)^2 + sdd^2) 1840 tmp_theta_y = atan(tmp_dist) //this is theta_y 1841 tmp_theta_x = atan( (data_realDistX/10  xctr/10)/sdd ) 1842 1843 else 1844 // horizontal orientation (T/B panels) 1845 // this is now a different tmp_dist 1846 // convert everything to cm first! 1847 // sdd is in [cm], everything else is in [mm] 1848 tmp_dist = (data_realDistX/10  xctr/10)/sqrt((data_realDistY/10  yctr/10)^2 + sdd^2) 1849 tmp_theta_y = atan(tmp_dist) //this is theta_y, along tube direction 1850 tmp_theta_x = atan( (data_realDistY/10  yctr/10)/sdd ) // this is laterally across tubes 1851 endif 1852 1853 1854 // identify if the 2D efficiency wave has been generated for the data wavelength 1855 // 1856 // if so, declare 1857 // if not, generate 1858 1859 if(WaveExists($"root:Packages:NIST:VSANS:Globals:Efficiency:eff") == 0) 1860 // generate the proper efficiency wave, at lambda 1861 NewDataFolder/O root:Packages:NIST:VSANS:Globals:Efficiency 1862 Print "recalculating efficiency table ..." 1863 V_TubeShadowEfficiencyTables_oneLam(lambda) 1864 // declare the wave 1865 Wave/Z effW = root:Packages:NIST:VSANS:Globals:Efficiency:eff 1866 else 1867 Wave/Z effW = root:Packages:NIST:VSANS:Globals:Efficiency:eff 1868 //is the efficiency at the correct wavelength? 1869 string str=note(effW) 1870 // Print "Note = ",str 1871 1872 if(V_CloseEnough(lambda,NumberByKey("LAMBDA", str,"="),0.1)) //absolute difference of < 0.1 A 1873 // yes, proceed, no need to do anything 1874 else 1875 // no, regenerate the efficiency and then proceed (wave already declared) 1876 Print "recalculating efficiency table ..." 1877 V_TubeShadowEfficiencyTables_oneLam(lambda) 1878 endif 1879 endif 1880 1881 1882 Variable ii,jj,numx,numy,xAngle,yAngle 1883 numx = DimSize(w,0) 1884 numy = DimSize(w,1) 1885 1886 // loop over all of the pixels of the panel and find the interpolated correction (save as a wave) 1887 // 1888 for(ii=0 ;ii<numx;ii+=1) 1889 for(jj=0;jj<numy;jj+=1) 1890 1891 // from the angles, find the (x,y) point to interpolate to get the efficiency 1892 1893 xAngle = tmp_theta_x[ii][jj] 1894 yAngle = tmp_theta_y[ii][jj] 1895 1896 xAngle = abs(xAngle) 1897 yAngle = abs(yAngle) 1898 1899 // the x and y scaling of the eff wave (2D) was set when it was generated (in radians) 1900 // simply reading the scaled xy value does not interpolate!! 1901 // tmp_corr[ii][jj] = effW(xAngle)(yAngle) // NO, returns "stepped" values 1902 tmp_corr[ii][jj] = Interp2D(effW,xAngle,yAngle) 1903 1904 endfor 1905 endfor 1906 // 1907 // 1908 // apply the correction and calculate the error 1909 // 1910 // Here it is! Apply the correction to the intensity (divide  to get the proper correction) 1911 w /= tmp_corr 1912 // 1913 // relative errors add in quadrature to the current 2D error 1914 // assume that this numerical calculation of efficiency is exact 1915 // 1916 // tmp_err = (w_err/tmp_corr)^2 + (lat_err/lat_corr)^2*w*w/lat_corr^2 1917 // tmp_err = sqrt(tmp_err) 1918 // 1919 // w_err = tmp_err 1920 // 1921 1922 // TODO 1923 //  clean up after I'm satisfied computations are correct 1924 // KillWaves/Z tmp_theta_x,tmp_theta_y,tmp_dist,tmp_err,tmp_corr 1925 1926 return(0) 1927 end 1928 1929 1930 1931 // the actual integration of the efficiency for an individual pixel 1932 Function V_Efficiency_Integral(pWave,in_u) 1933 Wave pWave 1934 Variable in_u 1935 1936 Variable lambda,th_x,th_y,u_p,integrand,T_sh,max_x,d_ss,d_He 1937 1938 lambda = pWave[0] 1939 th_x = pWave[1] 1940 th_y = pWave[2] 1941 1942 u_p = in_u + kTube_cc * cos(th_x) 1943 1944 // calculate shadow if th_x > 23.727 deg. th_x is input in radians 1945 max_x = 23.727 / 360 * 2*pi 1946 if(th_x < max_x) 1947 T_sh = 1 1948 else 1949 1950 // get d_ss 1951 if(abs(u_p) < kTube_ri) 1952 d_ss = sqrt( (kTube_ri + kTube_ss)^2  in_u^2 )  sqrt( kTube_ri^2  in_u^2 ) 1953 elseif (abs(u_p) < (kTube_ri + kTube_ss)) 1954 d_ss = sqrt( (kTube_ri + kTube_ss)^2  in_u^2 ) 1955 else 1956 d_ss = 0 1957 endif 1958 1959 // get d_He 1960 if(abs(u_p) < kTube_ri) 1961 d_He = 2 * sqrt( kTube_ri^2  in_u^2 ) 1962 else 1963 d_He = 0 1964 endif 1965 1966 //calculate T_sh 1967 T_sh = exp(2*kSig_ss*lambda*d_ss/cos(th_y)) * exp(kSig_8b_He*lambda*d_He/cos(th_y)) 1968 1969 endif 1970 1971 1972 // calculate the integrand 1973 1974 //note that the in_u value is used here to find d_ss and d_he (not u_p) 1975 // get d_ss 1976 if(abs(in_u) < kTube_ri) 1977 d_ss = sqrt( (kTube_ri + kTube_ss)^2  in_u^2 )  sqrt( kTube_ri^2  in_u^2 ) 1978 elseif (abs(in_u) < (kTube_ri + kTube_ss)) 1979 d_ss = sqrt( (kTube_ri + kTube_ss)^2  in_u^2 ) 1980 else 1981 d_ss = 0 1982 endif 1983 1984 // get d_He 1985 if(abs(in_u) < kTube_ri) 1986 d_He = 2 * sqrt( kTube_ri^2  in_u^2 ) 1987 else 1988 d_He = 0 1989 endif 1990 1991 integrand = T_sh*exp(kSig_ss*lambda*d_ss/cos(th_y))*( 1exp(kSig_8b_He*lambda*d_He/cos(th_y)) ) 1992 1993 return(integrand) 1994 end 1995 1996 // 1997 // Tube efficiency + shadowing 1998 // 1999 // function to generate the table for interpolation 2000 // 2001 // table is generated for a specific wavelength and normalized to eff(lam,0,0) 2002 // 2003 // below 24 deg (theta_x), there is no shadowing, so the table rows are all identical 2004 // 2005 // Only one table is stored, and the wavelength of that table is stored in the wave note 2006 //  detector correction checks the note, and recalculates the table if needed 2007 // (calculation takes approx 5 seconds) 2008 // 2009 Function V_TubeShadowEfficiencyTables_oneLam(lambda) 2010 Variable lambda 2011 2012 // storage location for tables 2013 SetDataFolder root:Packages:NIST:VSANS:Globals:Efficiency 2014 2015 //make waves that will be filed with the scattering angles and the result of the calculation 2016 // 2017 2018 //// fill arrays with the scattering angles theta_x and theta_y 2019 // 0 < x < 50 2020 // 0 < y < 50 2021 2022 // *** the definitions of x and y for the T/B panels is flipped so that x is always lateral WRT the tubes 2023 // and y is always along the length of the tubes 2024 2025 Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp,normVal 2026 numx = 25 2027 numy = numx 2028 2029 Make/O/D/N=(numx,numy) eff 2030 Make/O/D/N=(numx) theta_x, theta_y,eff_with_shadow,lam_cos 2031 2032 SetScale x 0,(numx*2)/360*2*pi,"", eff 2033 SetScale y 0,(numy*2)/360*2*pi,"", eff 2034 2035 Note/K eff // clear the note 2036 Note eff "LAMBDA="+num2str(lambda) 2037 2038 // theta_x = p*2 2039 theta_y = p *2 // value range from 0>45, changes if you change numx 2040 2041 //convert degrees to radians 2042 // theta_x = theta_x/360*2*pi 2043 theta_y = theta_y/360*2*pi 2044 2045 // Make/O/D/N=12 lam_wave 2046 // lam_wave = {0.5,0.7,1,1.5,2,3,4,6,8,10,15,20} 2047 2048 // Make/O/D/N=(12*numx) eff_withX_to_interp,lam_cos_theta_y 2049 // eff_withX_to_interp=0 2050 // lam_cos_theta_y=0 2051 2052 Make/O/D/N=3 pWave 2053 pWave[0] = lambda 2054 2055 2056 for(ii=0 ;ii<numx;ii+=1) 2057 2058 for(jj=0;jj<numx;jj+=1) 2059 2060 pWave[1] = indexToScale(eff,ii,0) //set theta x 2061 pWave[2] = indexToScale(eff,jj,1) //set theta y 2062 2063 eff_with_shadow[jj] = Integrate1D(V_Efficiency_Integral,kTube_ri,kTube_ri,2,0,pWave) // adaptive Gaussian quadrature 2064 eff_with_shadow[jj] /= (2*kTube_ri) 2065 2066 eff[ii][jj] = eff_with_shadow[jj] 2067 endfor 2068 2069 //eff[ii][] = eff_with_shadow[q] 2070 endfor 2071 2072 lam_cos = lambda/cos(theta_y) 2073 2074 Sort lam_cos,eff_with_shadow,lam_cos 2075 2076 // 2077 ////// // value for normalization at current wavelength 2078 pWave[0] = lambda 2079 pWave[1] = 0 2080 pWave[2] = 0 2081 //// 2082 normVal = Integrate1D(V_Efficiency_Integral,kTube_ri,kTube_ri,2,0,pWave) 2083 normVal /= (2*kTube_ri) 2084 // 2085 // print normVal 2086 // 2087 eff_with_shadow /= normVal // eff(lam,th_x,th_y) / eff(lam,0,0) 2088 2089 eff /= normVal 2090 2091 // TODO 2092 //  clean up after I'm satisfied computations are correct 2093 // KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err 2094 2095 SetDataFolder root: 2096 return(0) 2097 end 2098 2099 2100 2101 // 2102 // Tube efficiency + shadowing 2103 // 2104 // 2105 // TESTING function to generate the tables for interpolation 2106 // and various combinations of the corrections for plotting 2107 // 2108 Function V_TubeShadowEfficiencyTables_withX() 2109 2110 2111 Variable lambda 2112 lambda = 6 2113 2114 Variable theta_val=3 //the single theta_x value that is used 2115 2116 // TODO 2117 //  better storage location for tables 2118 // bad place for now... 2119 SetDataFolder root: 2120 2121 //make waves that will be filed with the scattering angles and the result of the calculation 2122 // 2123 2124 //// fill arrays with the scattering angles theta_x and theta_y 2125 // 0 < x < 50 2126 // 0 < y < 50 2127 2128 // *** the definitions of x and y for the T/B panels is flipped so that x is always lateral WRT the tubes 2129 // and y is always along the length of the tubes 2130 2131 Variable ii,jj,numx,numy,dx,dy,cos_th,arg,tmp,normVal 2132 numx = 10 2133 numy = 10 2134 2135 // Make/O/D/N=(numx,numy) eff 2136 Make/O/D/N=(numx) theta_x, theta_y,eff_with_shadow,lam_cos 2137 2138 theta_x = p*5 2139 theta_y = p*5 // value range from 0>45, changes if you change numx 2140 2141 //convert degrees to radians 2142 theta_x = theta_x/360*2*pi 2143 theta_y = theta_y/360*2*pi 2144 2145 Make/O/D/N=12 lam_wave 2146 lam_wave = {0.5,0.7,1,1.5,2,3,4,6,8,10,15,20} 2147 2148 Make/O/D/N=(12*numx) eff_withX_to_interp,lam_cos_theta_y 2149 eff_withX_to_interp=0 2150 lam_cos_theta_y=0 2151 2152 Make/O/D/N=3 pWave 2153 2154 for(jj=0;jj<12;jj+=1) 2155 2156 pWave[0] = lam_wave[jj] 2157 2158 for(ii=0 ;ii<numx;ii+=1) 2159 2160 pWave[1] = theta_val/360*2*pi //set theta x to any value 2161 pWave[2] = theta_y[ii] 2162 2163 eff_with_shadow[ii] = Integrate1D(V_Efficiency_Integral,kTube_ri,kTube_ri,2,0,pWave) // adaptive Gaussian quadrature 2164 eff_with_shadow[ii] /= (2*kTube_ri) 2165 2166 endfor 2167 2168 lam_cos = lam_wave[jj]/cos(theta_y) 2169 2170 // messy indexing for the concatentation 2171 lam_cos_theta_y[jj*numx,(jj+1)*numx1] = lam_cos[pjj*numx] 2172 eff_withX_to_interp[jj*numx,(jj+1)*numx1] = eff_with_shadow[pjj*numx] 2173 2174 endfor 2175 2176 Sort lam_cos_theta_y,eff_withX_to_interp,lam_cos_theta_y 2177 2178 // 2179 //////// // value for normalization at what wavelength??? 2180 // pWave[0] = 6 2181 // pWave[1] = 0 2182 // pWave[2] = 0 2183 ////// 2184 // normVal = Integrate1D(V_Efficiency_Integral,kTube_ri,kTube_ri,2,0,pWave) 2185 // normVal /= (2*kTube_ri) 2186 //// 2187 // print normVal 2188 //// 2189 // eff_withX_to_interp /= normVal // eff(lam,th_x,th_y) / eff(lam,0,0) 2190 2191 // TODO 2192 //  clean up after I'm satisfied computations are correct 2193 // KillWaves/Z tmp_theta,tmp_dist,tmp_err,lat_err 2194 2195 return(0) 2196 end 2197 2198 2199 2200 ///////////// 2201
Note: See TracChangeset
for help on using the changeset viewer.