Ignore:
Timestamp:
Nov 20, 2008 10:54:27 AM (14 years ago)
Author:
srkline
Message:

Bug fixes in the Monte Carlo routines for SASCALC, and some improvements for functionality, including the calculation on smeared model functions. For various reasons, the MC simulation must use the unsmeared model.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf

    r436 r454  
    2121//   by allowing the data arrays to accumulate. First pass at the XOP is done. Not pretty, not the speediest (5.8x) 
    2222//   but it is functional. Add spinCursor() for long calculations. See WaveAccess XOP example. 
    23 // - the background parameter for the model MUST be zero, or the integration for scattering 
    24 //    power will be incorrect. 
    25 // - fully use the SASCALC input, most importantly, flux on sample. 
     23// X the background parameter for the model MUST be zero, or the integration for scattering 
     24//    power will be incorrect. (now the LAST point in a copy of the coef wave is set to zero, only for the rad_dev calculation 
     25// X fully use the SASCALC input, most importantly, flux on sample. 
    2626// X if no MC desired, still use the selected model 
    2727// X better display of MC results on panel 
    28 // - settings for "count for X seconds" or "how long to 1E6 cts on detector" (but 1E6 is typically too many counts...) 
     28// X settings for "count for X seconds" or "how long to 1E6 cts on detector" (but 1E6 is typically too many counts...) 
     29// X warn of projected simulation time 
    2930// - add quartz window scattering to the simulation somehow 
    30 // - do smeared models make any sense?? Yes, John agrees that they do, and may be used in a more realistic simulation 
     31// -?- do smeared models make any sense?? Yes, John agrees that they do, and may be used in a more realistic simulation 
     32//   -?- but the random deviate can't be properly calculated... 
    3133// - make sure that the ratio of scattering coherent/incoherent is properly adjusted for the sample composition 
    3234//   or the volume fraction of solvent. 
    … …  
    4547// X ask John how to verify what is going on 
    4648// - a number of models are now found to be ill-behaved when q=1e-10. Then the random deviate calculation blows up. 
    47 //   a warning has been added - but the models are better fixed with the limiting value. 
     49//   a warning has been added - but some models need a proper limiting value, and some (power-law) are simply unuseable 
     50//   unless something else can be done. 
    4851// 
    4952// 
    … …  
    8285                 
    8386                // ?? I need explicit wave references? 
     87                // maybe I need to have everything in separate data folders - bu tI haven't tried that. seems like a reach. 
     88                // more likely there is something bad going on in the XOP code. 
    8489                if(i==0) 
    8590                        WAVE inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0 
    … …  
    531536// returns the random deviate as a wave 
    532537// and the total SAS cross-section [1/cm] sig_sas 
    533 Function        CalculateRandomDeviate(func,coef,lam,outWave,SASxs) 
     538Function CalculateRandomDeviate(func,coef,lam,outWave,SASxs) 
    534539        FUNCREF SANSModelAAO_MCproto func 
    535540        WAVE coef 
    … …  
    541546        qu = 4*pi/lam            
    542547         
    543         Make/O/N=(nPts_ran)/D root:Packages:NIST:SAS:Gq,root:Packages:NIST:SAS:xw               // if these waves are 1000 pts, the results are "pixelated" 
    544         WAVE Gq = root:Packages:NIST:SAS:gQ 
    545         WAVE xw = root:Packages:NIST:SAS:xw 
     548//      Make/O/N=(nPts_ran)/D root:Packages:NIST:SAS:Gq,root:Packages:NIST:SAS:xw               // if these waves are 1000 pts, the results are "pixelated" 
     549//      WAVE Gq = root:Packages:NIST:SAS:gQ 
     550//      WAVE xw = root:Packages:NIST:SAS:xw 
     551 
     552// hard-wired into the Simulation directory rather than the SAS folder. 
     553// plotting resolution-smeared models won't work any other way 
     554        Make/O/N=(nPts_ran)/D root:Simulation:Gq,root:Simulation:xw             // if these waves are 1000 pts, the results are "pixelated" 
     555        WAVE Gq = root:Simulation:gQ 
     556        WAVE xw = root:Simulation:xw 
    546557        SetScale/I x (0+1e-6),qu*(1-1e-10),"", Gq,xw                    //don't start at zero or run up all the way to qu to avoid numerical errors 
     558 
     559/// 
     560/// if all of the coefficients are well-behaved, then the last point is the background 
     561// and I can set it to zero here (only for the calculation) 
     562        Duplicate/O coef,tmp_coef 
     563        Variable num=numpnts(coef) 
     564        tmp_coef[num-1] = 0 
     565         
    547566        xw=x                                                                                            //for the AAO 
    548         func(coef,Gq,xw)                                                                        //call as AAO 
     567        func(tmp_coef,Gq,xw)                                                                    //call as AAO 
    549568 
    550569//      Gq = x*Gq                                                                                                       // SAS approximation 
    … …  
    665684        list = RemoveFromList(tmp, list  ,";") 
    666685         
    667 //      tmp = FunctionList("*X",";","KIND:4")           //XOPs, but these shouldn't show up if KIND:10 is used initially 
    668 //      Print "X* = ",tmp 
    669 //      print " " 
    670 //      list = RemoveFromList(tmp, list  ,";") 
    671          
    672686        //non-fit functions that I can't seem to filter out 
    673687        list = RemoveFromList("BinaryHS_PSF11;BinaryHS_PSF12;BinaryHS_PSF22;EllipCyl_Integrand;PP_Inner;PP_Outer;Phi_EC;TaE_Inner;TaE_Outer;",list,";") 
     688//////////////// 
     689 
     690        //simplify the display, forcing smeared calculations behind the scenes 
     691        tmp = FunctionList("Smear*",";","NPARAMS:1")            //smeared dependency calculations 
     692        list = RemoveFromList(tmp, list  ,";") 
    674693 
    675694        if(strlen(list)==0) 
    … …  
    746765End 
    747766 
     767// globals are initialized in SASCALC.ipf 
    748768Window MC_SASCALC() : Panel 
    749769        PauseUpdate; Silent 1           // building window... 
    750         NewPanel /W=(787,44,1088,563)  /N=MC_SASCALC as "SANS Simulator" 
    751         CheckBox MC_check0,pos={11,11},size={98,14},title="Use MC Simulation" 
    752         CheckBox MC_check0,variable= root:Packages:NIST:SAS:gDoMonteCarlo 
    753         SetVariable MC_setvar0,pos={11,38},size={144,15},bodyWidth=80,title="# of neutrons" 
     770        NewPanel /W=(92,556,390,1028)/K=1 as "SANS Simulator" 
     771        SetVariable MC_setvar0,pos={28,73},size={144,15},bodyWidth=80,title="# of neutrons" 
     772        SetVariable MC_setvar0,format="%5.4g" 
    754773        SetVariable MC_setvar0,limits={-inf,inf,100},value= root:Packages:NIST:SAS:gImon 
    755         SetVariable MC_setvar0_1,pos={11,121},size={131,15},bodyWidth=60,title="Thickness (cm)" 
     774        SetVariable MC_setvar0_1,pos={28,119},size={131,15},bodyWidth=60,title="Thickness (cm)" 
    756775        SetVariable MC_setvar0_1,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gThick 
    757         SetVariable MC_setvar0_2,pos={11,93},size={149,15},bodyWidth=60,title="Incoherent XS (cm)" 
     776        SetVariable MC_setvar0_2,pos={28,96},size={149,15},bodyWidth=60,title="Incoherent XS (cm)" 
    758777        SetVariable MC_setvar0_2,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gSig_incoh 
    759         SetVariable MC_setvar0_3,pos={11,149},size={150,15},bodyWidth=60,title="Sample Radius (cm)" 
     778        SetVariable MC_setvar0_3,pos={28,142},size={150,15},bodyWidth=60,title="Sample Radius (cm)" 
    760779        SetVariable MC_setvar0_3,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gR2 
    761         PopupMenu MC_popup0,pos={11,63},size={162,20},proc=MC_ModelPopMenuProc,title="Model Function" 
     780        PopupMenu MC_popup0,pos={13,13},size={165,20},proc=MC_ModelPopMenuProc,title="Model Function" 
    762781        PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()" 
    763         Button MC_button0,pos={17,455},size={50,20},proc=MC_DoItButtonProc,title="Do It" 
    764         Button MC_button1,pos={15,484},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D" 
    765          
     782        Button MC_button0,pos={17,181},size={130,20},proc=MC_DoItButtonProc,title="Do MC Simulation" 
     783        Button MC_button1,pos={181,181},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D" 
     784        SetVariable setvar0_3,pos={105,484},size={50,20},disable=1 
     785        GroupBox group0,pos={15,42},size={267,130},title="Monte Carlo" 
     786        SetVariable cntVar,pos={190,73},size={80,15},proc=CountTimeSetVarProc,title="time(s)" 
     787        SetVariable cntVar,format="%d" 
     788        SetVariable cntVar,limits={1,10,1},value= root:Packages:NIST:SAS:gCntTime 
     789         
     790        String fldrSav0= GetDataFolder(1) 
    766791        SetDataFolder root:Packages:NIST:SAS: 
    767         Edit/W=(13,174,284,435)/HOST=# results_desc,results 
    768         ModifyTable width(Point)=0,width(results_desc)=150 
    769         SetDataFolder root: 
     792        Edit/W=(13,217,283,450)/HOST=# results_desc,results 
     793        ModifyTable format(Point)=1,width(Point)=0,width(results_desc)=150 
     794        SetDataFolder fldrSav0 
    770795        RenameWindow #,T_results 
    771796        SetActiveSubwindow ## 
    772797EndMacro 
     798 
     799 
     800Function CountTimeSetVarProc(sva) : SetVariableControl 
     801        STRUCT WMSetVariableAction &sva 
     802 
     803        switch( sva.eventCode ) 
     804                case 1: // mouse up 
     805                case 2: // Enter key 
     806                case 3: // Live update 
     807                        Variable dval = sva.dval 
     808 
     809                        // get the neutron flux, multiply, and reset the global for # neutrons 
     810                        NVAR imon=root:Packages:NIST:SAS:gImon 
     811                        imon = dval*beamIntensity() 
     812                         
     813                        break 
     814        endswitch 
     815 
     816        return 0 
     817End 
     818 
    773819 
    774820Function MC_ModelPopMenuProc(pa) : PopupMenuControl 
    … …  
    794840                case 2: // mouse up 
    795841                        // click code here 
     842                        NVAR doMC = root:Packages:NIST:SAS:gDoMonteCarlo 
     843                        doMC = 1 
    796844                        ReCalculateInten(1) 
     845                        doMC = 0                //so the next time won't be MC 
    797846                        break 
    798847        endswitch 
    … …  
    814863        return 0 
    815864End 
     865 
     866// after a 2d data image is averaged in the usual way, take the waves and generate a "fake" folder of the 1d 
     867// data, to appear as if it was loaded from a real data file. 
     868// 
     869// currently only works with SANS data, but can later be expanded to generate fake USANS data sets 
     870// 
     871Function        Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,dataFolder) 
     872        WAVE qval,aveint,sigave,sigmaQ,qbar,fSubs 
     873        String dataFolder 
     874 
     875        String baseStr=dataFolder 
     876        if(DataFolderExists("root:"+baseStr)) 
     877                SetDataFolder $("root:"+baseStr) 
     878        else 
     879                NewDataFolder/S $("root:"+baseStr) 
     880        endif 
     881 
     882        ////overwrite the existing data, if it exists 
     883        Duplicate/O qval, $(baseStr+"_q") 
     884        Duplicate/O aveint, $(baseStr+"_i") 
     885        Duplicate/O sigave, $(baseStr+"_s") 
     886//      Duplicate/O sigmaQ, $(baseStr+"sq") 
     887//      Duplicate/O qbar, $(baseStr+"qb") 
     888//      Duplicate/O fSubS, $(baseStr+"fs") 
     889 
     890        // need to switch based on SANS/USANS 
     891        if (isSANSResolution(sigave[0]))                //checks to see if the first point of the wave is <0] 
     892                // make a resolution matrix for SANS data 
     893                Variable np=numpnts(qval) 
     894                Make/D/O/N=(np,4) $(baseStr+"_res") 
     895                Wave res=$(baseStr+"_res") 
     896                 
     897                res[][0] = sigmaQ[p]            //sigQ 
     898                res[][1] = qBar[p]              //qBar 
     899                res[][2] = fSubS[p]             //fShad 
     900                res[][3] = qval[p]              //Qvalues 
     901                 
     902                // keep a copy of everything in SAS too... the smearing wrapper function looks for  
     903                // data in folders based on waves it is passed - an I lose control of that 
     904                Duplicate/O res, $("root:Packages:NIST:SAS:"+baseStr+"_res") 
     905                Duplicate/O qval,  $("root:Packages:NIST:SAS:"+baseStr+"_q") 
     906                Duplicate/O aveint,  $("root:Packages:NIST:SAS:"+baseStr+"_i") 
     907                Duplicate/O sigave,  $("root:Packages:NIST:SAS:"+baseStr+"_s") 
     908        else 
     909                //the data is USANS data 
     910                // nothing done here yet 
     911//              dQv = -$w3[0] 
     912                 
     913//              USANS_CalcWeights(baseStr,dQv) 
     914                 
     915        endif 
     916 
     917        //clean up               
     918        SetDataFolder root: 
     919         
     920End 
     921 
     922 
    816923 
    817924/////UNUSED, testing routines that have not been updated to work with SASCALC 
Note: See TracChangeset for help on using the changeset viewer.