Changeset 534 for sans/Dev


Ignore:
Timestamp:
Jul 7, 2009 6:34:44 PM (14 years ago)
Author:
srkline
Message:

changes to move 2D MonteCarlo? to a hidden feature and include 1D simulation. Skeleton of functionality only. Math still needs to be done as well as interface.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r515 r534  
    109109         
    110110        // for the MC simulation 
     111        Variable/G root:Packages:NIST:SAS:doSimulation  =0              // == 1 if 1D simulated data, 0 if other from the checkbox 
    111112        Variable/G root:Packages:NIST:SAS:gRanDateTime=datetime 
    112113        Variable/G root:Packages:NIST:SAS:gImon = 10000 
    … …  
    124125        Make/O/D/N=10 root:Packages:NIST:SAS:results = 0 
    125126        Make/O/T/N=10 root:Packages:NIST:SAS:results_desc = {"total X-section (1/cm)","SAS X-section (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         
    126128         
    127129        //tick labels for SDD slider 
    … …  
    383385        CheckBox checkLens,value=root:Packages:NIST:SAS:gUsingLenses 
    384386         
    385         CheckBox checkSim,pos={6,175},size={44,14},proc=SimCheckProc,title="MC Simulation?" 
     387        CheckBox checkSim,pos={6,175},size={44,14},proc=SimCheckProc,title="Simulation?" 
    386388        CheckBox checkSim,value=0 
    387389         
    … …  
    686688End 
    687689 
    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 
    689693Function SimCheckProc(ctrlName,checked) : CheckBoxControl 
    690694        String ctrlName 
    … …  
    692696 
    693697        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 
    700713                endif 
    701714        else 
    702                 //get rid of the controls 
     715                //get rid of the controls (just try both) 
    703716                DoWindow MC_SASCALC 
    704717                if(V_flag!=0) 
     718                        NVAR do2D = root:Packages:NIST:SAS:gDoMonteCarlo 
     719                        do2D = 0 
    705720                        KillWindow MC_SASCALC 
    706721                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 
    709729        endif 
    710730 
    … …  
    855875         
    856876        // do the simulation here, or not 
    857         Variable r1,xCtr,yCtr,sdd,pixSize,wavelength 
     877        Variable r1,xCtr,yCtr,sdd,pixSize,wavelength,trans 
    858878        String coefStr,abortStr,str 
    859879 
    860         NVAR doMonteCarlo = root:Packages:NIST:SAS:gDoMonteCarlo                // == 1 if MC, 0 if other from the checkbox 
     880        // now the cases are: simulation (0|1), 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 
    861883        SVAR funcStr = root:Packages:NIST:SAS:gFuncStr          //set by the popup 
    862884 
    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 well-behaved.",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 non-threaded 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 
    936892                else 
    937                         Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 
    938                 endif 
    939                  
    940                 t0 = (stopMSTimer(-2) - t0)*1e-6 
    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 "non-reentrant". 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)*1e-6 
    967                 Printf  "MC sim time = %g seconds\r",t0 
    968                  
    969                 trans = results[8]                      //(n1-n2)/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((p-xCtr)^2+(q-yCtr)^2) < rad) ? 0 : 1          //behind beamstop = 0, away = 1 
     893                        //1D simulation 
    989894                         
    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                 // re-average 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 hard-wired - 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 
    1025933                                Abort "The coefficients and function type do not match. Please correct the selections in the popup menus." 
    1026934                        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 
    1034937                endif 
     938        else 
     939                //no simulation 
    1035940                 
     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 
    1036944                 
    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 
    1039950        endif 
    1040951         
Note: See TracChangeset for help on using the changeset viewer.