 Timestamp:
 Oct 27, 2009 2:07:14 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf
r572 r583 907 907 SetVariable cntVar,pos={185,73},size={90,15},proc=CountTimeSetVarProc,title="time(s)" 908 908 SetVariable cntVar,format="%d" 909 SetVariable cntVar,limits={1, 600,1},value= root:Packages:NIST:SAS:gCntTime909 SetVariable cntVar,limits={1,3600,1},value= root:Packages:NIST:SAS:gCntTime 910 910 Button MC_button2,pos={17,234},size={100,20},proc=SaveAsVAXButtonProc,title="Save 2D VAX" 911 911 CheckBox check0,pos={216,180},size={68,14},title="Raw counts",variable = root:Packages:NIST:SAS:gRawCounts … … 1037 1037 // will prompt for the sample label 1038 1038 // 1039 // currently hardwired for SAS data folder 1040 // 1039 1041 Function SaveAsVAXButtonProc(ctrlName) : ButtonControl 1040 1042 String ctrlName 1041 1043 1042 Write_RawData_File("SAS","",0) 1044 String fullpath="",destStr="" 1045 Variable refnum 1046 1047 fullpath = Write_RawData_File("SAS","",0) 1048 1049 // write out the results into a text file 1050 destStr = "root:Packages:NIST:SAS:" 1051 SetDataFolder $destStr 1052 1053 WAVE results=results 1054 WAVE/T results_desc=results_desc 1055 1056 //check each wave 1057 If(!(WaveExists(results))) 1058 Abort "results DNExist WriteVAXData()" 1059 Endif 1060 If(!(WaveExists(results_desc))) 1061 Abort "results_desc DNExist WriteVAXData()" 1062 Endif 1063 1064 Open refNum as fullpath+".txt" 1065 wfprintf refNum, "%30s\t\t%g\r",results_desc,results 1066 FStatus refNum 1067 FSetPos refNum,V_logEOF 1068 Close refNum 1069 1070 /////////////////////////////// 1071 1072 // could also automatically do the average here, but probably not worth the the effort... 1073 1074 SetDataFolder root: 1075 1076 return(0) 1043 1077 End 1044 1078 … … 1306 1340 ValDisplay valdisp0_3,pos={326,121},size={220,13},title="Estimated transmission" 1307 1341 ValDisplay valdisp0_3,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_1DEstTrans 1342 ValDisplay valdisp0_4,pos={326,145},size={220,13},title="Multiple Coherent Scattering" 1343 ValDisplay valdisp0_4,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_MultScattFraction 1308 1344 // set the flags here  do the simulation, but not 2D 1309 1345 … … 1425 1461 NVAR detCR = root:Packages:NIST:SAS:g_1DEstDetCR // estimated detector count rate 1426 1462 NVAR fractScat = root:Packages:NIST:SAS:g_1DFracScatt 1463 NVAR mScat = root:Packages:NIST:SAS:g_MultScattFraction 1427 1464 1428 1465 SIMProtocol[0] = junk … … 1431 1468 SIMProtocol[3] = "\tDetector countrate (1/s) = "+num2str(detCR) 1432 1469 SIMProtocol[4] = "\tFraction of beam scattered coherently = "+num2str(fractScat) 1433 SIMProtocol[5] = junk1470 SIMProtocol[5] = "\tFraction of multiple coherent scattering = "+num2str(mScat) 1434 1471 SIMProtocol[6] = "" 1435 1472 SIMProtocol[7] = "" … … 1565 1602 1566 1603 End 1604 1605 1606 /// called in SASCALC:ReCalculateInten() 1607 Function Simulate_1D(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs) 1608 String funcStr 1609 WAVE aveint,qval,sigave,sigmaq,qbar,fsubs 1610 1611 Variable r1,xCtr,yCtr,sdd,pixSize,wavelength 1612 String coefStr,abortStr,str 1613 1614 FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr) //a wrapper for the structure version 1615 FUNCREF SANSModelAAO_MCproto funcUnsmeared=$(funcStr) //unsmeared 1616 coefStr = MC_getFunctionCoef(funcStr) 1617 1618 if(!MC_CheckFunctionAndCoef(funcStr,coefStr)) 1619 Abort "Function and coefficients do not match. You must plot the unsmeared function before simulation." 1620 endif 1621 1622 Wave inten=$"root:Simulation:Simulation_i" // this will exist and send the smeared calculation to the corect DF 1623 1624 // the resolutionsmeared intensity is calculated, including the incoherent background 1625 func($coefStr,inten,qval) 1626 1627 NVAR imon = root:Packages:NIST:SAS:gImon 1628 NVAR ctTime = root:Packages:NIST:SAS:gCntTime 1629 NVAR thick = root:Packages:NIST:SAS:gThick 1630 NVAR trans = root:Packages:NIST:SAS:gSamTrans 1631 NVAR SimDetCts = root:Packages:NIST:SAS:g_1DTotCts //summed counts (simulated) 1632 NVAR estDetCR = root:Packages:NIST:SAS:g_1DEstDetCR // estimated detector count rate 1633 NVAR fracScat = root:Packages:NIST:SAS:g_1DFracScatt // fraction of beam captured on detector 1634 NVAR estTrans = root:Packages:NIST:SAS:g_1DEstTrans // estimated transmission of sample 1635 NVAR mScat = root:Packages:NIST:SAS:g_MultScattFraction 1636 1637 WAVE rw=root:Packages:NIST:SAS:realsRead 1638 WAVE nCells=root:Packages:NIST:SAS:nCells 1639 1640 pixSize = rw[10]/10 // convert pix size in mm to cm 1641 sdd = rw[18]*100 //convert header of [m] to [cm] 1642 wavelength = rw[26] // in 1/A 1643 1644 imon = beamIntensity()*ctTime 1645 1646 // calculate the scattering cross section simply to be able to estimate the transmission 1647 Variable sig_sas=0 1648 1649 // remember that the random deviate is the coherent portion ONLY  the incoherent background is 1650 // subtracted before the calculation. 1651 CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas) 1652 1653 // if(sig_sas > 100) 1654 // sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is wellbehaved.",sig_sas 1655 // endif 1656 1657 // calculate the multiple scattering fraction for display (10/2009) 1658 Variable ii,nMax=10,tau 1659 mScat=0 1660 tau = thick*sig_sas 1661 // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered 1662 // neutrons out of those that were scattered 1663 for(ii=2;ii<nMax;ii+=1) 1664 mScat += tau^(ii)/factorial(ii) 1665 // print tau^(ii)/factorial(ii) 1666 endfor 1667 estTrans = exp(1*thick*sig_sas) //thickness and sigma both in units of cm 1668 mscat *= (estTrans)/(1estTrans) 1669 1670 if(mScat > 0.1) // /Z to supress error if this was a 1D calc with the 2D panel open 1671 ValDisplay/Z valdisp0_4 win=Sim_1D_Panel,labelBack=(65535,32768,32768) 1672 else 1673 ValDisplay/Z valdisp0_4 win=Sim_1D_Panel,labelBack=0 1674 endif 1675 1676 1677 Print "Sig_sas = ",sig_sas 1678 1679 Duplicate/O qval prob_i,countsInAnnulus 1680 1681 // not needed  nCells takes care of this when the error is correctly calculated 1682 // Duplicate/O qval circle_fraction,rval,nCells_expected 1683 // rval = sdd*tan(2*asin(qval*wavelength/4/pi)) //radial distance in cm 1684 // nCells_expected = 2*pi*rval/pixSize //does this need to be an integer? 1685 // circle_fraction = nCells / nCells_expected 1686 1687 1688 // prob_i = trans*thick*nCells*(pixSize/sdd)^2*inten //probability of a neutron in qbin(i) that has nCells 1689 prob_i = trans*thick*(pixSize/sdd)^2*inten //probability of a neutron in qbin(i) 1690 1691 Variable P_on = sum(prob_i,inf,inf) 1692 Print "P_on = ",P_on 1693 1694 // fracScat = P_on 1695 fracScat = 1estTrans 1696 1697 // aveint = (Imon)*prob_i / circle_fraction / nCells_expected 1698 aveint = (Imon)*prob_i 1699 1700 countsInAnnulus = aveint*nCells 1701 SimDetCts = sum(countsInAnnulus,inf,inf) 1702 estDetCR = SimDetCts/ctTime 1703 1704 1705 NVAR doABS = root:Packages:NIST:SAS:g_1D_DoABS 1706 NVAR addNoise = root:Packages:NIST:SAS:g_1D_AddNoise 1707 1708 // this is where the number of cells comes in  the calculation of the error bars 1709 // sigma[i] = SUM(sigma[ij]^2) / nCells^2 1710 // and since in the simulation, SUM(sigma[ij]^2) = nCells*sigma[ij]^2 = nCells*Inten 1711 // then... 1712 sigave = sqrt(aveint/nCells) // corrected based on John's memo, from 8/9/99 1713 1714 // add in random error in aveint based on the sigave 1715 if(addNoise) 1716 aveint += gnoise(sigave) 1717 endif 1718 1719 // signature in the standard deviation, do this after the noise is added 1720 // start at 10 to be out of the beamstop (makes for nicer plotting) 1721 // end at 50 to leave the natural statistics at the end of the set (may have a total of 80+ points if no offset) 1722 sigave[10,50;10] = 10*sigave[p] 1723 1724 // convert to absolute scale 1725 if(doABS) 1726 Variable kappa = thick*(pixSize/sdd)^2*trans*iMon 1727 aveint /= kappa 1728 sigave /= kappa 1729 endif 1730 1731 1732 return(0) 1733 End 1734 1735 /// for testing only 1736 Function testProbability(sas,thick) 1737 Variable sas,thick 1738 1739 Variable tau,trans,p2p,p3p,p4p 1740 1741 tau = sas*thick 1742 trans = exp(tau) 1743 1744 Print "tau = ",tau 1745 Print "trans = ",trans 1746 1747 p2p = tau^2/factorial(2)*trans/(1trans) 1748 p3p = tau^3/factorial(3)*trans/(1trans) 1749 p4p = tau^4/factorial(4)*trans/(1trans) 1750 1751 Print "double scattering = ",p2p 1752 Print "triple scattering = ",p3p 1753 Print "quadruple scattering = ",p4p 1754 1755 Variable ii,nMax=10,mScat=0 1756 for(ii=2;ii<nMax;ii+=1) 1757 mScat += tau^(ii)/factorial(ii) 1758 // print tau^(ii)/factorial(ii) 1759 endfor 1760 mscat *= (Trans)/(1Trans) 1761 1762 Print "Total fraction of multiple scattering = ",mScat 1763 1764 return(mScat) 1765 End
Note: See TracChangeset
for help on using the changeset viewer.