Ignore:
Timestamp:
Oct 27, 2009 2:07:14 PM (13 years ago)
Author:
srkline
Message:

Added a calculation of the total fraction of multiple coherent scattering. This value is now displayed on the 1D sim panel and on the UCALC panel. If the total fraction of multiple coherent scattering (relative to the fraction coherently scattered) is greater than 0.10, then the field turns red as a warning.

Changed the numbering and suffix scheme for saving simulated VAX raw data files.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf

    r572 r583  
    907907        SetVariable cntVar,pos={185,73},size={90,15},proc=CountTimeSetVarProc,title="time(s)" 
    908908        SetVariable cntVar,format="%d" 
    909         SetVariable cntVar,limits={1,600,1},value= root:Packages:NIST:SAS:gCntTime 
     909        SetVariable cntVar,limits={1,3600,1},value= root:Packages:NIST:SAS:gCntTime 
    910910        Button MC_button2,pos={17,234},size={100,20},proc=SaveAsVAXButtonProc,title="Save 2D VAX" 
    911911        CheckBox check0,pos={216,180},size={68,14},title="Raw counts",variable = root:Packages:NIST:SAS:gRawCounts 
     
    10371037// will prompt for the sample label 
    10381038// 
     1039// currently hard-wired for SAS data folder 
     1040// 
    10391041Function SaveAsVAXButtonProc(ctrlName) : ButtonControl 
    10401042        String ctrlName 
    10411043 
    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) 
    10431077End 
    10441078 
     
    13061340        ValDisplay valdisp0_3,pos={326,121},size={220,13},title="Estimated transmission" 
    13071341        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 
    13081344        // set the flags here -- do the simulation, but not 2D 
    13091345         
     
    14251461                NVAR detCR = root:Packages:NIST:SAS:g_1DEstDetCR                // estimated detector count rate 
    14261462                NVAR fractScat = root:Packages:NIST:SAS:g_1DFracScatt 
     1463                NVAR mScat = root:Packages:NIST:SAS:g_MultScattFraction 
    14271464         
    14281465                SIMProtocol[0] = junk 
     
    14311468                SIMProtocol[3] = "\tDetector countrate (1/s) = "+num2str(detCR) 
    14321469                SIMProtocol[4] = "\tFraction of beam scattered coherently = "+num2str(fractScat) 
    1433                 SIMProtocol[5] = junk 
     1470                SIMProtocol[5] = "\tFraction of multiple coherent scattering = "+num2str(mScat) 
    14341471                SIMProtocol[6] = "" 
    14351472                SIMProtocol[7] = "" 
     
    15651602         
    15661603End 
     1604 
     1605 
     1606/// called in SASCALC:ReCalculateInten() 
     1607Function 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 resolution-smeared 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 well-behaved.",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)/(1-estTrans) 
     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 q-bin(i) that has nCells 
     1689        prob_i = trans*thick*(pixSize/sdd)^2*inten                      //probability of a neutron in q-bin(i)  
     1690         
     1691        Variable P_on = sum(prob_i,-inf,inf) 
     1692        Print "P_on = ",P_on 
     1693         
     1694//                              fracScat = P_on 
     1695        fracScat = 1-estTrans 
     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) 
     1733End 
     1734 
     1735/// for testing only 
     1736Function 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/(1-trans) 
     1748        p3p = tau^3/factorial(3)*trans/(1-trans) 
     1749        p4p = tau^4/factorial(4)*trans/(1-trans) 
     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)/(1-Trans) 
     1761 
     1762        Print "Total fraction of multiple scattering = ",mScat 
     1763 
     1764        return(mScat) 
     1765End 
Note: See TracChangeset for help on using the changeset viewer.