Ignore:
Timestamp:
Aug 13, 2009 10:37:11 AM (13 years ago)
Author:
srkline
Message:

more changes to UCALC, display vs. countrate and add a display of experimental empty measurements along with the simulation. empty beam and empty call measurements are so similar - they have not been differentiated, especially since they are only for comparison.

Location:
sans/Dev/trunk/NCNR_User_Procedures/Reduction/USANS
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/USANS/USANS_Includes.ipf

    r546 r547  
    2323// USANS simulation and required procedures 
    2424#include "U_CALC" 
     25#include "USANS_EmptyWaves" 
    2526#include "MultScatter_MonteCarlo_2D" 
    2627#include "SASCALC" 
  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/USANS/U_CALC.ipf

    r546 r547  
    55 
    66// to simulate the intensity from a USANS experiment for planning 
    7  
     7// see John's instrument paper: 
     8// 
     9// J. Appl. Cryst. (2005) 38 1004-1011. 
     10// 
    811// SRK JUL 2009 
    912 
     
    1518// - fast ways to increase/decrease counting time 
    1619// 
    17 // - default model of q-4 power law or spheres 
     20// - plot as countrate, not absolute scale 
     21// - 3e-5 cutoff 
     22// - ? don't plot lowest angle range (but needs to be in the count time) 
     23// - need empty beam and empty cell count rate vs. aperture (Cd vs. Gd?) 
    1824// 
    1925 
     
    2430// model the direct beam?? currently the "red" region from -1 to 0.6 is almost entirely 
    2531// the primary beam, so it's a bit artificial (qmin is really ~ 3e-5) 
     32// 
     33 
    2634// 
    2735// Need T_wide, T_rock, I peak for proper absolute scaling 
     
    104112        Variable/G gSamTrans=0.8 
    105113        Variable/G g_1D_DoABS = 0               //=1 for abs scale, 0 for just counts 
     114        Variable/G g_1D_PlotCR = 1              //=1 to plot countrate 
    106115        Variable/G g_1D_AddNoise = 1            // add in appropriate noise to simulation 
    107116         
     
    125134// so that the subwindow syntax doesn't break all of the other functionality 
    126135// 
    127  
    128136Window UCALC_Panel() : Graph 
    129137        PauseUpdate; Silent 1           // building window... 
     
    132140        DoWindow/C UCALC 
    133141        DoWindow/T UCALC,"USANS Simulation" 
    134 //      ShowTools/A 
    135142        ControlBar 320 
    136143         
     
    140147         
    141148        PopupMenu popup0,pos={17,18},size={165,20},title="Sample Aperture Diam (in)" 
    142         PopupMenu popup0,mode=2,popvalue="0.625",value="0.5;0.625;0.75;1.0;1.75;" 
    143 //      PopupMenu popup0,proc=Sim_USANS_SamplAperPopMenuProc 
     149        PopupMenu popup0,mode=3,popvalue="0.625",value="0.25;0.50;0.625;0.75;1.0;1.75;2.0;" 
    144150        PopupMenu popup2,pos={220,18},size={165,20},title="Presets" 
    145151        PopupMenu popup2,mode=3,popvalue="Long Count",value="Short Count;Medium Count;Long Count;" 
     
    233239        SetVariable setvar6e proc=CtTimeSetVarProc,limits={-1,50000,100} 
    234240         
    235         Button button0,pos={260,213},size={50,20},proc=U_SimPlotButtonProc,title="Plot" 
     241        Button button0,pos={255,180},size={60,20},fColor=(65535,65535,0),proc=U_SimPlotButtonProc,title="Plot" 
    236242        Button button1,pos={260,286},size={50,20},proc=U_SaveButtonProc,title="Save" 
    237243 
     
    252258        SetVariable setvar0_1,limits={0,10,0.01},value= root:Packages:NIST:USANS:Globals:U_Sim:g_BkgLevel 
    253259         
    254         CheckBox check0_2,pos={253,239},size={60,14},title="Abs scale?",variable= root:Packages:NIST:USANS:Globals:U_Sim:g_1D_DoABS 
     260        CheckBox check0_4 title="Show EMP?",pos={160,260},proc=ShowEMPCheckProc,value=0 
     261         
     262        CheckBox check0_2,pos={253,239},size={60,14},title="CountRate?",variable= root:Packages:NIST:USANS:Globals:U_Sim:g_1D_PlotCR 
    255263        CheckBox check0_3,pos={262,264},size={60,14},title="Noise?",variable= root:Packages:NIST:USANS:Globals:U_Sim:g_1D_AddNoise 
    256264         
    257265// a box for the results 
    258266        SetVariable totalTime,pos={338,185},size={150,15},title="Count time (h:m)",value= gTotTimeStr 
    259         ValDisplay valdisp0,pos={338,210},size={220,13},title="Total detector counts" 
    260         ValDisplay valdisp0,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:USANS:Globals:U_Sim:g_1DTotCts 
     267//      ValDisplay valdisp0,pos={338,210},size={220,13},title="Total detector counts" 
     268//      ValDisplay valdisp0,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:USANS:Globals:U_Sim:g_1DTotCts 
    261269        ValDisplay valdisp0_2,pos={338,234},size={220,13},title="Fraction of beam scattered" 
    262270        ValDisplay valdisp0_2,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:USANS:Globals:U_Sim:g_1DFracScatt 
     
    468476 
    469477 
     478Function ShowEMPCheckProc(cba) : CheckBoxControl 
     479        STRUCT WMCheckboxAction &cba 
     480 
     481        switch( cba.eventCode ) 
     482                case 2: // mouse up 
     483                        Variable checked = cba.checked 
     484 
     485                        String list,item,popStr,qval,CR 
     486                        Variable OK=1 
     487 
     488                        if(exists("root:Packages:NIST:USANS:Globals:Q_2p0")==0)                 // 
     489                                MakeUSANSEmptyWaves() 
     490                        endif 
     491                         
     492                         
     493                        if(checked) 
     494                                // put it on the graph 
     495                                SetDataFolder root:Packages:NIST:USANS:Globals 
     496                                ControlInfo/W=UCALC popup0 
     497                                popStr = S_Value 
     498                                strswitch(popStr)       // string switch 
     499                                        case "0.25":            // execute if case matches expression 
     500                                                qval = "Q_0p25" 
     501                                                CR = "CR_0p25" 
     502                                                break            
     503                                        case "0.50":     
     504                                                qval = "Q_0p50" 
     505                                                CR = "CR_0p50" 
     506                                                break    
     507                                        case "0.625":            
     508                                                qval = "Q_0p625" 
     509                                                CR = "CR_0p625" 
     510                                                break    
     511                                        case "1.75":             
     512                                                qval = "Q_1p75" 
     513                                                CR = "CR_1p75" 
     514                                                break    
     515                                        case "2.0":              
     516                                                qval = "Q_2p0" 
     517                                                CR = "CR_2p0" 
     518                                                break    
     519                                        default:                                                        // optional default expression executed 
     520                                                OK=0 
     521                                endswitch 
     522                                 
     523                                if(OK) 
     524                                        AppendToGraph/W=UCALC $CR vs $qval 
     525                                        ModifyGraph marker=19,mode($CR)=4,msize($CR)=3,rgb($CR)=(0,0,0) 
     526                                endif 
     527                                 
     528                        else 
     529                                //take it off of the graph 
     530                                SetDataFolder root:Packages:NIST:USANS:Globals 
     531                                list=WaveList("CR*", ";", "WIN:UCALC") 
     532                                item=StringFromList(0, list ,";")               //should be one item 
     533                                if(strlen(item) != 0) 
     534                                        RemoveFromGraph/W=UCALC $item 
     535                                endif 
     536                        endif 
     537                         
     538                        break 
     539        endswitch 
     540 
     541        SetDataFolder root: 
     542        return 0 
     543End 
     544 
     545 
     546 
    470547// based on the angle ranges above (with non-zero count times) 
    471548// plot where the data points would be 
     
    612689                prob_i = trans*thick*omega*Smeared_inten                        //probability of a neutron in q-bin(i) 
    613690                 
    614                 Variable P_on = sum(prob_i,-inf,inf) 
    615                 Print "P_on = ",P_on 
    616                 fracScat = P_on 
     691//              Variable P_on = sum(prob_i,-inf,inf) 
     692//              Print "P_on = ",P_on 
     693                fracScat = 1-estTrans 
    617694                 
    618695                inten = (Imon*countingTime)*prob_i 
     
    626703                 
    627704                NVAR doABS = root:Packages:NIST:USANS:Globals:U_Sim:g_1D_DoABS 
     705                NVAR plotCR = root:Packages:NIST:USANS:Globals:U_Sim:g_1D_PlotCR 
    628706                NVAR addNoise = root:Packages:NIST:USANS:Globals:U_Sim:g_1D_AddNoise 
    629707                                         
     
    638716                endif 
    639717 
    640                 // convert to absolute scale 
     718                // convert to absolute scale? Maybe not needed 
    641719                // does nothing yet - need Ipeak, Twide 
    642720//              if(doABS) 
     
    645723//                      inten /= kappa 
    646724//              endif 
     725 
     726                // plot as countrate - maybe easier to visualize, and all of the data overlaps 
     727                if(plotCR) 
     728                        inten /= countingTime 
     729                        sigave /= countingTime 
     730                endif 
    647731                 
    648732                GraphSIM() 
     
    650734        else 
    651735                //no function plotted, no simulation can be done 
    652                 DoAlert 0,"No function is selected or plotted, so no simulation is done. The default Sphere function is used." 
    653  
    654                 inten = U_SphereForm(1,9000,6e-6,0,qvals)                
     736                DoAlert 0,"No function is selected or plotted, so no simulation is done. The default power law function is used." 
     737 
     738                inten = U_Power_Law_Model(1e-6,3,0,qvals) 
     739//              inten = U_SphereForm(1,9000,6e-6,0,qvals)                
    655740         
    656741                GraphSIM() 
     
    746831End 
    747832 
     833// better default function 
     834Function U_Power_Law_Model(A,m,bgd,x) : FitFunc 
     835        Variable A, m,bgd,x 
     836//       Input (fitting) variables are: 
     837        //[0] Coefficient 
     838        //[1] (-) Power 
     839        //[2] incoherent background 
     840         
     841//      local variables 
     842        Variable inten, qval 
     843//      x is the q-value for the calculation 
     844        qval = x 
     845//      do the calculation and return the function value 
     846         
     847        inten = A*qval^-m + bgd 
     848        Return (inten) 
     849End 
     850 
     851 
    748852// mimics LoadBT5File 
    749853// creates two other waves to identify the set and the counting time for that set 
     
    843947                ModifyGraph mirror(left)=1 
    844948                ModifyGraph grid=2 
    845                  
    846                 // to make sure that the scales are the same 
    847                 SetAxis bottom 2e-06,0.003 
    848                 SetAxis top 2e-06/5.55e-5,0.003/5.55e-5 
    849                  
     949                ModifyGraph standoff=0 
     950                 
     951                // to make sure that the scales are the same (but fails on zoom) 
     952//              SetAxis bottom 2e-06,0.003 
     953//              SetAxis top 2e-06/5.55e-5,0.003/5.55e-5 
     954                 
     955                SetDrawEnv linefgc= (39321,1,1),dash= 3,linethick= 3 
     956                SetDrawEnv xcoord= bottom 
     957                DrawLine 3e-05,0.01,3e-05,0.99 
     958 
    850959                Label top "Angle" 
    851960                Label bottom "Q (1/A)" 
    852                 Label left "Intensity (1/cm) or Counts" 
     961                Label left "Counts or Count Rate" 
    853962                 
    854963                Legend 
     
    11431252// return the beam intensity based on the sample aperture diameter 
    11441253// 
    1145 /// based on only TWO known values at 0.625 (=25000) and 1.75 in (=57000) diam 
    1146 // 
    1147 // so this is a BAD interpolation!!! 
    1148 // get proper numbers from John 
     1254// based on the equation in John's instrument paper 
    11491255// 
    11501256Function GetUSANSBeamIntensity()         
    11511257 
    11521258        String popStr 
    1153         Variable flux,diam 
     1259        Variable flux,diam,rad 
    11541260         
    11551261        ControlInfo/W=UCALC popup0 
    11561262        popStr = S_Value 
    1157         diam=str2num(popStr) 
    1158          
    1159         // a switch would be easier, but cases need to be integer 
    1160         strswitch(popStr)       // string switch 
    1161                 case "0.625":           // execute if case matches expression 
    1162                         flux=25000 
    1163                         break                                           // exit from switch 
    1164                 case "1.75":            // execute if case matches expression 
    1165                         flux=57000 
    1166                         break 
    1167                 default:                                                        // optional default expression executed 
    1168                         flux=7200+28400*diam 
    1169         endswitch 
     1263        diam=str2num(popStr)            //in inches 
     1264        rad = diam/2*25.4                       //radius in mm 
     1265         
     1266        flux = 662*rad*rad-39.9*rad*rad*rad+0.70*rad*rad*rad*rad 
    11701267         
    11711268        return(flux) 
Note: See TracChangeset for help on using the changeset viewer.