Changeset 710 for sans/Dev/trunk/NCNR_User_Procedures
- Timestamp:
- Jul 1, 2010 12:10:31 PM (13 years ago)
- Location:
- sans/Dev/trunk/NCNR_User_Procedures
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
sans/Dev/trunk/NCNR_User_Procedures/Common/NIST_XML_v40.ipf
r709 r710 1143 1143 Variable fileref 1144 1144 1145 Open/R fileref as filestr1145 Open/R/Z fileref as filestr 1146 1146 FReadLine fileref, line 1147 1147 Close fileref -
sans/Dev/trunk/NCNR_User_Procedures/Common/Packages/Invariant/Invariant_v40.ipf
r639 r710 374 374 endif 375 375 376 Variable yesGuinier=0,nume,inv 376 Variable yesGuinier=0,nume,inv,scale 377 377 // do the extrapolation of the correct type 378 378 ControlInfo/W=Invariant_Panel check_0 //the Guinier box … … 397 397 398 398 if(yesGuinier) 399 Make/O/D G_coef={1000,-1000} //input 399 // Make/O/D G_coef={1000,-1000} //input 400 WaveStats/M=1/Q/R=[0,(nbeg-1)] iw 401 scale = V_avg 402 Make/O/D G_coef={(scale),-1000} //input -- with a better guess for the overall scale factor, 400 403 FuncFit Guinier_Fit G_coef iw[0,(nbeg-1)] /I=1 /X=qw /W=sw /D 401 404 extr_lqi= Guinier_Fit(G_coef,extr_lqq) -
sans/Dev/trunk/NCNR_User_Procedures/Common/Packages/PlotManager/PlotUtilsMacro_v40.ipf
r698 r710 1095 1095 1096 1096 // SANS Reduction bits 1097 tmp = "ASStandardFunction;Ann_1D_Graph;Avg_1D_Graph;BStandardFunction;CStandardFunction;Draw_Plot1D;MyMat2XYZ;NewDirection;SANSModelAAO_MCproto;Monte_SANS_Threaded;Monte_SANS_NotThreaded;Monte_SANS_W1;Monte_SANS_W2;Monte_SANS_W3;Monte_SANS_W4;Monte_SANS;FractionReachingDetector;" 1098 list = RemoveFromList(tmp, list ,";") 1097 tmp = "ASStandardFunction;Ann_1D_Graph;Avg_1D_Graph;BStandardFunction;CStandardFunction;Draw_Plot1D;MyMat2XYZ;" 1098 list = RemoveFromList(tmp, list ,";") 1099 tmp = "NewDirection;SANSModelAAO_MCproto;Monte_SANS_Threaded;Monte_SANS_NotThreaded;Monte_SANS_W1;Monte_SANS_W2;Monte_SANS_W3;Monte_SANS_W4;Monte_SANS;FractionReachingDetector;TwoLevel_EC;SmearedTwoLevel_EC;" 1100 list = RemoveFromList(tmp, list ,";") 1101 1099 1102 1100 1103 // USANS Reduction bits -
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/HFIR_Utils.ipf
r699 r710 432 432 // -- better to check the physical location every time 433 433 // 434 // tol changed to 401 per Gernot's experience 6/24/10 435 // 434 436 Function isTransFile(fName) /// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 435 437 String fname 436 438 437 Variable beamtrap_1y=0,beamtrap_2y=0,beamtrap_3y=0,beamtrap_4y=0,tol=4 51439 Variable beamtrap_1y=0,beamtrap_2y=0,beamtrap_3y=0,beamtrap_4y=0,tol=401 438 440 //Check by key "transsmission" 439 441 // if (stringmatch( getIsTrans(fName),"True")>0) -
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf
r708 r710 1 1 #pragma rtGlobals=1 // Use modern global access method. 2 2 #pragma IgorVersion=6.1 3 3 4 4 5 // … … 18 19 // the RNG issue is really not worth the effort. multiple copies with different RNG is as good as I need. Plus, 19 20 // whatever XOP crashing was happining during threading is really unlikely to be from the RNG 21 // 22 // -- June 2010 - calls from different threads to the same RNG really seems to cause a crash. Probably as soon 23 // as the different threads try to call at the same time. Found this out by accident doing the 24 // wavelength spread. Each thread called ran3 at that point, and the crash came quickly. Went 25 // away immediately when I kept the ran calls consistent and isolated within threads. 20 26 // 21 27 // *** look into erand48() as the (pseudo) random number generator (it's a standard c-lib function, at least on unix) … … 77 83 // 78 84 79 // - to add --- 80 // -- wavelength distribution = another RNG to select the wavelength 81 // -- quartz windows (an empirical model?? or measure some real data - power Law + background) 82 // -- blocked beam (measure this too, and have some empirical model for this too - Broad Peak) 83 // 85 // --- TO ADD --- 86 // X- wavelength distribution = another RNG to select the wavelength 87 // ---- done Jun 2010, approximating the wavelength distribution as a Gaussian, based on the triangular 88 // FWHM. Wavelength distribution added to XOP too, and now very accurately matches the shape of the 1D 89 // simulation. 90 // 91 // 92 // X- quartz windows (an empirical model?? or measure some real data - power Law + background) 93 // X- blocked beam (measure this too, and have some empirical model for this too - Broad Peak) 94 // --- Done (mostly). quartz cell and blocked beam have been added empirically, giving the count rate and predicted 95 // scattering. Count time for the simulated scattering is the same as the sample. The simulated EC 96 // data can be plotted, but only by hand right now. EC and blocked beam are combined. 97 // 98 // -- divergence / size of the incoming beam. Currently everything is parallel, and anything that is transmitted 99 // simply ends up in (xCtr,yCtr), and the "real" profile of the beam is not captured. 84 100 85 101 … … 116 132 //OK 117 133 if(nthreads>4) //only support 4 processors until I can figure out how to properly thread the XOP and to loop it 134 //AND - just use 4 threads rather than the 8 (4 + 4 hyperthread?) my quad-core reports. 118 135 nthreads=4 119 136 endif … … 364 381 Variable NDoubleCoherent,NMultipleScatter,countIt,detEfficiency 365 382 Variable NMultipleCoherent,NCoherentEvents 366 383 Variable deltaLam,v1,v2,currWavelength,rsq,fac //for simulating wavelength distribution 384 385 // don't set to other than one here. Detector efficiency is handled outside, only passing the number of 386 // countable neutrons to any of the simulation functions (n=imon*eff) 367 387 detEfficiency = 1.0 //70% counting efficiency = 0.7 368 388 … … 378 398 sig_incoh = inputWave[9] 379 399 sig_sas = inputWave[10] 400 deltaLam = inputWave[11] 380 401 381 402 // SetRandomSeed 0.1 //to get a reproduceable sequence … … 462 483 while(rr>r1) 463 484 485 //pick the wavelength out of the wavelength spread, approximate as a gaussian 486 // from NR - pg 288. Needs random # from [0,1]. del is deltaLam/lam (as FWHM) and the 487 // 2.35 converts to a gaussian std dev. 488 do 489 v1=2.0*abs(enoise(1))-1.0 490 v2=2.0*abs(enoise(1))-1.0 491 rsq=v1*v1+v2*v2 492 while (rsq >= 1.0 || rsq == 0.0) 493 fac=sqrt(-2.0*log(rsq)/rsq) 494 495 // gset=v1*fac //technically, I'm throwing away one of the two values 496 497 currWavelength = (v2*fac)*deltaLam*wavelength/2.35 + wavelength 498 499 464 500 do //Scattering Loop, will exit when "done" == 1 465 501 // keep scattering multiple times until the neutron exits the sample … … 504 540 // so get it from the wave scaling instead 505 541 Q0 =left + binarysearchinterp(ran_dev,abs(enoise(1)))*delta 506 theta = Q0/2/Pi* wavelength //SAS approximation. 1% error at theta=30 deg (theta/2=15deg)542 theta = Q0/2/Pi*currWavelength //SAS approximation. 1% error at theta=30 deg (theta/2=15deg) 507 543 508 544 //Print "q0, theta = ",q0,theta … … 545 581 546 582 Theta_z = acos(Vz) // Angle WITH respect to z axis. 547 testQ = 2*pi*sin(theta_z)/ wavelength583 testQ = 2*pi*sin(theta_z)/currWavelength 548 584 549 585 // pick a random phi angle, and see if it lands on the detector … … 554 590 555 591 // is it on the detector? 556 FindPixel(testQ,testPhi, wavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel)592 FindPixel(testQ,testPhi,currWavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 557 593 558 594 if(xPixel != -1 && yPixel != -1) … … 1195 1231 1196 1232 /// called in SASCALC:ReCalculateInten() 1197 Function 1233 Function Simulate_2D_MC(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs) 1198 1234 String funcStr 1199 1235 WAVE aveint,qval,sigave,sigmaq,qbar,fsubs … … 1225 1261 1226 1262 // do the simulation here, or not 1227 Variable r1,xCtr,yCtr,sdd,pixSize,wavelength 1263 Variable r1,xCtr,yCtr,sdd,pixSize,wavelength,deltaLam 1228 1264 String coefStr,abortStr,str 1229 1265 … … 1234 1270 pixSize = rw[10]/10 // convert pix size in mm to cm 1235 1271 wavelength = rw[26] 1272 deltaLam = rw[27] 1236 1273 coefStr = MC_getFunctionCoef(funcStr) 1237 1274 … … 1261 1298 Make/O/D/N=5000 root:Packages:NIST:SAS:nt=0,root:Packages:NIST:SAS:j1=0,root:Packages:NIST:SAS:j2=0 1262 1299 Make/O/D/N=100 root:Packages:NIST:SAS:nn=0 1263 Make/O/D/N=1 1root:Packages:NIST:SAS:inputWave=01300 Make/O/D/N=15 root:Packages:NIST:SAS:inputWave=0 1264 1301 1265 1302 WAVE nt = root:Packages:NIST:SAS:nt … … 1280 1317 inputWave[9] = sig_incoh 1281 1318 inputWave[10] = sig_sas 1319 inputWave[11] = deltaLam 1320 // inputWave[] 12-14 are currently unused 1282 1321 1283 1322 linear_data = 0 //initialize … … 1328 1367 Monte_SANS_Threaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 1329 1368 else 1330 Monte_SANS _NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)1369 Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 1331 1370 endif 1332 1371 … … 1340 1379 endif 1341 1380 1342 Print "counts on detector, including transmitted = ",sum(linear_data,-inf,inf)1381 // Print "counts on detector, including transmitted = ",sum(linear_data,-inf,inf) 1343 1382 1344 1383 // linear_data[xCtr][yCtr] = 0 //snip out the transmitted spike … … 1381 1420 linear_data[0][2] = 1 1382 1421 linear_data[1][1] = 1 1383 //linear_data[2][2] = 11422 linear_data[2][2] = 1 1384 1423 linear_data[1][0] = 0 1385 //linear_data[2][1] = 01424 linear_data[2][1] = 0 1386 1425 linear_data[0][1] = 0 1387 // linear_data[1][2] = 0 1426 linear_data[1][2] = 0 1427 1428 linear_data[0][3] = 0 1429 linear_data[1][3] = 0 1430 linear_data[2][3] = 0 1431 linear_data[3][3] = 0 1432 linear_data[3][2] = 0 1433 linear_data[3][1] = 0 1434 linear_data[3][0] = 0 1435 1388 1436 1389 1437 data = linear_data … … 1394 1442 Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,"Simulation") 1395 1443 1444 // simulate the empty cell scattering, only in 1D 1445 Simulate_1D_EmptyCell("TwoLevel_EC",aveint,qval,sigave,sigmaq,qbar,fsubs) 1446 NVAR ctTime = root:Packages:NIST:SAS:gCntTime 1447 Print "Sample Simulation (2D) CR = ",results[9]/ctTime 1448 if(WinType("SANS_Data") ==1) 1449 Execute "ChangeDisplay(\"SAS\")" //equivalent to pressing "Show 2D" 1450 endif 1396 1451 1397 1452 return(0) … … 1870 1925 endif 1871 1926 1872 1927 1928 Simulate_1D_EmptyCell("TwoLevel_EC",aveint,qval,sigave,sigmaq,qbar,fsubs) 1929 Print "Sample Simulation (1D) CR = ",estDetCR 1930 1873 1931 return(0) 1874 1932 End … … 1905 1963 return(mScat) 1906 1964 End 1965 1966 1967 1968 // 1969 // -- empirical simulation of the scattering from an empty quartz cell + background (combined) 1970 // - there is little difference vs. the empty cell alone. 1971 // 1972 // - data was fit to the TwoLevel model, which fits rather nicely 1973 // 1974 Function Simulate_1D_EmptyCell(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs) 1975 String funcStr 1976 WAVE aveint,qval,sigave,sigmaq,qbar,fsubs 1977 1978 Variable r1,xCtr,yCtr,sdd,pixSize,wavelength 1979 String coefStr,abortStr,str 1980 1981 FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr) //a wrapper for the structure version 1982 FUNCREF SANSModelAAO_MCproto funcUnsmeared=$(funcStr) //unsmeared 1983 1984 Make/O/D root:Packages:NIST:SAS:coef_Empty = {1,1.84594,714.625,5e-08,2.63775,0.0223493,3.94009,0.0153754,1.72127,0} 1985 WAVE coefW = root:Packages:NIST:SAS:coef_Empty 1986 1987 Wave samInten=$"root:Simulation:Simulation_i" // this will exist and send the smeared calculation to the corect DF 1988 Duplicate samInten, root:Simulation:Simulation_EC_i 1989 Wave inten_EC=$"root:Simulation:Simulation_EC_i" 1990 1991 // the resolution-smeared intensity of the empty cell 1992 func(coefW,inten_EC,qval) 1993 1994 NVAR imon = root:Packages:NIST:SAS:gImon 1995 NVAR ctTime = root:Packages:NIST:SAS:gCntTime 1996 // NVAR thick = root:Packages:NIST:SAS:gThick 1997 NVAR trans = root:Packages:NIST:SAS:gSamTrans 1998 // NVAR SimDetCts = root:Packages:NIST:SAS:g_1DTotCts //summed counts (simulated) 1999 // NVAR estDetCR = root:Packages:NIST:SAS:g_1DEstDetCR // estimated detector count rate 2000 // NVAR fracScat = root:Packages:NIST:SAS:g_1DFracScatt // fraction of beam captured on detector 2001 // NVAR estTrans = root:Packages:NIST:SAS:g_1DEstTrans // estimated transmission of sample 2002 // NVAR mScat = root:Packages:NIST:SAS:g_MultScattFraction 2003 NVAR detectorEff = root:Packages:NIST:SAS:g_detectorEff 2004 2005 // use local variables here for the Empty cell - maybe use globals later, if I really want to save them 2006 // - here, just print them out for now 2007 Variable SimDetCts,estDetCR,fracScat,estTrans,mScat,thick 2008 2009 // for two 1/16" quartz windows, thick = 0.32 cm 2010 thick = 0.32 2011 2012 WAVE rw=root:Packages:NIST:SAS:realsRead 2013 WAVE nCells=root:Packages:NIST:SAS:nCells 2014 2015 pixSize = rw[10]/10 // convert pix size in mm to cm 2016 sdd = rw[18]*100 //convert header of [m] to [cm] 2017 wavelength = rw[26] // in 1/A 2018 2019 imon = beamIntensity()*ctTime 2020 2021 // calculate the scattering cross section simply to be able to estimate the transmission 2022 Variable sig_sas=0 2023 2024 // remember that the random deviate is the coherent portion ONLY - the incoherent background is 2025 // subtracted before the calculation. 2026 CalculateRandomDeviate(funcUnsmeared,coefW,wavelength,"root:Packages:NIST:SAS:ran_dev_EC",sig_sas) 2027 2028 // calculate the multiple scattering fraction for display (10/2009) 2029 Variable ii,nMax=10,tau 2030 mScat=0 2031 tau = thick*sig_sas 2032 // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered 2033 // neutrons out of those that were scattered 2034 for(ii=2;ii<nMax;ii+=1) 2035 mScat += tau^(ii)/factorial(ii) 2036 // print tau^(ii)/factorial(ii) 2037 endfor 2038 estTrans = exp(-1*thick*sig_sas) //thickness and sigma both in units of cm 2039 mscat *= (estTrans)/(1-estTrans) 2040 2041 2042 Duplicate/O qval prob_i_EC,countsInAnnulus_EC 2043 2044 prob_i_EC = trans*thick*(pixSize/sdd)^2*inten_EC //probability of a neutron in q-bin(i) 2045 2046 Variable P_on = sum(prob_i_EC,-inf,inf) 2047 // Print "P_on = ",P_on 2048 2049 fracScat = 1-estTrans 2050 2051 // added correction for detector efficiency, since SASCALC is flux on sample 2052 Duplicate/O aveint root:Packages:NIST:SAS:aveint_EC,root:Packages:NIST:SAS:sigave_EC 2053 WAVE aveint_EC = root:Packages:NIST:SAS:aveint_EC 2054 WAVE sigave_EC = root:Packages:NIST:SAS:sigave_EC 2055 aveint_EC = (Imon)*prob_i_EC*detectorEff 2056 2057 countsInAnnulus_EC = aveint_EC*nCells 2058 SimDetCts = sum(countsInAnnulus_EC,-inf,inf) 2059 estDetCR = SimDetCts/ctTime 2060 2061 // Print "Empty Cell Sig_sas = ",sig_sas 2062 Print "Empty Cell Count Rate : ",estDetCR 2063 2064 NVAR doABS = root:Packages:NIST:SAS:g_1D_DoABS 2065 NVAR addNoise = root:Packages:NIST:SAS:g_1D_AddNoise 2066 2067 // this is where the number of cells comes in - the calculation of the error bars 2068 // sigma[i] = SUM(sigma[ij]^2) / nCells^2 2069 // and since in the simulation, SUM(sigma[ij]^2) = nCells*sigma[ij]^2 = nCells*Inten 2070 // then... 2071 sigave_EC = sqrt(aveint_EC/nCells) // corrected based on John's memo, from 8/9/99 2072 2073 // add in random error in aveint based on the sigave 2074 if(addNoise) 2075 aveint_EC += gnoise(sigave_EC) 2076 endif 2077 2078 // signature in the standard deviation, do this after the noise is added 2079 // start at 10 to be out of the beamstop (makes for nicer plotting) 2080 // end at 50 to leave the natural statistics at the end of the set (may have a total of 80+ points if no offset) 2081 sigave_EC[10,50;10] = 10*sigave_EC[p] 2082 2083 // convert to absolute scale, remembering to un-correct for the detector efficiency 2084 if(doABS) 2085 Variable kappa = thick*(pixSize/sdd)^2*trans*iMon 2086 aveint_EC /= kappa 2087 sigave_EC /= kappa 2088 aveint_EC /= detectorEff 2089 sigave_EC /= detectorEff 2090 endif 2091 2092 2093 return(0) 2094 End 2095 2096 2097 // instead of including the Beaucage model in everything, keep a local copy here 2098 2099 //AAO version, uses XOP if available 2100 // simply calls the original single point calculation with 2101 // a wave assignment (this will behave nicely if given point ranges) 2102 Function TwoLevel_EC(cw,yw,xw) 2103 Wave cw,yw,xw 2104 2105 #if exists("TwoLevelX") 2106 yw = TwoLevelX(cw,xw) 2107 #else 2108 yw = fTwoLevel_EC(cw,xw) 2109 #endif 2110 return(0) 2111 End 2112 2113 Function fTwoLevel_EC(w,x) 2114 Wave w 2115 Variable x 2116 2117 Variable ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg 2118 Variable erf1,erf2,prec=1e-15,scale 2119 2120 //Rsub = Rs 2121 scale = w[0] 2122 G1 = w[1] //equivalent to I(0) 2123 Rg1 = w[2] 2124 B1 = w[3] 2125 Pow1 = w[4] 2126 G2 = w[5] 2127 Rg2 = w[6] 2128 B2 = w[7] 2129 Pow2 = w[8] 2130 bkg = w[9] 2131 2132 erf1 = erf( (x*Rg1/sqrt(6)) ,prec) 2133 erf2 = erf( (x*Rg2/sqrt(6)) ,prec) 2134 //Print erf1 2135 2136 ans = G1*exp(-x*x*Rg1*Rg1/3) 2137 ans += B1*exp(-x*x*Rg2*Rg2/3)*(erf1^3/x)^Pow1 2138 ans += G2*exp(-x*x*Rg2*Rg2/3) 2139 ans += B2*(erf2^3/x)^Pow2 2140 2141 if(x == 0) 2142 ans = G1 + G2 2143 endif 2144 2145 ans *= scale 2146 ans += bkg 2147 2148 Return(ans) 2149 End 2150 2151 2152 Function SmearedTwoLevel_EC(s) 2153 Struct ResSmearAAOStruct &s 2154 2155 // the name of your unsmeared model (AAO) is the first argument 2156 Smear_Model_20(TwoLevel_EC,s.coefW,s.xW,s.yW,s.resW) 2157 2158 return(0) 2159 End 2160 2161 //wrapper to calculate the smeared model as an AAO-Struct 2162 // fills the struct and calls the ususal function with the STRUCT parameter 2163 // 2164 // used only for the dependency, not for fitting 2165 // 2166 Function fSmearedTwoLevel_EC(coefW,yW,xW) 2167 Wave coefW,yW,xW 2168 2169 String str = getWavesDataFolder(yW,0) 2170 String DF="root:"+str+":" 2171 2172 WAVE resW = $(DF+str+"_res") 2173 2174 STRUCT ResSmearAAOStruct fs 2175 WAVE fs.coefW = coefW 2176 WAVE fs.yW = yW 2177 WAVE fs.xW = xW 2178 WAVE fs.resW = resW 2179 2180 Variable err 2181 err = SmearedTwoLevel_EC(fs) 2182 2183 return (0) 2184 End 2185 2186 1907 2187 1908 2188 -
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/SASCALC.ipf
r665 r710 745 745 // 746 746 // depending on the state of the 2D flag, open the 1d or 2d control panel 747 Function SimCheckProc( ctrlName,checked) : CheckBoxControl748 S tring ctrlName749 Variable checked 750 751 if( checked)747 Function SimCheckProc(CB_Struct) : CheckBoxControl 748 STRUCT WMCheckboxAction &CB_Struct 749 750 751 if(CB_Struct.checked) 752 752 NVAR do2D = root:Packages:NIST:SAS:gDoMonteCarlo 753 754 if(CB_Struct.eventMod == 2) //if the shift key is down - go to 2D mode 755 do2D = 1 756 endif 753 757 754 758 if(do2D)
Note: See TracChangeset
for help on using the changeset viewer.