Changeset 534
 Timestamp:
 Jul 7, 2009 6:34:44 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/SASCALC.ipf
r515 r534 109 109 110 110 // for the MC simulation 111 Variable/G root:Packages:NIST:SAS:doSimulation =0 // == 1 if 1D simulated data, 0 if other from the checkbox 111 112 Variable/G root:Packages:NIST:SAS:gRanDateTime=datetime 112 113 Variable/G root:Packages:NIST:SAS:gImon = 10000 … … 124 125 Make/O/D/N=10 root:Packages:NIST:SAS:results = 0 125 126 Make/O/T/N=10 root:Packages:NIST:SAS:results_desc = {"total Xsection (1/cm)","SAS Xsection (1/cm)","number that scatter","number that reach detector","avg # times scattered","fraction single coherent","fraction double coherent","fraction multiple scattered","fraction transmitted","detector counts w/o beamstop"} 127 126 128 127 129 //tick labels for SDD slider … … 383 385 CheckBox checkLens,value=root:Packages:NIST:SAS:gUsingLenses 384 386 385 CheckBox checkSim,pos={6,175},size={44,14},proc=SimCheckProc,title=" MCSimulation?"387 CheckBox checkSim,pos={6,175},size={44,14},proc=SimCheckProc,title="Simulation?" 386 388 CheckBox checkSim,value=0 387 389 … … 686 688 End 687 689 688 //simulation ccontrols as a control bar that toggles on/off to the right 690 //simulation controls as a control bar that toggles on/off to the right 691 // 692 // depending on the state of the 2D flag, open the 1d or 2d control panel 689 693 Function SimCheckProc(ctrlName,checked) : CheckBoxControl 690 694 String ctrlName … … 692 696 693 697 if(checked) 694 //MoveWindow/W=SASCALC 5,44,763,570 //to resize 695 // draw the controls 696 DoWindow/F MC_SASCALC 697 if(V_flag==0) 698 Execute "MC_SASCALC()" 699 AutoPositionWindow/M=1/R=SASCALC MC_SASCALC 698 NVAR do2D = root:Packages:NIST:SAS:gDoMonteCarlo 699 700 if(do2D) 701 DoWindow/F MC_SASCALC 702 if(V_flag==0) 703 Execute "MC_SASCALC()" //sets the variable 704 AutoPositionWindow/M=1/R=SASCALC MC_SASCALC 705 endif 706 else 707 //draw 1D panel 708 DoWindow/F Sim_1D_Panel 709 if(V_flag==0) 710 Execute "Sim_1D_Panel()" 711 AutoPositionWindow/M=1/R=SASCALC Sim_1D_Panel 712 endif 700 713 endif 701 714 else 702 //get rid of the controls 715 //get rid of the controls (just try both) 703 716 DoWindow MC_SASCALC 704 717 if(V_flag!=0) 718 NVAR do2D = root:Packages:NIST:SAS:gDoMonteCarlo 719 do2D = 0 705 720 KillWindow MC_SASCALC 706 721 endif 707 //MoveWindow/W=SASCALC 5,44,463,570 //to resize 708 //KillWindow SASCALC#T_results 722 DoWindow Sim_1D_Panel 723 if(V_flag!=0) 724 KillWindow Sim_1D_Panel 725 endif 726 727 NVAR doSim = root:Packages:NIST:SAS:doSimulation 728 doSim=0 709 729 endif 710 730 … … 855 875 856 876 // do the simulation here, or not 857 Variable r1,xCtr,yCtr,sdd,pixSize,wavelength 877 Variable r1,xCtr,yCtr,sdd,pixSize,wavelength,trans 858 878 String coefStr,abortStr,str 859 879 860 NVAR doMonteCarlo = root:Packages:NIST:SAS:gDoMonteCarlo // == 1 if MC, 0 if other from the checkbox 880 // now the cases are: simulation (01), 1D (default) or 2D (hidden) 881 NVAR doMonteCarlo = root:Packages:NIST:SAS:gDoMonteCarlo // == 1 if 2D MonteCarlo set by hidden flag 882 NVAR doSimulation = root:Packages:NIST:SAS:doSimulation // == 1 if 1D simulated data, 0 if other from the checkbox 861 883 SVAR funcStr = root:Packages:NIST:SAS:gFuncStr //set by the popup 862 884 863 if(doMonteCarlo == 1) 864 WAVE rw=root:Packages:NIST:SAS:realsRead 865 866 NVAR imon = root:Packages:NIST:SAS:gImon 867 NVAR thick = root:Packages:NIST:SAS:gThick 868 NVAR sig_incoh = root:Packages:NIST:SAS:gSig_incoh 869 NVAR r2 = root:Packages:NIST:SAS:gR2 870 871 r1 = rw[24]/2/10 // sample diameter convert diam in [mm] to radius in cm 872 xCtr = rw[16] 873 yCtr = rw[17] 874 sdd = rw[18]*100 //conver header of [m] to [cm] 875 pixSize = rw[10]/10 // convert pix size in mm to cm 876 wavelength = rw[26] 877 coefStr = MC_getFunctionCoef(funcStr) 878 879 if(!MC_CheckFunctionAndCoef(funcStr,coefStr)) 880 doMonteCarlo = 0 //we're getting out now, reset the flag so we don't continually end up here 881 Abort "The coefficients and function type do not match. Please correct the selections in the popup menus." 882 endif 883 884 Variable sig_sas 885 // FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr) //a wrapper for the structure version 886 FUNCREF SANSModelAAO_MCproto func=$(funcStr) //unsmeared 887 WAVE results = root:Packages:NIST:SAS:results 888 WAVE linear_data = root:Packages:NIST:SAS:linear_data 889 WAVE data = root:Packages:NIST:SAS:data 890 891 results = 0 892 linear_data = 0 893 894 CalculateRandomDeviate(func,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",SIG_SAS) 895 if(sig_sas > 100) 896 sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is wellbehaved.",sig_sas 897 Abort abortStr 898 endif 899 900 WAVE ran_dev=$"root:Packages:NIST:SAS:ran_dev" 901 902 Make/O/D/N=5000 root:Packages:NIST:SAS:nt=0,root:Packages:NIST:SAS:j1=0,root:Packages:NIST:SAS:j2=0 903 Make/O/D/N=100 root:Packages:NIST:SAS:nn=0 904 Make/O/D/N=11 root:Packages:NIST:SAS:inputWave=0 905 906 WAVE nt = root:Packages:NIST:SAS:nt 907 WAVE j1 = root:Packages:NIST:SAS:j1 908 WAVE j2 = root:Packages:NIST:SAS:j2 909 WAVE nn = root:Packages:NIST:SAS:nn 910 WAVE inputWave = root:Packages:NIST:SAS:inputWave 911 912 inputWave[0] = imon 913 inputWave[1] = r1 914 inputWave[2] = r2 915 inputWave[3] = xCtr 916 inputWave[4] = yCtr 917 inputWave[5] = sdd 918 inputWave[6] = pixSize 919 inputWave[7] = thick 920 inputWave[8] = wavelength 921 inputWave[9] = sig_incoh 922 inputWave[10] = sig_sas 923 924 linear_data = 0 //initialize 925 926 Variable t0,trans 927 928 // get a time estimate, and give the user a chance to exit if they're unsure. 929 t0 = stopMStimer(2) 930 inputWave[0] = 1000 931 NVAR useXOP = root:Packages:NIST:SAS:gUse_MC_XOP //if zero, will use nonthreaded Igor code 932 933 if(useXOP) 934 //use a single thread, otherwise time is dominated by overhead 935 Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 885 if(doSimulation == 1) 886 if(doMonteCarlo == 1) 887 //2D simulation (in MultiScatter_MonteCarlo_2D.ipf 888 889 Simulate_2D_MC(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs) 890 891 //end 2D simulation 936 892 else 937 Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 938 endif 939 940 t0 = (stopMSTimer(2)  t0)*1e6 941 t0 *= imon/1000/ThreadProcessorCount //projected time, in seconds (using threads for the calculation) 942 inputWave[0] = imon //reset 943 944 if(t0>10) 945 sprintf str,"The simulation will take approximately %d seconds.\r Proceed?",t0 946 DoAlert 1,str 947 if(V_flag == 2) 948 doMonteCarlo = 0 949 reCalculateInten(1) //come back in and do the smeared calculation 950 return(0) 951 endif 952 endif 953 954 linear_data = 0 //initialize 955 // threading crashes!!  there must be some operation in the XOP that is not threadSafe. What, I don't know... 956 // I think it's the ran() calls, being "nonreentrant". So the XOP now defines two separate functions, that each 957 // use a different rng. This works. 1.75x speedup. 958 t0 = stopMStimer(2) 959 960 if(useXOP) 961 Monte_SANS_Threaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 962 else 963 Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 964 endif 965 966 t0 = (stopMSTimer(2)  t0)*1e6 967 Printf "MC sim time = %g seconds\r",t0 968 969 trans = results[8] //(n1n2)/n1 970 if(trans == 0) 971 trans = 1 972 endif 973 974 Print "counts on detector, including transmitted = ",sum(linear_data,inf,inf) 975 976 // linear_data[xCtr][yCtr] = 0 //snip out the transmitted spike 977 // Print "counts on detector not transmitted = ",sum(linear_data,inf,inf) 978 979 // or simulate a beamstop 980 NVAR MC_BS_in = root:Packages:NIST:SAS:gBeamStopIn //if zero, beam stop is "out", as in a transmission measurement 981 982 Variable rad=beamstopDiam()/2 //beamstop radius in cm 983 if(MC_BS_in) 984 rad /= 0.5 //convert cm to pixels 985 rad += 0. // (no  it cuts off the low Q artificially) add an extra pixel to each side to account for edge 986 Duplicate/O linear_data,root:Packages:NIST:SAS:tmp_mask//,root:Packages:NIST:SAS:MC_linear_data 987 WAVE tmp_mask = root:Packages:NIST:SAS:tmp_mask 988 tmp_mask = (sqrt((pxCtr)^2+(qyCtr)^2) < rad) ? 0 : 1 //behind beamstop = 0, away = 1 893 //1D simulation 989 894 990 linear_data *= tmp_mask 991 endif 992 993 results[9] = sum(linear_data,inf,inf) 994 // Print "counts on detector not behind beamstop = ",results[9] 995 996 // convert to absolute scale 997 Variable kappa //,beaminten = beamIntensity() 998 // kappa = beamInten*pi*r1*r1*thick*(pixSize/sdd)^2*trans*(iMon/beaminten) 999 kappa = thick*(pixSize/sdd)^2*trans*iMon 1000 1001 //use kappa to get back to counts => linear_data = round(linear_data*kappa) 1002 Note/K linear_data ,"KAPPA="+num2str(kappa)+";" 1003 1004 NVAR rawCts = root:Packages:NIST:SAS:gRawCounts 1005 if(!rawCts) //go ahead and do the linear scaling 1006 linear_data = linear_data / kappa 1007 endif 1008 data = linear_data 1009 1010 // reaverage the 2D data 1011 S_CircularAverageTo1D("SAS") 1012 1013 // put the new result into the simulation folder 1014 Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,"Simulation") 1015 endif 1016 1017 1018 if(doMonteCarlo != 1) 1019 if(exists(funcStr) != 0) 1020 FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr) //a wrapper for the structure version 1021 // FUNCREF SANSModelAAO_MCproto func=$(funcStr) //unsmeared 1022 coefStr = MC_getFunctionCoef(funcStr) 1023 1024 if(!MC_CheckFunctionAndCoef(funcStr,coefStr)) 895 if(exists(funcStr) != 0) 896 FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr) //a wrapper for the structure version 897 // FUNCREF SANSModelAAO_MCproto func=$(funcStr) //unsmeared 898 coefStr = MC_getFunctionCoef(funcStr) 899 900 if(!MC_CheckFunctionAndCoef(funcStr,coefStr)) 901 Abort "The coefficients and function type do not match. Please correct the selections in the popup menus." 902 endif 903 Wave inten=$"root:Simulation:Simulation_i" // this will exist and send the smeared calculation to the corect DF 904 func($coefStr,inten,qval) 905 906 NVAR imon = root:Packages:NIST:SAS:gImon 907 NVAR thick = root:Packages:NIST:SAS:gThick 908 909 WAVE rw=root:Packages:NIST:SAS:realsRead 910 WAVE nCells=root:Packages:NIST:SAS:nCells 911 912 ///////Trans is hardwired  need to input from somewhere!!! 913 trans = 0.8 914 915 916 pixSize = rw[10]/10 // convert pix size in mm to cm 917 sdd = rw[18]*100 //conver header of [m] to [cm] 918 919 duplicate/O qval prob_i 920 prob_i = trans*thick*nCells*(pixSize/sdd)^2*inten //probability of a neutron in bin(i) 921 922 Variable P_on = sum(prob_i,inf,inf) 923 Print "P_on = ",P_on 924 925 aveint = Imon*P_on*prob_i 926 927 928 929 aveint *= fSubS 930 931 sigave = sqrt(aveint) //reset for model calculation 932 else 1025 933 Abort "The coefficients and function type do not match. Please correct the selections in the popup menus." 1026 934 endif 1027 Wave inten=$"root:Simulation:Simulation_i" // this will exist and send the smeared calculation to the corect DF 1028 func($coefStr,inten,qval) 1029 inten *= fSubS 1030 aveint = inten //but aveint is displayed 1031 else 1032 aveint = S_Debye(1000,100,0.0,qval) 1033 aveint *= fSubS // multiply either estimate by beamstop shadowing 935 936 // end 1D simulation 1034 937 endif 938 else 939 //no simulation 1035 940 941 aveint = S_Debye(1000,100,0.0,qval) 942 aveint *= fSubS // multiply either estimate by beamstop shadowing 943 sigave = 0 //reset for model calculation 1036 944 1037 WAVE sigave=root:Packages:NIST:SAS:sigave 1038 sigave = 0 //reset for model calculation 945 //end no simulation 946 endif 947 948 if(doMonteCarlo == 1) 949 1039 950 endif 1040 951
Note: See TracChangeset
for help on using the changeset viewer.